Causal Effects of Female Reproductive and Hormonal Factors on Osteoporosis, Bone Mineral Density, and Osteoarthritis: A Two-Sample Mendelian Randomization Study

Introduction

Osteoporosis and its consequent fragility fractures are major public health concerns, particularly with a globally aging population. These conditions lead to substantial morbidity, mortality, and socioeconomic burden.^1,2 The global burden of disability-adjusted life-years (DALYs) and deaths attributable to low bone mineral density (BMD) and related fractures has increased significantly over the past three decades.^3–5 This burden is particularly substantial in low- and middle-income countries, where rapid population aging and limited screening resources exacerbate the impact of osteoporosis and osteoarthritis,^6–8 underscoring the urgent need for effective prevention and management strategies. The implementation of preventive measures necessitates a comprehensive understanding of the factors influencing bone metabolism. Bone metabolism is regulated by various factors, among which sex hormones play a pivotal role.^9–11 Of the many sex hormones that contribute to bone health and disease, estrogens have been the most widely researched due to the clear association between estrogen deficiency during menopause and menopause-related bone diseases.^12,13 Mechanistically, multiple studies have demonstrated that estrogens and their receptors can influence the sensitivity of bone cells to mechanical loading and the subsequent bone cell mechanotransduction process, highlighting the critical role of estrogens in bone health.^14,15 Estrogen deficiency, most notably following menopause, disrupts the balance between bone resorption and formation, leading to accelerated bone loss and a decline in bone strength, thereby increasing fracture risk.^16–18

Alongside osteoporosis, osteoarthritis (OA) represents another major cause of disability in aging populations. OA exhibits a marked sex disparity, with higher prevalence and greater severity in women, particularly after menopause. Estrogen receptors (ER), including ERα and ERβ, are expressed in articular cartilage and subchondral bone, supporting a role for estrogen signaling in chondrocyte metabolism and extracellular matrix homeostasis.¹⁹ Experimental evidence suggests that estrogen modulates inflammatory pathways and cartilage turnover; however, epidemiological studies report inconsistent associations between estrogen exposure, hormone replacement therapy (HRT), and OA risk.^20,21 These discrepancies highlight the need for causal inference approaches to clarify hormonal effects in joint disease. Investigating both osteoporosis and OA simultaneously is crucial, as they share risk factors like age and hormonal changes, yet their underlying pathophysiology and response to hormones may differ, creating a complex clinical challenge.

Throughout a woman’s lifespan, reproductive and hormonal Factors significantly influence endogenous estrogen levels. Observational studies have explored the association between various female reproductive and hormonal factors and bone health. For instance, later age at menarche, which corresponds to a shorter duration of endogenous estrogen exposure, has been linked to lower BMD in postmenopausal women.²² Similarly, the use of HRT is well-documented to prevent postmenopausal bone loss and reduce fracture rates.²³ However, the evidence landscape for other common exposures is fraught with controversy. For instance, the net effect of oral contraceptive (OCP) use on peak bone mass remains debated, while the lifelong skeletal impact of parity is still unclear. Crucially, the role of estrogen in OA (OA) is particularly enigmatic, with observational studies reporting conflicting evidence of both protective and harmful effects.^24–27 These observational findings are susceptible to biases from residual confounding (eg., lifestyle, socioeconomic status) and reverse causation, which obscure true causal relationships. Therefore, a robust etiological investigation is urgently needed to disentangle these complex associations.

Mendelian randomization (MR) is a powerful genetic epidemiological method that uses genetic variants as instrumental variables (IVs) to infer the causal effect of an exposure on an outcome. By leveraging the random allocation of genetic variants at conception, MR minimizes confounding and is not prone to reverse causation, thus emulating a natural randomized controlled trial.^28,29 This approach can provide more robust evidence on the causal nature of associations between hormonal exposures and bone health.

