key: cord-0898723-ni7ck5zp
authors: Medini, Hadar; Zirman, Amit; Mishmar, Dan
title: Immune system cells from COVID-19 patients display compromised mitochondrial-nuclear expression co-regulation and rewiring towards glycolysis
date: 2021-11-18
journal: iScience
DOI: 10.1016/j.isci.2021.103471
sha: d8a26e99ef55cc0889e2cf6311c311c0df0f830a
doc_id: 898723
cord_uid: ni7ck5zp

Mitochondria are pivotal for bioenergetics, as well as in cellular response to viral infections. Nevertheless, their role in COVID-19 was largely overlooked. Here, we analyzed available bulk RNA-seq datasets from COVID-19 patients and corresponding healthy controls (three blood datasets, N=48 healthy, 119 patients; two respiratory tract datasets, N=157 healthy, 524 patients). We found significantly reduced mtDNA gene expression in blood, but not in respiratory tract samples from patient. Next, analysis of eight single-cells RNA-seq datasets from peripheral blood mononuclear cells, nasopharyngeal samples and Broncho alveolar lavage fluid (N=1,192,243 cells), revealed significantly reduced mtDNA gene expression especially in immune system cells from patients. This associated with elevated expression of nuclear DNA-encoded OXPHOS subunits, suggesting compromised mitochondrial-nuclear co-regulation. This, together with elevated expression of ROS-response genes and glycolysis enzymes in patients, suggest rewiring towards glycolysis, thus generating beneficial conditions for SARS-CoV-2 replication. Our findings underline the centrality of mitochondrial dysfunction in COVID-19.

A significant number of COVID-19 patients develop severe fatal consequences of the respiratory tract infection, attributed to catastrophic inflammatory events, described jointly as a 'cytokine storm'. This inflammatory state associates with deleterious oxidative stress, and elevated production of reactive oxygen species (ROS)-both point towards altered mitochondrial activity (Singh et al. 2020) . This is further supported by evidence suggesting a central role for the mitochondria in modulating the activity of the innate immune system (Bordon 2018 ; Zuo and Wan 2019; Ganji and Reddy 2020).

Specifically, once viral genomes are recognized by the innate immune system, the retinoic acidinducible gene I (RIG-I)-like receptors (RLRs), RIG-I and MDA5, interact with mitochondrial antiviralsignalling protein (MAVS). MAVS, in turn, activates the intracellular signalling cascade that induce the transcription of genes encoding type I interferons (Rehwinkel and Gack 2020) .

Several lines of evidence underline the mitochondria as a pivotal cellular target during SARS-CoV-2 infection: Firstly, the discovery of high sequence similarity between SARS-CoV-2-encoded proteins and mitochondrial-localized proteins encoded by SARS-CoV-1 (3/29 protein coding genes) (Shi et al. 2014 ; Gordon et al. 2020) , which were shown to impact mitochondrial gene expression (Shao et al. 2006 ) and components of the innate immunity (Burtscher et al. 2020) . Secondly, recent extensive computational analyses predicted preferential mitochondrial localization of SARS-CoV-2 RNAs and sgRNAs (Wu et al. 2020 ) and experimental evidence revealed interaction of SARS-CoV-2-encoded proteins with the mitochondria (Gordon et al. 2020) . Third, recent evidence revealed alteration of mitochondrial dynamics in COVID-19 patients, which associated with increased fusion, and inhibition of mitochondrial fission (Holder and Reddy 2021) . It is thus likely, that mitochondrial function is hijacked, and altered, upon SARS-CoV-2 infection (Gatti et al. 2020; Singh et al. 2020) . Despite the known regulatory involvement of the mitochondria in the immune system, the regulatory response of the mitochondria to SARS-CoV-2 infection has been largely overlooked. Recent RNA-seq analysis in airway clinical samples and in lung-associated cell lines suggested little or no response of mitochondrial gene expression, either in the mitochondrial genome (mtDNA) or in the nucleus, to SARS-CoV-2 infection (Miller et al. 2021 ). This finding, which was reported during the course of our work, may lead to the conclusion that mitochondrial gene expression and subsequent mitochondrial activities are not significantly affected in COVID-19 patients. However, such findings do not dismiss possible tissue/cell type specific regulatory response of the mitochondria to SARS-CoV-2 infection. J o u r n a l P r e -p r o o f 4 Here, by analysing both bulk and single cells RNA-seq data in COVID-19 patients as compared to controls we discovered significant reduction in mtDNA gene expression levels in patient samples from blood, in contrast to inconsistent changes in samples from the respiratory tract. In consistence with this finding single cells RNA-seq analysis revealed that such altered mtDNA gene expression preferentially occurred in immune system cells, whereas respiratory tract epithelial cells tended to display less prominent alterations. The response of mitochondrial-related pathways in the nucleus (such as ROS generation, TCA, glycolysis) is discussed in light of the apparent dependence of SARS-CoV-2 replication on glycolysis (Codo et al. 2020) , likely on the expense of the mitochondrial oxidative phosphorylation. Our findings underline a tissue-and cell-type dependent negative response of mitochondrial gene expression regulation in COVID-19 patients.

