Inflammatory and interferon gene expression signatures in patients with mitochondrial disease

Background: People with mitochondrial disease (MtD) are susceptible to metabolic decompensation and neurological symptom progression in response to an infection. Increasing evidence suggests that mitochondrial dysfunction may cause chronic inflammation, which may promote hyperresponsiveness to pathogens and neurodegeneration. Methods: We collected whole blood from a cohort of MtD patients and healthy controls and performed RNAseq to examine transcriptomic differences. We performed GSEA analyses to compare our findings against existing studies to identify commonly dysregulated pathways. Results: Gene sets involved in inflammatory signaling, including type I interferons, interleukin-1β and antiviral responses, are enriched in MtD patients compared to controls. Monocyte and dendritic cell gene clusters are also enriched in MtD patients, while T cell and B cell gene sets are negatively enriched. The enrichment of antiviral response corresponds with an independent set of MELAS patients, and two mouse models of mtDNA dysfunction. Conclusions: Through the convergence of our results, we demonstrate translational evidence of systemic peripheral inflammation arising from MtD, predominantly through antiviral response gene sets. This provides key evidence linking mitochondrial dysfunction to inflammation, which may contribute to the pathogenesis of primary MtD and other chronic inflammatory disorders associated with mitochondrial dysfunction.


Background
In mitochondrial disease (MtD), a bidirectional relationship between MtD and systemic in ammation emerges, wherein mitochondrial dysfunction may trigger in ammatory cascades, which may then reciprocally contribute to the pathogenesis of MtD. Mouse models have linked primary mitochondrial dysfunction and systemic in ammation. Polg D257A/D257A mutator mice (hereafter: Polg mice), which accumulate mtDNA mutations causing a gradual reduction in mitochondrial respiration (1), display aberrant type I interferon (IFN-I) responses in the innate immune axis leading to immunometabolic dysfunction, accelerated aging, and reduced lifespan (2). The Ndufs4 −/− mouse, a model of neurodegenerative MtD, is also marked by widespread in ammation (3), including increases in serum levels of in ammatory cytokines (IFN-and IL-6), in ammatory markers in the skin and liver, and numbers of activated microglia (4). These responses may be initiated in part through mitochondrial components acting as damage-associated molecular patterns (DAMPs) to activate pattern recognition receptor (PRR) signaling, for example mtDNA activation of the cGAS/STING antiviral response or the NLRP3 in ammasome (5)(6)(7), which can trigger the production and release of IFN-I and interleukin-1β (IL-1β). In primary MtD, mitochondrial dysfunction may cause these pathways may be continuously activated, leading to chronic in ammation.
Chronic in ammation contributes to numerous disorders, including cardiovascular and metabolic disease, cancer, and neurodegeneration (8,9). Neurodegenerative diseases, such as many forms of MtD, may present a uniquely damaging intersection between in ammation and mitochondrial dysfunction.
Pro-in ammatory cytokines, released during systemic in ammation, reach the central nervous system (CNS) via multiple pathways including the blood brain barrier, the choroid plexus, and the vagus nerve (10), and may modulate region-speci c immune cell activation in the brain (11,12), leading to microglial activation, cytotoxicity, and immune dysregulation (10,13). Microglial activation releases cytokines, such as TNF , which impair neuronal mitochondria, causing both oxidative stress and activating additional in ammatory signaling (14). Recent studies have demonstrated that depletion of leukocytes, including microglia, abrogates neuronal death in the Ndufs4 −/− mouse (15,16). Consequently, primary mitochondrial defects may initiate systemic and CNS in ammation, which may contribute to neuronal damage observed in patients with MtD, which manifests clinically as seizures, developmental regression, and degeneration.
To date, most studies on systemic in ammation in MtD have been performed in model organisms. We performed RNAseq on peripheral blood mononuclear cells (PBMCs) in a heterogeneous group of patients with MtD and controls, using gene set enrichment analysis (GSEA) to identify positively and negatively enriched transcriptional signatures in MtD. GSEA has been extensively validated as a method to identify patterns of gene expression with robust biological relevance (17). We compared those RNA signatures with those from transcriptomic studies of mitochondrial encephalomyopathy, lactic acidosis, and strokelike episode (MELAS) patients and two mouse models of mitochondrial dysfunction. Across all four studies we observed enrichment of immune activation and in ammatory gene sets, particularly in antiviral pathways.

