Regulatory SNP of RREB1 is Associated With Bone Mineral Density in Chinese Postmenopausal Osteoporosis Patients

Postmenopausal osteoporosis (PMO) is the most common bone disorder in elderly Chinese women. Although genetic factors have been shown to have a pivotal role in PMO, studies on genetic loci associated with PMO in Chinese individuals are still lacking. We aimed to identify SNPs that contribute to PMO in Chinese individuals by conducting a genome-wide association study (GWAS). Bone mineral density (BMD) of postmenopausal Chinese women was assessed. Participants with T-score < −2.5 standard deviations (n = 341) were recruited and divided into a discovery group (n = 150) and a replication group (n = 191). GWAS was performed, with T-score as the quantitative trait, using linear regression. Our results revealed that an SNP cluster upstream of RREB1 showed a trend of association with BMD in Chinese PMO patients. The leading SNP of the cluster was rs475011 (p_combined = 1.15 × 10⁻⁶, beta = 0.51), which is a splicing quantitative trait locus (sQTL) of RREB1. This association was further supported by data from the UK Biobank (UKBB; p = 9.56 × 10⁻¹²). The high BMD-associated allele G of rs475011 is related to a high intron excision ratio. This SNP may increase BMD by upregulating mature RREB1 mRNA, based on data from the Genotype-Tissue Expression (GTEx) database. We identified BMD-associated SNPs that regulate RREB1 in Chinese PMO patients. Future functional experiments are needed to further link rs475011, RREB1, and PMO in Chinese individuals.

Introduction

Postmenopausal osteoporosis (PMO) is the most common bone disorder in elderly women worldwide (Black and Rosen, 2016; Eastell et al., 2016). In China, its estimated prevalence is 34.65% (Chen et al., 2016). It is characterized by altered bone microstructure and decreased bone mineral density (BMD), and it is the leading cause of nonstress fractures in the elderly (Zhuang et al., 2020; Porter and Varacallo, 2021). The risk factors for PMO include decreased postmenopausal estrogen levels, smoking, decreased physical activity, low body mass index (BMI), and genetic susceptibility (Eastell et al., 2016; Bijelic et al., 2017). The heritability of BMD is 0.6–0.8 (Peacock et al., 2002) and the heritability of osteoporotic fracture is 0.5–0.7 (Deng et al., 2002). These findings indicate that genetics play a pivotal role in PMO.

Genome-wide association studies (GWASs) are useful for identifying genetic susceptibility factors underlying diseases and phenotypes (Visscher et al., 2017). Early single-sample GWASs identified associations of SNPs in RANKL, OPG, ESR1, and TNFRSF11B with BMD in individuals of European descent (Richards et al., 2008; Styrkarsdottir et al., 2008). By integrating data from multiple studies, the GEnetic Factors for OSteoporosis Consortium (GEFOS) conducted a GWAS meta-analysis of >580,000 individuals. They identified 15 loci associated with BMD and fractures, including SNPs in SOST, WNT16, and ESR1 (Trajanoska et al., 2018). The susceptibility SNPs identified in the GEFOS study explain 6–10% of the phenotypic variance in BMD (Medina-Gomez et al., 2018; Yang et al., 2020). Furthermore, data on ∼420,000 individuals from the UK Biobank (UKBB) were used to link 518 genetic loci with BMD, explaining 20% of the phenotypic variance (Morris et al., 2019). The susceptibility genes are related to osteoblast and osteoclast differentiation, which further supports the importance of genetic factors in PMO (Yang et al., 2020).

The majority of the genetic susceptibility loci regarding BMD and osteoporosis have been obtained from individuals of European descent. There is a lack of GWAS data from Chinese populations. Kung et al. reported an association between rs2273061 of JAG1 with BMD among 800 Han Chinese women using a case-control GWAS (Kung et al., 2010). Guo et al. linked rs13182402 of ALDH7A with osteoporotic fractures in 700 elderly Han Chinese subjects (Guo et al., 2010). Tan et al. identified SNPs upstream of ATP6V1G1 that were associated with both age at menarche and BMD in 800 Chinese participants, using bivariate GWAS (Tan et al., 2015). The fact that the susceptibility loci identified in Chinese subjects exhibit little overlap with those identified in European subjects indicates differences in genetic risk factors for PMO between Chinese and European populations. More research on the genetic contributors to PMO in Chinese individuals is crucial.

We aimed to conduct a GWAS using the T-score (which is a critical measurement of BMD and osteoporosis) as the quantitative trait in Han Chinese PMO patients to identify genetic factors influencing the severity of the disease. We found that SNPs upstream of RREB1 were significantly associated with BMD. This finding provides new information on the genetic basis of PMO in Chinese individuals.

Materials and Methods

Ethics Statement

This study was conducted according to the principles of the Declaration of Helsinki and approved by the ethics committees of Beijing Jishuitan Hospital. All of the patients provided written informed consent, and the ethics committees approved the consent procedure.