Mitochondrial gene expression in COVID-19 patients is significantly reduced in peripheral blood, but not in the upper respiratory tract. To examine whether SARS-CoV-2 infection associates with altered mitochondrial regulation, we assessed mtDNA-encoded genes' expression in thirteen publicly available datasets from healthy and COVID-19 patients: five bulk RNA-seq datasets from nasooropharyngeal (NP/OP) swab (Lieberman et al. 2020; Mick et al. 2020) , peripheral blood and whole blood (Bernardes et al. 2020; Levy et al. 2021; Thair et al. 2021) , as well as eight scRNA-seq datasets from peripheral blood mononuclear cells (PBMC) (Bernardes et al. 2020; Wilk et al. 2020; Ren et al. 2021; Stephenson et al. 2021; Wilk et al. 2021) , Broncho alveolar lavage fluid samples (BALF) (Liao et al. 2020 ) and nasopharyngeal (NP) samples (Chua et al. 2020; Ziegler et al. 2021) (Fig 1, Table 1 -notice dataset numbering, used below).

In consistence with previous findings, differential expression analysis of mtDNA genes in bulk RNA-Seq from the upper airway in healthy individuals as compared to patients who were positively diagnosed for SARS-CoV-2 (Datasets I, IV) did not show any significant difference in mtDNA gene expression in the patients (Miller et al. 2021 ) (N=524 SARS-CoV-2 patients, N=157 healthy control; see Table 1 in Methods) (Fig S1) . In contrast, analysis of three independent bulk RNA-seq datasets from peripheral blood of healthy individuals and COVID-19 patients (N=34 healthy individuals and N=106 COVID-19 J o u r n a l P r e -p r o o f 5 patients; see Table 1 in Methods), revealed significant reduction in the expression of 8 / 13 and 11 / 13 mtDNA-encoded protein-coding transcripts, respectively (Fig 2, Fig S2, Table S1 ; Datasets II, V ).

We next analysed bulk RNA-seq data (Dataset III) collected from 14 healthy individuals, samples from 13 hospitalized COVID-19 patients, each having up to five samples collected during disease progression, termed pseudotime points (Bernardes et al. 2020 ). Since the sample size of pseudotime 7 (recovery long-term follow-up) was very small (N=2), it was excluded from further analysis. The pseudotime points were originally defined based on inflammatory markers and ventilation requirements as follows: pseudotimes 1,2 and 3 were considered 'severe', pseudotimes 4 and 5 were considered early/moderate and late convalesces, respectively, and pseudotime 6 were from recovered individuals (see Methods). Notably, to avoid minute sample sizes in bulk RNA-seq of Dataset III, we grouped pseudotimes 1, 2 and 3 (all originally defined as 'severe'). Pairwise comparison of gene expression between control and each of the remaining disease trajectory phases revealed significant reduction in the expression of all mtDNA protein-coding transcripts in the severe, moderate/early convalescent as well as in late convalescent individuals ( Fig 2B, Table S1 ).

Interestingly, no significant change was observed in mtDNA protein-coding transcript levels while comparing healthy controls to recovered individuals (pseudotime 6; Fig 2B) . Hence, our results suggest that mtDNA gene expression is consistently reduced in patients' blood, but not in the respiratory tract, and that such altered expression is likely reversible upon recovery.

Expression of nuclear DNA-encoded OXPHOS, ROS formation, TCA, glycolysis and mitochondriarelated viral response genes is generally elevated in COVID-19 patients. As positive co-expression of mitochondrial and nuclear DNA-encoded OXPHOS genes has been observed across many human tissues (Barshad et al. 2018) we asked whether the reduced mtDNA gene expression in COVID-19 patients correlate with changes in nuclear DNA-encoded OXPHOS genes and related pathways. To this end, we subjected Datasets II, III and V to assessment of differentially expressed genes involved in the following biochemical pathways: OXPHOS structural and assembly genes, mitochondrial ribosome (Wolf and Mootha 2014; Barshad et al. 2018) . We also analysed the following mitochondria-related Gene Ontology (GO) terms: tricarboxylic acid (TCA), glycolysis enzymes including L-lactate dehydrogenase, response to type I interferons, mitochondrial antiviral signalling , including genes involve in RLR-MAVS pathway (hereby annotated as MAVS pathway; Table S2 ) (Burtscher et al. 2020; Rehwinkel and Gack 2020) , generation and response to reactive oxygen species (ROS) formation, J o u r n a l P r e -p r o o f 6 along with ROS scavenging enzymes. Notably, since measurement of the average expression of a given pathway may mask the impact of SARS-CoV-2 infection on individual genes, we performed a gene-by-gene differential expression analysis (while stringently correcting for multiple testing) ( Fig   3A) . We screened for differentially expressed genes that were consistently significant in Datasets II, III and V, and considered only genes whose responses were consistent with at least one disease condition from Dataset III. Notably, while considering Dataset III, we focused on conditions that displayed significantly reduced mtDNA genes expression analysis. The analysis revealed a tendency towards elevated expression among OXPHOS subunits (7 elevated genes), TCA (6 elevated genes) and glycolysis (19 elevated genes, 6 reduced genes) in Datasets II and V and in severe, moderate/early and late convalescent individuals from Dataset III ( Fig 3A; Table S3 ). Specifically, while considering glycolysis, expression levels of a key glycolysis enzyme -glyceraldehyde-3-phosphate dehydrogenase (GAPDH), were elevated in Datasets II and V as well as in severe and moderate/early convalescent patients of Dataset III ( Fig 3B) . It is worth noting that gene expression analysis of mitochondrial ribosome subunits revealed a mixed response in patients, with only few statistically significant values (Table S3) .

