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 .
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  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 . 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
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
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% CO2.
ACP-based differential display RT-PCR was carried out using the predesigned arbitrary ACP (Seegene, Seoul, Korea) . First-strand cDNA synthesis was performed using the primer dT-ACT1 (5’-CTGTGAATGCTGCGACTACGATIIIII(T)18-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)15-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.
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
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.
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 CO2 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.
The subjects from the Korean Association Resource (KARE) study which were used in this study have been described in the previous report . 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 . 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.
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 . 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 . 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).
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
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
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
To identify the DEGs in the
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 (
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
To evaluate the accuracy of the identified DEGs
The mRNA expression levels of the identified DEGs were examined using the
We finally selected the 7 genes,
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
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
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
Notably, 8 SNPs in the 5 genes (1 SNP in
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 . 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 . 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 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,
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.
Basic characteristics of the women subjects in the KARE study cohort
|Characteristics||Quantitative analysis for bone density (n=3,570)||Case-control analysis for osteoporosis|
|Control (n=1,711)||Case (n=651)|
|Distal radius T-score||0.20±1.55||0.99±1.14||–1.26±1.64||<0.0001|
|Midshaft tibia T-score||–0.81±1.55||0.31±0.93||–3.11±0.99||<0.0001|
Values are presented as mean±standard deviation.
aSignificant differences in characteristics between the control and case were determined by the two-tailed Student’s
KARE, Korean Association REsource, BMI, body mass index.
List of the significantly DEGs in the dexamethasone-treated mouse MC3T3-E1 cells
|DEG no.||Expression level||Gene symbol||GenBank Accession no.||Gene definition|
|3||+||++||NM_015734||Collagen, type V, alpha 1|
|4||++||+++||NM_146007||Collagen, type VI, alpha 2|
|7||+||++||NM_008687||Nuclear factor I/B|
|8||+||+++||NM_028903||Scavenger receptor class A, member 5|
|9||+||+++||NM_009112||S100 calcium binding protein A10|
DEG, differentially expressed gene; Con, control; Dex, dexamethasone.
The results of association analysis between the SNPs in the 5 genes and bone density in the KARE women subjects
|Gene||SNP||Minor allele||MAF||Women (n=3,570)|
|BD-RT (T-score at distal radius)|
|BD-TT (T-score at midshaft tibia)|
SNP, single nucleotide polymorphism; KARE, Korean Association Resource; MAF, minor allele frequency; SEM, standard error; Add
The results of case-control association analysis between the SNPs in the 5 genes and osteoporosis in the KARE women subjects
|Gene||SNP||Minor allele||MAF||Women subjects (1,711 controls, 651 cases)|
|OR (95% CI)||Add
||OR (95% CI)||Dom
||OR (95% CI)||Rec
|rs17728338||T||0.095||1.39 (1.06-1.83)||0.018a||1.36 (1.01-1.82)||0.041a||3.53 (1.07-11.62)||0.038a|
|rs4335205||G||0.430||1.20 (1.02-1.43)||0.031a||1.40 (1.08-1.81)||0.010a||1.13 (0.83-1.53)||0.434|
|rs4319175||G||0.369||0.88 (0.74-1.05)||0.154||0.78 (0.61-0.99)||0.039a||1.02 (0.72-1.43)||0.931|
|rs9409917||G||0.065||1.16 (0.83-1.61)||0.377||1.23 (0.87-1.74)||0.243a||NA||0.998|
|rs12005720||G||0.127||0.85 (0.67-1.10)||0.214||0.90 (0.68-1.19)||0.448||0.39 (0.15-0.99)||0.048a|
|rs11103535||T||0.144||1.21 (0.96-1.53)||0.106||1.16 (0.89-1.52)||0.264||2.17 (1.04-4.51)||0.039a|
|rs10858284||A||0.148||0.73 (0.58-0.93)||0.012a||0.72 (0.55-0.95)||0.022a||0.49 (0.22-1.11)||0.087|
|rs11121247||T||0.067||1.41 (1.02-1.97)||0.041a||1.44 (1.01-2.04)||0.042a||1.68 (0.28-10.07)||0.573|
|rs2797581||C||0.108||0.81 (0.61-1.06)||0.121||0.86 (0.63-1.15)||0.307||0.20 (0.05-0.80)||0.022a|
|rs787665||G||0.109||0.83 (0.63-1.09)||0.186||0.87 (0.65-1.18)||0.374||0.25 (0.06-0.99)||0.048a|
|rs787667||T||0.109||0.84 (0.64-1.10)||0.202||0.89 (0.67-1.20)||0.461||0.21 (0.05-0.81)||0.024a|
|rs787668||A||0.109||0.84 (0.64-1.10)||0.192||0.89 (0.66-1.20)||0.442||0.21 (0.05-0.81)||0.024a|
|rs2726959||T||0.298||1.19 (0.99-1.43)||0.064||1.11 (0.88-1.41)||0.390||1.79 (1.18-2.71)||6.1E-03a|
|rs7002829||A||0.122||1.18 (0.92-1.52)||0.185||1.09 (0.82-1.45)||0.534||2.93 (1.32-6.50)||8.2E-03a|
SNP, single nucleotide polymorphism; KARE, Korean Association Resource; MAF, minor allele frequency; OR, odds ratio; CI, confidence interval; Add
The summary of the results from each step of the experiments in the cell line, mouse model, and humans
|Mouse Gene symbol||DD||qRT-PCR||Humans Gene symbol||Number of tested SNPs||Number of associated SNPs (lowest
|Cell line||Cell model||Mouse model||Human (woman)|
|+||++||1.0||2.8±0.1a||1.0||1.5±0.2c||18||2 (0.044)||4 (1.6×10-3)e||1 (0.018)|
|+||++||1.0||2.1±0.5b||1.0||1.9±0.5c||29||3 (2.5×10-3)||6 (5.3×10-4)e||6 (0.010)|
|+++||++||1.0||0.8±0.1b||1.0||0.7±0.1c||3||0||3 (6.3×10-3)||1 (0.041)|
|+||+++||1.0||3.2±0.2a||1.0||2.3±0.4d||29||6 (0.038)||9 (0.017)||4 (0.022)|
|+||+++||1.0||3.0±0.2a||1.0||2.2±0.3d||25||1 (0.044)||6 (6.4×10-3)||2 (6.1×10-3)|
Values are presented as mean±standard deviation.
DD, differential display; qRT-PCR, quantitative reverse transcriptase polymerase chain reaction; SNP, single nucleotide polymorphism; Con, control; Dex, dexamethasone; OVX, ovariectomy; BD-RT, bone density estimated by T-score at the distal radius; BD-TT, bone density estimated by T-score at the midshaft tibia; - , not tested.