Osteoporosis (OP) is a complex skeletal disease characterized by a generalized reduction in bone mineral density (BMD) and microarchitectural deterioration of bone tissue. It is clinically important due to its strong association with fragility fractures, particularly of the hip (HFx). As a multifactorial disease, it is well recognized that the interplay of genetic, environmental, and clinical factors significantly contributes to its development. The primary risk factors associated with OP include aging, low body mass index, prior fractures, and muscle weakness. In addition, numerous studies have demonstrated that genetic factors play an important role in OP susceptibility.

Among the various mechanisms involved in OP, decreased serum estrogen concentrations are a well-established contributor to disease pathogenesis. After menopause, the aromatization of androgenic precursors becomes the principal source of estrogen in women. Therefore, the estrogen receptor alpha (ESR1) gene and the aromatase (CYP19A1) gene, which encodes the enzyme responsible for the conversion of androgens to estrogens, are considered strong candidate genes for OP. Associations between single-nucleotide variants (SNVs) in ESR1 and CYP19A1 and osteoporosis have been extensively reported, although findings have varied across populations. For example, a meta-analysis found associations between ESR1 SNVs and OP that differed by ethnicity (Zhu et al., 2018; Shu et al., 2020). Similarly, one haplotype in CYP19A1, formed by SNVs rs700518, rs767199, rs4775936, and rs1062033, which are highly frequent in the Spanish population, was associated with increased OP risk and low CYP19A1 expression (Riancho et al., 2007). Individually, rs700518 and rs1062033 were associated with OP in a meta-analysis (Ma et al., 2020) and in Spanish women (Riancho et al., 2009), but not in Polish women (Kamiński et al., 2019).

Inflammatory cytokines also play an important role in OP. Interleukin 6 (IL6), interleukin 1 beta (IL1β), and their respective receptors exert significant effects on osteoclasts genesis and maturation, and bone resorption. Genetic variants in IL6 and IL6R have been associated with OP risk. One of the most studied SNVs in IL6 is rs1800795 (G174C), which has been described as a protective factor against osteoporosis in Caucasian and Asiatic populations (Ferrari et al., 2003; Chung et al., 2003). As part of a haplotype with other variants, rs1800795 was also associated with a favorable response to alendronate treatment in Mexican women (Villagómez Vega et al., 2023). Another SNV in IL6, rs1800796 (G572C), has been linked to BMD variations in Japanese women (Oishi et al., 2012) and a higher risk of vertebral compression fractures in Chinese women (Xiong et al., 2022). However, in Iranian patients, this SNV was not associated with BMD, but was linked to serum calcium levels (Eftekhari et al., 2018). The IL1β polymorphism rs1143634 has been associated with low BMD and OP risk in Taiwanese women (Chen et al., 2003) but not in Polish women (Czerny et al., 2010).

Disorders in bone metabolism are central to the pathogenesis of osteoporosis. The signaling system comprising the receptor activator of nuclear factor-κB (RANK), its ligand (RANKL), and osteoprotegerin (OPG) play a crucial role in coupling osteoblast and osteoclast activity to regulate bone remodeling. Conflicting results have been reported regarding the association of these genes with OP. The A allele of RANK rs3018362 has been associated with low BMD and an increased risk of fracture in European and Chinese populations (Paternoster et al., 2010; Liu et al., 2010). Conversely, in another study in Chinese women, the G allele of this SNV was associated with low BMD (Shang et al., 2013).

A range of lifestyle factors, including physical activity, smoking, alcohol consumption, and calcium and protein intake (Zhu and Prince, 2015), as well as the presence of comorbidities, may interact with genetic predispositions and influence OP risk. Most studies have focused on single-gene effects and have not considered how combinations of genetic variants and environmental exposures jointly influence osteoporosis risk, particularly in Mexican women, an underrepresented population in genetic osteoporosis research. This approach could also help explain the inconsistencies in results reported across population-based studies on osteoporosis risk.

Given the inconsistent associations of CYP19A1, ESR1, IL6, IL6R, IL1β, RANK and RANKL genes across populations and the limited research examining interactions effects, we aimed to evaluate both individual SNVs and gene–gene and gene–environment interactions in relation to OP and HFx risk in Mexican women.

