Osteoporosis is a calcium and metabolic disorder characterized by decreased bone mass, enhanced risk of bone fragility and susceptibility to fracture [1,2] caused by a failure of bone homeostasis, which is due to both an increase in osteoclastic bone resorption and a decrease in osteoblastic bone formation [3,4]. Osteoporotic fractures are an important cause of morbidity and mortality, particularly in elderly women and men [5].

Osteoporosis is involved in the interactions of multiple genetic and environmental risk factors [6,7]. Recently, genetic factors have attracted much attention of many investigators due to their high importance in the pathogenesis of osteoporosis [8,9]. Severe osteoporosis may be related to mutation in a single gene, otherwise, bone mineral density (BMD) or bone mineral mass can be accounted for by the common genetic variations in multiple genes with relatively small effects or by the rare genetic variations in specific genes with relatively large effects [10,11]. Heritability studies in twins and families have demonstrated that between 50% and 85% of the variance in peak BMD is genetically determined [11,12].

The identification of genetic variants that contribute to osteoporosis and BMD phenotypes can be helpful not only for elucidation of the molecular mechanisms of osteoporosis, but also for the development of effective treatment for osteoporosis. So far, many osteoporosis susceptibility loci and genes have been identified by genome-wide linkage analysis [13,14] , candidate gene association studies [15] and genome-wide association studies (GWAS) [16-18], and recent important and representative findings through molecular genetic studies of gene identification for osteoporosis have been well summarized [19]. However, the vast majority of the monogenic and polygenic genetic factors in osteoporosis still remain to be discovered.

In this study, we aimed to identify the novel genes for susceptibility to osteoporosis that influence the pathogenesis of osteoporosis in humans in a particular way, by using the in vitro and in vivo models for osteoporosis. As an in vitro model for osteoporosis, we used the glucocorticoid (GC)-treated mouse osteoblastic MC3T3-E1 cell line in this study. GC causes osteoporosis-like bone loss due to decreasing bone formation and increasing bone resorption [20,21]. GC-induced apoptosis in osteoblastic cells has been reported [22]. The mouse MC3T3-E1 cell line has been demonstrated to be a suitable in vitro model of osteoblast development due to typical osteoblast differentiation and formation of a bone-like mineralized extracellular matrix [23,24]. Therefore, GC-treated osteoblastic MC3T3-E1 cells may provide a useful in vitro system for the first screening step of the differentially expressed genes (DEGs) between controls and osteoporosis models. In case of the in vivo models for osteoporosis, many animal models have been established for osteoporosis research [25,26]. Among the mouse models for osteoporosis, ovariectomized mice and GC-induced mice represent a reliable in vivo model to investigate bone loss in osteoporosis [27].

We first screened and selected the candidate genes that have a higher possibility to be linked with osteoporosis by comparative analysis of gene expression between controls and in vitro osteoporosis model. Subsequently, we carried out evaluation of the candidate genes using the in vivo osteoporosis model. Finally, we performed an association analysis between the genetic variants in the selected genes and bone density and osteoporosis in human subjects. This approach may provide an accurate identification of novel genes for susceptibility to osteoporosis.

1. Cell culture

Mouse osteoblastic MC3T3-E1 cells were purchased from the RIKEN cell bank (Tsukuba, Japan) and grown in α-MEM medium supplemented with 10% FBS, penicillin (100 U/mL) and streptomycin (100 µg/mL). Cultured cells were incubated at 37°C in a humidified atmosphere containing 5% CO₂.

2. Annealing control primer (ACP)-based differential display reverse-transcription polymerase chain reaction (RT-PCR)

ACP-based differential display RT-PCR was carried out using the predesigned arbitrary ACP (Seegene, Seoul, Korea) [28]. First-strand cDNA synthesis was performed using the primer dT-ACT1 (5’-CTGTGAATGCTGCGACTACGATIIIII(T)₁₈-3’: ‘I’ is inosine) and 1 µL of M-MLV reverse transcriptase (200 U/µL) (Fermentas, Burlington, Canada). PCR amplification was conducted using GeneFishing DEG kit (Seegene) in 20 µL reaction volumes containing 10 µL of 2× SeeAmp ACP Master Mix; 2 µL of 5 µM each arbitrary ACP; 1 µL of 10 µM dT-ACP2 (5’-CTGTGAATGCTGCGACTACGATIIIII(T)₁₅-3’,); and 3 µL of diluted first-strand cDNA. Each kit comprises 120 different arbitrary annealing control primers. The amplified PCR products were separated in a 2% agarose gel and stained with ethidium bromide.