Observational studies have linked individual reproductive factors with bone outcomes, but these associations may be confounded by lifestyle and reverse causation.^22,30,31 Recent MR studies have evaluated age at menarche,³² age at menopause,³³ and selected reproductive traits³⁴ in relation to osteoporosis or specific BMD measures, such as total or site-specific BMD, with some including age-stratified BMD but without assessing OA. To date, no MR study has systematically integrated multiple reproductive and hormonal factors while jointly evaluating osteoporosis, OA, and age-stratified BMD within a single analytical framework.

Therefore, this study aimed to conduct a comprehensive two-sample MR analysis to systematically evaluate the causal effects of a wide range of female reproductive and hormonal factors, categorized as reproductive timing (age at menarche, age at menopause, menstrual cycle length), reproductive history (age at first and last live birth, parity-related traits), and hormonal interventions (OCP use, HRT, and bilateral oophorectomy), on the risk of osteoporosis, BMD at different life stages, fragility fractures, and OA. By providing evidence on these causal pathways, we aim to inform risk stratification, identify potential targets for intervention, and ultimately contribute to the development of personalized prevention strategies for these debilitating bone and joint diseases.

Material and Methods

Study Design

This study employed a two-sample MR framework to assess the causal relationships between 12 reproductive and hormonal exposures and 7 bone-related outcomes. The design relies on three core assumptions for the genetic variants used as instrumental variables (IVs): (i) they are robustly associated with the exposure of interest; (ii) they are not associated with any confounders of the exposure-outcome relationship; and (iii) they affect the outcome only through the exposure. Because genetic variants are fixed at conception, MR estimates reflect the cumulative lifetime effect of genetically predicted exposures. An overview of the study design is presented in Figure 1. This study is reported in accordance with the STROBE guidelines (Supplementary File 1).

Figure 1 Schematic overview of the two-sample Mendelian randomization study design. Instrumental variables (SNPs) were selected based on predefined criteria and used to estimate the associations between reproductive/hormonal exposures and skeletal outcomes. The dashed lines represent the three core MR assumptions. Red cross symbols indicate potential violations of MR assumptions, including confounding pathways or horizontal pleiotropy.

Data Sources

All analyses were based on publicly available summary-level data from large-scale genome-wide association studies (GWAS), primarily involving individuals of European ancestry. The 12 exposures included factors related to the reproductive cycle and fertility history (age at menarche, age at menopause, length of menstrual cycle, age at first live birth, age at last live birth) and hormonal interventions (ever used HRT, age started HRT, age last used HRT, age started OCP, age when last used OCP, bilateral oophorectomy, and age at bilateral oophorectomy). The 7 outcomes comprised osteoporosis (GCST90086118), total body bone mineral density (BMD) (ebi-a-GCST005348) and its age-stratified measures (for ages >60, 45–60, and 30–45; ebi-a-GCST005349, ebi-a-GCST005350, ebi-a-GCST005346, respectively), fragility fractures (finn-b-OSTEOPOROSIS_FRACTURE_FG), and OA (GCST90566795).

Because exposure data were derived from UK Biobank-based GWAS, potential sample overlap with outcome datasets was evaluated. Osteoporosis (GCST90086118) was based on a U.S. population and did not include UK Biobank participants. For the BMD GWAS,³⁵ UK Biobank contributed 1,553 individuals (2.76% of 56,284), indicating minimal overlap. For OA (GCST90566795; PMID: 40205036),³⁶ UK Biobank contributed 425,643 participants (21.69% of 1,962,069), indicating partial overlap. Detailed information on the GWAS consortia, sample sizes, and populations for each trait is provided in Supplementary Table 1.

Instrumental Variable Selection