Notably, we found that the expression of two subunits of lactate dehydrogenase (LDH), a key enzyme that regulates pyruvate generation by glycolysis, responded differently in the peripheral blood samples of patients versus controls, as follows: The expression of LDHA, a subunit which is more abundant in glycolytic tissues (Read et al. 2001; Porporato et al. 2011) , was elevated in patients from Datasets II and V and in severe, moderate/early and late convalescent of Dataset III ( Fig 3B) . In contrast, the expression of LDHB, a subunit that is more abundant in high OXPHOS tissues (Read et al. 2001; Porporato et al. 2011) , was reduced in patients from Datasets II and V and in severe, moderate/early convalescent samples of Dataset III ( Fig 3B) . As LDH activity governs the choice between generation of pyruvate (which promotes OXPHOS function) versus lactate, these results further support the interpretation that OXPHOS function is likely compromised and rewired to glycolysis in patients versus controls.

Lastly, genes that participate in type I interferon pathway (likely a reflection of the so-called mitochondria-related cytokine storm), genes involved in MAVS pathway, ROS formation and ROS response, mostly displayed elevated expression in COVID-19 patients in Dataset II, V and in severe, moderate/early and late convalescent samples of Dataset III (Fig 3A) . Specifically, there were 10 J o u r n a l P r e -p r o o f 7 elevated genes in IFN1 pathway, and three genes involved in MAVS pathway. Moreover, the expression of retinoic acid-inducible gene I (RIG-I), melanoma differentiation-associated protein 5 (MDA5) (DDX58 and IFIH1, respectively) and TANK-binding kinase 1 (TBK1), was consistently elevated in COVID-19 patients ( Fig 3A, Table S3 ), supporting the increase in type I interferon pathway. Next, most of the significant genes from ROS formation pathway (30 out of 34) displayed elevated expression in patients; similarly, most ROS response pathway genes (49 out of 66) displayed elevated expression as well. Notably, the expression of certain ROS scavenging enzymes (e.g., Peroxiredoxin 3 (PRDX3) and superoxide dismutase 2 (SOD2) (Yoboue et al. 2018) ) was also significantly elevated in COVID-19 patients in 2/3 Datasets (Table S3) . Taken together, these findings suggest that along with the apparent rewiring towards glycolysis, COVID-19 patients also suffer from possible oxidative stress.

We hypothesized that the reduction of mtDNA gene expression levels in COVID-19 patients could be regulated, rather than spontaneous. To identify candidate regulatory factors that potentially explain such putative downregulation, we assessed differential expression of a set of known and candidate mtDNA regulatory factors of transcription, regulatory factors of mtDNA replication, and nuclear DNAencoded factors with known mitochondrial RNA-binding activity (Wolf and Mootha 2014; Cohen et al. 2016 ), as well as RNA and DNA-binding proteins that were recently identified in human mitochondria (Ardail et al. 1993; Fernandez-Vizarra et al. 2008; She et al. 2011; Blumberg et al. 2014; Lambertini et al. 2015; Chatterjee et al. 2016) (Table S2 ). Our analysis revealed that in all peripheral and whole blood bulk RNA-seq datasets (Datasets II, III and V), the expression of JUN (c-Jun) was significantly elevated while the expression levels of JunD and of mitochondrial RNA polymerase (POLRMT) were reduced especially in peripheral blood datasets (Datasets II, III) ( Fig 3C, Fig S2, Table S3 ). The reduction in POLRMT expression in patients, the only known mitochondrial RNA polymerase, provides an attractive explanation to the reduction in mtDNA genes expression ( Fig 3C, Fig S2) . The altered expression of c-Jun and JunD is intriguing -these are known nuclear gene expression regulators that were also found to bind human mtDNA in vivo (Blumberg et al. 2014) . Our results suggest that these are candidate regulators of gene expression in both the nuclear and mitochondrial genomes (see Discussion). patients, especially in immune system cells. Our above-described analysis of bulk RNA-seq data suggest that the response of mitochondrial genes' expression in patients' peripheral blood cells is more consistent than cells coming from the respiratory tract. As blood cells are enriched by immune system cell types, it is possible that cells belonging to the immune system went through a more prominent response of mtDNA gene expression to COVID-19, than other cell types. To address this possibility, we analysed eight single cells RNA-seq (scRNA-seq) datasets including a total of 1,192,243 high-quality and annotated cells (Table 1, Table S4 ). In total, we analysed cells from 77 critical/severe COVID-19 patients, 65 COVID-19 patients with a moderate phenotype, 65 healthy individuals and scRNA-seq data from the patients studied in Dataset III (above analysed also for bulk RNA-seq in peripheral blood) (see Table 1 in Methods). To avoid technical noise due to known elevated proportions of zero read counts in scRNA-seq data (Vallejos et al. 2017) , we included in our analysis only cells whose sequencing read counts were higher than zero for each of our analysed mtDNA genes (Table S4) . Notably, since scRNA-seq libraries were prepared while enriching for polyA mRNAs, we focused our analysis on nine mtDNA genes with known longer 3' polyA (ND1-4, CO1-3, CytB and ATP6) (Slomovic et al. 2005) , and hence reduced the false discovery rates of lower mtDNA gene expression in certain genes. Finally, we considered a finding biologically meaningful if it was observed in at least two independent datasets.