3. Quantitative reverse transcription-PCR (qRT-PCR)

Total RNAs from cells were treated with RNase-free DNase I (Invitrogen) at room temperature for 15 minutes to avoid amplification of genomic DNA; denatured at 70°C for 10 minutes; and subsequently reverse transcribed by Superscript II reverse transcriptase (Invitrogen) with 0.5 µg of oligo(dT)_15-18 primer in a volume of 20 µL according to manufacturer instructions.

The specific primers used for qRT-PCR were as follows: 5’-CACAGGGTGCCATGTACCG-3’ and 5’-GAGGTCCTTGCCATACAGGG-3’ for mouse AnxA6; 5’-GCTCTCCGCCGAAGTTAAGAA-3’ and 5’-TTCGCACAATATGATGCCGTC-3’ for mouse Cnn3; 5’-CAATTTGCCCTCAGGGGTAAC-3’ and 5’-TCCTCGGGAAAACCAGACTCA-3’ for mouse Col5a1; 5’-GCTCCTGATTGGGGGACTCT-3’ and 5’-CCAACACGAAATACACGTTGAC-3’ for mouse Col6a2; 5’-ACACTGGGCTTCATCATGCC-3’ and 5’- ACTGCGAAGATCATCCTCAGG-3’ for mouse Gper; 5’-CCCTGAAGACTCGGGCCTA-3’ and 5’-CAATTACAAGCGAAATGAGAGCC-3’ for mouse Kitl; 5’-TGATCGAGGGCCGTCAGTTAT-3’ and 5’-CTGTCACTCACCTTAAATTCCCC-3’ for mouse Myof; 5’-AGCACCATCATCCCGGAATAC-3’ and 5’-GTACCAGGACTGGCTCGTTTG-3’ for mouse Nfib; 5’-ACCGGACAGCTCGTTTTGG-3’ and 5’-AGGGGACAGTACAAGTCACCC-3’ for mouse Scara5; 5’-TGGAAACCATGATGCTTACGTT-3’ and 5’-GAAGCCCACTTTGCCATCTC-3’ for mouse S100a10; 5’-GCTGCCTCCGAGTTCTACAG-3’ and 5’- GCAGGGATTCGGTCACAGAG-3’ for mouse Eno1. To normalize the efficiency of qRT-PCR reactions, the mouse Gadph gene was used as an internal standard with the following primers: 5’-TGACCACAGTCCATGCCATC-3’ and 5’-GACGGACACATTGGGGGTAG-3’. qRT-PCR was performed using SYBR Green PCR premix (Takara, Shiga, Japan). All measurements were performed in triplicate.

4. Ovariectomized (OVX) mouse model

The OVX (n=10) and sham-operated (Sham, n=10) 8-weeks-old female ddY mice were purchased from Shizuoka Laboratory Center Inc. (Hamamatsu, Japan). Mice were maintained on a diet (5.0 g/day) of Formula-M07 (Feedlab Co., Ltd., Hanam, Korea) and tap water (15 mL/day). All mice were housed individually in clear plastic cages under controlled temperature (23±2°C), humidity (55±5%), and illumination (12-hour light/dark cycle). After 8 weeks of feeding, the BMD between the two groups of mice was measured. The animal research protocol was approved by the Animal Care and Use Committee of the Ajou University School of Medicine, and all experiments were conducted in accordance with the institutional guidelines established by the Committee.

5. Measurement of bone mineral density and tissue sample preparation

Whole body BMD of mice was measured using a PIXI-mus bone densitometer (GE Lunar, Madison, WI, USA). After anesthetization using tiletamine/zolazepam (Zoletil; Virbac Laboratories, Carros, France), the mice were placed on the specimen tray for measurements. All mice were placed carefully in the same position. After measurement of BMD, the mice were killed by CO₂ asphyxiation and cervical dislocation. Mice femurs were excised, and the isolated femur bones and skeletal muscles were then frozen by liquid nitrogen and deep-freeze, respectively. The frozen samples were homogenized using a porcelain mortar and pestle and then lysed using RIPA lysis buffer and used for Western blot analysis.

6. Human subjects

The subjects from the Korean Association Resource (KARE) study which were used in this study have been described in the previous report [18]. The participants were recruited from two community-based epidemiological cohorts, the rural community of Ansung and the urban community of Ansan cities. A total of 3,570 women subjects were investigated in this study. The basic characteristics of the study subjects are described in Table 1.

