Osteoarthritis (OA) is the most common form of joint disease around the world, and it is estimated that approximately 302 million individuals suffer from OA (1). Numerous studies have shown that OA has been ranked as the leading cause of disability and accounts for a heavy burden of disease (2, 3). As the aging population and obesity pandemic increase, OA is more prevalent and has become a major public health concern (4). However, the exact mechanism underlying the pathogenesis of OA has not been fully clarified.

Epidemiological studies indicate a distinct sexual disparity in the incidence of OA in which females develop the disease more frequently and more severely than males, particularly after menopausal age (5). Individuals who had undergone postmenopausal hormone therapy exhibited a higher risk of developing OA than those who did not, revealing that hormonal reproductive factors may play a critical role in the initiation and progression of OA. There is a wealth of data concerning the associations of hormonal reproductive factors, including age at menarche (AAM), age at nature menopause (ANM) and age at first birth (AFB), with incident OA. For example, findings from a prospective cohort study suggested that an early AAM increased the risk of hip and knee replacement for OA, while the ANM was not associated with the risk (6). Another cross-sectional study demonstrated that women with ANM < 45 years were associated with a 2.60-fold risk of developing OA compared with those with an ANM ≥ 45 years (7). These discrepancies may be due to potential confounding factors (such as obesity and age) and reverse causality. Reproductive behavior is shaped by biology and environment, while AFB represents an accurate measure of complex reproductive outcomes, are frequently recorded and consistently measured. It is indicated that the heritability of AFB shifted from 9% for women born in 1940 to 22% in 1965 (8).

Mendelian randomization (MR) is an epidemiological method that utilizes genetic variants robustly associated with exposure as instrumental variables (IVs) to estimate the causal effect of exposure on an outcome (9, 10). Because of the bias from confounders, reverse causation and measurement error, even there was a significant statistical association between the exposure and outcome, traditional analyses have limitations for the assessment of causality, while MR offers an alternative way to probe it. It is considered important to target the management of relevant health problems in OA susceptible individuals, which requires us to first clarify the root causes of OA. In this study, confounders including environmental factors and BMI could be excluded due to the application of MR. To the best of our knowledge, the MR study focusing on this topic do not exist. Therefore, it is necessary to do the MR design to investigate it.

In total, after implementing the pruning strategy previously described, there were 349 SNPs achieved genome-wide significance for AAM (11) and 121, 54 and 10 IVs for AAM (12), ANM and AFB, respectively. F-statistic values for individual instrumental SNPs were all above the threshold 10, with means of 64.27, 58.00, 68.14 and 36.49 for AAM (11), AAM (12), ANM and AFB, respectively (Table 1). SNPs were excluded or substituted with highly correlated (r² > 0.8) proxy SNPs due to unavailability in outcome datasets or palindromic with ambiguous A/T or G/C. Detailed information for incorporated instrumental SNPs is presented in Supplementary Tables 2–4.

The MR analysis indicated that genetically determined each additional year in AAM (11) was associated with a decreasing risk of knee OA after correcting for multiple testing (IVW: OR = 0.86, 95% CI: 0.76–0.95, P = 1.58 × 10⁻³). The causality between AAM (11) and overall OA (P = 5.95 × 10⁻³) and hip and/or knee OA (P = 3.64 × 10⁻²) was only normally significantly positive, with IVW OR_{per − SDincrement(95%CI)} of 0.91 (0.85–0.98) and 0.92 (0.84–1.00), respectively. To the best of our knowledge, such causality was not observed in another AAM (12) exposures, with P values above the threshold and OR of 0.94, 1.16, 1.02 and 1.02 for OA and three subtypes, respectively. The available evidence also made it difficult to explain the causality of ANM to OA and subtypes. Furthermore, the IVW method indicated that AFB exerted a causal effect on OA and all subtypes in an unfavorable manner, with OR of 0.80 and 0.79 for overall (95% CI: 0.68–0.92, P = 1.80 × 10⁻⁴) and hip and/or knee OA (95% CI: 0.64–0.93, P = 1.33 × 10⁻³), respectively. However, as shown in Table 1 and Figure 1, some nominally significant positive causality correlations with OA did not pass multiple-testing correction for the Bonferroni method (hip OA: OR = 0.76, 95% CI: 0.51–1.00, P = 2.65 × 10⁻²; knee OA: OR = 0.81, 95% CI: 0.63–0.98, P = 1.57 × 10⁻²). Additional methods, including weighted median, weighted mode and MR-Egger section, validate the uniformity conclusion (Supplementary Table 5).

To effectively control pleiotropy, we next investigated the MR-PRESSO and p value for the MR-Egger intercept test. The sensitivity analysis revealed the absence of outliers (Supplementary Table 7) and horizontal pleiotropy (Supplementary Figure 2). As shown in Supplementary Figure 2, symmetry in funnel plots did not show any evidence of publication bias. In addition, the leave-one-out analysis showed that none of the individual genetic markers drove the majority of the association signal. The scatter plots for effect sizes of SNPs for three hormonal related exposures and those for OA and subtypes are shown in Supplementary Figure 1.