For each exposure, we selected single-nucleotide polymorphisms (SNPs) as IVs following a rigorous protocol. Initially, SNPs strongly associated with the exposure were identified at a genome-wide significance threshold (P < 5 × 10⁻⁸). For exposures with fewer than five significant SNPs at this level (Age last used HRT, Age at bilateral oophorectomy, Age started OCP, Age when last used OCP), the threshold was relaxed to P < 5 × 10⁻⁶ to ensure sufficient statistical power. This relaxation was applied to maintain adequate instrument number and estimation precision for traits with limited genome-wide significant variants. To ensure independence between IVs, we performed linkage disequilibrium (LD) clumping using a strict threshold (r² < 0.001, window size = 10,000 kb) based on the European 1000 Genomes Project reference panel. SNPs with a minor allele frequency (MAF) ≤ 0.01 were excluded. When an IV was not available in the outcome GWAS, we searched for a proxy SNP in high LD (r² > 0.8). We harmonized the exposure and outcome datasets to ensure that the effect of each SNP on the exposure and outcome corresponded to the same allele. Palindromic SNPs with ambiguous strand orientation were removed. The strength of the selected IVs was assessed using the F-statistic, calculated as F = R²(N-2)/(1-R²), where R² is the proportion of variance in the exposure explained by the IV and N is the sample size of the exposure GWAS. An F-statistic > 10 was considered indicative of a sufficiently strong instrument. In our analysis, all retained IVs had F-statistics > 10, suggesting that weak instrument bias was unlikely.

Mendelian Randomization Analysis

The primary analysis was conducted using the random-effects inverse-variance weighted (IVW) method, which combines the Wald ratios of each SNP to provide a pooled causal estimate. The IVW method assumes either the absence of horizontal pleiotropy or that pleiotropy is balanced. To assess the robustness of our findings, we employed several complementary MR methods. The MR-Egger regression provides a causal estimate that is robust to directional pleiotropy, although it has lower statistical power. The weighted median method can provide a consistent estimate even if up to 50% of the IVs are invalid. The weighted mode method provides an estimate based on the largest cluster of SNPs with similar causal effects.

For BMD outcomes, the original GWAS has adjusted for age, height, weight, principal components, and study-specific covariates, and applied inverse normal transformation to residuals. Thus, SNP effect sizes represent per 1 standard deviation (SD) change in standardized BMD. In MR analyses, causal estimates (β) were exponentiated (odds ratio [OR] = e^β), and reported ORs therefore reflect the change in outcome risk per 1 SD genetically predicted difference in BMD.

Sensitivity Analysis

We performed a comprehensive suite of sensitivity analyses to evaluate the validity of the MR assumptions. Heterogeneity among the IVs was assessed using Cochran’s Q statistic for the IVW method. Significant heterogeneity (P < 0.05) may indicate the presence of pleiotropic effects. The intercept term from the MR-Egger regression was used to test for directional pleiotropy, with a P-value < 0.05 suggesting its presence. The MR-PRESSO (Pleiotropy RESidual Sum and Outlier) test was used to detect and correct for horizontal pleiotropy by identifying and removing outlier SNPs. A leave-one-out sensitivity analysis was conducted by systematically removing one SNP at a time and recalculating the IVW estimate to assess whether the causal association was driven by a single influential SNP. Finally, we applied the Steiger directionality test to verify that the genetic variants used as IVs had a stronger effect on the exposure than on the outcome, thus supporting the assumed causal direction.

Statistical Analysis

All statistical analyses were performed using the “TwoSampleMR” package in R (version 4.0.5). To account for multiple comparisons across the numerous exposure-outcome pairs tested, we applied the False Discovery Rate (FDR) correction. A causal association was considered significant if the FDR-corrected P-value was less than 0.05.

Results

Instrumental Variable Selection and Characteristics

After a rigorous selection and quality control process, a set of independent and strong IVs was identified for each of the 12 exposures. The number of SNPs per exposure ranged from 6 for Bilateral oophorectomy to 201 for Age at menarche. The F-statistic for each IV was substantially greater than the conventional threshold of 10 (mean F-statistics ranged from 22.69 to 82.84), indicating that weak instrument bias was unlikely to affect our results. Exposure-specific instrumental variables are provided in Supplementary Table 2, and the harmonized SNP-level datasets used in the MR analyses are presented in Supplementary Table 3.

Main Mendelian Randomization Results

After applying the FDR correction, our primary MR analysis, conducted after removing pleiotropic outliers, revealed several significant associations. A complete summary of the estimates from all MR methods is provided in Table 1. Detailed multi-method forest plots are provided in Figure 2. For consistency, all effect estimates reported here are ORs from the primary IVW analysis.

