Impact of genetically predicted characterization of mitochondrial DNA quantity and quality on osteoarthritis

Background: Existing studies have indicated that mitochondrial dysfunction may contribute to osteoarthritis (OA) development. However, the causal association between mitochondrial DNA (mtDNA) characterization and OA has not been extensively explored. Methods: Two-sample Mendelian randomization was performed to calculate the impact of mitochondrial genomic variations on overall OA as well as site-specific OA, with multiple analytical methods inverse variance weighted (IVW), weighted median (WM), MR-Egger and MR-robust adjusted profile score (MR-RAPS). Results: Genetically determined mitochondrial heteroplasmy (MtHz) and mtDNA abundance were not causally associated with overall OA. In site-specific OA analyses, the causal effect of mtDNA abundance on other OA sites, including hip, knee, thumb, hand, and finger, had not been discovered. There was a suggestively protective effect of MtHz on knee OA IVW OR = 0.632, 95% CI: 0.425–0.939, p-value = 0.023. No causal association between MtHz and other different OA phenotypes was found. Conclusion: MtHz shows potential to be a novel therapeutic target and biomarker on knee OA development. However, the variation of mtDNA abundance was measured from leukocyte in blood and the levels of MtHz were from saliva samples rather than cartilage or synovial tissues. Genotyping samples from synovial and cartilage can be a focus to further exploration.


Introduction
Osteoarthritis (OA) is a common joint disease with degenerative changes in articular cartilage and remains major impact on public health worldwide (Nelson, 2018). According to the estimation of the Global Burden of Disease Study (GBD), nearly 300 million individuals throughout the world are suffering from OA, and the number is still increasing. OA is one of the leading causes of pain, disability, and great socioeconomic burden in developed countries (Glyn-Jones et al., 2015). However, the mechanism of OA remains unclear. Identifying risk factors such as a disorder-related biomarker is essential to understand OA pathogenesis, decrease the incidence of the disease and develop efficiency prevention strategies.
Mitochondria are complex multifunctional organelles involved in various cell functions, e.g., heat regulation, calcium homeostasis, reactive oxygen species (ROS) production, and proinflammatory cytokines production, the damage of which may affect chondrocyte health at an extent (Blanco et al., 2011;Fernandez-Moreno et al., 2017;Blanco et al., 2018). They have their own genome (mitochondrial DNA, mtDNA) which contain 37 genes and encode 13 proteins, 22 transfer RNAs (tRNAs) and two ribosomal RNAs (rRNAs) (Anderson et al., 1981). As multicopy genome, sequence mutation and copy number variation of mtDNA are prevalent in population (McArdle et al., 2013;McGuire, 2019). As one of the mtDNA quantity characteristics, mtDNA abundance is recognized as a rough estimation for the number of mitochondria (Castellani et al., 2020). Previous studies provided supportive evidence for extensive pathogenicity of mtDNA abundance (Frahm et al., 2005;Hu et al., 2016). As one of the mtDNA quality characteristics, mitochondrial herteroplasmy (MtHz), where mtDNA with distinct sequences coexist, has been found in a large spectrum of human disease, including classical mitochondrial diseases and complex disorders (Ye et al., 2014). MtDNA abundance and MtHz can integrate many aspects of mitochondrial function and serve as promising biomarkers in probing interactions between mitochondria and disorders (Castellani et al., 2020;Tian et al., 2021).
In a recent review, OA was recognized as a potential mitochondrial disease considering the impact of mitochondrial dysfunction on cartilage degradation (Fernandez-Moreno et al., 2022). However, the causality between mitochondria and OA has rarely been investigated. Therefore, the goal of our study was to probe causal effects of characterization of mtDNA quantity and quality on OA development which is crucial to understand the role of mitochondria in OA etiology. A two-sample Mendelian randomization (MR) analysis was performed to investigate potential causal associations between mtDNA abundance and MtHz with OA using summary statistics from large-scale genome-wide association study (GWAS). MR-Steiger test was applied to ascertain whether variation of mtDNA characterization is a cause or a consequence of OA development.

Study design overview
The design of our research is displayed in Figure 1. We adopted a two-sample MR design to compute the causal effect of the characterization of mtDNA quantity and quality on overall OA and site-specific OA separately (Boer et al., 2021). Genetic association estimates for mtDNA abundance were derived from the United Kingdom biobank study (N = 291,950) (Hagg et al., 2021). Single Nucleotide Polymorphisms (SNPs) associated with MtHz were obtained from the 23andME research program (N = 982,072), a personal genomics and biotechnology company (Nandakumar et al., 2021). The participants included in both GWASs were of European ancestry.

