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.

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.

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).

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⁻⁵.

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”.

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.

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.

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).

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 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).

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.