A total of 609 women were included in this study: 169 with OP, 205 with HFx, and 235 controls. Characteristics of the study population are summarized in Table 2. The mean age was 57.78 ± 8.44 years in the control group, 69.62 ± 9.88 years in the OP group, and 77.7 ± 10.89 years in the HFx group. The mean body mass index (BMI, kg/m²) was 28.76 ± 4.01 in the control group, 25.05 ± 3.99 in the OP group, and 24.87 ± 4.6 in the HFx group. Variables such as age at menarche, number of pregnancies, and estrogen intake showed significant differences between groups, while age at menopause and smoking habits did not. Regarding comorbidities, hypertension and diabetes differed between the OP and HFx groups compared with controls. The incidence of arthritis also differed significantly between the OP and control groups (Table 2).

Regarding fracture risk predictions using the FRAX tool, the 10-year probability of a major fracture without BMD was significantly higher in both the OP and HFx groups compared with controls (p = 0.0001 in both cases). Although the 10-year probability of a major fracture based on BMD was not calculated for the HFx group, significant differences were observed between the OP and control groups (p = 0.0001) (Table 2).

All SNVs were in HWE (p > 0.05), except for rs17523880 in CYP19A1 and rs1800796 in IL6. The A allele of RANK rs3018362 was significantly more frequent in OP cases than in controls (51.0% vs. 42.0%) and was associated with increased OP risk (OR = 1.5, 95% CI: 1.1–1.9, p = 0.007). No other alleles showed significant differences in distribution between the OP or HFx groups and controls (Table 3).

Genotypic analysis revealed that RANK rs3018362 was the most significantly associated SNV with OP. Under the codominant model, both GA and AA genotypes were associated with increased risk in unadjusted analyses, and these associations remained significant after adjustment (Adjusted OR = 2.08, 95% CI: 1.08–3.98 for GA; and 2.76, 95% CI: 1.21–6.30 for AA). The dominant model also yielded a significant association (Adjusted OR = 2.23, 95% CI: 1.19–4.17). For ESR1 rs2234693, the CC genotype was associated with a reduced OP risk under both codominant and recessive models, but only in the adjusted analysis (Adjusted OR = 0.28, 95% CI: 0.11–0.69 and 0.25, 95% CI: 0.10–0.61, respectively) (Supplementary Table S1A).

In the HFx group, the CC genotype of CYP19A1 rs700518 was associated with lower fracture risk in both the codominant and recessive models (OR = 0.45, 95% CI: 0.21–0.96 and OR = 0.47, 95% CI: 0.22–0.98, respectively), although these associations lost significance after adjustment. The CC genotype of ESR1 rs2234693 remained significantly associated with reduced fracture risk in the adjusted recessive model (OR = 0.30, 95% CI: 0.12–0.75). Conversely, the AA genotype of RANK rs3018362 was associated with increased fracture risk in the adjusted codominant model (OR = 2.4, 95% CI: 1.01–6.06) (Supplementary Table S1B).

Haplotypes of CYP19A1, ESR1, and IL6R were constructed separately for each gene to assess their associations with OP and HFx. However, no statistically significant associations were observed (Table 4).

In the MDR analysis, a significant gene–gene interaction was identified between ESR1 rs2228480 and RANK rs3018362 in OP, with a testing accuracy of 0.5656 and a CVC of 10/10. This interaction was also associated with increased OP risk (OR = 2.1, 95% CI: 1.4–3.2) (Table 5; Figure 1, panel A). In the HFx group, the top gene–gene interaction model included ESR1 (rs2228480), IL6R (rs4845617), IL1β (rs16944), RANK (rs3018362), and RANK-L (rs12585014), with a testing accuracy of 0.5577 and a CVC of 10/10 (Table 6; Figure 2). A significant gene–environment interaction was observed in OP, involving CYP19A1 SNVs (rs700518, rs1062033, rs4775936, and rs767199), IL1β (rs16944), and 10-year probability of hip fracture (FRAXH-BMD). This model yielded a testing accuracy of 0.8488 and a CVC of 10/10 (Table 5; Figure 1, panel B). No significant gene–environment interactions were detected in the HFx group (Table 6).