Bone density was estimated by T-score by dividing the difference of measured speed of sound (SOS) from mean SOS in healthy young adult population by the standard deviation of SOS in young adult population. Bone SOS was measured by quantitative ultrasound at the distal radius and mid-shaft tibia, using the Omnisense 7000P QUS (Sunlight Medical Ltd, Tel-Aviv, Israel). For the case-control analysis of osteoporosis, the subjects whose bone density T-scores at either the distal radius or mid-shaft tibia were less than –2.5 standard deviation (SD) were allocated to case and the subjects whose bone densities T-scores at both the distal radius and mid-shaft tibia were more than –1 SD were allocated to control, according to the general diagnostic categories to be established for adult women [29]. This study was approved by the institutional review board of the Korean National Institute of Health (KBN-2017-046). Written informed consent was obtained from all subjects.

7. Genotyping and selection of SNPs

The genotype data were provided by the Center for Genome Science, the Korea National Institute of Health. The detailed genotyping and quality control processes have been described in the previous report [18]. Briefly, most DNA samples were isolated from the peripheral blood of participants and genotyped using the Affymetix Genome-Wide Human SNP array 5.0 (Affymetrix, Santa Clara, CA, USA). The accuracy of the genotyping was calculated by Bayesian Robust Linear Modeling using the Mahalanobis Distance genotyping algorithm [30]. The SNPs in the 7 genes that we analyzed were selected based on their locations within the gene boundary (5 kb upstream and downstream of the first and last exons, respectively) according to NCBI human genome build 36. The locations of the SNPs were validated with the Ensemble BioMart database (http://www.ensembl.org/biomart).

8. Statistical analysis

In the qRT-PCR analysis, all experiments were repeated at least 3 times unless stated otherwise and results were presented as the means±SD as indicated. Statistical significance between groups was calculated by a two-tailed Student’s t-test. Probability values less than 0.05 (P<0.05) were considered statistically significant.

Statistical analyses for association studies were performed using the PLINK version 1.07 (http://pngu.mgh.harvard.edu/~purcell/plink) and IBM SPSS Statistics for Windows, Version 25.0 (IBM Co., Armonk, NY, USA). Linear regression was used for quantitative analysis of bone density, controlling for cohort and age as covariates. Logistic regression was used for case-control analysis of osteoporosis. All association tests were performed under the additive, dominant and recessive models, and P-values were adjusted for multiple tests by using the Bonferroni-corrected significance level (P<0.00185). The Haploview version 4.2 program (Whitehead Institute for Biomedical Research, Cambridge, MA, USA) was used to examine the structure of the linkage disequilibrium (LD) block [31] using the KARE genotype data and the HapMap database (International HapMap Project, http://www.hapmap.org/). We examined the LD coefficient r2 between all pairs of biallelic loci [32].

1. Study design

The flow chart of the study design is shown in Fig. 1. We carried out a series of experiments on cell line model, mouse model and humans step-by-step for identifying novel genes for susceptibility to osteoporosis. The first experiment was the screening and identification of the DEGs in Dex-treated osteoblastic MC3T3-E1 cell line, using a RT-PCR-based gene expression differential display approach, the ACP-based PCR GeneFishing DEG screening method. Next, the identified DEGs were validated by quantitative real-time PCR with the gene-specific primers in the Dex-treated cells. In the next step, we tried to evaluate the accuracy of the identified DEGs in vivo and carried out quantitative real-time PCR with the gene-specific primers in the ovariectomized mice. Lastly, to determine whether the genetic variations of the selected DEGs were associated with bone density and osteoporosis, we performed association analysis in a large Korean Women’s Cohort (n=3,570).

2. Screening and identification of the DEGs in Dex-treated osteoblastic cell line model

To identify the DEGs in the in vitro osteoporosis model, we performed a whole-genome comparative expression study using a RT-PCR based gene expression differential display, the ACP-based PCR GeneFishing DEG screening method [28]. To establish an in vitro osteoporosis model, mouse osteoblastic cell line, MC3T3-E1, was treated with 1 µM of synthetic GC, Dex for 2 days. Total RNAs were isolated from the cells and used for first-strand cDNA synthesis. The first-strand cDNAs were subjected to gene expression differential display.

Using 120 arbitrary ACP primers, GeneFishing DEG screening was performed and a total of 10 DEGs that showed clear differences between the two treatment groups were found. The gel images for the 10 DEGs were shown in Fig. 2: 9 DEGs had increased mRNA expression levels in the Dex-treated cells compared with the controls and 1 DEG showed decreased mRNA expression level in the Dex-treated cells.

To identify the DEGs, the RT-PCR bands were extracted, re-amplified, and PCR fragments were isolated from gels, cloned, and sequenced. BLASTN and BLASTX searches in the NCBI GenBank revealed that all the 10 DEGs were known genes as listed in Table 2. The expression levels of Annexin A6 (AnxA6), Calponin 3 (Cnn3), Collagen type V alpha 1 (Col5a1), Collagen type VI alpha 2 (Col6a2), Kit ligand (Kitl), Myoferlin (Myof), Nuclear factor I/B (Nfib), Scavenger receptor class A member 5 (Scara5), and S100 calcium binding protein A10 (S100a10), were increased in the Dex-treated cells, however the expression level of Enolase 1 (Eno1) was decreased in the Dex-treated cells.

3. Validation of the identified DEGs in the cell line model by quantitative real-time RT-PCR

To confirm the efficacy and accuracy of screening by the ACP-based differential display RT-PCR, fluorescence-monitored quantitative real-time RT-PCR analysis was employed for the 10 genes. Gene-specific primers were designed to amplify RT-PCR products ranging from 100 to 250 bp. Quantitative real-time RT-PCR results of the 10 genes are shown in Fig. 3 and are presented as relative ratios compared with the mouse Gapdh gene (internal control) with a value of 1.0. The mRNA expression levels of 8 genes, AnxA6, Col5a1, Col6a2, Kitl, Myof, Nfib, Scara5 and S100a10, were increased significantly in the Dex-treated cells compared with the controls, and Eno1 gene expression was decreased significantly, thereby indicating that these results are consistent with the results of ACP-based differential display shown in Fig. 2 and Table 2. The expression level of the Cnn3 gene, however, was not different between the Dex-treated and untreated cells.

4. Evaluation of the identified DEGs in the ovariectomized mouse model by quantitative real-time RT-PCR

To evaluate the accuracy of the identified DEGs in vivo, we carried out comparative analysis of gene expression levels in the identified DEGs between the OVX mice group and Sham control mice group. Each of the ten 8-week-old female OVX and sham ddY mice was purchased and maintained in our laboratory for 4 weeks. Five OVX mice died during maintenance. At 16 weeks after ovariectomy, whole-body BMD of the 10 sham mice and 5 OVX mice was calculated using a PIXI-mus bone densitometer (Fig. 4A). The OVX mice group demonstrated a significant reduction in BMD (approximately 1.5-fold) compared with the sham-operated mice group. After scarifying the mice, the femur bones were excised and frozen in liquid nitrogen. Total RNAs were isolated from the frozen samples.

The mRNA expression levels of the identified DEGs were examined using the in vivo samples by quantitative RT-PCR analysis (Fig. 4B). The expression levels of 7 genes, AnxA6, Col5a1, Col6a2, Eno1, Kitl, Myof, and Scara5, among the 10 identified DEGs were significantly altered in the OVX group compared with the sham group, but the expression level of Nfib and S100a10 genes were not significantly different between the groups. These results suggested that these 7 genes may be involved in the development of osteoporosis in the OVX mice model.

5. Association analysis of the genetic variation in the identified DEGs with bone density and osteoporosis in humans

We finally selected the 7 genes, AnxA6, Col5a1, Col6a2, Eno1, Kitl, Myof, and Scara5, as the target genes for further study in human subjects. To investigate whether the genetic variations in these 7 selected genes influenced the bone density and susceptibility to osteoporosis, we performed the quantitative trait analysis for bone density and osteoporosis case-control analysis for the SNPs of the 7 genes in the KARE Women’s Study Cohort (3,570 subjects). The basic characteristics of the study subjects are shown in Table 1. The mean age of the women subjects was 51.02 years, the mean bone density estimated by T-score at the distal radius (BD-RT) was 0.20±1.55, and the mean bone density estimated by T-score at the midshaft tibia (BD-TT) was –0.81±1.55.

Linear regression analysis was used to associate the genotypes with bone density traits, controlling for age and cohort as covariates. The 116 SNPs were genotyped in the 7 genes (Supplementary Table 1). The genotyped 116 SNPs of the 7 genes were partitioned into a total 27 LD blocks, which was demonstrated by the Haplotype and PLINK program using the KARE data. Therefore, the Bonferroni-corrected significance P-value threshold was calculated as 0.00185 (0.05/27 LD blocks).

The results of association analysis between the 116 SNPs in the 7 genes and bone density in the 3,570 KARE women subjects are summarized in Table 3. Total 12 SNPs in the 4 genes (2 SNPs in ANXA6, 3 SNPs in COL5A1, 6 SNPs in MYOF and 1 SNP in SCARA5) were significantly associated with BD-RT trait, and total 28 SNPs in the 5 genes (4 SNPs in ANXA6, 6 SNPs in COL5A1, 3 SNPs in ENO1, 9 SNPs in MYOF and 6 SNPs in SCARA5) were significantly associated with BD-TT trait. Particularly, 1 SNP, rs868641 in the ANXA6 gene and 2 SNPs, rs7875570 and rs6537942 in the COL5A1 gene showed a highly significant association with BD-TT trait and their P-values satisfied the Bonferroni-corrected significance level (P<0.001852).

For osteoporosis case-control association analysis, the control subjects (n=1,711) and osteoporosis case subjects (n=651) were analyzed. The results of case-control association analysis between the 116 SNPs in the 7 genes and osteoporosis in the KARE women subjects are summarized in Table 4. Total 14 SNPs in 5 genes (1 SNP in ANXA6, 6 SNPs in COL5A1, 1 SNP in ENO1, 4 SNPs in MYOF and 2 SNPs in SCARA5) were significantly associated with osteoporosis.

Notably, 8 SNPs in the 5 genes (1 SNP in ANXA6, 1 SNP in COL5A1, 1 SNP in ENO1, 4 SNPs in MYOF and 1 SNP in SCARA5) were significantly associated with both the bone density and osteoporosis traits (Tables 3, 4). In all the 8 SNPs, their β-values in BD-RT and BD-TT traits were in the same direction and showed consistent trends with the odds ratios of osteoporosis. The location and basic LD of the analyzed SNPs in the ANXA6, COL5A1, MYOF and SCARA5 genes are shown in Supplementary Fig. 1. The SNPs that were significantly associated with bone density and/or osteoporosis in the KARE women subjects are also indicated. Interestingly, the 4 SNPs in the MYOF gene showing a significant association with both bone density and osteoporosis were located in the same LD block of the gene (Supplementary Fig. 1C).

Many approaches for identifying the genetic factors contributing to the pathogenesis of osteoporosis have been studied and have contributed to the detection of numerous genes for susceptibility to osteoporosis [11]. Among them, GWAS has been extensively executed for identifying the loci and genes that are significantly associated with bone density and osteoporosis. A major advantage of GWAS is that it offers the ranking for significance in multiple association signals across the genome. Since the statistical significance thresholds are very stringent due to the analysis of a large number of SNPs, many polymorphisms having a true association with osteoporosis but with a relatively small effect size can be missed [11]. This may lead to missing an opportunity to identify the novel osteoporosis susceptibility genes. In this study, to identify the novel genes more accurately as well as more effectively, we combined two methods, i.e whole genome expression profiling for screening of candidate genes and candidate gene association study.

In the series of experiments in the in vitro and in vivo osteoporosis models and eventually in human subjects, we identified 5 novel osteoporosis susceptibility genes, ANXA6, COL5A1, ENO1, MYOF, and SCARA5. The results from each step of the experiments in the cell line model, mouse model and humans are summarized in Table 5. In the screening step of the DEGs in the cell line model, 10 candidate genes were screened. During the validation and evaluation steps in the cell line and mouse models, respectively, the 5 genes showing false positive results were ruled out, and the 5 genes that passed all the steps of the experiment were finally selected.

The ANXA6 gene encodes Annexin A6 which belongs to a family of calcium-dependent membrane and phospholipid binding proteins and is involved in matrix vesicle calcification [33,34]. Annexin A6 binds to phospholipids in cellular membranes in a dynamic and reversible fashion and is implicated in membrane-related events along the exocytosis and endocytosis pathways [35]. Previous reports have documented that Annexin A6 is involved in cell proliferation, growth, and apoptosis [36,37]. Annexin A6 participates in the regulation of EGFR/Ras signaling pathway and cholesterol homeostasis [35,38,39]. The SNPs in the ANXA6 gene have been reported to be associated with osteonecrosis of the femoral head in the Korean population [40]. Although the phenotypes in this disease differ from those in osteoporosis, since the polymorphisms in this gene are associated with the phenotypes related with bone loss, it is very likely that the ANXA6 gene plays an important role in the pathogenesis of osteoporosis.

The COL5A1 gene encodes the alpha 1 chain of type V collagen, one of the low abundance fibrillar collagens. The COL1A1 gene encoding the alpha 1 chain of type I collagen has long been implicated in the pathogenesis of osteoporosis because type I collagen is the main protein in bone. Many studies on the association between polymorphisms in the COL1A1 gene and osteoporosis have been published [41,42], however, none of association results between the COL5A1 gene and osteoporosis have been reported. In the meantime, it has been reported that the COL5A1 gene is associated with various diseases including chronic Achilles tendinopathy [43], Achilles tendon injuries [44,45] and anterior cruciate ligament rupture in female participants [46]. In addition, mutations within the COL5A1 gene have been implicated in Ehlers-Danlos syndrome which is a multisystemic disorder that primarily affects the soft connective tissues [47]. These results suggested the COL5A1 gene can be a candidate high risk factor for osteoporosis and we found a significant association between the SNPs in the COL5A1 gene and bone density and osteoporosis phenotypes.

The ENO1 gene encodes a key glycolytic enzyme alpha-enolase that acts as a 2-phospho-D-glycerate hydrolase. It is also involved in various processes such as growth control, hypoxia tolerance and autoimmune responses [48,49]. Previous report showed that the ENO1 gene was significantly down-regulated in postmenopausal women compared with premenopausal women [50]. Our results also showed that the expression level of the ENO1 gene was decreased in both the cell line and mouse models for osteoporosis (Figs. 3 and 4). Based on the fact that deficiency of estrogen level due to menopause is closely related with an increase in osteoclast life span and a concomitant decrease in osteoblast life span and further osteoporosis [51,52], these results suggest that the ENO1 gene may be involved in bone metabolism.

The MYOF gene encodes Myoferlin which is a member of the ferlin family of proteins that promotes endomembrane fusion with the plasma membrane in muscle cells and endothelial cells [53]. Since Myoferlin was identified as a protein highly homologous to Dysferlin, the gene product of the limb girdle muscular dystrophy (LGMD) 2B locus, MYOF has been suggested as a candidate gene and potential modifier for muscular dystrophy [54]. The SCARA5 gene encodes scavenger receptor class A member 5 which is involved in the host defense properties of populations of human epithelial cells [55]. Scara5 is also known to be a ferritin receptor mediating non-transferrin iron delivery [56]. Class A scavenger receptor promotes osteoclast differentiation [57].

By in silico analysis of transcription factor binding of the significant SNPs using the TRANSFAC database (http://www.cbrc.jp/research/db/TFSEARCH.html), we found the binding sites of transcription factors in several SNPs. The sequence region in the minor allele of the SNP, rs12005720 in the COL5A1 gene contained the binding site for ETS1 (93.1 scoring point). ETS1 has been reported to be associated with systemic lupus erythematosus (SLE) [58]. SLE often accompanies osteoporosis. The sequence region in the minor allele of the SNP, rs787667 in the MYOF gene contained the binding site for GATA1 (93.0 scoring point). In addition, POU2F could bind to the sequence regions of the SNP, rs6660137 in the ENO1 gene (90.0 scoring point) and the SNP, rs2726941 in the SCARA5 gene (94.5 scoring point).

In conclusion, we identified 5 novel osteoporosis susceptibility genes through candidate gene selection in cells, evaluation of the DEGs in cells and mice, and association analysis in humans. There was a significant association between the SNPs in the 5 genes, ANXA6, COL5A1, ENO1, MYOF, and SCARA5, and bone density and/or osteoporosis. Nevertheless, further replication studies in other ethnic populations and functional studies on these 5 genes are needed. Notably, the SNPs in the ANXA6 and COL5A1 genes were highly significantly associated with bone density. These results indicate that these two genes may play an important role in regulation of bone metabolism and further the pathogenesis of osteoporosis.

The authors declare that they do not have any conflicts of interest.

This research was supported by a grant of the Korea Health Technology R&D Project through the Korea Health Industry Development Institution (KHIDI), funded by the Ministry of Health & Welfare, Republic of Korea (grant no. HR22C1734), the Korea Initiative for fostering University of Research and Innovation Program of the National Research Foundation (NRF) funded by the Korean government (MSIT) (No. NRF2021M3H1A104892211), and an intramural research grant from the Korea National Institute of Health (2020-NG-021).

Conception and design: HSJ, SYJ. Acquisition of data: BYK, EP, JK, CGL. Analysis and interpretation of data: DWK, BYK, EP. Drafting the article: BYK, DWK, EP, HSJ, SYJ. Final approval of the version to be published: HSJ, SYJ.