The effect of systemic iron status on osteoarthritis: A mendelian randomization study

Objective: To assess the causal effect of systemic iron status by using four biomarkers (serum iron; transferrin saturation; ferritin; total iron-binding capacity) on knee osteoarthritis (OA), hip OA, total knee replacement, and total hip replacement using 2-sample Mendelian randomization (MR) design. Methods: Three instrument sets were used to construct the genetic instruments for the iron status: Liberal instruments (variants associated with one of the iron biomarkers), sensitivity instruments (liberal instruments exclude variants associated with potential confounders), and conservative instruments (variants associated with all four iron biomarkers). Summary-level data for four OA phenotypes, including knee OA, hip OA, total knee replacement, and total hip replacement were obtained from the largest genome-wide meta-analysis with 826,690 individuals. Inverse-variance weighted based on the random-effect model as the main approach was conducted. Weighted median, MR-Egger, and Mendelian randomization pleiotropy residual sum and outlier methods were used as sensitivity MR approaches. Results: Based on liberal instruments, genetically predicted serum iron and transferrin saturation were significantly associated with hip OA and total hip replacement, but not with knee OA and total knee replacement. Statistical evidence of heterogeneity across the MR estimates indicated that mutation rs1800562 was the SNP significantly associated with hip OA in serum iron (odds ratio, OR = 1.48), transferrin saturation (OR = 1.57), ferritin (OR = 2.24), and total-iron binding capacity (OR = 0.79), and hip replacement in serum iron (OR = 1.45), transferrin saturation (OR = 1.25), ferritin (OR = 1.37), and total-iron binding capacity (OR = 0.80). Conclusion: Our study suggests that high iron status might be a causal factor of hip OA and total hip replacement where rs1800562 is the main contributor.


Introduction
Osteoarthritis (OA) is a common joint disease characterized by loss of articular cartilage, remodeling of the synovitis, and alterations of periarticular structures (Castañeda and Vicente, 2017). Any joint can develop OA, but symptoms linked to OA most commonly affect the knees, hips, hands, and feet (Katz et al., 2021). It is reported that persons with knee or hip OA have excess mortality compared with age-matched controls (Nüesch et al., 2011;Katz et al., 2021). Currently, there is no curative drug for OA. According to the estimation of the World Health Organization, there are about 300 million OA patients worldwide, and the prevalence of OA can reach up to 10%-20% due to the increasing obese and longevous population (GBD, 2017 Disease andInjury Incidence andPrevalence Collaborators, 2018). Given the high health economic burden, a better understanding of the risk factors associated with the occurrence of OA, especially modifiable factors, is needed.
Iron is an essential mineral to various biochemical processes, including DNA synthesis, ATP generation, and oxygen transport, where its absence or depletion can cause abnormal metabolization (Abbaspour et al., 2014). However, excess iron can also be toxic as it would be deposited into organs forming free radicals (Abbaspour et al., 2014). Therefore, disorders of iron homeostasis are involved in a wide scope of diseases (Bogdan et al., 2016). In humans, systemic iron status can be measured by clinical biomarkers: Serum iron, transferrin saturation, ferritin, and total iron-binding capacity (Bell et al., 2021). High serum iron, transferrin saturation, and ferritin signify high iron, while high total iron-binding capacity signifies low iron (Winter et al., 2014). Increasing evidence has found that systemic iron status is associated with OA. For example, high synovial iron was associated with a faster progression of OA in a murine model (Camacho et al., 2016). In vitro, iron overload was found to have detrimental effects on various joint components, such as synovium, cartilage, and subchondral bone, leading to synovial hyperplasia and inflammation, abnormal osteoblast function, and chondrocyte apoptosis (van Vulpen et al., 2015). Besides, some epidemiological evidence has shown the relationship between iron and OA risk. In a 2-year longitudinal study with 127 OA patients, radiographic findings showed that higher ferritin was significantly associated with narrower baseline joint space width and higher risks in the prediction of Kellgren-Lawrence grade severity (Kennish et al., 2014). Among 40 subjects with symptomatic knee OA, elevated serum ferritin levels were found to be associated with faster progression of cartilage damage assessed by arthroscopy (Nugzar et al., 2018). However, whether high iron status is a causal factor of OA is unclear because of the inherent defects of observational studies, such as residual confounding and reverse causality.
Mendelian randomization (MR) is a study design that can strengthen the causal inference on exposure-outcome relationships by minimizing the effect of confounding and excluding the potential reverse causality based on the use of genetic variants as instruments (Sheehan et al., 2008). The rationale of the causality assessment in MR is that genetic variants solely associated with the exposure are randomly assorted so that genetic effects on the outcome cannot be affected by potential confounders. Also, genetic variants are not modified by the occurrence or development of any diseases where reverse causality is impossible. For these reasons, MR represents robust indirect evidence of a causal relationship between exposure and outcome if any effect of the selected instruments on diseases is entirely mediated through the exposure (Sheehan et al., 2008). Therefore, we conducted a 2-sample MR (Lawlor, 2016) study to examine the causal effect of four clinical iron biomarkers on four OA phenotypes, including knee OA, hip OA, total knee replacement, and total hip replacement.