Firstly, in consistence with our analysis of bulk RNA-seq from peripheral blood, significantly reduced mtDNA-encoded genes' expression was observed in several T cells sub-types (Table S1 , Data S1-S4 http://dx.doi.org/10.17632/8kd3xjfrh4.1.). Specifically, reduced mtDNA gene expression was consistently observed in CD8+ T cells from patients in BALF samples (Dataset VI) and from CD8 T cells, including CD8 memory T cells, from PBMC (Datasets VII, XI) ( Fig 4A, Data S4 , Table S1 http://dx.doi.org/10.17632/8kd3xjfrh4.1.). A smaller effect on mtDNA gene expression (measured by the fraction of mtDNA genes which displayed a significant change in patients), yet still reduction in expression, was observed in CD8 effector cells from PBMC (Datasets VII, XII) (Fig 4B, Data S1, Table S1 http://dx.doi.org/10.17632/8kd3xjfrh4.1.). Similarly, significantly reduced mtDNA gene expression in patients was also observed in CD4 T cells from 2 / 5 PBMC datasets, CD4 memory T cells (Dataset VII) ( Fig 4B) and additional sub-types of T cells that were represented only in one dataset each, such as CCR7+ T cells and Follicular helper CD4 T (Dataset VI and Dataset XII, respectively) (Table S1 ).

J o u r n a l P r e -p r o o f 9 Significantly reduced expression of most mtDNA-encoded genes was also observed in patients' NK cells in BALF (Fig 4A) , as well as in NK cells from PBMC samples (3 / 4 datasets) regardless of their severity level (Data S1, Fig 4B, Table S1 http://dx.doi.org/10.17632/8kd3xjfrh4.1.). It is worth noting that in both CD4 T and NK cells we observed a significant reduction in mtDNA gene expression patients was significantly elevated only in one out of two datasets (Data S3, Data S5, Data S7, Table S1 http://dx.doi.org/10.17632/8kd3xjfrh4.1.) and proliferating T cells and Macrophages displayed a mixed tendency in mtDNA gene expression in patients (Data S2, Table S1 http://dx.doi.org/10.17632/8kd3xjfrh4.1.). As currently we consider significance only for results that were replicated by at least two datasets, the biological importance of the mixed response should be evaluated once additional RNA-seq datasets become available. Taken together our analysis strongly suggest a consistent reduction in mtDNA gene expression in COVID-19, especially in immune system cells.

Disease severity does not clearly impact the reduction in mtDNA gene expression in patients. While considering Datasets IV and IX (NP), division into individual cell types belonging to the immune system led to cell sample sizes below our threshold for further analysis. To overcome this problem, we grouped all cell types belonging to the immune system for further comparisons of differential expression between patients and controls; this analysis revealed consistent and significant reduction in mtDNA genes' expression in patients (Data S8, Data S9 http://dx.doi.org/10.17632/8kd3xjfrh4.1.).

