Mitochondrial Genetics Reinforces Multiple Layers of Interaction in Alzheimer’s Disease

Simple Summary Nuclear DNA remains the main source of genome-wide loci association in neurodegenerative diseases, only partially accounting for the heritability of Alzheimer’s Disease (AD). In this context, mitochondrial DNA (mtDNA) is gaining more attention. Here, we investigated mitochondrial genes and genetic variants that may influence mild cognitive impairment and AD, through an integrative analysis including both differential gene expression and mitochondrial genome-wide epistasis analysis. Our results highlight important layers of interactions involving mitochondrial genetics and suggest specific molecular alterations as potential biomarkers for AD. Abstract Nuclear DNA has been the main source of genome-wide loci association in neurodegenerative diseases, only partially accounting for the heritability of Alzheimer’s Disease (AD). In this context, mitochondrial DNA (mtDNA) is gaining more attention. Here, we investigated mitochondrial genes and genetic variants that may influence mild cognitive impairment and AD, through an integrative analysis including differential gene expression and mitochondrial genome-wide epistasis. We assessed the expression of mitochondrial genes in different brain tissues from two public RNA-Seq databases (GEO and GTEx). Then, we analyzed mtDNA from the ADNI Cohort and investigated epistasis regarding mitochondrial variants and levels of Aβ1−42, TAU, and Phosphorylated TAU (PTAU) from cognitively healthy controls, and both mild cognitive impairment (MCI) and AD cases. We identified multiple differentially expressed mitochondrial genes in the comparisons between cognitively healthy individuals and AD patients. We also found increased protein levels in MCI and AD patients when compared to healthy controls, as well as novel candidate networks of mtDNA epistasis, which included variants in all mitochondrially-encoded oxidative phosphorylation complexes, 12S rRNA and MT-DLOOP. Our results highlight layers of potential interactions involving mitochondrial genetics and suggest specific molecular alterations as potential biomarkers for AD.


Introduction
Alzheimer's Disease (AD) is a progressive neurodegenerative condition with a complex origin that leads to a myriad of symptoms, such as severe memory loss, confusion, multiple cognitive deficiencies and personality changes. AD is the most common cause of dementia, responsible for 60-80% of cases worldwide [1]. AD progression is slowgradually worsening over years or even decades-and a final diagnosis is commonly

Materials and Methods
Our methodology consisted of different approaches to analyzing gene regulation and genomics from mtDNA data. First, differential expression (DE) analyses of mitochondrial genes were performed. DE analyses were based on RNA-seq data from brain tissues sampled from healthy young adults (HYA), healthy elderly adults (HEA) and Alzheimer's Disease cases (AD). Next, we investigated epistasis of mitochondrial genome variants genotyped for AD, MCI and healthy individuals from the ADNI Cohort. Then, genome-wide association analyses were performed based on mitochondrial data. SNP-SNP interaction analyses were performed for mtDNA variants in a quantitative design to investigate candidate associations with cerebrospinal levels of Aβ 1−42 , TAU and PTAU.

RNA-seq Transcriptome Data from Brain Tissues
First, we processed two available RNA-seq experiments found in the Gene Expression Omnibus (GEO) for brain tissues. Transcriptome data were generated for lateral temporal lobe and fusiform gyrus tissues by third parties and both experiments are stored in GEO under accession numbers GSE104704 and GSE125583, respectively. Sample numbers and age ranges of the individuals included in these experiments are presented in Table 1. In addition to the aforementioned brain tissues, we extracted RNA-seq data from the Genotype Tissue Expression database (GTEx, v.8), which provides genome-wide data and a summary of expression quantitative trait loci analysis for a large set of tissue donors [18]. This allowed us to extract transcriptome data from 54 tissues, from 948 donors and a total of 17,382 samples. Around 50% of the donors presented ages ranging from 60 to 70 years old. Out of GTEx large data, we extracted expression data for 13 different brain tissues. For further analyses, these data were used to explore gene expression patterns in healthy brains, which allowed us to investigate how different genes and pseudogenes are expressed at baseline.