Study Subjects

A total of 341 unrelated Han Chinese postmenopausal women aged >54 years were recruited from Beijing Jishuitan Hospital. They were diagnosed with PMO using World Health Organization criteria (WHO Scientific Group on the Prevention and Management of Osteoporosis, 2003). In brief, the lumbar vertebra BMD of each participant was measured by dual-energy X-ray absorptiometry (DXA; Prodigy; GE Healthcare, Boston, MA) or quantitative computed tomography (qCT; uCT710; United-imaging, Shanghai, China). T-score was obtained by comparison with the BMD of healthy young women. Patients with a T-score less than −2.5 SD were diagnosed with PMO. Data on PMO-related characteristics and environmental exposures (such as BMI, age of menopause, age of menarche, hypertension, type 2 diabetes, and smoking and drinking habits) were also collected. No participants had other diseases that may affect bone mass or were taking hormone replacement therapy.

Pearson’s correlation test was adopted to test the relationship between pairs of all collected phenotypes and environmental exposures in our samples. To calculate the correlation coefficient between binary traits and quantitative traits, binary traits or exposures such as hypertension, type 2 diabetes (T2D), and smoking and drinking habits were coded as 1 for affected or exposed and 0 for unaffected or unexposed. The relationship between BMD and hypertension or T2D was further tested by Student’s t-test. The statistical analysis was performed on R 4.0 software.

DNA Extraction and Genotyping

A whole-blood sample was obtained from each participant using an ethylenediaminetetraacetic acid (EDTA) blood collection tube (BD, Franklin Lakes, NJ). The DNA was extracted using a QIAamp DNA Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer’s protocol. The blood samples and extracted DNA were stored at -80 °C prior to genotyping.

In the discovery stage of GWAS, whole-genome genotyping was performed using an Illumina Infinium Asian Screening Array (Illumina, San Diego, CA) according to the manufacturer’s protocol. Data management and experimental quality control were conducted using GenomeStudio 2.0 software (Illumina, San Diego, CA).

