Epigenome-wide association study of nuclear DNA methylation in relation to mitochondrial heteroplasmy

Abstract

We analyze 10,986 participants (mean age 77; 63% women; 54% non-White) across seven U.S. cohorts to study the relationship between mitochondrial DNA (mtDNA) heteroplasmy and nuclear DNA methylation. We identify 597 CpGs associated with heteroplasmy burden, generally showing lower methylation. These CpGs are enriched in dynamically regulated island shores and depleted in CpG islands, indicating involvement in context-specific rather than constitutive gene regulation. In HEK293T cells, we introduce a truncating mtDNA mutation (MT-COX3, mt.9979) and observe a positive correlation between variant allele fraction and methylation at cg04569152, supporting a direct mtDNA-nDNA epigenetic link. Many heteroplasmy-associated CpGs overlap with known methylation-trait associations for metabolic and behavioral traits. Composite CpG scores predict all-cause mortality and incident CVD, with one-unit increases associated with 1.27-fold and 1.12-fold higher hazards, respectively. These findings suggest an mtDNA-nDNA epigenetic connection in aging and disease, though its direction and mechanisms remain to be studied.

Fig. 1. Study flowchart.

Whole-genome sequencing was conducted to identify mitochondrial DNA (mtDNA) heteroplasmy in 10,986 individuals across seven epidemiological cohorts. Cohort- and ancestry-specific association analyses were performed to investigate the relationship between mtDNA heteroplasmy and nuclear DNA (nDNA) methylation levels at CpG sites. Random-effects meta-analyses were subsequently conducted in both pooled and race-stratified datasets. CpG sites demonstrating statistically significant associations (false discovery rate [FDR] < 0.05) were subjected to downstream bioinformatic analyses, including gene ontology enrichment and transcriptomic integration. Experimental validation using targeted genome editing in cell lines, followed by pyrosequencing, provided functional evidence supporting a causal role of mtDNA heteroplasmic variation in modulating nuclear epigenetic states. Integrative analyses further linked heteroplasmy-associated CpGs to complex traits through methylation quantitative trait loci (mQTL) mapping, interrogation of the Epigenome-Wide Association Studies (EWAS) Catalog and Genome-Wide Association Studies (GWAS) Catalog, and Mendelian randomization. Finally, composite methylation scores derived from heteroplasmy-associated CpG sites were evaluated for association with all-cause mortality and cardiovascular disease risk.

Fig. 2. Distribution of heteroplasmy.

A Distribution of heteroplasmy counts (MHC) in study participants (n = 10,986) across the double-stranded mitochondrial DNA (mtDNA) molecule, following quality control procedures (see Methods). Heteroplasmic variants were identified at 2264 unique sites across the mtDNA genome, with each site occurring in less than 1% of participants in the combined sample (Supplementary Data 3). Pink regions indicate the D-loop regions. Blue regions indicate the 13 protein-coding genes. B Distribution of MHC and mitochondrial local constraint score sum (MSS) score across the study participants. MHC represents the total number of heteroplasmic variants in an individual, while MSS reflects the predicted deleteriousness and functional impact of these variants. Approximately 32.4% of participants carried at least one heteroplasmic variant. Among all participants, 12.4% had MSS values between 0.01 and 0.25, 4.4% had values between 0.251 and 0.5, and 7.9% had MSS values greater than 0.5. C The distribution of MHC across age groups indicates a trend of increasing MHC levels with advancing age.

Fig. 3. Association and meta-analyses of mtDNA heteroplasmy with DNA methylation.

A Cohort- and race-specific association analyses using linear regression were followed by meta-analyses in pooled samples (n = 10,986) to identify CpG sites associated with MSS (top) and MHC (bottom). We found 479 MSS- and 166 MHC-associated CpGs at a false discovery rate (FDR) < 0.05 in a two-sided test. The x-axis indicates chromosomes in different colors. The gray dotted line represents the FDR-adjusted significance threshold (FDR = 0.05), and the black dotted line indicates the Bonferroni threshold (p = 1 × 10^-7). B Comparison of effect sizes (beta values) for MSS vs. MHC meta-analyses, showing a Pearson correlation r of 0.75. Red dots indicate CpGs with positive betas and blue dots with negative betas. C Volcano plot of DNA methylation with MSS in meta-analysis. Most significant CpGs (94%, FDR < 0.05) showed negative associations (blue dots), while a minority showed positive associations (red points). D Forest plot of the top CpG, cg15233611, from linear regression analyses of MSS in African American (AA) and White American (EA) participants, and in the meta-analysis. All tests were two-sided t-tests. Among AA participants, beta-estimates (95% CI) were 0.00087 (–0.0097, 0.011) in ARIC, –0.0077 (–0.032, 0.016) in CARDIA, –0.018 (–0.038, 0.0021) in GENOA, –0.022 (–0.033, –0.0097) in JHS, 0.0059 (–0.027, 0.039) in MESA, and –0.021 (–0.066, 0.024) in WHI, with a meta-analysis estimate of –0.010 (–0.021, 0.00063) in all AA participants. For White American (EA) participants, the β (95% CI) values were –0.0087 (–0.022, 0.0041) in ARIC, –0.0094 (–0.031, 0.012) in CARDIA, –0.015 (–0.021, –0.0092) in FHS, –0.013 (–0.035, 0.0086) in MESA, and –0.039 (–0.058, –0.0019) in WHI. The meta-analysis of EA participants yielded a beta of −0.015 (–0.020, –0.010), and the overall meta-analysis combining AA and EA participants produced a beta of –0.014 (–0.019, –0.0099). CI confidence interval, ARIC Atherosclerosis Risk in Communities, CARDIA Coronary Artery Risk Development in Young Adults, FHS Framingham Heart Study, GENOA Genetic Epidemiology Network of Arteriopathy, JHS Jackson Heart Study, MESA Multi-Ethnic Study of Atherosclerosis Study, WHI Women’s Health Initiative.