Participants
All 81 participants were consented and enrolled in an IRB approved longitudinal natural history study of viral infection and immunity in children with MtD (NIH MINI Study, NCT01780168, www.clinicaltrials.gov) and evaluated at the NIH Clinical Center. Participant characteristics are shown in Table S1. The diagnosis of MtD was made by the referring provider (i.e., neurologist, clinical geneticist), and Modi ed Walker criteria score of "probable" (P) or "de nite" (D) was assigned. The mean age of the control and MtD cohorts were 14.2 (Std dev = 10.6) and 18.4 (Std dev = 16.1) years of age (P = 0.15), respectively. Molecular testing was available for 30 out of 32 patients (94%) with MtD.
Transcriptomic analysis PBMC were isolated from 5mL whole blood using LeucoSep tubes (Greiner Bio-one) and Ficoll-Paque Plus (GE Healthcare) for density gradient centrifugation, before lysis in TRIzol (Thermo Fisher, Waltham, MA). For RNA extraction samples were batched according to their age, gender and case/control status, and two reference samples were simultaneously processed with each batch. Total RNA was isolated and puri ed with miRNeasy kit (Qiagen, Hilden, Germany), with RNA quality and quantity estimated using Nanodrop (Thermo Scienti c, Wilmington, DE) and Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA). Stranded cDNA sequencing libraries were generated with TruSeq Stranded mRNA Library Prep Kit (Illumina, San Diego, CA) following the manufacturer's instructions. Brie y 500 ng of total RNA was used for mRNA selection. After the reverse transcription to 1st strand cDNA, strand info was reserved with dUTPs during 2nd strand synthesis. The dsDNA fragments then had the addition of a single 'A' base and subsequent ligation of the adapter. The products were then puri ed and enriched with PCR to create the nal cDNA library. The library was quali ed with Agilent Bioanalyzer and quanti ed with Qubit 2.0 uorometer. The cluster generation and paired-end (2x75 bp) sequencing was run on Illumina HiSeq 3000 at NHLBI Sequencing Core. All 81 barcoded samples were pooled for one single run, which yielded at least 25M passed lter paired reads per sample.