Our results demonstrate an association between variants in ESR1 and RANK and the risk of both OP and HFx. The ESR1 SNVs have been implicated in numerous conditions including cartilage and bone diseases. Variants rs9340799 and rs2228480 are associated with risk of osteoarthritis, but rs2234693 was not in a meta-analysis from many populations throughout the world (Ma et al., 2015). In Mexican population rs2234693 and rs9340799 SNVs, were associated with a reduced risk of osteoarthritis (Borgonio-Cuadra et al., 2012).

With regard to OP and Fx, the variant rs2234693 has previously been associated with higher T-scores in postmenopausal Slovak women (Mondockova et al., 2018), and with increased femoral neck BMD and Z-scores in Asian populations (Zhu et al., 2018).

Particularly in Mexican population this variant has been linked to a reduced risk of OP and HFx in Mexican women, (García-Rojas et al., 2022), and of distal radius fracture (Farías et al., 2020). The variants rs3020331, rs3020404 and rs4870044 showed association with a reduced risk of osteoporosis and/or hip fracture in a cohort of Mexican women (Hidalgo-Bravo et al., 2019). These results are consistent with our actual findings showing the protective association of ESR1 SNVs in independent-effect analyses. There is one study with Mexican postmenopausal women, were the rs2234693 did not show association with BMD variations (Rojano-Mejía et al., 2014), However, in Mexican population with only limited exceptions ESR1 variants are usually associated as protective factors against low BMD and osteoporosis.

Variants of the RANK/RANKL/OPG genes have showed association with several bone conditions. For example, RANK rs3826620 and RANKL rs9594738 were associated with periodontitis (Petean et al., 2019). Genetic variants of RANK/RANKL/OPG increase the risk of rheumatoid arthritis in Mexican individuals (Nava-Valdivia et al., 2024) and Chinese populations (Pei et al., 2025) and the OPG rs2073618 is related with increased risk of Legg-Calvé-Perthes disease in Mexican children (Cruz-Ortíz et al., 2024).

With regard to osteoporosis and fragility fractures, previous reports associated the RANK rs3018362 SNV with reduced cortical BMD, increased bone resorption, with low-trauma and atraumatic fractures in several Caucasian populations (Chen et al., 2003; Kemp et al., 2014; Styrkarsdottir et al., 2008). In Mexican women, it was reported association between this variant and elevated OP risk (Casas-Avila et al., 2019). In the present study, was associated with increased risk of OP and HFx in Mexican women, consistent with previous reports that support the role of RANK rs3018362 as a risk variant across diverse populations.

The variants in IL6 gene explored in this work were not associated with OP or Fx in the individual analysis. As occurred in previous investigations with Mexican population, where the rs1800795 variant of IL6 did not showed association with risk of fracture (Ponce de León-Suárez et al., 2018) nor with BMD variations in postmenopausal obese women (Méndez et al., 2013). However, it showed an important association with severe radiographic damage of the hands in patients with rheumatoid arthritis in Mexican adult patients (Zavaleta-Muñiz et al., 2016). In the other hand, the rs1800796 (GG genotype) of IL6 has been associated with reduced risk of hip fracture in Mexican women (Ponce de León-Suárez et al., 2018) but it was not associated in this study. Few studies examined the effects of IL6R variants on OP in Mexican population, but the SNVs rs4845617 and rs2228145 were not associated with hip fragility fracture in postmenopausal Mexican women (Ponce de León-Suárez et al., 2018) just as happened in this study.

Aromatase contributes importantly to bone metabolism. Variants in CYP19A1 gene have been studied with regard to the risk of vertebral fractures (Koudu et al., 2012) and with regard to the BMD in women with hormonal replacement therapy (Masi et al., 2014). In Mexican population the rs700518 was associated to response to pharmacological treatment for OP (Villagómez Vega et al., 2023) and (TTTA)n microsatellite combined with a TCT deletion were associated to HFx risk (Casas-Avila et al., 2015), but they have not been widely studied in Mexican. As well as IL6 and IL6R variants, CYP19A1 SNVs were not associated directly with OP or HFx in this study; however, they associate when interact with other factors, emphasizing the value of interaction analyses for detecting risk-modifying effects that are not apparent in single-variant models.