Fig. 4. Functional characterization of heteroplasmy-associated CpGs.

A A total of 370 genes annotated to the heteroplasmy-associated CpGs were analyzed using DAVID, a comprehensive tool for functional annotation, disease association, and Gene Ontology (GO) analysis. Enrichment was assessed using DAVID’s modified Fisher’s exact test (EASE score), a two-sided test, detecting both over- and under-representation of terms relative to the background set. Enrichment analysis using the Genetic Association Database (GAD) revealed that these genes are enriched in metabolic traits and processes related to multicellular organization (a false discovery rate, FDR < 0.05). B GO biological pathway analysis of these 370 genes indicated enrichment in pathways related to signaling and neurodevelopment (FDR < 0.05). C Correlation between variant allele fraction (VAF) and CpG methylation at cg04569152. Functional validation using HEK293T cell line editing and pyrosequencing demonstrated a significant positive correlation between VAF of an mtDNA nonsense mutation (mt.9979 in MT-COX3) and methylation levels at cg04569152 (p = 0.003 from a two-sided t-test). The mutation was introduced using the FusX TALE-based editing system (FusXTBE), and methylation was quantified by pyrosequencing. Linear correlation was assessed using Pearson’s method, supporting a potential epigenetic influence of mtDNA heteroplasmic variation on nuclear CpG methylation (see Supplementary Fig. 16). D Integrative analysis of cg15233611(SETD1B) illustrates its associations with mtDNA heteroplasmy burden, measured by the mitochondrial local constraint score sum (MSS), as well as with smoking, HDL cholesterol, and triglycerides, supported by methylation-trait associations and Mendelian randomization (MR). Pink circles with the positive sign (+) indicate a positive association, while blue circles with the negative sign (-) indicate a negative association. Association analyses were performed using linear regression, with significance assessed by two-sided t-tests for regression coefficients. Mendelian randomization (MR) analyses used the inverse-variance weighted (IVW) method, with significance assessed by two-sided Wald z-tests.

Fig. 5. Forest plot: association analysis of heteroplasmy count (MHC)-associated DNA methylation score with mortality and incident CVD.

We applied elastic net–regularized Cox proportional hazards regression in the Framingham Heart Study (FHS) training cohort to identify heteroplasmy count (MHC)-associated CpGs, yielding 57 associated with all-cause mortality in panel A and 18 associated with cardiovascular disease (CVD) in panel B. Statistical significance of regression coefficients was assessed using two-sided Wald z-tests. The forest plot presents hazard ratios (HRs) and 95% confidence intervals (CIs) for associations between weighted MHC-CpG scores and outcomes (mortality and incident CVD) across training and testing cohorts, adjusted for age, self-reported sex, and smoking status (never, former, current). In FHS, the 57- and 18-CpG scores yielded HRs (95% CI) of 1.61 (1.50–1.72) for all-cause mortality and 1.29 (1.23–1.36) for CVD, respectively. Internal validation in the Jackson Heart Study (JHS) and the Multi-Ethnic Study of Atherosclerosis (MESA) showed, for all-cause mortality, HRs (95% CI) of 0.99 (0.77–1.28) in MESA White American (EA) participants, 1.18 (1.06–1.31) in JHS, and 1.02 (0.70–1.48) in MESA African American (AA) participants for the 57-CpG score. For CVD, the 18-CpG score yielded HRs (95% CI) of 0.97 (0.78–1.21) in MESA EA, 1.19 (1.05–1.36) in JHS, and 1.41 (0.97–2.04) in MESA AA. Meta-analysis of JHS and MESA participants produced HRs (95% CI) of 1.12 (0.99–1.27) for all-cause mortality and 1.15 (0.97–1.36) for CVD. External validation in the Health and Retirement Study (HRS), which was not used for CpG selection, found that for all-cause mortality, the 57-CpG score yielded HRs (95% CI) of 1.28 (1.16–1.42) in EA, 1.62 (1.31–2.01) in AA, and 1.78 (1.18–2.70) in Hispanic American (HA) participants. For CVD, the 18-CpG score yielded HRs of 1.07 (0.97–1.19) in EA, 0.83 (0.59–1.16) in AA, and 1.48 (1.06–2.05) in HA. Meta-analysis of MESA, JHS, and HRS participants produced HRs (95% CI) of 1.27 (1.09–1.48) for all-cause mortality and 1.12 (0.99–1.26) for CVD.