BMI is a known modifiable risk factor that plays an important role in the etiology of OA and shapes reproductive exposures. A total of 82 SNPs were found to be associated with BMI (P < 5 × 10⁻⁸), including 51 in AAM (11), 29 in AAM (12) and 2 in AFB. SNPs associated with BMI have been annotated with ^* in the Supplementary Tables 2–4. As shown in Table 1 and Figure 2, no casualty between AAM (11) no BMI and keen OA were found when we performed MR again (IVW: OR = 1.01, 95% CI: 0.88–1.16, P = 0.89), which indicates that BMI-associated SNPs confounded causality in our initial calculation. The negative result was confirmed in AAM (12) no BMI. Similarly, after deleting two BMI-associated SNPs, there was not enough evidence to indicate the causal relationship between AFB no BMI and overall OA (IVW: OR = 0.83, 95% CI: 0.72–0.94, P = 4.91 × 10⁻³). However, a similar causal effect association between AFB and hip and/or knee OA could still be measured (IVW: OR = 0.78, 95% CI: 0.66–0.92, P = 3.22 × 10⁻³). Additional methods section validate the uniformity conclusion (Supplementary Table 6).

In the current study, we performed two-sample MR to investigate the causality between three hormonal reproductive factors and the risk of OA. The concern was that 1 year later in AFB was associated with a reduced risk of OA. These findings support the hypothesis that AFB may play a causal role in the pathway of developing OA. The adverse causal effects were robust in our MR after pruning potential confounder BMI.

Several observational studies confirm the causal association of hormonal reproductive factors with OA, while current results on this topic from conventional epidemiological studies remain controversial. A large cohort study with 30.727 cases is consistent with our conclusions, affirming the role of AFB in OA etiology (15). Other studies also point toward some positive association, and a prospective cohort study containing over 30,000 women found that older age at menarche was associated with a decreased risk of total knee replacement (TKR) due to primary OA (28). Researchers proposed that one possible explanation could be that lower AAM may be a marker of other factors, such as higher BMI when young (6). High BMI is known as a risk factor for OA, but it is unlikely that BMI would explain the casualty found here, as we adjusted our analyses for BMI, and our findings were consistently observed within subgroups of no BMI. However, the assumed relationship between the female hormonal aspects and OA was not clinically significant in another cohort study (29). Since hormonal reproductive factors are prone to bias due to interference from potential confounders which difficult to be excluded by traditional epidemiology, MR estimates reflect the causality at the genetic level.

Given the complexity of these confounders, the underlying mechanism of hormonal reproductive factors in the development of OA remains to be elucidated. Estrogen is considered to strongly associate with the female hormonal reproductive cycle, in which receptors are found on bone and chondrocyte cells (30). Several studies show evidence of the associations between radiographic changes in OA and high bone density since considering estrogen could prevent bone loss (31, 32). Consequently, greater exposure to estrogen, while preventing bone loss, may plausibly promote OA. Gao et al. studied estrogen and estrogen metabolites in Chinese women with OA. Compared to the controls (healthy and rheumatoid arthritis women), premenopausal women with OA had a significantly lower concentration of 2-hydroxyestrone and free estrogen in serum. In postmenopausal women, the serum concentration of 2-hydroxyestradiol was increased compared to that in controls, while free and total estrogen were significantly decreased. Apart from estrogen deficiency, rapidly elevated serum levels of 2-hydroxyestrone in the perimenopausal period may correlate with the pathogenesis of OA (33, 34). Furthermore, it has been reported that estrogen may have different effects on the initiation and progression of OA. Hence, it is difficult to ratiocinate the biological mechanisms that underlie this study due to these heterogeneous effects. However, the effects of female hormones on OA can be further explored in animal models and in vitro studies.

Our study has several strengths. The large sample size and richness of the data set for reproductive variables of interest led the estimated effects to be close to the truth. In addition, three different reproductive traits (AAM, ANM, and AFB) were incorporated to reflect the length of the reproductive period and complementing each other well. To reduce the interference of potential factors, we examined OA directly rather than proxies of OA, such as hospitalization or joint replacement. Moreover, we were able to adjust for not only hormonal reproductive factors but also reported confounders of OA, such as BMI. Notwithstanding, we must acknowledge several limitations. Firstly, there was no stratification of sex in the existing GWAS data set, while our selection of hormonal reproductive factors was female specific. However, since most cases in the GWAS dataset were from females (63.7% in overall OA), we thought the estimated effects would be close to the truth. Secondly, to diminish population stratification, our samples were restricted to the European population, which leads our findings to be applicable for European populations. Finally, the design of our study precluded us from considering other factors, such as environmental effects and hormone use, in addition to the only confounder BMI regarded in the current study.

In summary, our findings add to a growing body of evidence surrounding the unfavorable effects of early age at first birth on OA risk, suggesting the essential for relevant health problem management in susceptible populations. Further large-scale studies or longitudinal studies are required to validate our findings.

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

Each participating study obtained written informed consent from all participants and received approval from the appropriate local institutional review boards.

ZY was responsible for the conception, design of the study, full access to all the data in the study, took responsibility for the integrity of the data, and the accuracy of the data analysis. BW and JW Performed the statistical analysis and drafted the manuscript. All authors were involved in drafting the article or revising it critically for important intellectual content, and approved the final version to be published.

This study was supported by the University Natural Foundation of Anhui Educational Committee (KJ2020A0199) and Nature Science Foundation of Anhui Medical University (2021xkj156).

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

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.