Finally, we took advantage of available expression data for B-cells and monocytes from scRNA-seq in the frame of Dataset III. Such analysis revealed significantly reduced expression of most mtDNA analysed genes in CD14+ monocytes (pseudotime points 1,5), in CD16+ monocytes (mainly in pseudotime points 1,4 and 5), and in B cells (mainly in pseudotime points 1,4) (Data S6 http://dx.doi.org/10.17632/8kd3xjfrh4.1.). Notably, while comparing mtDNA gene expression of the immune system cell types between the pseudo time points (Dataset III-scRNA-seq), and between the severity levels (e.g. moderate and critical) in Datasets VI-XII (Table S1) (Table S1) Table S1 http://dx.doi.org/10.17632/8kd3xjfrh4.1.). Finally, IFNG-responsive epithelial cells (IRC), displayed a significant reduction in the expression of only a single mtDNA gene (ATP6) (Data S8 http://dx.doi.org/10.17632/8kd3xjfrh4.1.). These results suggest that whereas cells belonging to the immune system consistently show reduction in mtDNA gene expression, whether isolated from the peripheral blood or from the respiratory tract, such response is weaker, and notably varies among respiratory tract epithelial cells.

Analysis of scRNA-seq databases identify candidate regulators that explain reduced mtDNA gene expression in COVID-19 patients. We next asked whether reduction in mtDNA gene expression in certain cell types from patients associated with expression changes in genes with mitochondrial function encoded by the nuclear genome. To address this question, we performed differential expression analysis of OXPHOS genes and mitochondrial-related pathways that were analysed in the bulk RNA-seq datasets (see above). We found that in cell types which displayed reduction in the J o u r n a l P r e -p r o o f 11 expression of most mtDNA genes, the significant expression of most nuclear DNA-encoded OXPHOS genes showed elevated expression (in all tested Datasets) (Table S3 , Data S10 http://dx.doi.org/10.17632/8kd3xjfrh4.1.). Similar to the analysis of bulk RNA-seq, our analysis of mitochondria-related pathways revealed that such cells displayed mostly elevated expression of genes belonging to glycolysis, response to type I interferon and MAVS pathway genes (Data S10 http://dx.doi.org/10.17632/8kd3xjfrh4.1.).

To further test our hypothesis that the reduced mtDNA gene expression in certain cell types from patient samples reflect downregulation of mtDNA gene expression upon SARS-CoV-2 infection, we assessed changes in regulatory factors of mtDNA genes' expression (Table S3) 

Our results indicate a tissue-and cell-type-dependent reduction in mtDNA-encoded gene expression in response to SARS-CoV-2 infection. Firstly, bulk RNA-seq analysis of samples coming from the peripheral and whole blood showed such reduction, in contrast to samples coming from the respiratory tract that did not. Secondly, analysis of scRNA-seq revealed reduced mtDNA gene expression levels in several tested immune system cells, as compared to little or no reduction of mtDNA gene expression in respiratory tract epithelial cells. These findings support cell-type and tissue-dependent differences in the magnitude of mitochondrial gene expression regulatory effect in COVID-19 patients. We are tempted to interpret these results, as partially explaining the apparent J o u r n a l P r e -p r o o f 12 lack of change in mtDNA gene expression in patients' respiratory tract bulk RNA-seq samples: the weaker mtDNA gene expression differences in epithelial cells from patients (as well as in some immune system cells isolated from the respiratory tract) might have masked the difference in mtDNA gene expression in the respiratory tract, but not in blood. Since two out of the five analyzed bulk RNAseq datasets stem from different sequencing platforms (Datasets IV, V), the differences in the impact on mtDNA gene expression could be attributed to the different tissues of origin, although variation in sample collection methods and/or sequencing platforms cannot be excluded. Another parameter that may affect gene expression differences in the mitochondria, is mtDNA genetic backgrounds (haplogroups) (Cohen et al. 2016) . It would therefore be of interest to study the contribution of haplogroups to gene expression differences between healthy and COVID-19 samples once larger sample sizes become available.

We hypothesized that reduction in mtDNA gene expression in samples from COVID-19 patients stem from a general alteration in OXPHOS gene co-regulation. The general tendency towards increased expression in nuclear DNA-encoded OXPHOS genes in COVID-19 patients supports our hypothesis, as by-and-large healthy human tissues display positive co-expression of mitochondrial and nuclear DNAencoded members of the OXPHOS machinery (Barshad et al. 2018 ). This leads to two possible interpretations: the altered coordination of mito-nuclear gene expression in COVID-19 patients may (A) reflect a compromised mito-nuclear co-regulation, and/or (B) controlled rewiring of the OXPHOS machinery to glycolysis. Indeed, we found that COVID-19 patients exhibited with a general increase in the expression of genes encoding glycolysis enzymes, reduced expression of LDHB and increased expression of LDHA. This is consistent the observed dependence of SARS-CoV-2 replication in monocytes on active glycolysis (Codo et al. 2020) , suggesting that OXPHOS malfunction in COVID-19 patients primarily enables generation of a glycolytic environment, which in turn promotes the cytokine storm. This interpretation is supported by our observed increased expression of genes involved in the IFN pathway, as well as increased expression of genes involved in RIG-I-like receptors and MAVS pathway as previously shown (Gordon et al. 2020; Rui et al. 2021) . The latter might lead to activation of a signaling cascade that positively induces the genes encoding type I interferons (Rehwinkel and Gack 2020) . We speculate that our observed consistent alteration in mtDNA gene expression and related pathways enabled the viral induced cytokine storm, and hence is fundamental to the disease etiology. Testing for this possibility requires a future time-course controlled experiment J o u r n a l P r e -p r o o f 13 in cells (preferably from the immune system) in which gene expression is assessed in time intervals following SARS-CoV-2 infection.

We noticed, that cells and tissues that displayed a consistently significant difference in mitochondrial gene expression between patients and controls (in at least two independent datasets), displayed reduction in mtDNA gene expression in patients (both in bulk RNA-seq and scRNA-seq), regardless of difference between sample collection sites and sequencing platforms used. In contrast, elevated mtDNA gene expression in patients was inconsistent among databases (such as the observation of Platelets(. This finding strongly suggests that SARS-CoV-2 affects mitochondrial regulation in a similar manner in different cells. Nevertheless, although the direction of the effect is consistent (e.g. reduction in mtDNA gene expression in COVID-19 patients) the regulatory phenomenon is influenced by cell-types, which raises an interesting question as to the role of the affected cell types in the disease etiology.

As mentioned above, our findings indicate that bulk RNA-seq from blood, but not from the respiratory tract, showed significant reduction in mtDNA gene expression. Apart from the apparent cell type composition, and differences between these two types of samples, extracellular intact cell-free active mitochondria have been found in blood (Al Amir Dache et al. 2020). It will therefore be of interest to measure the number of extracellular mitochondria in both blood and respiratory tract samples from COVID-19 patients as compared to control, and assess their transcriptional signatures. This may add another layer to our understanding of the alteration in mitochondrial gene expression patterns in COVID-19 patients.

To identify first clues for the candidate mechanism underlying reduced mtDNA gene expression in COVID-19 patients we took a candidate gene approach and sought for association between altered mtDNA gene expression in patients with known regulators of mtDNA transcription, post transcription and replication. Our findings indicate that in patient peripheral and whole blood samples, reduction in mtDNA gene expression associated with reduced expression of the mitochondrial RNA-polymerase (POLRMT), along with elevated expression of c-Jun, and reduced expression of JunD. The reduction in POLRMT expression provides a simple explanation for the reduced mtDNA gene expression in the patients' peripheral blood, but this observation was not clearly identified in single cells. Interestingly, unlike POLRMT, the altered expression of c-Jun and JunD was consistently found in both the bulk J o u r n a l P r e -p r o o f 14 RNA-seq and sc-RNA-seq datasets, e.g. the expression of c-Jun increased and JunD decreased in patients in cell types that showed significantly reduced mtDNA gene expression. This particularly attracted our attention, as c-Jun and JunD binding in human mtDNA was identified in certain cell lines (Blumberg et al. 2014) , thus suggesting their possible involvement in mtDNA transcriptional regulation. Our identified association of altered c-Jun and JunD expression in patient samples suggests that they likely act as a repressor and activator, respectively, of mtDNA gene expression, and that such impact is induced in patients. As c-Jun and JunD are well-known for their regulatory targets in the nucleus it will be of great interest to experimentally assess their mechanism of action while in the mitochondria, their possible involvement in mito-nucelar co-regulation and whether their mitochondrial localization is altered in response to SARS-CoV-2 infection.

While further considering the mechanism by which mtDNA gene expression was down regulated in COVID-19 patients, one should consider additional factors, such as microRNAs. Indeed, recent microRNA-seq GO-term analysis of peripheral blood ) revealed altered expression of microRNAs in COVID-19 patients that target genes which associate with the mitochondrial matrix.

Furthermore, microRNA-seq differential expression analysis from blood plasma (Farr et al. 2021) revealed that the expression levels of hsa-miR-542-5p, a miRNA which likely co-localize with human mitochondria (Borralho et al. 2015) and was previously associated with mitochondrial dysfunction (Garros et al. 2017) , was significantly reduced in COVID-19 patients. Third, the expression of hsa-miR-483-5p which binds and represses the activity of FIS1 (Purohit and Saini 2021), a key mitochondrial fission component (James et al. 2003) , was significantly elevated in COVID-19 patients who were treated by supplemental oxygen. Hence, there is room for future investigation of the role of mitochondrial-targeted microRNAs in mtDNA regulation in general, and in COVID-19 in particular.

Along with the observed reduced expression of mitochondrial genes we identified elevated expression of several genes involved in the IFN1 response pathway. This is interesting, as it was previously shown that treatment of mammalian cells with IFN leads to reduction in mtDNA gene expression, thus suggesting a regulatory impact of IFN on mtDNA transcripts (Shan et al. 1990 ). This offers an attractive explanation for the connection between the cellular response to SARS-CoV-2 infection and changes in mitochondrial regulation. As our results indicate that such changes occur in a cell type-dependent manner, mostly apparent in immune system cells, it is not surprising that previous analyses of (lung) epithelial cells following SARS-CoV-2 infection displayed a rather low J o u r n a l P r e -p r o o f 15 response of type I and III interferon (Blanco-Melo et al. 2020) , and inconsistent changes in mtDNA genes expression (Miller et al. 2021) . Hence, it will be of great interest to carefully assess in the future the underlying mechanism by which IFN modulates mitochondrial gene expression regulation.

In summary, we observed reduced levels of mtDNA gene expression in multiple cell types, yet preferentially in cells belonging to the immune system, regardless of collection from blood or from the respiratory tract. This finding, along with apparent opposite gene expression changes in nuclearencoded OXPHOS genes in COVID-19 patients suggest departure from co-expression regulation of the mitochondrial and nuclear genomes. This interpretation may explain the elevated expression of genes involved ROS production which likely reflect cellular response to mitochondrial dysfunction. This change was also accompanied by elevated expression of glycolytic enzymes, especially LDHA, and reduction of LDHB expression. Such findings suggest that upon SARS-CoV-2 infection, cells which particularly belong to the immune system, rewire to glycolysis. Analysis of mtDNA gene expression in the unique set of patients who provided samples during COVID-19 disease progression suggests that the reduction of mtDNA gene expression is reversible upon recovery. It is thus possible that recovery of mitochondrial function predicts better health conditions of COVID-19 patients, thus underlining the mitochondria as an important drug target to ameliorate patients' health conditions.

Our work is based on analysis of gene expression at the RNA level. Nevertheless, expression variability in COVID-19 patients may be also observed at the protein level, which would be interesting to assess, yet obviously cannot be observed in analysis of RNA-seq. Secondly, RNA-seq analysis, which corresponds to the steady-state RNA levels, cannot decipher whether the observed changes in mitochondrial gene expression in COVID-19 patients, were due to transcriptional regulation, RNA decay or both. Clues to such were identified by observing correlation between the changes in mitochondrial gene expression in patients and potential regulators of mitochondrial transcription (POLRMT, c-Jun and JunD). Such associations should be considered with caution until subsequent experimental validations are performed.

Author contributions: HM analyzed the data and participated in writing the manuscript; AZ Analyzed bulk RNA-seq data; DM conceived the study, wrote the paper and supervised the analysis. Table 1 for resources). For scRNA-seq data quality control, cells were filtered to retain only those with read counts in mtDNA encoded-protein genes. Finally, differential expression analysis of mitochondrial genes was performed. *Dataset III was used both for bulk and scRNA RNA-seq analyses, and hence was considered two separate datasets.

Box plot of Bulk RNA-seq analysis in peripheral blood display lower mtDNA gene expression in COVID-19 patients as compared to healthy controls (Dataset II). (B) Box plot displaying mtDNA gene expression across COVID-19 pseudotimes as compared to control (Dataset III -bulk RNA-seq). X-axis -mtDNA genes, Y-axis -normalized read counts, which account for expression levels. Significance: * -p<0.05, ** -p<0.005, *** -p<0.0005. See also Figures S1, S2 and Table S1 . Fig S2 for Dataset V). X-axis -gene names; Y-axis -normalized read counts as in Fig. 2 . Significance: * -p<0.05, ** -p<0.005, *** -p<0.0005. See also Figure S2 and Table S2 , S3 and Data S1-10(http://dx.doi.org/10.17632/8kd3xjfrh4.1.).

). X-axis -gene names; Y-axis -normalized UMI (unique molecular identifier) counts. (B) Bar plot presenting the fraction of significantly altered mtDNA genes' expression per cell type. X-axis: cell types present in at least two datasets, Y-axis: fraction of mtDNA genes with significantly altered expression. Orange -Epithelial cells, Green -Immune system cells. Dataset numbers are indicated in parenthesis. The dashed line represents a threshold of significance in half of the analyzed mtDNA genes. See also Data S1-S10 and Table S1, S4 (http://dx.doi.org/10.17632/8kd3xjfrh4.1.). 

Further information and requests for resources should be directed to and will be fulfilled by the Lead Contact, Dan Mishmar (dmishmar@bgu.ac.il).

This study did not generate new unique reagents.

J o u r n a l P r e -p r o o f 20

Links of the source data used in this paper are available in the Key Resources The bulk and scRNA-seq codes used for data analysis are available from GitHub https://github.com/medinih/COVID_mishmarlab.

Bulk RNA-seq raw data (read) counts were downloaded from GEO and Github (Table 1) scRNA-seq data: scRNA-seq Seurat processed objects and sparse count matrices were obtained from each respective study (Stuart et al. 2019 and Cambridge but not from London in Dataset XII. Furthermore, for Dataset VI, a filtered integrated matrix was created from the filtered expression matrices using Seurat v3 pipeline as previously described to remove batch effects across all donors (Butler et al. 2018) . For all datasets Seurat objects/h5ad objects were analysed with the same cell type annotations. For the sake of consistency,

Dataset XII that was downloaded as h5ad objects was converted into Seurat objects using "sceasy" R package (Cakir et al. 2020 ) and h5 objects (Dataset VI) were converted into Seurat objects using "Read10X_h5" Seurat function and was processed as previously described (Liao et al. 2020) .

Additionally, expression matrices (Dataset IX, XI) were converted into Seurat objects using high quality annotated cells as performed previously (Ren et al. 2021; Ziegler et al. 2021 ).

Five datasets (Datasets III, VI, VIII, XI and XII) used 10X Genomics sequencing platform while datasets VII, IX, X, utilized the seq-well sequencing platform. We used the following Gene Ontology (GO) Terms to download and perform differential expression analysis: Genes under the "TCA" ( To avoid technical noise, cells with at least one mtDNA gene with zero reads were removed from further analysis for all scRNA-seq datasets. Notably, each scRNA-seq dataset were separately analysed.

Bulk RNA-seq raw expression counts matrices were normalized using DESeq normalization method available from the DESeq2 R package (Love et al. 2014 (Vallejos et al. 2017) , we included in our analysis only cells whose sequencing read counts were higher than zero in our analysis of mtDNA genes. For each respective Seurat object, the analysis was performed per cell type. Permutation of cells from each cell type was applied to control for variation in cell numbers collected from healthy individuals and patients. Equal number of cells was sampled (n=50 cells, for each group compared) to generate subset of Seurat object 1000 times per cell type. To avoid small sample sizes, we removed cell types represented by less than 50 cells per group of individuals (i.e., control healthy individuals, moderate and severe/critical patients or in each of the pseudotime points in Dataset III) and with at least 300 cells in total. Two comparisons were performed: in the first analysis we compared gene expression levels between cells from healthy individuals and COVID-19 patients without considering severity levels (sampling 1000 times n=50 cells for each group). The second analysis was performed using pairwise comparison of gene expression levels between cells from different conditions. In each iteration "FindAllMarkers" function of Seurat was used (performing Wilcoxon signed rank test), followed by BH (false discovery rate) correction. P-values for all iterations were computed by dividing the number of insignificant p-values (p>0.05) by 1000. For gene-by-gene nuclear differential expression analysis "FindAllMarkers" function was applied to all cells per cell type with min.pct=0.25 parameter. Genes from each group of genes were extracted from the resulted list of genes and were corrected with (BH) correction for multiple testing. Heatmaps were generate using the 'ComplexHeatmap' R package (Gu et al. 2016 ) using consistent genes with log fold change greater than 0.2 or less than −0.2. Box plot and bar plot graphs were generated using "ggplot2" (Wickham 2009 ) and RColorBrewer R packages.

Tables were generated using "openxlsx", "tidyr" and "reshape2" (Zhang 2016) R packages.

Descriptive titles -Supplemental Tables S1-S4 Table S3 . Differential expression of nuclear-mitochondrial related genes in Bulk RNA-seq -Healthy vs COVID-19 patients. Yellow background -genes with log fold change higher than 0.2 or less than -0.2, J o u r n a l P r e -p r o o f 23 Bold: genes with log fold change higher than 0.5 or less than -0.5. Related to Figure 3 , Figure S2 , Data S10 (http://dx.doi.org/10.17632/8kd3xjfrh4.1.). Table S4 . Single cell dataset numbers (before and after zero-mtDNA gene expression filtering).

Related to Figure 4 , and Data S1-S10 (http://dx.doi.org/10.17632/8kd3xjfrh4.1.). V)  OGDHL  IDH1  DLD  IDH2  ME2  SDHB  GCKR  IGF1  NUPR1  MLXIPL  PRKAG3  HK3  PFKFB3  PFKFB2  OGDHL  LDHA  NUP58  HDAC4  GAPDH  HK2  EP300  HIF1A  PKM  PGM1  PGAM1  NUP93  NUP88  ENO2  ALDOC Cell type Dataset *** *** *** *** *** *** *** *** *** *** *** *** *** * *** * *** *** *** *** *** *** * *** * *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** * * *** * * Ciliated (NP) *** *** *** *** *** *** *** ** ** ** ** ** *** ***

Blood contains circulating cell-free respiratory competent mitochondria

Evidence for the presence of alpha and beta-related T3 receptors in rat liver mitochondria

Human primitive brain displays negative mitochondrialnuclear expression correlation of respiratory genes

Longitudinal Multi-omics Analyses Identify Responses of Megakaryocytes, Erythroid Cells, and Plasmablasts as Hallmarks of Severe COVID-19

Imbalanced Host Response to SARS-CoV-2 Drives Development of COVID-19

Transcription factors bind negatively-selected sites within human mtDNA genes

mtDNA synthesis ignites the inflammasome

microRNAs in Mitochondria: An Unexplored Niche

Mitochondria: In the Cross Fire of SARS-CoV-2 and Immunity

Integrating single-cell transcriptomic data across different conditions, technologies, and species

Comparison of visualization tools for singlecell RNAseq data

COVID-19 severity correlates with airway epithelium-immune cell interactions identified by single-cell analysis

Elevated Glucose Levels Favor SARS-CoV-2 Infection and Monocyte Response through a HIF-1alpha/Glycolysis-Dependent Axis

Ancient Out-of-Africa Mitochondrial DNA Variants Associate with Distinct Mitochondrial Gene Expression Patterns

Altered microRNA expression in COVID-19 patients enables identification of SARS-CoV-2 infection

Mitochondrial gene expression is regulated at multiple levels and differentially in the heart and liver by thyroid hormones

Impact of COVID-19 on Mitochondrial-Based Immunity in Aging and Age-Related Diseases

MicroRNA-542 Promotes Mitochondrial Dysfunction and SMAD Activity and Is Elevated in Intensive Care Unit-acquired Weakness

Mitochondria Targeted Viral Replication and Survival Strategies-Prospective on SARS-CoV-2

A SARS-CoV-2 protein interaction map reveals targets for drug repurposing

Complex heatmaps reveal patterns and correlations in multidimensional genomic data

The COVID-19 Effect on the Immune System and Mitochondrial Dynamics in Diabetes, Obesity, and Dementia. The Neuroscientist : a review journal bringing neurobiology

hFis1, a novel component of the mammalian mitochondrial fission machinery

Osteogenic differentiation of human MSCs: Specific occupancy of the mitochondrial DNA by NFATc1 transcription factor

CD177, a specific marker of neutrophil activation, is associated with coronavirus disease 2019 severity and death

Differential microRNA expression in the peripheral blood from human patients with COVID-19

Single-cell landscape of bronchoalveolar immune cells in patients with COVID-19

In vivo antiviral host transcriptional response to SARS-CoV-2 by viral load, sex, and age

Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2

Upper airway gene expression differentiates COVID-19 from other acute respiratory illnesses and reveals suppression of innate immune responses by SARS-CoV-2. medRxiv : the preprint server for health sciences

Host mitochondrial transcriptome response to SARS-CoV-2 in multiple cell models and clinical samples

Anticancer targets in the glycolytic metabolism of tumors: a comprehensive review

Mitochondrial microRNA (MitomiRs) in cancer and complex mitochondrial diseases: current status and future perspectives

Structural basis for altered activity of M-and Hisozyme forms of human lactate dehydrogenase

RIG-I-like receptors: their regulation and roles in RNA sensing

COVID-19 immune features revealed by a large-scale single-cell transcriptome atlas

Unique and complementary suppression of cGAS-STING and RNA sensing-triggered innate immune responses by SARS-CoV-2 proteins

Interferon selectively inhibits the expression of mitochondrial genes: a novel pathway for interferon-mediated responses

Upregulation of mitochondrial gene expression in PBMC from convalescent SARS patients

Direct regulation of complex I by mitochondrial MEF2D is disrupted in a mouse model of Parkinson disease and in human patients

SARS-coronavirus open reading frame-9b suppresses innate immunity by targeting mitochondria and the MAVS/TRAF3/TRAF6 signalosome

Decoding SARS-CoV-2 hijacking of host mitochondria in COVID-19 pathogenesis

Polyadenylation and degradation of human mitochondrial RNA: the prokaryotic past leaves its mark

Single-cell multi-omics analysis of the immune response in COVID-19

Transcriptomic similarities and differences in host response between SARS-CoV-2 and other viral infections

Normalizing single-cell RNA sequencing data: challenges and opportunities

Ggplot2 : elegant graphics for data analysis

Multi-omic profiling reveals widespread dysregulation of innate immunity and hematopoiesis in COVID-19

A single-cell atlas of the peripheral immune response in patients with severe COVID-19

Functional genomic analysis of human mitochondrial RNA processing

RNA-GPS Predicts SARS-CoV-2 RNA Residency to Host Mitochondria and Nucleolus

Redox crosstalk at endoplasmic reticulum (ER) membrane contact sites (MCS) uses toxic waste to deliver messages

Reshaping and aggregating data: an introduction to reshape package