ADNI Cohort Data: mtDNA Sequence Data and Biomarkers
Mitochondrial genomes were sequenced for 809 individuals according to the protocol previously described in [19]. According to [20,21], the individuals included in this cohort are mainly of mitochondrial European ancestry, mostly belonging to haplogroups H, I, J and K, widely present throughout the European continent; by protocol, individuals were clinically analyzed using medical resonance imaging, cognition measurement scores, APOE genotyping and a series of biomarkers extracted from plasma and CSF.

Differential Gene Expression Analysis in Brain Tissues
In order to process the RNA-seq experiments, we downloaded .sra files and converted them to .fastq files using the SRA-Toolkit [23]. QC was assessed by FastQC and MultiQC [24,25]. Trimmomatic was performed for read cleaning and a second round of QC [26].
After QC, reads were aligned to the UCSC reference genome (University of California Santa Cruz), version hg19 (http://hgdownload.soe.ucsc.edu/downloads.html, accessed on 4 August 2021), using the STAR tool [27], and annotated following the coordinate regions of the Gencode genome (v19) (https://www.gencodegenes.org/human/release _19.html, accessed on 4 August 2021). Finally, read counting was performed using the HTseq library [28], implemented in the Python programming language, using the function scripts.count, thus generating files with the abundance of each gene per sample.
DE analyses of the three experiments were performed with the edgeR package in R version 3.5.2 [29]. In this analysis, the gene count matrix was filtered for reads with more than 10 counts. The resulting matrices were normalized, thereby converting expression values to a scale that allows adequate data comparisons. Lastly, we adjusted the gene counts to a generalized negative binomial log-linear model and controlled for multiple testing with False Discovery Rate (FDR ≤ 0.05).

Gene-Gene Interaction Network Analysis
GeneMANIA (www.genemania.org, accessed on 4 August 2021) is a user-friendly web tool for generating hypotheses about gene function and gene prioritization support [30]. GeneMANIA finds similar genes using functional genomics, transcriptomics and proteomics data given a seed gene list. The tool receives as input a list of genes to build a gene-gene interaction network, which is based on extensive biological datasets for functional similarity analysis, co-expression or even interactions at protein level. Here, we employed GeneMANIA to analyze whether seed mtDNA genes and their interactions could be involved in the development of AD or interact with other nuclear genes related to neurogenerative processes.

Epistasis between mtDNA Variants with an Impact on CSF Biomarker Levels
Epistasis analyses were performed with multiple linear regression for quantitative traits. The linear model is based on the allele dosage for each mtDNA variant, A and B, and fits the following model, in which the interaction test is based on the b 3 coefficient: The mtDNA variants were filtered for minor allele frequency (MAF) greater than 0.01 and 99% genotyping rate; monomorfic SNPs were removed. MAF threshold was used to exclude rare SNPs from our analysis. Then, epistasis analyses were performed in the context of AD considering the ADNI groups (AD, MCI and healthy subjects, N = 475) and association analyses with AB, PTAU and TAU levels. Next, we conducted a mitochondrial genome-wide epistasis analysis using a linear regression model for quantitative trait analysis association. The linear regression was implemented in PLINK software and was executed with the -epistasis command, release v1.9 [31]. Association hits with a p-value ≤ 5.0 ×10 −5 were considered statistically significant.

Screening for Differential Expression of mtDNA Genes and Pseudogenes
Three comparative analyses of gene expression were performed for the two datasets: HYA vs. AD (GSE104704), HEA vs. AD (GSE104704) and HEA vs. AD (GSE125583). When comparing HYA vs. AD (lateral temporal lobe), 12 genes were found to be DE, while three were DE in HEA vs. AD comparisons (lateral temporal lobe and fusiform gyrus). Considering them individually, there were 15 DE genes in total: five tRNA genes, eight pseudogenes and two isoforms.
As seen in Table 2, the following genes were differentially expressed in the investi- Notably, MTRNR2L1 and MTND1P23 appeared in the results of more than one analysis: MTRNR2L1 was DE in two (with downregulated expression in both lateral temporal lobe tissues) and MTND1P23 was DE in all three analyses (with upregulated expression in both lateral temporal lobe tissues, while downregulated in the fusiform gyrus tissue). This pattern suggests a tissue-specific regulation of such genes in AD. The other DE genes presented individual patterns in the investigated tissues, i.e., only appeared once.
In addition to Table 2, we plotted the results from the previous analysis for all three experiments (HYA vs. AD in lateral temporal lobe, HEA vs. AD in lateral temporal lobe and HEA vs. AD in fusiform gyrus), as a clear representation of the DE genes, as shown in Figure 1A. To gain insight into the expression of such genes, we performed a comparative expression analysis of all 15 genes in 13 healthy brain tissues with data extracted from the GTEx portal, which did not include fusiform gyrus or lateral temporal lobe tissues specifically ( Figure 1B). From these results, interaction networks were searched for each gene in GeneMANIA, but were only found for MTRNR2L1 ( Figure 1C). As seen in Figure 1B, out of the 15 genes, five were highly overexpressed in all of the investigated tissues (MTND2P28, MTND1P23, MT-TM, MT-TV and MT-TL1), with two exceptions (MT-TM was mildly overexpressed in the cerebellum and cerebellar hemisphere). It should be noted that MTND2P28 presented much higher levels of expression than others. Other mildly overexpressed genes in all tissues included MT-TH, MTND5P11 and MT-TS2. All other genes showed similar underexpression patterns and, to some extent, were grouped together. As for the tissues, it is noteworthy that, based on gene expression patterns, the group formed by the amygdala and hippocampus was similar to nucleus accumbens (basal ganglia) and that this cluster was also similar to another, formed by the anterior cingulate cortex and hypothalamus; these tissues also presented great similarity with the caudate and putamen, both basal ganglia regions and, to a lesser extent, substantia nigra (also part of the basal ganglia). Cortex and frontal cortex were also grouped together, and these tissues were close to the spinal cord. Lastly, the cerebellum and cerebellar hemisphere formed a cluster that, although still similar, was the most separated from the others.
In the interaction network for MTRNR2L1 shown in Figure 1C, we found 10 strongly linked genes, of which five seemed to present physical interactions (MPHOSPH8, PPA1, TRIM2, TRIM11 and TSFM) and the other five shared protein domains not only with MTRNR2L1 but also with each other (MTRNR2L3, MTRNR2L4, MTRNR2L5, MTRNR2L8 and MTRNR2L10). In addition, PPA1 (Inorganic Pyrophosphatase 1) presented genetic interactions with TSFM (Ts Translation Elongation Factor, Mitochondrial) and TRIM2 (Tripartite Motif Containing 2).

Epistasis between mtDNA and CSF PTAU and TAU Levels
Then, we hypothesized that some patterns of differential expression observed in the previous analysis could result from the effects caused by the presence of mtDNA variants, so we searched the ADNI data for epistasis between mtDNA genes related to Aβ 1−42 , PTAU and TAU levels extracted from CSF. Interestingly, linear regression results revealed epistatic interactions among mtDNA variants and PTAU and TAU levels, but not for Aβ 1−42 .
For PTAU levels, we identified a significant gene network between MT-RNR1(709) and two variants, namely MT-ATP6(9632) and MT-ND4(12083) ( Figure 2); however, we did not find statistically significant differences between groups of genotype pairs when analyzing PTAU levels within each sample group (see Figure 2A,B).
When analyzing CSF TAU levels among sample groups in relation to the genotypes of these three variants, we found statistically significant associations for MT-DLOOP(194) and MT-ND5(13135), but not for MT-DLOOP(152) (Figure 3A-C). For instance, for MT-DLOOP(194), we found that individuals with MCI carrying CC/AA genotypes in this variant and MT-ATP8(8701), respectively, had significantly higher CSF TAU levels when compared to CC/GG genotype carriers ( Figure 3D (Figure 3H), but we found higher levels of CSF TAU in healthy and AD individuals with the TT/GG genotype for MT-COX1(7476) and MT-ND5(13135) in comparison to the analyzed pair for these variants ( Figure 3I). Carriers of GG/GG for MT-ND3(10172) and MT-ND5(13135) with MCI had lower TAU levels and this same pattern was also seen for carriers of GG/GG in MT-ND5(13135) and MT-CYB(15257) among healthy and AD individuals ( Figure 3J,K). No statistical significance was observed in the comparisons for MT-ND5(13135) and MT-CYB(15812) ( Figure 3L).

Baseline -TAU Levels Epistasis Network
Genotype influence on TAU levels

Discussion
Mitochondria are organelles crucial for cellular balance and function, being responsible for several pathways and standing out for their role in energy generation through aerobic respiration-the tricarboxylic acid (TCA) cycle and, particularly, oxidative phosphorylation (OXPHOS) [32]. This function is so vital that the human mitochondrial genome (mtgenome or mitogenome) is highly specialized: it contains only 37 genes: 13 encoding proteins, all of which are subunits of the electron transport chain (ETC) complexes where OXPHOS occurs, and the remaining are part of the RNA machinery (22 tRNAs and 2 rRNAs), in addition to non-coding regions, so that all other mitochondrial proteins (over 1000) are encoded by the nucleus [33]. Hence, alterations in these mitochondrial or nuclear genes lead to mitochondrial dysfunction, affecting cellular homeostasis to promote the development of diseases. In this context, neurons would be specially affected, as they consume a large amount of energy, thus being especially dependent on mitochondrial metabolism for energy generation and other functions, such as neurotransmission and neuroplasticity [34,35].
In 2004, a mitochondrial cascade hypothesis was first proposed for late-onset sporadic AD considering that: (i) mitochondrial function decreases with age; (ii) mtDNA damage may be triggered by excessive reactive oxygen species (ROS), which are produced during OXPHOS; and (iii) alterations in mtDNA might reduce ETC efficiency from a basal inherited level [36]. Once the mitochondrial dysfunction reaches a certain threshold, a compensatory response would be induced and some histopathological characteristics of AD would arise as a consequence of this response [37]. Although AD is considered a multifactorial disease, this hypothesis stands out for emphasizing the importance of mtDNA to the development of this type of AD, placing mitochondria in a central position in this process. In addition, as genomic ancestry may play important roles in the development of different diseases, multiple mitochondrial haplogroups have been related to AD in independent studies [38,39].
Notably, in AD, atrophy may occur in the different areas of the brain, but atrophy of the hippocampus (located in the medial temporal lobe) has been traditionally considered a core feature of the disease [40]. In fact, researchers have observed an early divergence of the AD brain model from the normal aging trajectory in the hippocampus, lateral ventricles and amygdala [41]. Interestingly, it has also been demonstrated that volume reduction in AD may start in the medial temporal lobes (although not necessarily the hippocampus) and in the fusiform gyrus at least three years prior to AD progression, spreading to the other lobes before diagnosis [42]. Recently, a study highlighted differences across typical and atypical AD phenotypes in TAU accumulation and atrophy regions, including the lateral temporal lobe [43]. Therefore, the heterogeneous course of this disease should be further explored.
Here, we found DE mitochondrial genes when analyzing RNA-seq data from AD lateral temporal lobe and fusiform gyrus tissues. The lateral temporal lobe (or the lateral surface of the temporal lobe), delimited by superior and inferior temporal sulci and composed of three gyri (superior, middle and inferior temporal gyrus), is responsible for different visual and auditory functions, such as facial recognition, language comprehension and hearing [44,45]. The fusiform gyrus (or occipitotemporal gyrus) is located in the inferior region of the temporal and occipital lobes, being associated with high-level vision functions such as the recognition of faces, bodies and objects, as well as reading [46,47]. In AD, facename memory-the ability to recognize faces and recall names-is markedly impaired, which could be related to degeneration in such brain regions [48]. This, in turn, could be due to molecular alterations such as genetic mutations and different gene expression levels.
Among the DE genes in the studied tissues, five are transfer RNA (tRNA) genes (MT-TH, MT-TL1, MT-TM, MT-TS2 and MT-TV), which are crucial parts of the mitochondrial translational machinery [49]. Considering that the presence of mutations and the altered expression of mitochondrial genes may affect the functioning of mitochondria, especially in muscular and neuronal tissues, and lead to the onset and progress of different diseases [50], the alteration of tRNA genes is of great interest. The downregulated expression of tRNA genes can reduce translation efficiency and is likely to be associated with protein deficiency, which may be linked to the pathogenesis of AD.
Moreover, in our analyses, 10 of the 15 genes found to be DE in the studied tissues are classified as pseudogenes or isoforms of the following mitochondrial genes: MT-ND1 (MTND1P20, MTND1P21 and MTND1P23), MT-ND2 (MTND2P12 and MTND2P28), MT-ND4 (MTND4P9), MT-ND5 (MTND5P11), MT-ND6 (MTND6P3) and MT-RNR2 (MTRNR2L1 and MTRNR2L2). To date, there are not many studies investigating these specific pseudogenes and isoforms in the global literature. In fact, to the best of our knowledge, this is the first study to explore their expression, with the exception of the MT-RNR2 isoforms MTRNR2L1 and MTRNR2L2-and no previous studies were found on the involvement of MTRNR2L1 with any of these genes or the other genetic interactions shown here. Curiously, the isoform MTRNR2L12 has been suggested as a potential biomarker for early AD-like dementia in individuals with Down Syndrome [51], but it was not DE in our study. Regardless, this highlights the potential role mitochondrial that isoforms might have in AD and other types of dementia.
The MT-RNR2 gene encodes the mitochondrial 16S ribosomal RNA (rRNA), and it is also associated with the production of humanin (HN), a peptide that has been shown to be a neuroprotective factor for AD through the suppression of apoptotic cell death when discovered [52] and that, since then, has been associated with different processes and age-related diseases [53]. Humanin presents 13 isoforms encoded by nuclear MT-RNR2-like genes, such as MTRNR2L1 (HN1) and MTRNR2L2 (HN2) [54]. In the last decade, both of these isoforms have been investigated in different multifactorial diseases, although there are few studies so far.
For instance, a genome-wide association study showed a statistically significant relation of MTRNR2L2 with the progression of Huntington's Disease, so the role of this humanin isoform may vary among various diseases [55]. In fact, a recent review on humanin by Hazafa et al. [56] reinforced the importance of this mitochondrial-derived peptide and its isoforms in cytoprotection through the regulation of different mechanisms, including mitochondrial pathways, with a potential influence on the development and treatment of multifactorial diseases related to oxidative stress and apoptosis, which comprise neurodegenerative diseases such as AD.
Here, we found MTRNR2L1 to be underexpressed in both lateral temporal lobe experiments and MTRNR2L2 to be overexpressed in the fusiform gyrus of AD patients. Considering the established neuroprotective function of these humanin isoforms, these results suggest a progression of the disease in the lateral temporal lobe, but not in the fusiform gyrus. Notably, we found MTND1P23 to be DE in all three experiments, being overexpressed in both lateral temporal lobe experiments and underexpressed in the fusiform gyrus, which suggests that MTND1P23 might play an important role in the evolution of AD. This is particularly interesting given that all other genes in the lateral temporal lobe of AD patients in the HYA vs. AD analysis were underexpressed, possibly being protective factors for AD. Similarly, it is also possible to hypothesize, for instance, that the overexpression of MTND2P12 may be indicative of AD progression and that the underexpression of MTND6P3 may be a protective factor against neurodegeneration in AD. However, in the GTEx analysis with healthy brain tissues that did not include the lateral temporal lobe or fusiform gyrus, MTND1P23 was overexpressed in all of the explored tissues. Hence, more studies are needed to clarify the possibilities involving these genes.
To further explore the influence of mitochondrial genetics for AD progression, particularly the occurrence of epistasis, we analyzed the ADNI database with whole mtgenome sequencing of AD, MCI and cognitively healthy subjects, as well as AB, PTAU and TAU levels in CSF from these groups.
Interestingly, we found differences in the CSF levels for the three proteins and epistatic interactions between multiple mtDNA variants for PTAU and TAU. For PTAU, the variant is located in the 12S rRNA gene (MT-RNR1) and interacts with variants in Complex I (MT-ND4) and Complex V (MT-ATP6) genes; however, these interactions do not seem to affect PTAU levels in CSF.
For TAU, two variants are located in MT-DLOOP and one variant is located in a Complex I gene (MT-ND5). Of the MT-DLOOP variants, one (152) interacts with variants in genes of Complexes I, III and IV (MT-ND4, MT-CYB and MT-COX1, respectively), but these interactions do not seem to affect TAU levels in CSF. The other  showed interactions with variants in genes of Complexes I, III, IV and V (MT-ND4, MT-CYB, MT-COX3 and MT-ATP8, respectively); the joint presence of certain genotypes of these variants may affect the levels of TAU in CSF for individuals with MCI, increasing these levels, a trend seen for the investigated AD cases. As for the variant located in the Complex I gene, MT-ND5(13135), it presented interactions with variants in genes encoding Complexes I (one variant in MT-ND2 and one in MT-ND3), III (two variants in MT-CYB) and IV (one variant in MT-COX1; notably, specific genotypes of MT-COX1(7476) and MT-CYB(15257) jointly with the GG genotype for MT-ND5(13135) may affect TAU levels in both cognitively healthy and AD individuals, and the same pattern was observed for MT-ND3(10172) in MCI patients. These findings reflect the intricate network of mtDNA epistasis in the progression of AD.
Currently, only a few studies are found in the global literature on mtDNA epistasis in diseases, and most of them focus on mitonuclear interactions, reinforcing the need for a closer look at this phenomenon specifically within the mitochondrial genome. For instance, a recent study by Duarte-Guterman et al. [57] found sex differences in hippocampal volume and greater memory decline in females compared to males, due to CSF tau-pathology being elevated in female carriers of APOE 4 in the ADNI cohort, corroborating our findings in the sex analysis for the mtDNA epistasis.
Furthermore, a study by Andrews et al. [58] investigated mitonuclear interactions in AD, analyzing associations between mtDNA haplogroups and nuclear-encoded mitochondrial genes for polygenic risk scores of this neurodegenerative disease, and reported both positive and negative epistatic interactions, indicating that epistasis between nuclear and mitochondrial genomes may influence the risk and the age of onset of AD. In this context, to the best of our knowledge, our study is the first to explore epistasis within the mitochondrial genome in AD.

Limitations
We considered p ≤ 0.05 as statistically significant after FDR correction, but we acknowledge the recent movement in the scientific community to lower this traditional threshold to 0.005. Thus, we encourage readers to take this into consideration when interpreting our results. In addition, future studies with an increased sample size and/or functional approaches are needed to further validate the associations suggested here.

Conclusions
In this exploratory study, by employing an integrative analysis with mitochondrial gene expression and genome-wide epistasis approaches, we identified differentially expressed genes in brain tissues from AD patients and epistatic interactions within the mitochondrial genome with a potential influence on the CSF levels of AD-related proteins, revealing important layers of interactions involving the mitochondrial genetics and molecular alterations with a potential impact on the development and progression of AD. Future investigations with larger cohorts or with functional approaches are encouraged to strengthen these findings.  ADNI is a global research effort that actively supports the investigation and development of treatments that slow or stop the progression of AD and subjects have been recruited from over 50 sites across the US and Canada. The overall goal of ADNI is to determine biomarkers for use in Alzheimer's disease clinical treatment trials. To date, it has three phases: ADNI-1, ADNI-GO and ADNI-2, consisting of cognitively normal (CN) individuals, early mild cognitive impairment (EMCI) to late mild cognitive impairment (LMCI) and dementia or AD. For more information, see http://www.adni-info.org, accessed on 4 August 2021. Institutional review board approval was conducted at each ADNI site. Written informed consent was obtained from all participants or authorized representative.