Materials and methods
A 2-sample MR study was conducted to investigate the potential causal relationship between systemic iron status and knee OA, hip OA, total knee replacement, and total hip replacement by using summary-level data. Systemic iron status was comprehensively represented by serum iron, transferrin saturation, ferritin, and total iron-binding capacity. Figure 1 showed the overview of the exposures, the outcomes, and the three assumptions of the genetic instruments for the MR design. Since this study is based on existing publications and public databases, both ethical approval and participant consent have been received by each relevant institutional review committee.

Genetic instruments for four clinical iron biomarkers
We constructed three sets of genetic instruments for the clinical iron biomarkers: Liberal instruments, sensitivity instruments, and conservative instruments. From the recent meta-analysis of three genome-wide association studies for iron homeostasis biomarkers: serum iron (N = 163,511), transferrin saturation (N = 131,471), ferritin (N = 246,139), and total iron-binding capacity (N = 135,430) (Bell et al., 2021), we selected single nucleotide polymorphisms (SNPs) at the genome-wide significance threshold (p < 5 × 10 −8 ) and in low linkage disequilibrium (r 2 < 0.01) as the potential instruments. After excluding five SNPs that are not in our outcome database because their minor allele frequency within the 1000G dataset is 0 (Supplementary Tables S1-S4), 14 SNPs for serum iron, 10 SNPs for transferrin saturation, 37 SNPs for ferritin, and 15 SNPs for total iron-binding capacity were left as the liberal instruments. To avoid the violation of assumption two of the instruments (Figure 1), we searched the Phenoscanner database (http://www.phenoscanner.medschl.cam.ac.uk/) to identify whether these liberal SNPs were associated with potential confounding factors. As interleukin 6, C-reactive protein, glycosylated hemoglobin, cholesterol, and body mass index are potential risk factors of OA (Livshits et al., 2009;Louati et al., 2015;Zheng and Chen, 2015;Farnaghi et al., 2017;Kozijn et al., 2019), we further excluded SNPs that were strongly associated with these phenotypes (p < 5 × 10 −8 ) leaving 12 SNPs for serum iron, nine SNPs for transferrin saturation, 27 SNPs for ferritin, and 14 SNPs for total iron-binding capacity as the sensitivity instruments (Supplementary Tables S1-S4). We did not exclude rs855791 from TMPRSS6, rs1799945 from HFE, and rs1800562 from HFE here since they were missense mutations of TMPRSS6 or HFE. While HFE and TMPRSS6 involve the signaling cascade of iron homeostasis directly Frontiers in Genetics frontiersin.org (Bell et al., 2021), these three SNPs may have a vertical effect from iron to the potential confounding factors, namely glycosylated hemoglobin, low-density lipoprotein, and total cholesterol mentioned in the Supplementary Tables. Thus, these three SNPs may have vertical pleiotropy with the possible confounding factors which does not violate assumption two of MR (Davies et al., 2018) ( Figure 1). Similar instrument selection criteria are applied to evade excluding vertical SNPs (Chong et al., 2022). Conservative instruments were constructed where only SNPs associated with all four clinical iron biomarkers were eligible; therefore, rs1800562, rs1799945, rs855791, and rs57659670 were used as conservative instruments. Because the MR estimates were for each iron biomarker, we set the liberal instruments that have greater power as the main instrument. The strength of the instruments was evaluated by the F statistic using the formula F = R 2 (N-2)/(1-R 2 ) (Palmer et al., 2012), where R 2 is the proportion of the variance of the four clinical iron biomarkers explained by the genetic variant and N is the sample size of the gene-each of the iron biomarker association. Thus, the computed F statistics range from 5136 to 6388 for serum iron, from 7012 to 8489 for transferrin saturation, from 1644 to 5844 for ferritin, and from 4206 to 9724 for total iron-binding capacity which was much greater than 10 demonstrating sufficient statistical strength (Burgess et al., 2013).

Genetic associations with four OA phenotypes
To obtain the associations of the genetic instruments with OA, summary-level data were extracted from the largest genome-wide meta-analysis to date across 826,690 individuals with 177,517 cases and 649,173 controls. The mean age (standard deviation) is 62.4 (11.9) for the cases and 52.4 (17.4) for the controls. The percentage of females is 62% for cases and 52% for controls. The meta-analysis was a combination of data from 13 cohorts with a population of >97% European descent (Boer et al., 2021).
In this study, data from knee OA (n = 62,497), hip OA (n = 36,445), total knee replacement (n = 18,200), total hip replacement (n = 23,021), and a max of healthy controls (n = 333,557) were used (Boer et al., 2021). Knee (or hip) OA was defined as OA or joint replacement at the knee joint (or hip joint). Total knee (or hip) replacement was defined as having undergone total knee (or hip) replacement due to OA in the knee joint (or hip joint). All the definitions of the cases were self-reported, clinically diagnosed, ICD10, codes, or radiographic depending on the data available in the cohort. Controls were osteoarthritis-free or population-based with or without ICD code exclusions. Detailed association estimates for the genetic instruments with the four OA phenotypes were provided in Supplementary Tables S1-S4.

MR estimates
To derive MR estimates, the inverse-variance weighted method (IVW) was used as the main approach. This method first assesses the effect of each SNP on the outcome by calculating the Wald ratio and then uses the inverse variance of SNPs as weights to obtain a combined causal effect (Lee et al., 2016). Random-effects model was selected for this method because the fixed-effects models may be at risk of yielding artificially precise estimates in the presence of heterogeneity . Furthermore, heterogeneity across the instrumental SNPs was evaluated using I 2 statistics and Cochran's Q test (Higgins et al., 2003). For each SNP, the MR estimates were obtained from the Wald ratio method with standard error derived using the Delta method (Thompson et al., 2016). The MR estimates are presented using odds ratios which are scaled to one standard deviation increment of genetically predicted clinical iron biomarkers.
To investigate the robustness of the findings and assess the possible horizontal pleiotropy, we further performed the weighted

FIGURE 1
An overview of the Mendelian randomization design. Assumption 1: The instrumental variables must be strongly associated with the exposures; Assumption 2: The instrumental variables must be independent of the potential confounders of the association between the exposure and outcome; Assumption 3: The instrumental variables should not be associated with the outcomes directly.
Frontiers in Genetics frontiersin.org median, MR-Egger, and Mendelian randomization pleiotropy residual sum and outlier (MR-PRESSO) methods as sensitivity MR approaches. The weighted median method first orders the ratio estimate of each instrument SNP on the outcome by its magnitude of weight and then produces an overall MR estimate based on the median value with standard error derived using the parametric bootstrap method (Bowden et al., 2016). The MR-Egger method can detect the overall horizontal pleiotropy of genetic instruments by examining whether the intercept of the association between exposure and outcome differs from zero . However, the estimates of MR-Egger are generally small especially when the number of instruments is small , and thus, we did not use MR-Egger when using conservative instruments. MR-PRESSO identifies possible pleiotropic outliers and generates estimates after the outliers are removed (Verbanck et al., 2018). Moreover, a leave-one-out analysis was computed where we excluded one SNP at a time and conducted IVW on the rest SNPs to test whether any single variant was driving the causal association between exposure and outcome . Because four OA outcomes and four biomarkers were involved in the analyses, we set a Bonferroni corrected p-value of 0.003 (0.05/ 16) as the statistical significance level. All statistical data analyses were conducted with R software. The MendelianRandomization and the forestplot packages were used to facilitate the MR analyses and display the results (Yavorska and Burgess, 2017).

Repeated analysis in the second-largest meta-analysis of OA GWAS
We also repeated our analyses in the second-largest metaanalysis of OA GWAS (Tachmazidou et al., 2019). This metaanalysis was a combination of data from UK Biobank and Arthritis Research UK Osteoarthritis Genetics (arcOGEN) with 455,221 individuals of European descent only (Tachmazidou et al., 2019). Due to the large number of OA participants in these two cohorts, there are some overlapping participants between the largest and the second-largest meta-analysis of OA GWAS. As the summary data of total knee replacement and total hip replacement were not available in this meta-analysis, we only repeated our analyses in the knee OA (n = 24,955) and hip OA (n = 15,704) with healthy controls (n = 378,169) (Tachmazidou et al., 2019). The OA definitions were from self-report and hospital records in UK Biobank and were from total joint replacement records or radiographic evidence of disease in arcOGEN.

Results
Tables 1-4 show the associations of the four systemic iron status biomarkers and knee OA, hip OA, knee replacement, and hip replacement using the three instrument sets. In the IVW analyses using liberal instruments, genetically predicted serum iron and transferrin saturation were found to be significantly associated with hip OA and hip replacement, and the other associations were not evident (Tables 1-4; Figure 2). The odds ratios of hip OA were 1.18 (95% CI: 1.07-1.30) per one standard deviation increment in serum iron and 1.15 (95% CI: 1.07-1.24) in transferrin saturation (Tables 1-4; Figure 2).
Similar results were seen in analyses using sensitivity instruments. In the analyses using conservative instruments, genetically predicted total iron-binding capacity was negatively associated with hip OA and hip replacement with ORs ranging from 0.80 to 0.82 (Tables 2, 4). However, statistical evidence of heterogeneity across the MR estimates was found especially for hip OA and hip replacement (Tables 2, 4; Figure 2). For example, in the MR estimates of ferritin on hip OA and hip replacement, the I 2 of the heterogeneity was 64.8%-91.1%. When we used the Wald ratio method to assess the MR estimates of four individual SNPs (rs57659670, rs855791, rs1799945, and rs1800562) with hip OA and hip replacement, rs1800562 was the SNP that was significantly associated with hip OA in serum iron (OR = 1.48; 95% CI: 1.29-1.68), transferrin saturation (OR = 1.57; 95% CI: 1.17-1.37), ferritin (OR = 2.24; 95% CI: 1.71-2.94), and total-iron binding capacity (OR = 0.79; 95% CI: 0.73-0.86), and hip replacement in serum iron (OR = 1.45; 95% CI: 1.25-1.69), transferrin saturation (OR = 1.25; 95% CI: 1.14-1.37), ferritin (OR = 1.37; 95% CI: 1.58-2.99), and total-iron binding capacity (OR = 0.80; 95% CI: 0.73-0.88). The other three SNPs did not show evident associations with hip OA or hip replacement (data not shown). Also, the leave-one-out analyses demonstrated that the causal effect of systemic iron status on hip OA and total hip replacement was mainly resulted from rs1800562 (Supplementary Table S5). We repeated the MR estimates of the four systemic iron status biomarkers on knee OA and hip OA in the second-largest meta-analysis of OA GWAS, and similar results were shown in Supplementary Table S6. That is, serum iron, transferrin saturation, and total iron-binding capacity were associated with hip OA, and the other associations were not evident.

Discussion
The present 2-sample MR study used the largest public dataset to date to comprehensively evaluate the causal role of systemic iron status for knee OA, hip OA, total knee replacement, and total hip replacement. In a pattern concordant with an effect on systemic iron status, we found that genetically predicted clinical iron biomarkers were significantly associated with hip OA and total hip replacement. Given the heterogeneity across the genetic instrumental variables, these associations were mainly resulted from rs1800562, a missense mutation of the HFE gene.
Our findings are in line with previous conclusions from observational and experimental studies that iron is involved in the development of OA. In an observational study, synovial fluid iron concentrations determined by the colorimetric method were significantly higher in OA patients compared to the controls (Yazar et al., 2005). Another 2-year follow-up study reported that higher iron status was associated with more severe radiographic progression in knee OA patients (Kennish et al., 2014). In murine models, cellular iron accumulation in the knee joint resulted from systemic iron overload could induce the early onset of OA and accelerate OA progression via compromising chondrocyte metabolism and over-expressing local inflammatory mediators (Camacho et al., 2016;Simão et al., 2019;Burton et al., 2020). Several studies using the MR method also assessed the Frontiers in Genetics frontiersin.org    instrumented transferrin saturation as the iron load, Pilling et al. (2019) found that homozygotes of rs1800562 or compound heterozygotes of rs1800562 & rs1799945 were positively associated with incident OA in both men and women. Using rs800562, rs1799945, and rs855791 as the instruments of serum iron, no significant association was observed between iron levels and overall OA in UK Biobank (Zhou et al., 2020;Zhou et al., 2021), but per standard deviation increment in iron was associated with

FIGURE 2
Forest plots of Mendelian randomization estimates for the association between genetically predicted systemic iron biomarkers with liberal instruments and knee osteoarthritis, hip osteoarthritis, total knee replacement, and total hip replacement. OR: Odds ratio; CI: Compatibility/confidence interval; IVW: Inverse-variance weighted. Points indicate the odds ratio per standard deviation increment increase in the iron biomarkers. Error bars indicate 95% compatibility/confidence intervals. (A) for knee osteoarthritis, (B) for hip osteoarthritis, (C) for total knee replacement, and (D) for total hip replacement.
Frontiers in Genetics frontiersin.org increased risk of OA in males (Zhou et al., 2021). Using rs800562, rs1799945, rs855791, and rs8177240 as the instruments of nutritional iron and summary data from the UK Biobank and arcOGEN cohorts, no significant association was observed between the iron and knee OA, hip OA, or overall OA, but the positive causal effect of iron levels on overall OA was observed in females (Qu et al., 2020). Using the instruments from the Genetics of Iron Status Consortium and the summary-level data of outcomes from the UK Biobank and arcOGEN cohorts, Xu et al. (2022) found that transferrin saturation was positively associated with both knee OA and hip OA, transferrin was negatively associated with hip OA, but serum iron and ferritin did not show a prominent effect on OA outcomes. In our study, we assessed the iron status using four biomarkers with three instrument sets while the previous studies used one biomarker with three or four instruments. Our study used the largest meta-analysis data for OA while the others included the relatively smaller size of OA patients. The inconsistencies between the previous studies with the present study can be partly explained by the different biomarkers with non-specific instruments and underpowered samples. Thus, the conclusion from our study may be more representative of the association between iron status and OA. Our findings only suggest causal associations among hip OA and total hip replacement, but not knee OA or total knee replacement. This could be explained by the fact that genetic heritability contributes more to hip OA than knee OA [14.7% for knee OA and 51.9% for hip OA (Tachmazidou et al., 2019)], and the exposures in the MR study are genetically predicted. In the main analyses and the repeated analyses, we only found causal associations between hip OA & total hip replacement and serum iron, transferrin saturation, & total iron-binding capacity, but not ferritin. This can be partly explained by the heterogeneity among the variants in the analyses (I 2 = 91% in hip OA and I 2 = 83% in the total hip replacement). Another possible explanation for this is that ferritin is a representative biomarker of body iron stores in a non-inflammatory state (Adams, 2015;Bell et al., 2021), while OA is regarded as an inflammatory disease (Dainese et al., 2022), and thus, there might be a weak relationship between ferritin and OA. A potential concern for MR conclusions relates to the horizontal pleiotropy of the instruments. Our MR study does not have a such concern, because the findings were robust to analyses using sensitivity instruments that exclude variants associated with potential confounders. Though rs1800562 from the HFE gene is associated with potential confounders glycosylated hemoglobin, low-density lipoprotein, and total cholesterol, we did not exclude this variant in our analyses. This is because the HFE protein encoded by the HFE gene is the central regulator of systemic iron homeostasis and has no other well-established roles (Barton et al., 2015), any effect of the second phenotype may probably be acting downstream of iron status, rather than independent of it. Therefore, rs1800562 may have a vertical pleiotropy with the outcome which does not violate the assumption of MR (Davies et al., 2018). However, there is still a possibility that rs1800562 has other horizontal effects on the outcomes independent of iron.
Besides the causal evidence of the MR relationships, the biological plausibility adds further support for the possibility of causality. Hereditary hemochromatosis is a common disease in Europeans characterized by iron overload in multiple organs where the mutation of rs1800562 is the main contributor (Husar-Memmer et al., 2014). It is reported that people with hemochromatosis are at increased risk of OA and OA is one of the most frequent complications in people with hemochromatosis (Elmberg et al., 2013;Husar-Memmer et al., 2014). The missense mutation rs1800562 from the HFE gene affects the synthesis of HFE protein which is associated with the regulation of the circulating iron by regulating the interaction of the transferrin receptor with transferrin. Studies showed that excess iron could act as a catalyst to produce a large number of reactive oxygen species which are deleterious agents involved in OA cartilage degradation (Henrotin et al., 2003;Yang et al., 2018). Also, evidence showed that iron homeostasis dysregulation could lead to the secretion of pro-inflammatory cytokines which promote chondrocytes apoptosis and subsequent OA cartilage degeneration (Simão et al., 2019;Jing et al., 2021).
Since high iron status is a treatable condition, our findings have important clinical and public health implications. It is reported that phlebotomy, oral chelation, and dietary changes are options for people with iron overload (Palmer et al., 2018). There may be a benefit to the careful reduction of iron indices for reducing OA risk in people with high iron status, especially among hemochromatosis patients or people with the mutation rs1800562. However, the interpretation of the MR finding that iron decrement might be a benefit for OA patients requires justification, because MRs are the effects of lifelong exposures while interventions consider short-term pharmacological treatments. Taken together, iron decrement could be an attractive new target for hip OA treatment which needs to be validated in randomized controlled trials in the future.
The major strength of the present study is the 2-sample MR design, which minimizes the confounding and reverse causality seen in observational studies. We comprehensively investigated the causation of the systemic iron status in four representative biomarkers to knee OA, hip OA, total knee replacement, and total hip replacement using the largest summary-level data. Also, we used three instrument sets to obtain robust conclusions. Furthermore, we repeated our analyses in the second largest OA GWAS dataset with European descent, thereby diminishing bias from population stratification. However, we must acknowledge some potential limitations. First, rs1800562 was the only SNP associated with OA which needs to be verified in another dataset independent from the UK Biobank study and the arcOGEN consortium. Second, the findings of our study relate to patterns of iron status largely within the normal range, and thus, cannot be used to make inferences on the effect of abnormally high or low serum iron levels. Third, our study does not offer insight into whether the estimates are equally applicable to both men and women. Despite these limitations, the results of this work show consistent and biologically plausible effects.
In summary, the present MR findings suggest that high iron status might be a causal factor of hip OA and total hip replacement where rs1800562 is the main contributor. Given the modifiable nature of the iron status, further clinical trials are warranted to validate the therapeutic role of iron decrement in people with hip OA, especially among those with the mutation rs1800562.
Frontiers in Genetics frontiersin.org

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/yzhang666666/ IronMR_23022023.

Ethics statement
Since this study is based on existing publications and public databases, both ethical approval and participant consent have been received by each relevant institutional review committee.