To the best of our knowledge, no prior study has simultaneously evaluated all these genes and comprehensively examined both gene–gene and gene–environment interactions in the context of OP and HFx risk. Epistasis, defined as gene–gene interactions between two or more loci, may alter disease risk independently of individual gene effects (Perelygin et al., 2025). Recent research has explored the role of epistasis in OP. For example, an interaction between RMBS3 and ZNF516 has been shown to influence BMD in Caucasian and African populations (Yang et al., 2013), while rs3751143 of P2X7R interacting with ESR1 rs2234693 increases OP susceptibility in postmenopausal Chinese women (Wang et al., 2018). Additionally, a VDR–TNF-α interaction has been associated with OP risk in elderly women (Liao et al., 2019). In our study, a significant gene–gene interaction between ESR1 and RANK was observed in women with OP, as well as a gene–environment interaction involving CYP19A1, IL1β, and the 10-year probability of major fractures. Although in the graphic representation might seem that there is no gene-environment interaction in OP (Figure 1B), positive IG values between nodes and values in Table 5 for the best model, show synergistic interactions between these factors. For HFx, a gene–gene interaction was identified involving ESR1, IL6R, IL1β, RANK and RANKL, while no significant gene–environment interactions were detected in this group.

Bone cells, including osteocytes, osteoclasts, and osteoblasts, are major targets of estrogen. Estrogens, through their receptors (ESR), contribute to bone homeostasis by reducing bone resorption, partly through inhibition of RANKL expression in osteoblasts and increased production of osteoprotegerin (OPG). However, this protective mechanism declines after menopause, when bone resorption begins to exceed bone formation (Fischer and Haffner-Luntzer, 2022). The interaction observed between ESR1 and RANK in our study may indicate that in postmenopausal women, the role of ESR1 shifts from protective to a risk-modifying factor when interacts with RANK. Notably, this interaction involved SNV’s that were not significant when analyzed independently.

In HFx, the interaction among ESR1, IL6R, IL1β, RANK and RANKL genes conferred a substantially increased fracture risk. It can be explained because IL1β and IL6, and their respective receptors, are key mediators of inflammatory signaling and bone resorption in osteoporosis and other diseases leading bone damage (Tseng et al., 2022; Takeuchi et al., 2021). As estrogen deficiency also downregulates OPG and increases RANKL activity, the combined effect of proinflammatory and osteoclast-activating signals may account for the high fracture susceptibility observed in the presence of these interacting variants.

We also identified a gene–environment interaction in OP involving CYP19A1, and the 10-year probability of major fractures, in addition to IL1β (Table 6). This finding is biologically plausible, as aromatase encoded by CYP19A1 gene, is responsible for converting androgens to estrogens, a crucial source of estrogens in postmenopausal women, and the rs700518 and rs1062033 CYP19A1 variants have been previously associated with reduced expression (Riancho et al., 2007). All these facts likely may contribute to OP in our sample.

We recognized that our study could have some limitations. As a hospital-based case–control study, it carries an inherent risk of selection bias. To address this concern, we took comprehensive steps to minimize bias in our research design. To avoid selection bias, we carefully assessed cases with OP and HFx alongside their corresponding controls. This was accomplished by applying specific inclusion criteria for both groups, ensuring that comparisons were valid and reliable. Besides, we also control potential confounders. This was achieved through the use of exclusion criteria and by applying multivariate statistical analysis during data evaluation. These measures enhanced the accuracy of our findings by reducing the influence of variables outside the scope of our study. Population stratification is a known complication in association studies involving admixed populations. Given that Mexicans represent an admixed population, we took steps to limit stratification effects. Specifically, cases and controls were selected from particular regions of Mexico to minimize the impact of subpopulation differences. Additionally, we utilized a validated questionnaire to assess ethnic ancestry and anthropometric variations among participants. In summary, while we acknowledge the limitations inherent in our study, we are confident that the main strength of our work lies in the careful and appropriate selection of study groups. This approach enhances the overall validity and reliability of our research findings.

Our findings indicate that both individual effect of genetic variants as well as gene–gene and gene–environment interactions contribute to OP and HFx susceptibility in Mexican women. These results highlight the importance of evaluating polygenic and environmental interactions rather than isolated genetic effects, particularly in populations with distinct complex ancestral backgrounds.

This study contributes to a more comprehensive understanding of osteoporosis genetics in Mexican women. Future studies with larger and multi-ethnic cohorts are necessary to validate and expand upon these findings.