Table 1 Causal Estimates of Reproductive and Hormonal Factors on Bone-Related Outcomes from the Initial Mendelian Randomization Analysis

Figure 2 Continued.

Figure 2 Forest plots of Mendelian randomization estimates for reproductive factors and skeletal outcomes before and after outlier removal. Forest plots display odds ratios (ORs) with 95% confidence intervals (CIs) derived from inverse-variance weighted (IVW), MR-Egger, weighted median, and weighted mode methods. (A) Estimates before removal of pleiotropic outliers. (B) Estimates after removal of pleiotropic outliers identified by MR-PRESSO. The vertical dashed line represents the null value (OR = 1).

Reproductive Lifespan, Menstrual Cycle, and Bone Health

As shown in Figure 3, factors related to a shorter endogenous estrogen exposure window were associated with poorer bone health. Genetically predicted later age at menarche was associated with an increased risk of osteoporosis (IVW OR: 1.59, 95% CI: 1.24–2.04; FDR = 0.008, Figure 3A). This was complemented by findings indicating that later menarche led to lower total body BMD (IVW OR: 0.88, 95% CI: 0.82–0.94; FDR = 0.004, Figure 3B) and specifically lower BMD in the 45–60 age group (IVW OR: 0.87, 95% CI: 0.79–0.95; FDR = 0.028, Figure 3C).

Figure 3 Continued.

Figure 3 Forest plots of significant causal associations between reproductive cycle factors and bone health outcomes from the initial analysis (before outlier removal). The plots display the Odds Ratios (ORs) and 95% confidence intervals (CIs) for significant causal associations identified in the initial analysis. The vertical dashed line represents the null effect (OR = 1). (A) Age at menarche and osteoporosis. (B) Age at menarche and total body bone mineral density (BMD). (C) Age at menarche and BMD (age 45–60). (D) Age at menopause and osteoporosis. (E) Age at menopause and BMD (age 45–60). (F) Length of menstrual cycle and total body BMD. (G) Length of menstrual cycle and BMD (age 45–60). Abbreviations: OR, Odds Ratio; CI, Confidence Interval; BMD, Bone Mineral Density.

In contrast, a longer reproductive lifespan due to later age at menopause showed protective effects, significantly reducing the risk of osteoporosis (IVW OR: 0.74, 95% CI: 0.62–0.88; FDR = 0.009, Figure 3D). Furthermore, later menopause was associated with higher BMD in women aged 45–60 (IVW OR: 1.13, 95% CI: 1.03–1.24; P=0.008, FDR = 0.057, Figure 3E).

We also found that a longer length of the menstrual cycle was associated with detrimental effects on bone health, specifically leading to lower total body BMD (IVW OR: 0.83, 95% CI: 0.73–0.94; FDR = 0.028, Figure 3F) and lower BMD in the 45–60 age group (IVW OR: 0.73, 95% CI: 0.61–0.89; FDR = 0.016, Figure 3G). Corresponding scatter plots are shown in Figure 4A–G.

Figure 4 Scatter plots for significant associations between reproductive cycle factors and bone health outcomes (before outlier removal). Each point represents a single nucleotide polymorphism (SNP). The x-axis shows the SNP’s effect on the exposure, and the y-axis shows its effect on the outcome. The slopes of the lines represent the causal estimate from different MR methods (IVW, MR-Egger, Weighted Median, Weighted Mode) in the initial analysis. (A) Age at menarche and osteoporosis. (B) Age at menarche and total body bone mineral density (BMD). (C) Age at menarche and BMD (age 45–60). (D) Age at menopause and osteoporosis. (E) Age at menopause and BMD (age 45–60). (F) Length of menstrual cycle and total body BMD. (G) Length of menstrual cycle and BMD (age 45–60).

Reproductive History, Hormonal Interventions, and Bone Outcomes

Our analysis identified significant protective effects of reproductive timing on OA. Later age at both first live birth (IVW OR: 0.81, 95% CI: 0.76–0.86; FDR < 0.001, Figure 5A) and last live birth (IVW OR: 0.75, 95% CI: 0.65–0.87; FDR = 0.006, Figure 5B) were associated with a reduced risk of OA.

