Multi-Omic Insight Into the Molecular Networks of Mitochondrial Dysfunction in the Pathogenesis of Primary Open-Angle Glaucoma

Abstract

Purpose: The primary objective of this study is to comprehensively investigate the genetic association between mitochondrial gene expression, methylation, and primary open-angle glaucoma (POAG) using Mendelian randomization (MR) analysis.

Methods: In this study, mitochondrial-related genes and methylation sites were extracted from two expression quantitative trait locus (eQTL) datasets (eQTLGen and GTEx V8) and one methylation quantitative trait loci dataset for summary data-based Mendelian randomization (SMR) analysis and heterogeneity in dependent instrument (HEIDI) testing. Further MR analysis was conducted to explore their genetic association with POAG. Additionally, colocalization analysis was used to investigate whether there are shared genetic variations between the them.

Results: Our results indicate a genetic association between the mitochondrial genes ISCA2 and ME3 and POAG. The findings were consistent across different eQTL datasets (eQTLGen and GTEx V8). Additionally, the methylation sites corresponding to ISCA2 (cg16374328, cg15981604, and cg0584866) and ME3 (cg03646605, cg03955932, cg12030550, cg14883291, cg15698545, and cg19918734) also exhibited a genetic association with POAG. Colocalization analysis further revealed that these genes and their associated methylation sites share genetic variations with POAG.

Conclusions: Our study suggests that the mitochondrial genes ISCA2 and ME3, along with their corresponding methylation sites, may play a significant role in the pathogenesis of POAG.

Translational Relevance: This study provides a potential idea for the development of mitochondrial-based drugs for the treatment of POAG.

Introduction

Glaucoma is one of the leading causes of blindness worldwide, particularly affecting individuals aged 40 and above.¹ Although numerous studies have identified risk factors for glaucoma, such as age, race, and intraocular pressure, its exact pathogenesis remains unclear.^2,3 With the aging population, the incidence of glaucoma is steadily increasing. It is estimated that by 2040, over 100 million people globally will be affected by glaucoma. Primary open-angle glaucoma (POAG) is the most common subtype of glaucoma. An important pathological feature of glaucoma is the progressive structural and functional loss of retinal ganglion cells (RGCs).⁴ Currently, although various treatment options are available for POAG, some patients still experience persistent RGCs damage despite treatment, ultimately leading to blindness.⁵ 

Mitochondria are vital organelles in cells, playing a crucial role in regulating cellular energy supply. Previous studies have indicated an association between mitochondrial DNA (mtDNA) and POAG.⁵ Research involving African American and European populations has also shown variations in mtDNA that are associated with risk genes for POAG.^6,7 Additionally, animal experiments have confirmed that alterations in mitochondrial DNA can lead to the development of POAG.^8,9 Gu et al. analyzed the retinal gene expression profiles in an animal glaucoma model and found that elevated intraocular pressure leads to oxidative stress in mDNA and mRNA, with stronger associations to mitochondrial RNA.⁸ In addition, Zeng et al. demonstrated that increased intraocular pressure alters the dynamics of dynamin-related protein 1 (Drp1), leading to mitochondrial dysfunction and contributing to the development of glaucoma.⁹ Although the critical role of mitochondria in glaucoma's pathogenesis and progression has been recognized in recent years, the relationship between mitochondrial genes and glaucoma, as well as the impact of mitochondrial gene methylation on glaucoma, remains underexplored. 

Mendelian randomization (MR) is a novel approach for inferring causal relationships between exposures and outcomes.¹⁰ By using genetic data, such as Genome-Wide Association Studies (GWAS) data and molecular quantitative trait loci (QTL) data, this method could infer causality while mitigating the effects of confounding factors.^11,12 Consequently, it has been widely applied in drug target discovery and in exploring gene-disease associations. In this study, we used expression quantitative trait locus (eQTL) and methylation quantitative trait loci (mQTL) data to investigate the association among mitochondrial gene expression, methylation and POAG. 

Methods

Study Design