Selection of genetic instruments
SNPs serving as instrumental variables (IVs) for MtHz and mtDNA abundance, all reach genome-wide significance (p < 5 × 10-8) (Hagg et al., 2021;Nandakumar et al., 2021). Twenty SNPs were related to MtHz and accounted for 32% observed SNPheritability, 64 SNPs were included to estimate the genetic liability of mtDNA abundance and explained approximately 8.3% SNP-heritability. We performed linkage disequilibrium (LD) pruning in PLINK with 1,000 Genomes Europeans as the reference panel to ascertain whether these genetic variants are independent of each other (r 2 < 0.001) . After LD test, five genetic variants linked to MtHz were removed, and 17 IVs related to mtDNA abundance were excluded. These SNPs were then matched in summary GWASs of OA in subsequent analysis and those not available in the outcome GWAS were removed or replaced by the proxy SNPs. F-statistic was calculated to filter the weak instruments. The threshold of F-statistic that was sufficient for identifying causal effect was 10 (Burgess et al., 2011). The selected genetic variants are demonstrated in Supplementary Table S1, S2.

Osteoarthritis and data sources
Summary statistics data on overall OA and its specific sites, comprising knee, hip, spine, thumb, and hand OA, were obtained from the latest publicly available GWAS of OA. This GWAS covered nine populations with up to 826,690 participants (177,517 OA patients) (Boer et al., 2021). The definition of OA satisfied the criteria of Genetics of Osteoarthritis (GO) including self-reported status, hospital diagnosed, ICD10 codes or radiographic as defined by the TREAT-OA consortium (Boer et al., 2021). All studies contributing data to our analyses were approved by the relevant ethics committees, and all study participants in these studies provided written, informed consent.

Statistical analysis
Two-sample MR analysis was executed in R software (R version 4.3.0) with the TwoSampleMR (version 0.5.6), MR-PRESSO (version 1.0) and mr. raps (version 0.2) R packages (Hemani et al., 2018;Verbanck et al., 2018). The methods applied in our research are presented in Figure 2. MR-steiger test was applied to infer causal direction between traits under investigation (Hemani et al., 2017). Four methodologies including inverse variance weighted (IVW), weighted median (WM), MR-Egger and MRrobust adjusted profile score (MR-RAPS) were employed to estimate causality between mtDNA abundance and MtHz and OA. IVW was taken as a primary MR analysis method in our research, which is based on the hypothesis that all selected genetic instruments are valid and give an overall causal estimate strengthening causal inference (Lawlor et al., 2008;Burgess et al., 2013). MR-Egger regression takes presence of directional pleiotropy into consideration and measures horizontal pleiotropy with regression intercept, whereas the results of this method are susceptible to outlying genetic variants (Bowden et al., 2015). Weighted Median based on hypothesis that at least 50% of the variants are valid, improves power of causal effect detection but reduces precision (Bowden et al., 2016). MR-RAPS applies robust adjusted profile scores to correct for pleiotropy and makes our results more reliable (Hartwig et al., 2017).
Several sensitivity analyses were applied to detect and correct for heterogeneity and pleiotropy. MR-PRESSO was conducted to detect the existence of horizontal pleiotropy and correct the causal estimate affected by possible pleiotropic outliers (Verbanck et al., 2018). Corresponding p-values were derived based on 1,000 simulations (Verbanck et al., 2018). And the estimation of MR Egger regression intercept was also employed to reflect presence of pleiotropy (p < 0.05) . IVW method was utilized to investigate heterogeneity. The level of heterogeneity was quantified by Cochran Q statistics (Carnegie et al., 2020). Leave-one-out sensitivity analysis was performed to identify possibly influential SNPs, which repeated MR analysis with each SNP excluded in turn (Carnegie et al., 2020). As five separate outcomes were tested in our study, main results had statistical significance at p-value<0.01 (0.05/ 5) after Bonferroni correction.

Causal effect of mitochondrial heteroplasmy on OA
Fifteen LD-independent genetic variants were taken for repeated MR analysis (Supplementary Table S1). We extracted data of abovementioned SNPs from summary GWAS of outcome traits (overall OA and specific-site OA) and one SNP (rs2286639) was removed in all outcomes except finger OA due to the effect of non-concordant alleles (e.g., A/G vs A/C). In finger OA, two SNP (rs2286639, rs3702096) were excluded for absence of proxy SNPs on the online platform SNiPA (https://snipa.helmholtzmuenchen.de/ snipa3/). The mean F-statistic was 36.357 which was above the threshold F value of 10. In total, there were 14 IVs for overall, knee, hip, thumb, and hand OA and 13 IVs for finger OA.
Results of the casual association between MtHz and OA are summarized in Table 1 Table S3). The results of leave-oneout sensitivity analysis and forest plots demonstrated that our study in genetically prediction was robust (Supplementary Figure S1-S3).