Figure 5 Continued.

Figure 5 Forest plots of significant causal associations between reproductive history, hormonal interventions, and bone health outcomes (before outlier removal). The plots display the Odds Ratios (ORs) and 95% confidence intervals (CIs). (A) Age at first live birth and osteoarthritis. (B) Age at last live birth and osteoarthritis. (C) Age at last use of oral contraceptive pill and bone mineral density (BMD) (age 45–60). (D) Age at start of oral contraceptive pill and osteoarthritis. (E) Bilateral oophorectomy and osteoarthritis. Abbreviations: OR, Odds Ratio; CI, Confidence Interval; BMD, Bone Mineral Density.

Hormonal interventions demonstrated complex and distinct associations with bone health. A later age at the last use of OCP was associated with lower BMD in the 45–60 age group (IVW OR: 0.65, 95% CI: 0.48–0.88; FDR = 0.048, Figure 5C). However, a later age at starting OCP use was linked to a reduced risk of OA (IVW OR: 0.88, 95% CI: 0.82–0.95; FDR = 0.009, Figure 5D). Interestingly, having a bilateral oophorectomy was associated with a decreased risk of OA (IVW OR: 0.93, 95% CI: 0.89–0.98; FDR = 0.052, Figure 5E). Given the borderline FDR value and the unexpected direction of effect, this association should be regarded as hypothesis-generating and interpreted with caution. Corresponding scatter plots are shown in Figure 6A–E. For all other exposure-outcome pairs tested, no significant associations were observed after FDR correction (Table 1). As age when last used OCP, age started OCP, and bilateral oophorectomy were instrumented using fewer SNPs selected under relaxed genome-wide thresholds, the corresponding effect estimates may be less precisely estimated and should be interpreted as supportive.

Figure 6 Scatter plots for significant associations between reproductive history, hormonal interventions, and bone health outcomes (before outlier removal). (A) Age at first live birth and osteoarthritis. (B) Age at last live birth and osteoarthritis. (C) Age at last use of oral contraceptive pill and bone mineral density (BMD) (age 45–60). (D) Age at start of oral contraceptive pill and osteoarthritis. (E) Bilateral oophorectomy and osteoarthritis. Each point represents a single nucleotide polymorphism (SNP). The x-axis shows the SNP’s effect on the exposure, and the y-axis shows its effect on the outcome. The slopes of the lines represent the causal estimate from different MR methods in the initial analysis.

Sensitivity Analysis

Sensitivity analyses demonstrated significant heterogeneity (Cochran’s Q P < 0.05) for several exposure-outcome relationships, including age at menopause with OA (Q = 249.801, P < 0.001). MR-Egger intercept tests revealed no significant directional pleiotropy except for age started HRT with osteoporosis (P = 0.026) (Table 2). MR-PRESSO identified influential outliers in key associations such as age at menarche and total body bone mineral density (Global P = 0.039), and age at last OCP use and total body BMD (age 45–60) (Global P = 0.045) (Table 3). Alternative methods (weighted median, MR-Egger) yielded directionally consistent results with IVW for primary findings like age at first birth with osteoporosis (all methods P < 0.05), while age at menarche with OA showed robust consistency across methods. Steiger tests confirmed correct directionality for all significant associations except 4 exposure-outcome pairs (Table 4).

Table 2 Heterogeneity and Pleiotropy Diagnostics from the Initial Analysis

Table 3 MR-PRESSO Test Results for Outlier Detection and Distortion Analysis

Table 4 Steiger Directionality Test Results

We re-analyzed the data after excluding these outliers. The results from the four MR methods are summarized in Supplementary Table 4. The forest plots for reproductive cycle factors are shown in Supplementary Figure 1A–H. Notably, for the association between age at menopause and BMD in the 45–60 age group, the FDR p-value decreased from 0.057 before outlier removal to 0.004 after exclusion (Supplementary Figure 1H). The forest plots for reproductive history and hormonal interventions are presented in Supplementary Figure 2A–F. For the association between ever use of HRT and total body BMD (age 45–60), the FDR p-value changed from 0.994 before outlier removal to 0.004 after exclusion (Supplementary Figures 2F).