In the replication stage of GWAS, genotyping of candidate SNPs (selected based on the results of the discovery stage) was conducted using a MassARRAY system (Agena Bioscience, San Diego, CA) according to the manufacturer’s protocol. The genotyping experiments were designed using online Agena tools (https://agenacx.com). Data management and experimental quality control were conducted using Typer 4.0 software (Agena Bioscience, San Diego, CA).

Discovery Stage

In the discovery stage of GWAS, whole-genome genotyping data were analyzed using a PLINK 1.07 software suite (Purcell et al., 2007). First, stringent filtering was conducted; individuals with genotyping call rate <95%, inbreeding coefficient >0.025, or pairwise identity-by-descent (IBD) > 0.05 were excluded. SNP loci with call rate <90%, minor allele frequency (MAF) < 5%, or significant deviation from Hardy–Weinberg equilibrium (HWE, p < 0.01) were also removed. Single-allele association analysis was performed by comparing SNP allele counts with T-score using linear regression, with age and BMI as covariates. Since T-score is a measurement of standard normalized BMD, we used this value as a quantitative trait for association directly (Lewiecki and Lane, 2008). The genome-wide trend of association threshold was set at p < 10⁻⁵.

Replication Stage

In the replication stage, the leading SNPs of each of the BMD-associated SNP clusters and 1 additional SNP in the SNP cluster upstream of RREB1 (which was the second leading SNP in this cluster) were selected for validation. Linear regression was used to verify the associations between candidate SNP allele counts and T-score. Age and BMI were used as the covariates. The significance threshold was set at p < 0.05.

To replicate our result on public data, summary statistics of GWAS results of heel T-score were downloaded from the UK Biobank (UKBB) pan-ancestry genetic analysis portal (https://pan.ukbb.broadinstitute.org). The selected phenotype was “Heel bone mineral density (BMD) T-score, automated”.

Meta-Analysis

The final association results of our samples were obtained by combining the results from the discovery and replication stages using the meta-analysis module of PLINK 1.07 (Purcell et al., 2007), using a random- or fixed-effects model for loci with Cochrane’s Q statistic <0.05 or >0.05, respectively. The meta-analysis between our data and UKBB data was also performed by PLINK 1.07.

Results

Participant Grouping and Characteristics

In total, 341 PMO patients were recruited. All the patients were Han Chinese and diagnosed with PMO using DXA or qCT, based on T-score less than -2.5 SD. They were divided into two groups: 150 were used in the discovery stage of GWAS and 191 in the replication stage. The majority of samples in the discovery stage were collected during 2019–2020, and samples in the replication stage were collected during 2020–2021. Data on age, age of menopause, age of menarche, BMI, BMD, and T-score were collected (Table 1). Data on whether the participants had hypertension or T2D were also collected (Table 1). By analyzing the relationships among these variables in our samples (Figure 1 and Supplementary Table S1), we found that age (p = 1.25 × 10⁻¹¹, Pearson’s r = −0.44) and years since menopause (YSM, p = 3.83 × 10⁻¹⁰, Pearson’s r = −0.41) were significantly negatively correlated with BMD, and BMI was significantly positively correlated with BMD (p = 4.85 × 10⁻⁵, Pearson’s r = 0.27), which is consistent with previous reports (Ravn et al., 1999; Warming et al., 2002). However, hypertension (p = 0.73, t-test) and T2D (p = 0.67, t-test) did not affect BMD. As only ∼70% of the participants had YSM data available (Table 1) and there was a high correlation between age and YSM (p < 10⁻¹⁶, Pearson’s r = 0.87), we finally used age and BMI as covariates in the GWAS.

TABLE 1

TABLE 1. PMO-related characteristics of the participants.

FIGURE 1

FIGURE 1. Correlations among different phenotypes in PMO patients. (A) Heatmap of correlation coefficients between pairs of phenotypes. Each box represents the Pearson’s r² between the column and row phenotypes. The shade of red indicates the value of r². (B) Heatmap of correlation p values between pairs of phenotypes. Each box represents the log-transformed p value of the correlation between the column and row phenotypes. BMD: bone mineral density; YSM: years since menopause; BMI: body mass index; T2D; type 2 diabetes.

The Discovery Stage of GWAS Identified Loci Associated With T-Score in PMO Patients

To assess the contribution of genetic polymorphisms to BMD in PMO patients, we conducted a GWAS using T-score as the quantitative trait. In the discovery stage, genome-wide genotyping was performed on 150 patients (Table 1) using an Illumina Infinium Asian Screening Array. In total, 705,335 autosomal SNPs were genotyped. As part of the stringent quality control process, SNPs with missing rate >10%, MAF <5%, or significant deviation from HWE (p < 0.01) were removed from further analysis, leaving 311,319 loci. We also performed patient-level quality control, which involved identifying individuals with call rate <95%, inbreeding coefficient >0.025, or pairwise IBD >0.05. No genotyped patients met these thresholds, and so all of them were included in the GWAS.

We carried out linear regression with age and BMI as covariates to assess the associations between genotypes and T-score. As shown in Figure 2A, the quantile-quantile plot of the genome-wide p value exhibited significant deviation from the null distribution at the tail of the distribution. The genomic inflation factor was calculated as 1.014. This indicates little population stratification and the existence of genetic factors related to BMD in the PMO patients. Using a lenient threshold of p < 10⁻⁵ as a trend of genome-wide association, we identified 6 SNP clusters associated with BMD in the PMO patients (Figure 2B; Table 2 and Supplementary Table S2).

FIGURE 2

FIGURE 2. GWAS of T-score in PMO patients. (A) Quantile-quantile plot of the GWAS p value. X-axis represents the expected p value. Y-axis represents the observed p value based on linear regression. Red line shows the theoretical distribution under the null hypothesis of no association in the data. (B) Manhattan plot of GWAS results. X-axis represents physical position of each tested SNP. Y-axis indicates log-transformed p value. Blue line shows the trend of association threshold of 10⁻⁵. The names of leading SNPs in each associated SNP cluster are shown in the panel.

TABLE 2

TABLE 2. SNPs showed trend of association with BMD in PMO patients.

The most significantly associated SNP cluster was located at 6p25.1 and its strongest signal came from rs473437 (p = 1.4 × 10⁻⁷). Its minor allele T led to an increase in the T-score of 0.84 (95% CI: 0.54–1.13) per allele in the discovery group. The SNP cluster was located 67 kb upstream of RREB1. RREB1 is a zinc finger transcription factor that binds to RAS-responsive elements in gene promoters. It can upregulate calcitonin (Thiagalingam et al., 1996). SNPs in the RREB1 gene have been reported to be associated with heel BMD in individuals of European descent (Kemp et al., 2017; Kim, 2018; Morris et al., 2019).

The Replication Stage of GWAS Verified Associations Between SNPs Upstream of RREB1 and BMD in PMO Patients

The leading SNPs in each of the 6 BMD-associated SNP clusters and 1 additional SNP in the SNP cluster upstream of RREB1 (which was the second leading SNP in this cluster; the leading SNP of this cluster had a missing rate of 8.6%, Supplementary Table S2) were selected for validation in the additional 191 PMO patients (Table 1). Using linear regression with age and BMI as covariates, three out of the seven tested SNPs were significantly associated with BMD. These comprised 2 upstream of RREB1 and 1 from NELL1 on chromosome 11 (Table 2, p < 0.05). After combining the results from the discovery and replication stages in a meta-analysis, only the two SNPs upstream of RREB1 (rs473437 and rs475011, Table 2) were associated with BMD in PMO patients. These two SNPs exhibit perfect linkage disequilibrium (LD; r² = 1) in Han Chinese individuals, according to the data from the 1000 Genome Project. And rs475911 showed the most significant association signal with minor allele G related to increased T-score (p = 1.15 × 10⁻⁶, beta = 0.51, Table 2).

As a relatively small sample size was used in our study, we used a publicly available data set to further validate our results. According to data from The Musculoskeletal Knowledge Portal, the A allele of rs475011 is significantly associated with decreased estimated BMD among ∼570,000 samples (p = 5.69 × 10⁻¹², beta = ⁻0.02) (Kiel et al., 2020). This result is consistent with our finding that the G allele of rs475011 is associated with an increased T-score (Table 2). We also analyzed the relationship between rs475011 and heel BMD T-score using summary statistics data from UKBB. As shown in Figure 3, the G allele of rs475011 was significantly associated with increased T-score in European samples (p = 6.05 × 10⁻¹¹, beta = 0.030). Except for samples of Admixed American and Middle Eastern descents, the G allele was related to increased T-score in the remaining populations (Figure 3). By combining all data sets assessed (discovery, replication, and UKBB), the association between rs475011 and T-score reached p = 9.56 × 10⁻¹², beta = 0.031 (Figure 3). Meta-analysis of our samples with East Asian descents of UKBB also indicated the contribution of the rs475011 G allele to T-score increasing (p = 0.0035, beta = 0.153).

FIGURE 3

FIGURE 3. Forest plot of meta-analysis result of rs475011. By combining result from our samples with UKBB summary statistics, a significantly association between rs475011 and T-score was obtained (p = 9.56 × 10⁻¹², beta = 0.0308). Sample size of each ancestry and analysis group is listed on the left of the figure and the association results are list on the right. X-axis represents the beta value and Y-axis represents the six ancestry groups from UKBB and the discovery and replication group of our study. UKBB ancestry groups are rank by the sample size. The position of the square is the point estimation of the beta value from each group. The diamond represents the combined results. The width of the diamond indicate the 95% confidence interval. UKBB: UK biobank; EUR: European ancestry; CSA: Central/South Asian ancestry; AFR: African ancestry; EAS: East Asian ancestry; AMR: Admixed American ancestry; MID: Middle Eastern ancestry.

BMD-Associated SNPs May Affect RREB1 Expression

As the BMD-associated SNP cluster was located 67 kb upstream of RREB1, we further checked the regulatory effects of these SNPs by querying the Genotype-Tissue Expression (GTEx) database (Consortium, 2020). As shown in Figure 4A, this region is known to act as a splicing quantitative trait locus (sQTL) of RREB1 in the skin (Garrido-Martin et al., 2021). The high BMD-associated allele G of rs475011, which was the leading SNP in this cluster in the meta-analysis, was significantly associated with an increased intron excision ratio of RREB1 (Figure 4B, p = 2.1 × 10⁻⁶). Moreover, this SNP is located in the binding site of two transcription factors, RFX3 and FOXB1, according to chromatin immunoprecipitation (ChIP)-seq data and prediction, respectively, based on information in the RegulomeDB database (Jolma et al., 2010; Boyle et al., 2012; Jolma et al., 2013). This indicates that rs475011 affects BMD in PMO patients by regulating RREB1 expression.

FIGURE 4

FIGURE 4. rs475011 acts as a splicing quantitative trait locus (sQTL) of RREB1. (A) Upper layer: physical position (based on human reference genome hg38) of the SNP cluster containing rs475011. Middle (blue) layer: major alleles of the labeled SNPs are negatively associated with the intron excision ratio of RREB1. Lower triangle plot: linkage disequilibrium (D′) of SNPs in this cluster. Data were obtained from GTEx database. (B) Allele G of rs475011 is associated with an increased intron excision ratio of RREB1. Y-axis represents the normalized intron excision ratio of RREB1. The X-axis represents different genotypes of rs475011. Data were obtained from GTEx database.

Discussion

In this study, we performed a GWAS of BMD in PMO patients and identified SNPs upstream of RREB1 that significantly influenced BMD. GWAS is a powerful tool for identifying genetic factors related to various phenotypes, and it has been used to identify osteoporosis susceptibility genes in many studies (Yang et al., 2020; Zhu et al., 2021). A GWAS of half a million European participants linked BMD with >500 genetic loci (Morris et al., 2019), but this kind of research on Chinese individuals is still lacking. Although we only included 341 participants, the two-step GWAS design, the replication on UKBB summary statistics, and significant results make our study valuable for understanding the genetic basis of PMO in Chinese individuals.

Sample size is a critical parameter that influences the statistical power of GWAS (Spencer et al., 2009). Based on the GWAS power calculator described by Visscher et al. (Visscher et al., 2017), along with the beta value (0.51) and MAF of rs475011 (0.108) and the sample size of 341, the estimated detection power of our study was 0.98. This further demonstrates the reliability of our results. As genetic research on Chinese PMO is limited, further research involving large samples is highly important for finding genetic loci with small effect sizes (Spencer et al., 2009; Yang et al., 2020).

PMO is a late-onset disease with high prevalence. As it is hard to define appropriate controls (Warming et al., 2002; Chen et al., 2016; Bijelic et al., 2017), we conducted a quantitative trait GWAS to avoid misclassification of controls. The genetic correlation coefficients among various BMD measurements assessed using different platforms (e. g., DXA, qCT, or ultrasound) varies from 0.505 to 0.917 (Karasik et al., 2017). Hence, to reduce the influence of heterogeneity caused by different measurement platforms, we used the T-score (based on standard normalized BMD data measured using different platforms for different subjects), which is a relative value of BMD, as the quantitative trait in the GWAS. Moreover, BMD measured at different skeleton sites may also introduce heterogeneity to the data (Yang et al., 2020). Therefore, we obtained lumbar vertebra BMD values for all participants to avoid this heterogeneity. These measures decreased the confounding factors in our study.

The most significant high BMD-associated SNP in our study, rs475011, was located upstream of RREB1. This region is known to act as an sQTL of RREB1, according to data from the GTEx database (Figure 3). The high BMD-associated allele G of rs475011 is also linked to an increased intron excision ratio of RREB1, based on GTEx data. Higher intron excision indicates upregulation of mature RREB1 mRNA, suggesting that RREB1 upregulation may help to prevent bone mineral loss in elderly women. Our results are consistent with previous GWASs of individuals of European descent (Kemp et al., 2017; Kim, 2018), which reported that the minor allele A of rs525678, upstream of RREB1, was associated with increased BMD. According to data from the 1000 Genomes Project and GTEx database, in individuals of European descent, rs525678 was in high LD (D’ = 0.998) with rs2842895, which is an expression quantitative trait locus (eQTL) of RREB1 (Genomes Project et al., 2015). The minor allele C of rs2842895 is linked to RREB1 upregulation (Consortium, 2020). Moreover, by analyzing data from UKBB, we found the rs475011 is also associated with BMD in individuals of European descent. There is moderate LD between rs475011 and rs2842895 (D’ = 0.719), which indicates that both of the loci contribute to BMD in individuals of European descent. In contrast, rs2842895 is a rare SNP in Han Chinese individuals (MAF = 0.005), which indicates that rs475011 is the major SNP that affects BMD in Han Chinese individuals. Taking the results together, high RREB1 expression may be a protective factor regarding the development of PMO.

RREB1 is a transcription factor that binds to RAS-responsive elements in gene promoters. The calcitonin gene promoter has a binding site for RREB1, and RREB1 binding upregulates calcitonin (Thiagalingam et al., 1996). Calcitonin is widely used to treat osteoporosis in menopausal women (Reginster, 1993). It may play a protective role regarding BMD by inhibiting the bone resorption function of osteoclasts (Keller et al., 2014). The high BMD-associated alleles identified in both our research (G allele of rs475011) and previous research (A allele of rs525678) may upregulate calcitonin via RREB1 upregulation (Kemp et al., 2017; Morris et al., 2019). In turn, upregulated calcitonin may lower osteoclast activity and prevent bone mineral loss in postmenopausal women.

To overcome the under statistical power caused by the small sample size of our study, we obtained genome-wide association summary statistics of rs475011 with T-score from UKBB pan-ancestry genetic analysis portal. Although a meta-analysis of all sample groups showed a significant result, among UKBB populations, a significant association was only observed in European descents (p = 6.05 × 10⁻¹¹, Figure 3). This might be caused by the relatively small sample size of the non-European populations, and the results in 242,926 European samples may be more robust and trustable. Moreover, UKBB samples were collected from individuals under various health conditions, of which the PMO only accounted for a small proportion (Sudlow et al., 2015). In contrast, our samples were all PMO patients with T-score < −2.5, and rs475011 may contribute to the protective against PMO, which can also be inferred by the role of RREB1 on regulation calcitonin. The protective role can explain the smaller beta value in the UKBB data in comparison with our data (Figure 3). Additionally, since East Asian descents of UKBB may include limited individuals of PMO, this can explain the absence of the association in them. This speculation was further supported by the fact that the protective effect of the rs475011 G allele was larger in our discovery group than the replication group (beta = 0.69 vs. 0.32, Table 2) since the discovery group had a smaller average T-score or severe PMO condition (p = 0.016, t-test, Table 1 and Figure 3).

The results of our study may help to understand PMO in Chinese individuals; however, several limitations should be considered. First, the sample size was relatively small. Although the estimated detection power of our study for loci with large effect sizes (beta >0.5) is 0.98, the sample size is insufficient for identifying loci with small effect sizes, which may act as major genetic contributors of common diseases (Boyle et al., 2017). Second, the information on the relationships of the BMD-associated SNPs with RREB1 was gleaned from databases. Although our study and previous publications (Kemp et al., 2017; Kim, 2018) obtained concordant results, there is a lack of direct functional evidence linking the reported SNPs with the expression of RREB1 and calcitonin.

In conclusion, we conducted a GWAS of BMD based on Chinese PMO patients and found that an RREB1 allele that increases the intron excision ratio is associated with increased BMD. The results are consistent with previous research on individuals of European descent. However, this is a preliminary study, and future association studies should include large numbers of Chinese PMO patients. eQTL, sQTL, and functional analyses of rs475011, RREB1, and calcitonin should be carried out.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Ethics Statement

The studies involving human participants were reviewed and approved by the ethics committees of Beijing Jishuitan Hospital. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

SF and WC designed the research and wrote the manuscript. SF, YY, and WC performed the data analysis. HW and JA collected the data. XS conducted the genotyping experiment. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by a grant from Beijing JST Research Funding (ZR-201922, to SF), a grant from the National Natural Science Foundation of China (31701081, to WC) and a grant from Beihang Biomedical Engineering Excellence and Innovation Program (JKF-YG-21-B008, to WC).

Conflict of Interests

Author XS is employed by Beijing GuardianHealth Technology Co., Ltd.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors gratefully thank all the participants for volunteering to participate in the study. The authors also thank Charlesworth Author Services for editing and proofreading the manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.756957/full#supplementary-material

References

Bijelic, R., Milicevic, S., and Balaban, J. (2017). Risk Factors for Osteoporosis in Postmenopausal Women. Med. Arch. 71 (1), 25–28. doi:10.5455/medarh.2017.71.25-28

PubMed Abstract | CrossRef Full Text | Google Scholar

Black, D. M., and Rosen, C. J. (2016). Postmenopausal Osteoporosis. N. Engl. J. Med. 374 (3), 254–262. doi:10.1056/NEJMcp1513724

CrossRef Full Text | Google Scholar

Boyle, A. P., Hong, E. L., Hariharan, M., Cheng, Y., Schaub, M. A., Kasowski, M., et al. (2012). Annotation of Functional Variation in Personal Genomes Using RegulomeDB. Genome Res. 22 (9), 1790–1797. doi:10.1101/gr.137323.112

PubMed Abstract | CrossRef Full Text | Google Scholar

Boyle, E. A., Li, Y. I., and Pritchard, J. K. (2017). An Expanded View of Complex Traits: From Polygenic to Omnigenic. Cell 169 (7), 1177–1186. doi:10.1016/j.cell.2017.05.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, P., Li, Z., and Hu, Y. (2016). Prevalence of Osteoporosis in China: a Meta-Analysis and Systematic Review. BMC Public Health 16 (1), 1039. doi:10.1186/s12889-016-3712-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Consortium, G. T. (2020). The GTEx Consortium Atlas of Genetic Regulatory Effects across Human Tissues. Science 369 (6509), 1318–1330. doi:10.1126/science.aaz1776

PubMed Abstract | CrossRef Full Text | Google Scholar

Deng, H.-W., Mahaney, M. C., Williams, J. T., Li, J., Conway, T., Davies, K. M., et al. (2002). Relevance of the Genes for Bone Mass Variation to Susceptibility to Osteoporotic Fractures and its Implications to Gene Search for Complex Human Diseases. Genet. Epidemiol. 22 (1), 12–25. doi:10.1002/gepi.1040

PubMed Abstract | CrossRef Full Text | Google Scholar

Eastell, R., O'Neill, T. W., Hofbauer, L. C., Langdahl, B., Reid, I. R., Gold, D. T., et al. (2016). Postmenopausal Osteoporosis. Nat. Rev. Dis. Primers 2, 16069. doi:10.1038/nrdp.2016.69

PubMed Abstract | CrossRef Full Text | Google Scholar

Garrido-Martín, D., Borsari, B., Calvo, M., Reverter, F., and Guigó, R. (2021). Identification and Analysis of Splicing Quantitative Trait Loci across Multiple Tissues in the Human Genome. Nat. Commun. 12 (1), 727. doi:10.1038/s41467-020-20578-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Genomes Project, C., Auton, A., Brooks, L. D., Durbin, R. M., Garrison, E. P., Kang, H. M., et al. (2015). A Global Reference for Human Genetic Variation. Nature 526 (7571), 68–74. doi:10.1038/nature15393

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, Y., Tan, L.-J., Lei, S.-F., Yang, T.-L., Chen, X.-D., Zhang, F., et al. (2010). Genome-wide Association Study Identifies ALDH7A1 as a Novel Susceptibility Gene for Osteoporosis. Plos Genet. 6 (1), e1000806. doi:10.1371/journal.pgen.1000806

PubMed Abstract | CrossRef Full Text | Google Scholar

Jolma, A., Kivioja, T., Toivonen, J., Cheng, L., Wei, G., Enge, M., et al. (2010). Multiplexed Massively Parallel SELEX for Characterization of Human Transcription Factor Binding Specificities. Genome Res. 20 (6), 861–873. doi:10.1101/gr.100552.109

PubMed Abstract | CrossRef Full Text | Google Scholar

Jolma, A., Yan, J., Whitington, T., Toivonen, J., Nitta, K. R., Rastas, P., et al. (2013). DNA-binding Specificities of Human Transcription Factors. Cell 152 (1-2), 327–339. doi:10.1016/j.cell.2012.12.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Karasik, D., Demissie, S., Zhou, Y., Lu, D., Broe, K. E., Bouxsein, M. L., et al. (2017). Heritability and Genetic Correlations for Bone Microarchitecture: The Framingham Study Families. J. Bone Miner Res. 32 (1), 106–114. doi:10.1002/jbmr.2915

PubMed Abstract | CrossRef Full Text | Google Scholar

Keller, J., Catala-Lehnen, P., Huebner, A. K., Jeschke, A., Heckt, T., Lueth, A., et al. (2014). Calcitonin Controls Bone Formation by Inhibiting the Release of Sphingosine 1-phosphate from Osteoclasts. Nat. Commun. 5, 5215. doi:10.1038/ncomms6215

PubMed Abstract | CrossRef Full Text | Google Scholar

Kemp, J. P., Morris, J. A., Medina-Gomez, C., Forgetta, V., Warrington, N. M., Youlten, S. E., et al. (2017). Identification of 153 New Loci Associated with Heel Bone mineral Density and Functional Involvement of GPC6 in Osteoporosis. Nat. Genet. 49 (10), 1468–1475. doi:10.1038/ng.3949

PubMed Abstract | CrossRef Full Text | Google Scholar

Kiel, D. P., Kemp, J. P., Rivadeneira, F., Westendorf, J. J., Karasik, D., Duncan, E. L., et al. (2020). The Musculoskeletal Knowledge Portal: Making Omics Data Useful to the Broader Scientific Community. J. Bone Miner Res. 35 (9), 1626–1633. doi:10.1002/jbmr.4147

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, S. K. (2018). Identification of 613 New Loci Associated with Heel Bone mineral Density and a Polygenic Risk Score for Bone mineral Density, Osteoporosis and Fracture. PLoS One 13 (7), e0200785. doi:10.1371/journal.pone.0200785

PubMed Abstract | CrossRef Full Text | Google Scholar

Kung, A. W. C., Xiao, S.-M., Cherny, S., Li, G. H. Y., Gao, Y., Tso, G., et al. (2010). Association of JAG1 with Bone mineral Density and Osteoporotic Fractures: a Genome-wide Association Study and Follow-Up Replication Studies. Am. J. Hum. Genet. 86 (2), 229–239. doi:10.1016/j.ajhg.2009.12.014

CrossRef Full Text | Google Scholar

Lewiecki, E. M., and Lane, N. E. (2008). Common Mistakes in the Clinical Use of Bone mineral Density Testing. Nat. Rev. Rheumatol. 4 (12), 667–674. doi:10.1038/ncprheum0928

PubMed Abstract | CrossRef Full Text | Google Scholar

Medina-Gomez, C., Kemp, J. P., Trajanoska, K., Luan, J., Chesi, A., Ahluwalia, T. S., et al. (2018). Life-Course Genome-wide Association Study Meta-Analysis of Total Body BMD and Assessment of Age-specific Effects. Am. J. Hum. Genet. 102 (1), 88–102. doi:10.1016/j.ajhg.2017.12.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Morris, J. A., Kemp, J. P., Youlten, S. E., Laurent, L., Logan, J. G., Chai, R. C., et al. (2019). An Atlas of Genetic Influences on Osteoporosis in Humans and Mice. Nat. Genet. 51 (2), 258–266. doi:10.1038/s41588-018-0302-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Peacock, M., Turner, C. H., Econs, M. J., and Foroud, T. (2002). Genetics of Osteoporosis. Endocr. Rev. 23 (3), 303–326. doi:10.1210/edrv.23.3.0464

PubMed Abstract | CrossRef Full Text | Google Scholar

Porter, J. L., and Varacallo, M. (2021). Osteoporosis. (Treasure Island (FL)): StatPearls.

Google Scholar

Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M. A. R., Bender, D., et al. (2007). PLINK: a Tool Set for Whole-Genome Association and Population-Based Linkage Analyses. Am. J. Hum. Genet. 81 (3), 559–575. doi:10.1086/519795

CrossRef Full Text | Google Scholar

Ravn, P., Cizza, G., Bjarnason, N. H., Thompson, D., Daley, M., Wasnich, R. D., et al. (1999). Low Body Mass Index Is an Important Risk Factor for Low Bone Mass and Increased Bone Loss in Early Postmenopausal Women. J. Bone Miner Res. 14 (9), 1622–1627. doi:10.1359/jbmr.1999.14.9.1622

CrossRef Full Text | Google Scholar

Reginster, J. Y. (1993). Calcitonin for Prevention and Treatment of Osteoporosis. Am. J. Med. 95 (5A), 44S–47S. doi:10.1016/0002-9343(93)90382-y

CrossRef Full Text | Google Scholar

Richards, J., Rivadeneira, F., Inouye, M., Pastinen, T., Soranzo, N., Wilson, S., et al. (2008). Bone mineral Density, Osteoporosis, and Osteoporotic Fractures: a Genome-wide Association Study. The Lancet 371 (9623), 1505–1512. doi:10.1016/S0140-6736(08)60599-1

CrossRef Full Text | Google Scholar

Spencer, C. C. A., Su, Z., Donnelly, P., and Marchini, J. (2009). Designing Genome-wide Association Studies: Sample Size, Power, Imputation, and the Choice of Genotyping Chip. Plos Genet. 5 (5), e1000477. doi:10.1371/journal.pgen.1000477

PubMed Abstract | CrossRef Full Text | Google Scholar

Styrkarsdottir, U., Halldorsson, B. V., Gretarsdottir, S., Gudbjartsson, D. F., Walters, G. B., Ingvarsson, T., et al. (2008). Multiple Genetic Loci for Bone mineral Density and Fractures. N. Engl. J. Med. 358 (22), 2355–2365. doi:10.1056/NEJMoa0801197

CrossRef Full Text | Google Scholar

Sudlow, C., Gallacher, J., Allen, N., Beral, V., Burton, P., Danesh, J., et al. (2015). UK Biobank: an Open Access Resource for Identifying the Causes of a Wide Range of Complex Diseases of Middle and Old Age. Plos Med. 12 (3), e1001779. doi:10.1371/journal.pmed.1001779

PubMed Abstract | CrossRef Full Text | Google Scholar

Tan, L.-J., Wang, Z.-E., Wu, K.-H., Chen, X.-D., Zhu, H., Lu, S., et al. (2015). Bivariate Genome-wide Association Study ImplicatesATP6V1G1as a Novel Pleiotropic Locus Underlying Osteoporosis and Age at Menarche. J. Clin. Endocrinol. Metab. 100 (11), E1457–E1466. doi:10.1210/jc.2015-2095

PubMed Abstract | CrossRef Full Text | Google Scholar

Thiagalingam, A., De Bustros, A., Borges, M., Jasti, R., Compton, D., Diamond, L., et al. (1996). RREB-1, a Novel Zinc finger Protein, Is Involved in the Differentiation Response to Ras in Human Medullary Thyroid Carcinomas. Mol. Cel Biol 16 (10), 5335–5345. doi:10.1128/mcb.16.10.5335

PubMed Abstract | CrossRef Full Text | Google Scholar

Trajanoska, K., Morris, J. A., Oei, L., Zheng, H.-F., Evans, D. M., Kiel, D. P., et al. (2018). Assessment of the Genetic and Clinical Determinants of Fracture Risk: Genome Wide Association and Mendelian Randomisation Study. BMJ 362, k3225. doi:10.1136/bmj.k3225

PubMed Abstract | CrossRef Full Text | Google Scholar

Visscher, P. M., Wray, N. R., Zhang, Q., Sklar, P., McCarthy, M. I., Brown, M. A., et al. (2017). 10 Years of GWAS Discovery: Biology, Function, and Translation. Am. J. Hum. Genet. 101 (1), 5–22. doi:10.1016/j.ajhg.2017.06.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Warming, L., Hassager, C., and Christiansen, C. (2002). Changes in Bone mineral Density with Age in Men and Women: a Longitudinal Study. Osteoporos. Int. 13 (2), 105–112. doi:10.1007/s001980200001

PubMed Abstract | CrossRef Full Text | Google Scholar

WHO Scientific Group on the Prevention and Management of Osteoporosis, (2003). Prevention and Management of Osteoporosis: Report of a WHO Scientific Group. Geneva: World Health Organization.

Google Scholar

Yang, T.-L., Shen, H., Liu, A., Dong, S.-S., Zhang, L., Deng, F.-Y., et al. (2020). A Road Map for Understanding Molecular and Genetic Determinants of Osteoporosis. Nat. Rev. Endocrinol. 16 (2), 91–103. doi:10.1038/s41574-019-0282-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, X., Bai, W., and Zheng, H. (2021). Twelve Years of GWAS Discoveries for Osteoporosis and Related Traits: Advances, Challenges and Applications. Bone Res. 9 (1), 23. doi:10.1038/s41413-021-00143-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhuang, H. f., Wang, P. w., Li, Y. z., Lin, J. k., Yao, X. d., and Xu, H. (2020). Analysis of Related Factors of Brittle Hip Fracture in Postmenopausal Women with Osteoporosis. Comparison of Three Types of Tha 12 (1), 194–198. doi:10.1111/os.12605

PubMed Abstract | CrossRef Full Text | Google Scholar