Computational analysis
Sequence alignment, quality control ltering, and count matrix generation were performed using STAR (18), QoRTs (19), and RSEM (20) on the NIH HPC Biowulf cluster. All subsequent statistical analysis and graphical presentations were performed in R (https://cran.r-project.org/). RSEM-corrected transcript counts were imported using tximport (21), and differentially expressed genes were identi ed using DESeq2 (22). Preliminary gene set enrichment analysis (GSEA) was performed using clusterPro ler (23) and visualized using GOplot (24). Subsequent GSEA analysis used the fgsea (25) package and blood transcription module (BTM) gene sets (26). Each BTM is a set of genes, which has been shown to show coherent expression across many biological samples (27,28). Gene set variation analysis (GSVA) was used to quantify participant level variation in signatures and test correlation with participant age (29) ( Figure S1C). GEO2R and GEOquery (30) were used to perform differential expression analysis from gene expression microarray data in LIMMA (31), whereas additional RNAseq data was analyzed with DESeq2.

Results
Mouse models of mitochondrial dysfunction have revealed immune signatures marked by elevated in ammation and interferon responses (2,3,7,15). To understand whether similar transcriptomic signatures occur in people with MtD (Appendix Table S1), we performed bulk RNAseq on PBMCs. All participants were in their usual state of health at the time of collection and did not display any symptoms or signs of infection. We examined sample variance using a sample distance correlation matrix of all control participants and MtD patients ( Figure S1A). Through this analysis, we identi ed two females (one control, one MtD) that segregated from other samples, and excluded them from subsequent analysis.
Using PCA, we observed that our samples do not cluster into diagnosis group ( Figure S1B). Differential expression analysis identi ed differentially expressed genes (DEGs) between control and MtD groups (Appendix Table S2). With a threshold of p < 0.1 and log2 fold change (log2FC) of < |0.5|), we detected few DEGs including CXCL2 (Fig. 1A). By ranking the top 50 DEGs by t statistic, we found that the differences in relative expression between these groups was able to improve clustering of control and MtD groups ( Figure S1C).
To identify patterns of differential gene expression, we performed GSEA using Gene Ontology (GO) categories. This analysis revealed signi cant positive and negative enrichment of 172 gene sets (adjusted p < 0.05) (Appendix Table S3). In the six most signi cant positive and negatively enriched categories, we found decreased expression of genes in ribosomal, mitochondrial, B cell and natural killer (NK cell) categories, and increased expression of genes related to stimulus response, ion channel and G protein-coupled receptor (GPCR) activation, pattern recognition receptors (PRR), and interleukin-1β (IL-1β) production (Fig. 1B). Positively enriched gene sets included many involved in immune activation and signaling, including IL-1 and IL-1β production (n = 6), PRR and Toll-like receptor (TLR) binding (n = 2), IFNβ production (n = 2) and viral response regulation (n = 2). Negatively enriched gene sets included mitochondrial proteins and complexes (n = 28), ribosomes and translation (n = 30), natural killer (NK) cell (n = 9), B cell (n = 3), and major histocompatibility complex (MHC) activity (n = 8) (Fig. 1C).
As these results suggested immune and translational dysregulation in our MtD patients, we performed GSEA using blood transcriptional modules (BTMs) to identify targeted pathway enrichment (26). We identi ed 62 signi cant gene sets (Appendix Table S4). Examining the consensus pathways from these modules, we observed a positive enrichment in monocyte, TLR, IFN, and immune activation clusters, and negative enrichment in mitochondria, transcription and translation, and NK and T cell clusters ( Figure S2). We also found positive enrichment of in ammatory response, type I IFN response, and innate antiviral response modules, and negative enrichment of plasma cell/immunoglobulins and B cell modules trending toward signi cance (p ≤ 0.076).
We compared our ndings against three others: whole blood from MELAS patients (32), bone marrowderived macrophages (BMDMs) from Polg mutator mice (2), and mouse embryonic broblasts from Tfam +/− mice (7). We performed BTM GSEA on these sample sets and compared the enrichment of the previously identi ed MtD-signi cant modules ( Fig. 2A). Examining mitochondrial clusters con rms previously observed trends -MELAS subjects and Polg mice have a positive enrichment, while MtD patients and Tfam +/− MEFs have a negative enrichment. We discovered key overlaps between MtD patients and the previous MELAS study, including negative enrichment in NK cell and T cell modules. Positive enrichment of monocyte, dendritic cells, and neutrophils are unique to MtD patients, though immune and in ammatory cluster enrichment is shared by MtD and MELAS patients. While these modules were not signi cant in MELAS patients, related viral sensing and IRF2 targets clusters were (Appendix Table S4). Putative PAX3 target gene sets were positively enriched across all four datasets. We examined the genes that compose antiviral and RIG-I like receptor (RLR) (IFN-related), IRF2 targets, and putative PAX3 target modules (Fig. 2B). In IFN gene sets, we observed increased expression of core cluster genes PML, RARA, IFIT1, OAS3, IFIH1, DHX58, IRF7, and TRIM25 in all four datasets. In IRF2 targets, which were signi cantly enriched in MELAS subjects, Polg mice, and Tfam +/− MEFs, USP18, TAP1, BST2, IFI35 were consistently increased. Finally, in putative PAX3 targets, NR4A1, GEM, ATF3, and DUSP1 were increased in all four datasets.
To nd explainable sources of variance within our dataset, we considered the sex and age of our participants, as PCA plots suggested that a proportion of variance could be explained by sex ( Figure  S1B). We assessed the contribution of age by performing gene set variation analysis (GSVA) by evaluating a subset of BTMs across our sample population against age ( Figure S3). We did not observe a signi cant trend between age and diagnosis group. We divided our samples based on sex and reperformed differential expression analysis. In male samples, we continued to observe poor withindiagnosis sample correlation ( Figure S4A) and observed that a signi cant proportion of variance remained unexplained ( Figure S4B). We observed similar trends in female samples ( Figure S4C and S4D).
In males, we observed that only four DEGs crossed a log2FC > |0.5| and p < 0.1 ( Figure S5A and Table S2). Nonetheless, two downregulated genes, FCRL6 and KLRC4-KLRK1, are linked to cytotoxic T and/or NK cells. Ranking DEGs by t-statistic and examining the relative expression of the top 25 genes improved clustering between control and MtD ( Figure S5B). In females, no genes were signi cantly differentially expressed ( Figure S5C and Table S2), but we observed within-diagnosis group clustering when examining the t-statistic ranked top 25 DEGs ( Figure S5D).
In males, GSEA of GO terms identi ed 87 signi cant gene sets, including negative enrichment of gene sets relating to ribosomes (n = 12), mitochondria (n = 23), MHCs (n = 2), NK cells (n = 9) and T cells (n = 3) ( Fig. 3A and Table S3). In females, signi cant GO GSEA results (n = 78) revealed similar negatively enriched gene sets, including ribosomes (n = 21), mitochondria (n = 3), and B cells (n = 5), but also positively enriched gene sets including IL-1β production (n = 6) and endothelial morphogenesis (n = 2) ( Fig. 3B and Table S3). Based on these disparate results, we compared male and female GSEA using BTMs. While there were many overlapping signi cantly enriched categories (Appendix Table S4), we observed a few intriguing differences between the two (Fig. 3C). TLR and in ammatory signaling and monocyte gene sets are more positively enriched in female MtD samples than males, and IRF2 targets are signi cantly enriched in females. In contrast, several B cell gene sets are signi cantly negatively enriched in females, whereas these sets are either nonsigni cant or positively enriched in males. Finally, we observed that while NK and T cell gene sets are negatively enriched in males and females, these sets were more signi cant in males.