After removing outliers identified by MR-PRESSO (Supplementary Table 5), heterogeneity remained in some associations, such as age at menarche with total body BMD (Cochran’s Q P < 0.001; Supplementary Table 6), but the direction and significance of the estimates were consistent across different MR methods (IVW, Weighted Median, MR-Egger). Importantly, the MR-Egger intercept test did not detect significant directional pleiotropy for any of the reported associations (Supplementary Table 6), and the MR-PRESSO global test identified no additional outliers among the reported associations (all P > 0.05; Supplementary Table 7). This conclusion was further supported by the general symmetry observed in the corresponding scatter plots (Supplementary Figure 3A–H; Supplementary Figure 4A–F) and funnel plots (Supplementary Figures 5A–H and 6A–F), suggesting that the primary IVW estimates are unlikely to be biased by directional pleiotropic effects. Leave-one-out analyses confirmed that no single SNP was disproportionately driving any of the significant associations (Supplementary Figures 7A–H and 8A–F). Finally, the Steiger directionality test supported the assumed direction (from exposure to outcome) for all reported significant findings (Supplementary Table 8), with the exception of one pair, which was excluded from the main results. Collectively, these sensitivity analyses support the overall robustness of our primary findings.

Discussion

This study employed a two-sample Mendelian randomization approach to systematically elucidate the associations between female reproductive and hormonal factors on osteoporosis and OA. Our primary findings indicate that later age at menarche was associated with higher osteoporosis risk and that later age at menopause and HRT use were associated with higher BMD. In addition, we observed that advanced age at first birth reduces OA risk. These results not only validate hypotheses from previous observational studies but also provide a more comprehensive genetic framework by incorporating 12 distinct exposures, thereby establishing a foundation for bone health risk stratification and personalized intervention strategies.

While estrogen’s protective role in bone mass is well-established, its effect on joint cartilage remains complex and debated, motivating our concurrent investigation of both diseases. Our finding that genetically predicted later age at menarche was associated with higher osteoporosis risk and decreased BMD aligns with and strengthens the conclusions from both observational studies and a recent MR study in an East Asian population.^22,37 This consistency across different study designs and populations underscores the importance of the total duration of endogenous estrogen exposure for achieving and maintaining peak bone mass.

The protective effect of later age at menopause on BMD and osteoporosis risk is also consistent with its role in extending the period of estrogen protection against accelerated bone turnover.³⁸ The underlying biological mechanism is well-established: estrogen inhibits osteoclast differentiation and activity, primarily through the RANKL/OPG signaling pathway, and may also promote osteoblast function, thereby maintaining a balance that favors bone formation over resorption.^16,17,39

The study also suggests a protective effect of HRT on BMD in mid-life (45–60 years), which is consistent with the known therapeutic benefits of estrogen replacement in preventing postmenopausal bone loss.^40,41 Interestingly, we found that later age at last use of OCPs was associated with lower BMD in the 45–60 age group. The relationship between OCP use and bone health has been controversial, with some studies suggesting a neutral or beneficial effect, particularly with modern formulations.^25,42 The result may reflect complex interactions between the timing of OCP use, dosage, and an individual’s underlying hormonal milieu that warrant further investigation. For instance, certain progestins—particularly early-generation formulations—may exhibit mild androgenic effects, as seen with levonorgestrel, whereas fourth-generation progestins such as drospirenone possess anti-androgenic properties. Given that androgens also influence bone mass, we hypothesize that the observed reduction in BMD among individuals aged 45–60 years may be associated with their historical use of OCP containing earlier progestin formulations.⁴³

In addition to the primary osteoporosis and BMD findings, we observed a protective effect of later age at first and last live birth against OA. The role of hormonal factors in OA is complex, with estrogen potentially influencing chondrocyte function through various pathways.²⁷ The protective association observed may be mediated by hormonal changes during pregnancy or by confounding factors linked to reproductive timing that were not fully captured by the genetic instruments. In contrast, a borderline association was noted between bilateral oophorectomy and reduced OA risk. Given its modest effect size and proximity to the FDR threshold, this finding should be interpreted cautiously. The underlying mechanisms remain uncertain and require further investigation.