First, the eQTL data were downloaded from eQTLGen consortium and Genotype Tissue Expression Project database (GTEx v8). Mitochondrial genes were subsequently filtered from this dataset and subjected to summary-data-based MR (SMR) analysis to investigate their association with POAG. Second, mQTL data were used to analyze the association between gene methylation and POAG, and the results of the SMR analyses were organized and mitochondrial genes were screened out. Subsequently, we performed MR analysis and colocalization analysis on the selected genes and methylation sites (cGp) to further clarify the genetic associations between these genes and cGp sites with POAG. The results of the analyses were integrated, and common genes were identified (Fig. 1). 

Figure 1.

 

Flow chart of the study. This study screened mitochondrial genes from two eQTL databases and one mQTL database. Subsequently, the filtered data were used to perform SMR, HEIDI test, MR analysis and colocalization analyses with POAG. Finally, the results of the analyses were integrated and categorized into three groups (Tier 1, Tier 2, and Tier 3).

Data Source

The eQTLs data were sourced from the eQTLGen¹³ and Genotype-Tissue Expression (GTEx) V8 databases (https://www.gtexportal.org/home). The eQTLGen database provided gene expression data specifically for whole blood, whereas the GTEx V8 database included eQTL data from various tissue types and whole blood data was selected for analysis in this study. The mQTL data were obtained from McRae et al.,¹⁴ and this data can be downloaded in SMR format from the Yang Laboratory (https://yanglab.westlake.edu.cn/software/smr/#DataResource). The POAG data were sourced from the FinnGen database, which comprises genetic data from nearly 500,000 individuals of Finnish descent.¹⁵ This study used the tenth release of the database and selected POAG for analysis.¹⁵ Mitochondria-related genes were obtained from the MitoCarta3.0 database, which includes a total of 1136 human mitochondrial genes.¹⁶ All data information could be found in Table 1. 

Table 1.

 

Detail Information of GWAS Data

 

Ethics Approval and Consent to Participate

The data used in this study were obtained from public databases and all data are publicly available on the corresponding websites. The corresponding ethical review has been approved in the original literature. 

SMR Analysis and HEIDI Test

Under large sample conditions, SMR analysis can achieve higher statistical power.¹⁷ In this study, we used the SMR software (Version 1.3.1) for analysis, using the default parameters during the execution of SMR analysis.¹⁷ To control for potential false positives, the Benjamini-Hochberg method was applied to adjust the P values obtained from the SMR analysis. Results with an adjusted P value (P_{FDR of SMR}) less than 0.05 were considered statistically significant, whereas results with P_{FDR of SMR} > 0.05 but P_SMR < 0.05 were deemed suggestively associated.¹⁸ In addition, heterogeneity in dependent instrument (HEIDI) testing was conducted as part of the SMR analysis to evaluate whether the observed associations were influenced by linkage disequilibrium. Based on previous studies, P_HEIDI > 0.01 indicates no pleiotropic effects, while results with P_HEIDI < 0.01 were discarded and excluded from further analysis.¹¹ 

MR Analysis and Colocalization Analysis

After obtaining the results from the SMR analysis, we extracted the corresponding single-nucleotide polymorphisms (SNPs) associated with genes and methylation sites (P_SMR < 0.05 and P_HEIDI > 0.01 in both three databases) from the SMR data as candidate instrumental variables (IVs). Using the method provided by the Yang Laboratory, this study extracted the corresponding IVs through the SMR software.¹⁷ Detailed operational procedures can be found at the following website: https://yanglab.westlake.edu.cn/software/smr/#Extract/removeasubsetofdata. 

Subsequently, to exclude the effects of linkage disequilibrium and weak instrument bias, we used PLINK (v1.9) software¹⁹ and referenced the 1000 Genomes Project²⁰ in European populations to perform clumping (window size = 1000 kb, r² = 0.3, P < 5e-8),²¹ excluding SNPs with F-statistics < 10.²² The calculation formula for the F-statistic is referenced from Zhao et al., utilizing the beta values and standard errors (SE) of the IVs to compute an approximate value of the F statistic (F = Beta²/SE²).²³ Additionally, this study excluded IVs with a minor allele frequency (MAF) < 0.01.²⁴ The selected IVs were then matched with POAG. When the number of matched IVs exceeded one, we used the inverse variance weighted method^25,26 and the constrained maximum likelihood and model averaging (cML-MA) method²⁷ for MR analysis. Cochran's Q test was used to assess the heterogeneity of the MR analysis results,²⁸ whereas MR-Egger intercept test²⁹ and MR-PRESSO³⁰ were used to evaluate horizontal pleiotropy. If the number of IVs equaled one, the Wald ratio method was applied for analysis, and heterogeneity and horizontal pleiotropy were not assessed. 

After conducting MR analysis and sensitivity analysis, the results of the MR analysis were further adjusted. The Benjamini-Hochberg method was used to correct the P values obtained from the MR analysis to reduce the likelihood of false-positives. Results with an adjusted P value (P_{FDR of MR}) < 0.05 were considered statistically significant, whereas results with P_{FDR of MR} > 0.05 but P_MR < 0.05 were deemed suggestively associated.¹⁸ 

To determine whether the identified genes share genetic variants with the disease, this study used the coloc R package to perform colocalization analysis for genes (P_SMR < 0.05 and P_HEIDI > 0.01 in both three databases) and methylation sites with POAG. The colocalization analysis reports five distinct posterior probabilities derived from different hypotheses (PPH): (1) no shared causal variants exist between POAG and the gene at the genetic locus (H0); (2) causal variants exist only for POAG (H1); (3) causal variants exist only for the gene (H2); (4) POAG and the gene have distinct causal variants (H3); and (5) POAG and the gene share causal variants (H4).³¹ Based on prior studies, genes with the combined posterior probability of H3 and H4 exceeding 0.8 were considered positive results (PPH3 + PPH4 > 0.80).³² This study analyzed regions extending 1000 kb upstream and downstream of the gene for eQTL and 500kb for mQTL, while keeping all other parameters at their default settings. 

Integrating Results From eQTLs and mQTLs

To comprehensively describe the association between mitochondrial gene expression, methylation, and POAG, this study integrated data from multiple levels. Since the colocalization results provide clearer insights into the association between POAG and mitochondrial genes, this study considers the findings from the colocalization analysis as the most important reference, and categorized into three tiers: (1) Tier 1: Genes that were significant associated (P_{FDR of SMR} and P_{FDR of MR} < 0.05) in both eQTLGen and GTEx v8 datasets as well as in the mQTLs datasets, with colocalization analysis showing PPH3 + PPH4 > 0.80 in both three datasets; (2) Tier 2: Genes that were significant associated (P_{FDR of SMR} and P_{FDR of MR}< 0.05) in eQTLGen and GTEx v8, and suggestive associated (P_SMR and P_{FDR of MR} < 0.05) in mQTL, and the corresponding colocalization analysis showing PPH3 + PPH4 > 0.80; (3) Tier 3: Genes that were suggestive associated (P_SMR and P_{FDR of MR} < 0.05) in eQTLGen, GTEx v8 and mQTLs, and the corresponding colocalization analysis showing PPH3 + PPH4 > 0.80. 

Statistical Software

This study was conducted using R software (version 4.2.3) and SMR software (version 1.3.1). For colocalization analysis, the coloc R package (version 5.2.2) was used.³¹ For MR analysis, the TwoSampleMR package (version 0.6.14) was used and cMR-MA analysis was performed by MRcML package (version 0.0.0.9). 

Results

SMR Results of Mitochondrial Gene Expression and Methylation

The results of the SMR analysis for the eQTLGen and GTEx V8 database are presented in the Supplementary Table S1. A total of 826 mitochondrial genes were matched in the eQTLGen, with suggestively associated genes (P_SMR < 0.05, P_{FDR of SMR} > 0.05 and P_HEIDI > 0.01) including 68 related to POAG (Supplementary Table S2). The significantly associated genes (P_{FDR of SMR} < 0.05 and P_HEIDI > 0.01) included 1 related to POAG (Supplementary Table S2, Fig. 2). The GTEx V8 database matched only 410 mitochondrial genes (Supplementary Table S1). After screening for suggestively associated data, 35 genes were found related to POAG (Supplementary Table S2). Among the significantly associated genes in the GTEx V8 database, there was one related to POAG (Supplementary Table S2 and Fig. 2). Methylation data matched 703 mitochondrial genes, including 2745 CpG sites, with detailed results provided in the Supplementary Table S1. After filtering for genes with P_SMR < 0.05 and P_HEIDI > 0.01, 112 genes (244 CpG sites) remained causally associated POAG (Supplementary Table S2). After performing multiple corrections, 1 CpG site (ISCA2: cg16374328) remained significantly associated with POAG (Supplementary Table S2, Fig. 3). Subsequently, we took the intersection of the suggestively associated genes from the three datasets (P_SMR < 0.05 and P_HEIDI > 0.01) and ultimately identified 14 genes and 68 methylation sites that exhibit genetic associations with POAG (Figs. 2, 3). 

Figure 2.

 

The result of SMR analysis and HEIDI test in eQTL data. The SMR analysis revealed that in the three datasets, the P_SMR for 14 genes in eQTLs databases (eQTLgen and GTEx V8) were less than 0.05, and the HEIDI test results were greater than 0.01.

Figure 3.

 

The results of SMR analysis and HEIDI test in mQTL database. After taking the intersection of the SMR analysis results from the mQTL dataset and the other two eQTL datasets, we ultimately identified 14 genes and 68 methylation sites that have a causal association with POGA (P_SMR < 0.05, P_HEIDI > 0.01).

MR Analysis and Colocalization Analysis

After taking the intersection of the suggestively associated genes from the three datasets, we identified 14 genes and 68 methylation sites that exhibited genetic associations across all three datasets. We subsequently extracted data corresponding to these genes and methylation sites from the SMR-format data for MR analysis. The selected IVs are presented in Supplementary Table S3. After taking the intersection of the MR analysis results (Supplementary Table S4, Figs. 4, 5) from the three datasets, we ultimately identified eight genes (ISCA2, ME3, MRPL21, MRPS10, NDUFS5, NSUN4, PISD, and SARDH) and 45 methylation sites that exhibit genetic associations with POAG (P_{FDR of SMR} < 0.05). Sensitivity analyses indicate that there is no horizontal pleiotropy or heterogeneity in the MR analysis (Supplementary Table S5). 

Figure 4.

 

The results of MR analysis in eQTL data. MR analysis on the genes and data associated with POGA identified in the SMR analysis, which included 14 relevant genes.

Figure 5.

 

The results of MR analysis in mQTL data. The results of the MR analysis between methylation sites in the mQTL and POGA included 14 genes and 68 methylation sites.

Colocalization analysis of these sites revealed that SARDH, ME3, and ISCA2 had PPH3 + PPH4 values greater than 0.8 in both the eQTLGen and GTEx v8 datasets. Among the methylation data, there were 20 methylation sites associated with four genes (including NDUFS5, ME3, ISCA2, and PISD) that shared genetic variation with POAG (Supplementary Table S6). After taking the intersection of the colocalization results for the three databases, we found that ISCA2 and ME3, along with their methylation sites, exhibited genetic associations with POAG (Figs. 6–8). 

Figure 6.

 

The results of colocalization analysis in eQTLGen and GTEx V8. (A, B) ISCA2 and ME3 in eQTLGen. (C, D) ME3 and ISCA2 in GTEx V8.

Figure 7.

 

The results of colocalization analysis in methylation sites of ISCA2. Three methylation sites of ISCA2, including cg16374328 (7A), cg05848660 (7B), cg15981604 (7C).

Figure 8.

 

The results of colocalization analysis in methylation sites of ME3. Six methylation sites of ME3, including cg03646605 (8A), cg03955932 (8B), cg12030550 (8C), cg14883291 (8D), cg15698545 (8E), and cg19918734 (8F).

Integrating Results

The ISCA2 gene (eQTLGen: odds ratio [OR] = 1.147; 95% confidence interval [CI], 1.083–1.215; GTEx V8: OR = 1.400; 95% CI, 1.209–1.622) exhibits a Tier 1 relationship with POAG, demonstrating significant associations in both eQTLGen and GTEx V8 database, as well as in the mQTL data (ISCA2 [cg16374328]: OR = 1.163; 95% CI, 1.089–1.242). The colocalization analyses from three datasets both indicate PPH3 + PPH4 > 0.8 (Table 4). The other two methylation sites of ISCA2 (cg15981604 and cg05848660) also meet the Tier 2 criteria (Tables 2, 3). 

Table 2.

 

The Result of SMR Analysis After Integrating Results

 

Table 3.

 

The Result of MR Analysis After Integrating Results

 

Table 4.

 

The Result of Colocalization Analysis After Integrating Results

 

In the results of the SMR analysis, cg15981604 (OR = 0.751; 95% CI, 0.634–0.890) and cg05848660 (OR = 1.380; 95% CI, 1.138–1.673) were identified as genetic protective and risk factors for POAG, respectively. In the MR analysis, cg15981604 (OR = 0.751; 95% CI, 0.649–0.870) and cg05848660 (OR = 1.380; 95% CI, 1.172–1.624) yielded consistent conclusions. Colocalization analysis indicated that both methylation sites share genetic variants with POAG (PPH3 + PPH4 > 0.8). These two methylation sites demonstrate significant genetic associations only in the MR analysis (P_{FDR of MR} < 0.05), while exhibiting suggestive genetic associations in the SMR analysis (P_{FDR of SMR} > 0.05 and P_SMR < 0.05). 

ME3 (eQTLGen: OR = 0.912; 95% CI, 0.851–0.978; GTEx V8: OR = 0.893; 95% CI, 0.819–0.975) and six methylations (cg03646605, cg03955932, cg12030550, cg14883291, cg15698545, cg19918734) met the Tier 3 criteria, showing a suggestive association in the SMR analysis, while exhibiting significant genetic associations in the MR analysis (Tables 2, 3). The SMR analysis revealed that cg03646605 (OR = 1.080; 95% CI, 1.018–1.147), cg03955932 (OR = 1.057; 95% CI, 1.013–1.104), cg12030550 (OR = 1.093; 95% CI, 1.02–1.171), cg14883291 (OR = 1.094; 95% CI, 1.02–1.173), and cg15698545 (OR = 1.064; 95% CI, 1.015–1.116) are genetically associated with the occurrence of POAG and are linked to an increased risk of developing POAG. Conversely, from a genetic association perspective, cg19918734 (OR = 0.836; 95% CI, 0.721–0.970) is associated with a decreased risk of POAG. Consistent results were obtained from the MR analysis, which showed cg03646605 (OR = 1.086; 95% CI, 1.048–1.125), cg03955932 (OR = 1.064; 95% CI, 1.040–1.088), cg12030550 (OR = 1.109; 95% CI, 1.069–1.150), cg14883291 (OR = 1.098; 95% CI, 1.054–1.143), cg15698545 (OR = 1.056; 95% CI, 1.029–1.083), and cg19918734 (OR = 0.854; 95% CI, 0.765–0.953). Additionally, colocalization analysis indicated that all six methylation sites share genetic variations associated with POAG (Supplementary Table S6). 

Discussion

In this study, we used SMR and MR analysis for the first time to comprehensively investigate the impact of mitochondrial-related gene expression and methylation changes on POAG. From a genetic perspective, our research predicts that increased expression of the ISCA2 gene may elevate the risk of developing POAG, with its CpG site (cg16374328) also contributing to the likelihood of developing this condition. This conclusion was validated across three datasets (Tier 1). Moreover, we found that two CpG sites of ISCA2(cg15981604 and cg05848660) exhibited secondary associations with POAG (Tier 2). Specifically, six CpG sites on the ME3 gene (cg03646605, cg03955932, cg12030550, cg14883291, cg15698545, cg19918734) were associated with the risk of POAG (Tier 3). 

Iron-sulfur cluster (ISC) assembly is a critical cellular process essential for the production of various ISC-containing proteins in the nucleus, mitochondria, and cytoplasm.³³ ISCA2 plays a pivotal role in this process. Mutations in human ISCA2 lead to multiple mitochondrial dysfunction syndrome type 4 (MMDS4), a condition in which affected infants are unable to survive.³⁴ The neurological manifestations of ISCA2 deficiency include neurodevelopmental dysfunction.³⁴ Our study is the first to predict the impact of altered ISCA2 gene expression on the risk of POAG from a genetic perspective. However, research on the relationship between ISCA2 and POAG remains limited. One possible mechanism involves ISCA2 regulating RGCs ferroptosis, thereby contributing to POAG pathogenesis. A recent study demonstrated that ISCA2 inhibitors reduce HIF production and induce ferroptosis.³⁵ The study demonstrates that the downregulation of ISCA2 in renal cancer cells significantly reduces the levels of HIF proteins.³⁵ Additionally, the administration of drugs or siRNA that lowers ISCA2 expression can induce iron overload within cells, leading to the occurrence of ferroptosis.³⁵ Previous literature has also reported that the stability of iron metabolism is crucial for the survival of RGCs.³⁶ Furthermore, the neurotoxicity of glutamate can similarly cause damage to RGCs. Inevitably, glutamate, as a neurotransmitter involved in phototransduction, can accumulate excessively in RGCs, resulting in neurodegenerative damage.³⁷ Additionally, clinical studies have indicated that serum iron levels in patients with glaucoma are higher than those in healthy individuals.³⁸ Interestingly, our study found a positive correlation between ISCA2 gene expression levels and POAG occurrence, which may be attributed to tissue specificity. Further investigation into the relationship among ISCA2, ferroptosis, and POAG may provide novel therapeutic approaches for POAG. 

The ME3 gene encodes mitochondrial malic enzyme, which primarily catalyzes the production of NADPH. In 2023, gene risk scoring based on the GTEx database revealed that ME3 may be involved in the onset of POAG.³⁹ Moreover, ME3 expression exhibited a negative correlation with POAG, consistent with our findings. The results from SMR analyses (eQTLGen: OR = 0.912; 95% CI, 0.851–0.978; GTEx V8: OR = 0.893; 95% CI, 0.819–0.975) and MR analyses (eQTLGen: OR = 0.919; 95% CI, 0.888–0.951; GTEx V8: OR = 0.892; 95% CI, 0.841–0.946) indicate a potential genetic association between the ME3 gene and POAG. The possible mechanism is that the downregulation of ME3 may lead to a reduction in NADPH production, which is a crucial molecule for cellular antioxidant defense, involved in glutathione reduction and the clearance of reactive oxygen species.^40,41 This may result in oxidative stress within RGCs, contributing to the development of POAG. Furthermore, this study also found that the six methylation sites of ME3 are involved in the regulation of POAG. Unfortunately, there is currently a lack of corresponding studies to demonstrate the biological functions of these six methylation sites. This research provides potential insights for future investigations in this area. 

Although this study analyzed the genetic correlation between POAG and mitochondrial genes, it has several limitations. First, due to the lack of corresponding protein quantitative trait loci (pQTL) data, the study did not explore the association between the identified genes and POAG at the protein expression level. Second, the analysis was limited to eQTL data from European populations due to the absence of eQTL data from other populations. Therefore the conclusions of this study are confined to European populations and cannot be generalized to other ethnic groups, necessitating further exploration in future research. Third, the study only examined the impact of gene expression and methylation on POAG risk from a genetic correlation perspective, without providing experimental validation through animal models. Thus the findings of this study may require further confirmation through subsequent research. 

Conclusions

This study is the first to systematically analyze the effects of gene expression and methylation sites on the risk of POAG using MR and SMR analysis. It identified two genes and their methylation sites that may potentially influence the risk of POAG development. 

Acknowledgments

Disclosure: L.-H. Chen, None; H.-J. Zhang, None; Q.-A. Xiao, None