Discussion
Our study is among the rst to examine the transcriptomic pro le of PBMCs from a cohort of MtD patients. We used GSEA to identify key pathway perturbations, and, using both GO terms and validated BTMs, we observed consistent set enrichments. We observed positive enrichment of IL-1β, IFN, and TLR/PRR GO terms, and monocyte, TLR, and IFN signaling BTMs, suggesting elevated basal levels of in ammation in MtD patients. Conversely, we observed negative enrichment of ribosomal proteins, cell killing, and mitochondrial GO terms, and T cell, NK cell, B cell, and transcription and translation BTMs, suggesting suppression of immune system activity coupled with de cits in translation. These data are a valuable resource to inform future study of the mechanisms of peripheral in ammation in patients with MtD.
While we observed few signi cant DEGs, our sample population is highly heterogenous. We included male and female, adult and pediatric patients and controls, and patients with multiple MtD diseasecausing variants, including both mtDNA and nDNA mutations. Few signi cant DEGs were detected in the MELAS study, despite restricting to adult patients with the same m.A3243G variant (32). Further, this distinction in study populations likely underlies the oppositional enrichment between MtD and MELAS groups in mitochondrial gene sets. The MELAS study found a positive enrichment of mitochondrial genes, and our con rmation of this nding with our BTM GSEA approach is an important validation of our technique. The signi cant negative enrichment of mitochondrial genes among our MtD patients may re ect the differing genetic origins of MtD.
Our dataset allowed us to perform a preliminary analysis of male and female patients with MtD. While sex did not drive all variance observed, we found interesting differences between males and females. In males, T cell and NK cell GO terms were signi cantly negatively enriched. In females, B cell and immunoglobulin GO terms were negatively enriched, while IL-1β and colin-1 gene sets were positively enriched. BTM terms furthered these ndings, demonstrating oppositional enrichment of B cell gene sets between males and females, stronger negative enrichment of T cell and NK cell sets in males, and stronger positive enrichment of monocyte and TLR sets in females. This intriguing segregation of enriched pathways suggests potential sex-speci c manifestations of immune dysfunction in MtD.
Across all four datasets, we observed a positive enrichment of antiviral IFN or viral sensing signaling modules, supporting previous ndings. In Tfam +/− MEFs, a pathway was identi ed by which mitochondrial stress promotes the release of mtDNA into the cytosol, activating the cGAS-STING-IRF3 antiviral response (7). In Polg BMDMs, IFN-I signaling genes are basally elevated (2). Together, these demonstrated increased IFN signaling and in ammatory signaling in multiple mouse models of mitochondrial dysfunction. We found signi cant enrichment of similar pathways in our MtD patients and in MELAS, suggesting basal IFN and antiviral pathway hyperactivation may be common mechanisms in MtD. We identi ed signi cant enrichment of an "innate activation by cytosolic DNA sensing" module in Polg and Tfam datasets, but this module was not enriched in either MtD or MELAS. Consequently, while we con rm enrichment of IFN and in ammatory pathways in human MtD, we emphasize that additional research will be needed to characterize underlying mechanisms.

Conclusions
Using multiple GSEA approaches, we demonstrate that MtD patients have positively enriched IFN, in ammatory signaling, and monocyte gene sets, and negatively enriched NK cell, T cell, and ribosomal gene sets. Our approach con rms the enrichment of type I IFN signaling in two mouse models of mitochondrial dysfunction and detects a previously unreported enrichment of viral sensing genes in a study of MELAS patients. MtD and mitochondrial dysfunction may induce basal peripheral in ammation, which likely contributes to the susceptibility to elevated in ammatory responses and sepsis in people with MtD. Further, peripheral in ammation exacerbates neuroin ammation and neurodegeneration in animal models (33,34). This may point to an important mechanism underlying the acceleration of neurodegeneration following infection in MtD, and the exacerbation of neurodegeneration in disorders with mitochondrial dysfunction. illustrates the GO functional annotation of each gene. (C) Bubble plot of signi cant GO terms following GSEA analysis and reduction to consolidate redundant terms. Bubble size re ects the size of the gene set, while bubble color re ects the GO category of each set. Z score is derived from the comparison of upregulated and downregulated genes in the set. Yellow line indicates the signi cance threshold of adjusted p < 0.05. For clarity, only a subset of positive and negative enrichment sets is labeled (see Table   S3 for a full list).  Yellow line indicates the signi cance threshold of adjusted p < 0.05. See Table S3 for a full list of signi cant GO terms. (C) NES plots from selected modules from males (M) and females (F). Bubble size re ects adjusted p value, bubble color re ects NES. Modules are grouped into categories shown at right. IFN; interferon, Inf/TLR; In ammation and Toll-like receptor signaling. See Table S4