Our research has significant implications for risk stratification, prevention, and management of osteoporosis and OA. First, reproductive characteristics such as later menarche (≥16 years) and early menopause (<45 years) could serve as valuable indicators for identifying high-risk women who may benefit from earlier BMD screening. Integrating these factors into existing risk assessment tools (eg., FRAX) could improve the precision of osteoporosis prediction.⁴⁴ For OA, reproductive history may help identify susceptible individuals, enabling targeted lifestyle or therapeutic interventions before structural joint damage occurs. Second, our results support personalized approaches to hormonal interventions. HRT may be considered for early postmenopausal women (aged 45–60), especially those with shorter reproductive lifespans, to mitigate accelerated bone loss. For OCP users, prolonged use or late discontinuation may adversely affect mid-life BMD, suggesting the need for alternative contraception strategies or adjunct bone-protective measures in select populations. Adolescents initiating OCPs shortly after menarche may require particular vigilance.⁴⁵

From the public health perspective, these findings highlight the need for updated clinical guidelines and awareness campaigns. Reproductive milestones should be incorporated into osteoporosis and OA management recommendations to emphasize their lifelong impact on skeletal health. Policymakers and researchers should prioritize studies in diverse populations to ensure equitable applicability, as well as investigations into hormonally modulated therapies for OA. These findings may inform future refinement of screening protocols and preventive strategies, ultimately reducing the global burden of bone diseases in women.

The major strengths of our study include the use of an MR design, which minimizes confounding and reverse causation, and the utilization of large-scale GWAS summary data, which provides substantial statistical power. The systematic assessment of a wide array of exposures and outcomes, coupled with extensive sensitivity analyses, enhances the robustness and reliability of our conclusions.

However, this study has several limitations. First, the GWAS data were predominantly from individuals of European ancestry, which may limit the generalizability of our findings to other populations. Second, although extensive sensitivity analyses did not detect significant pleiotropy for the main findings, residual or balanced pleiotropy inherent to the MR framework cannot be entirely excluded. Third, our analyses were based on summary-level data, which precluded stratified analyses (eg., by BMI or other lifestyle factors). Finally, partial sample overlap between UK Biobank-based exposure data and certain outcome GWAS, particularly OA (21.69%), whereas overlap for BMD was minimal (2.76%) and absent for osteoporosis. Such overlap may modestly inflate Type I error.

Future research should aim to replicate these findings in diverse ancestral populations. Furthermore, multivariable MR analyses could help disentangle the effects of correlated reproductive traits, and functional studies are needed to elucidate the molecular mechanisms underlying the observed associations.

Conclusion

In conclusion, this study provides genetic evidence supporting associations between female reproductive and hormonal factors and skeletal outcomes. Estrogen-related reproductive factors demonstrated relatively consistent protective associations with BMD and osteoporosis, whereas associations with OA were more heterogeneous and less well established, reflecting distinct biological pathways underlying bone remodeling and joint degeneration. A longer reproductive lifespan, characterized by earlier menarche and later menopause, was associated with more favorable bone outcomes. An unexpected association was observed between bilateral oophorectomy and reduced OA risk, warranting further mechanistic investigation. In contrast, several exposures, including selected parity-related traits and certain hormone-fracture associations, did not demonstrate significant associations after correction, suggesting that not all reproductive or hormonal factors exert measurable effects across skeletal outcomes and thereby informing both future mechanistic research priorities and clinical risk stratification. Findings supported by fewer genetic instruments or near the statistical threshold, such as the associations between bilateral oophorectomy and OA and between age at last OCP use and mid-life BMD, should be interpreted cautiously. Overall, these results contribute to improved risk stratification and hypothesis generation for future research; however, direct translation into clinical practice requires further validation. In addition, the timing of hormonal exposure, including HRT initiation and cessation, may represent a critical determinant of skeletal outcomes and warrants further investigation.