Causal effect of mitochondrial abundance on OA
After LD test, 47 independent genetic variants were chosen for twosample MR analysis. Several genetic variants could not be matched in summary statistic GWASs of knee (rs35734242), finger (rs1065853), and hand OA (rs12451555). We searched their proxy SNPs on SNiPA and those whose proxy SNPs were absent on SNiPA were excluded. The mean F-statistic for IVs of mtDNA abundance was 93.380, which was above the threshold 10 (Supplementary Table S2). Finally, the SNPS were selected as IVs for thumb, hip, overall OA was 47, and 46 SNPS for hand, finger, knee. Regarding overall OA, MR-PRESSO identified two outliers s (rs59488041, rs12924138) and therefore removed them for repeated MR analysis. In the presence of heterogeneity (Cochran's Q: p= 0.033) and absence of pleiotropy (Egger-intercept: p= 0.687 > 0.05), the results of WM analysis on overall OA were nominally significant (OR = 1.108, 95% CI: 1.001-1,226, p= 0.048, Table 2). IVW method in random effects model was utilized to correct results potentially impacted by heterogeneity, and the results suggested that genetically elevated mtDNA abundance was not casually associated with overall OA (OR = 1.040, 95% CI: 0.963-1.125, p= 0.318). MR-RAPS analysis agreed with the results of IVW analysis (OR = 1.041, 95% CI: 0.974-1.112, p= 0.228,  Supplementary Table S3. Scatter plots, forest plots and leave-on-out sensitivity analysis plots were displayed in Supplementary Figure S4-6. The results of leave-one-out analysis implicated that the selected genetic variants potentially impacted the pooled results, which suggested that careful interpretations for the results was crucial.

Discussion
The causal roles of MtHz and mtDNA abundance in OA pathogenesis were poorly studied in previous works. To our best knowledge, this is the first MR study to evaluate the causal association between mitochondrial genome traits and OA. Four MR methods were employed to estimate the causal association between mtDNA characterization and OA. The results of MR steiger test verified the causal direction of our research (MtHz and mtDNA abundance were exposure and OA were outcome). No effect of mtDNA abundance and MtHz on overall OA was detected. In subgroup analysis, a suggestively protective role of MtHz on knee OA was observed, but not on other sites. And we did not find that mtDNA abundance was causally associated with any site-specific OA.
Although mitochondria is widely recognized as an important factor of OA development (Blanco et al., 2011), the role of MtHz has not been thoroughly investigated (Suliman and Piantadosi, 2016). Indeed, heteroplasmic mutations in mtDNA are often pathogenetic (Ye et al., 2014).Ananimalstudyhasindicatedthatthestateofheteroplasmyitselfwas deleterious when the two mtDNA sequences contain no pathogenic variants (Sharpley et al., 2012). However, there is potentially a particular level threshold for MtHz (McCormick et al., 2020). For instance, when the A3243G mutation in mitochondrial DNA is present in more than 10%, patients can manifest Type 2 diabetes (Wallace, 2005). And low-frequency mtDNA variants (0.2%-2% heteroplasmy) are extensively presented in healthysubjects(Payneetal.,2013).Furthermore,MtHzcanbebeneficialin health promotion as an intermediate state in emergence of novel mtDNA haplogroups (Wallace, 2016). In previous study, MtHz was found to be significantly relevant to several haplogroups (haplogroup H, J, K, T, U and X)withdifferentcharacterizationsamonghaplogroups (Nandakumaretal., 2021). Mitochondrial haplogroup J has been extensively found to mediate the development of OA (5). A Spanish cohort-study (n case = 457; n control = 262) had reported that haplogroup J was associated with a decreased risk of knee OA (OR = 0.460, 95% CI: 0.282-0.748, p= 0.002) (Rego-Perez et al., 2008). Furthermore, a meta-analysis in European cohorts also suggested that haplogroup J was associated with a lower risk of knee OA (HR = 0.702, 95% CI: 0.541-0.912, p= 0.008) (Fernandez-Moreno et al., 2017). However, the association between mitochondrial DNA variants and OA has not been verified in a large sample observational study (Hudson et al., 2013). The study has only explored the causal relationship in terms of single mutation, but we take MtHz (variation at a wholemtDNA level) as an exposure which contribute to understand the causal role of mitochondria in OA development. In terms of biological mechanism, MtHz has been reported a regulation role of metabolic and epigenomic changes. Kopinski PK et al. had found the levels of mitochondrially drove acetyl-Frontiers in Genetics frontiersin.org CoA decreased at high heteroplasmy (Kopinski et al., 2019). And the reduction of acetyl-CoA levels might influence histone acetylation and activate AMP-activated protein kinase (AMPK) to protect from OA development and progression . Besides, MtHz could play a role in other processes involvedinOApathogenesistoreducetherisk ofOA, including regulating energy production, maintaining mitochondrial proteostasis, suppressing matrix metalloproteinase expression, reducing ROS generation, and promoting mitophagy (Blanco et al., 2018). In conjunction with prior observational findings and mechanism studies, our epidemiologic and genetic findings provided supportive evidence that MtHz may have therapeutic value in knee OA and can be regarded as a candidate biomarker after precisely identifying the relationship between MtHz and each haplogroup (Blanco et al., 2018). However, the reason that MtHz showed diverse effects on different sites of OA remains unclear. So far, we did not find relevant research that investigated the relationship between MtHz and hand, finger, and thumb OA. A case-control study suggested that mitochondrial haplogroups was associated with hip OA (Rego et al., 2010) (OR = 0.661, 95% CI: 0.440-0.993, p= 0.045), however, with no evidence on the causality. In addition, the data sources and sample size used in MR analysis varied in site-specific OA, which mainly included more knee OA cases and less other skeletal joints cases. More MR studies and additional GWAS covering more other site-specific OAs are needed to determine a role of MtHz in the risk of different OA sites.
MtDNA abundance has been regarded as a potential biomarker of mitochondrial function and plays a role in several human diseases (Castellani et al., 2020;Clyde, 2022). Existing literature has suggested that mtDNA abundance may involve in production of inflammatory mediators and regulation of immune function that can influence OA development (Blanco et al., 2018;Zhan et al., 2020). In addition, mtDNA abundance is also associated with sex, advanced age, and elevated BMI, which are also risk factors for OA (16). However, there are few studies to estimate the effects of mtDNA abundance on OA. Only two case-control studies in Asian population were found that reported an association between mtDNA copy number and OA, but their findings were inconsistent (Fang et al., 2014;Zhan et al., 2020). One in Thailand (n case = 204; n control = 169) had found mtDNA abundance in the OA group was significantly lower than that in the control group (p < 0.0001), whereas the results were not adjusted by sex and age (Zhan et al., 2020). Another casecontrol study carried out in southern Chinese (n case = 187; n control = 420) observed a general increase of mtDNA abundance in OA patients (p= 0.019), but obesity was not adjusted in the analysis (Fang et al., 2014). These findings from case-control studies could not determine causal relationships between mtDNA abundance and OA, and the potential confounding variables were not comprehensively considered. Our MR analysis overcame these shortcomings and suggested that genetically determined mtDNA abundance was unrelated to OA. Furthermore, previous study has implicated that the relative contribution of mtDNA abundance might differ between different ethnic groups (Ruiz-Pesini et al., 2004). Considering that our studies have been carried out in European population, the comparisons with results from other ethnic groups such as Asian ancestry require careful considerations.
Our research applies MR methods to investigate causal relationships between mitochondrial genome characterization and OA. However, there are some limitations in our analysis. Firstly, we did not estimate causal effects of mitochondrial genome traits on OA stratified by gender. Mitochondrial genome are maternally inherited, and females may have lower MtHz than males from the perspective of mitochondrial inheritance (Nandakumar et al., 2021). Therefore, the effects of MtHz and mtDNA abundance on OA may differ in gender. Secondly, only participants of European descent are included in the study, but the impact of specific mtDNA variants on diseases could vary in different ethnic groups (Ruiz-Pesini et al., 2004). Additional MR studies on other ethnic groups are needed to probe a causal association between mtDNA characterization and OA. Thirdly, the genetic instruments for mtDNA abundance explain a relatively small amount of phenotypic variance (8.3%), IVs that can account for more variance of mtDNA abundance are warranted to draw robust conclusions. Besides, considering the Bonferroni correction of multiple independent tests, our findings about the association between MtHz and knee OA are deemed suggestive evidence of possible associations (0.01 < p< 0.05). Furthermore, the variation of mtDNA abundance and the levels of MtHz varied in different tissues and cell types. And in original GWASs, mtDNA abundance was measured from leukocyte in blood and MtHz from saliva samples rather than cartilage or synovial tissues. Considering that OA is mainly characterized by progressive loss of cartilage and synovial hyperplasia, the application of data from blood samples and saliva samples would limit the explanation ability of our study to a certain extent (Sellam and Berenbaum, 2010). Observational epidemiology studies exploring an association between concrete levels of MtHz and a risk of OA are needed to improve causal inference.

Conclusion
In conclusion, our MR analyses elucidated that MtHz is a suggestively protective factor of knee OA, implying that MtHz could be a genetically prediction factor and a therapeutic target in the development of knee OA. No causal association was found between mtDNA abundance and OA. Additional MR analyses are warranted to probe causal relationships between mitochondrial genome traits and OA stratified by gender. Moreover, GWAS covering more than one ethnic population are needed to detect the effect of characterization of mtDNA quantity and quality in different ethnic groups.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.