SARS-CoV-2 induces a unique mitochondrial transcriptome signature

SARS-CoV-2 has challenged global healthcare systems in part because its clinical manifestations are heterogeneous. Variable symptoms of SARS-CoV-2 could be attributed to the virus’ ability to mildly induce an innate immune response, as prior transcriptomic data has suggested. Mitochondrial dynamics might partially mediate the effect of SARS-CoV-2 on innate immunity. Many proteins encoded by SARS-CoV have been shown to localize to mitochondria and inhibit Mitochondrial Antiviral Signaling protein (MAVS). We analyzed multiple publicly available RNASeq data in order to unravel the metabolic and mitochondrial transcriptome signature of SARS-CoV-2 in primary cells, cell lines, and clinical samples (i.e., BALF and lung). We report here that SARS-CoV-2 does not dramatically regulate (1) mitochondrial-gene expression or (2) MAVS expression, but (3) downregulates nuclear-encoded mitochondrial (NEM) genes related to cellular respiration and Complex 1 assembly. We also report cell-specic and tissue-specic effects of SARS-CoV-2 on the mitochondrial-encoded and NEM transcriptome that could inform future experimental paradigm selection.


Introduction
The emergence of SARS-CoV-2 in late 2019 (COVID19) has put enormous pressure on global economic and health systems. 1 SAS-CoV-2 patients present differing degrees of disease severity due to underlying comorobidites. 2 Diabetic patients hospitalized in New York City were more likely to receive severe intensive care interventions than patients without diabetes. 3 It is unclear how certain metabolic conditions increase mortality in SARS-COV-2 patients. However, SARS-CoV-2 induces a modest pro in ammatory response that is distinct from other respiratory viruses. 4 This SARS-CoV-2 muted proin ammatory response is a phenomenon that might explain why individuals without metabolic comorbidities typically present less severe symptoms compared to patients with underlying comorbidities.
The effects of SARS-CoV-2 on the host cell innate immune response could be mediated by mitochondrial dynamics. Independent analyses of the SARS-CoV-2 transcriptomic data set revealed enrichment for mitochondrial organization processes, and host cell innate immunity is regulated by the Mitochondrial Antiviral Signaling protein (MAVS). 5,6 MAVS normally interacts with MFN2 under resting conditions. 7 But after viral infection, mitochondria-associated ER membranes and nearby mitochondria become tethered by MFN2 and RIG-1, forming a complex that recruits TRIM25 and the molecular chaperone 14-3-3e into a translocon structure. 8 This translocon localizes to the mitochondrion where it binds MAVS, after which MAVS interacts with TANK binding kinase 1, IKKA, and IKKB. 9 The host cell's immune and apoptotic response is ampli ed when MAVS induces phosphorylation and nuclear translocation of IRF3 ( Figure 2B). SARS-CoV, which was identi ed in 2002 amid the international SARS outbreak, displayed MAVS-inhibiting functions. For example, RIG-1 and MDA5 were inhibited by SARS-CoV, and the SARS-CoVORF9b targeted the mitochondrial-associated adaptor molecule MAVS signalosome. 10 Notably, ORF9b of the novel SARS-CoV-2 has been characterized as one of the largest hubs in interactome analysis. 11 SARS-CoVOR3b also localized to the mitochondria in vitro, and mitochondrial-gene expression was upregulated in peripheral blood mononuclear cells (PBMC) of infected patients. 12 The interaction between SARS-CoV and mitochondria suggest that the mitochondrion also participates in the SARS-COV-2-speci c host cell response. It is plausible that SARS-CoV-2 has speci c effects on both the nuclear-encoded mitochondrial (NEM) transcriptome mitochondrial transcriptome and mitochondrial-encoded transcriptome. In the following analyses, we assessed the effect of SARS-CoV-2 on the mitochondrial genetic network transcriptome by reanalyzing publicly available RNASeq data.

Selection of COVID datasets
In order to examine the NEM and mtDNA expression signature in SARS-CoV-2 infection, we utilized data sets that were uploaded to GEO (GSE147507 and GSE110551) and the BIG Data Center (CRA002390). These RNASeq data sets were derived from A549, A549 (ACE2), Calu-3, and NHBE cells as well as from SARS-CoV-2 patients' lung autopsies and bronchoalveolar lavage fluid (BALF). A549 cells were infected with seasonal influenza A virus (IAV), human orthopneumovirus (respiratory syncytial virus; RSV), human parainfluenza virus 3 (HPIV3), and SARS-CoV-2. ACE2-expressing A549 cells, Calu-3 cells, and NHBE cells were infected with SARS-CoV-2. NHBE cells were also infected with IAV. The original authors who curated the in vitro data infected SARS-CoV-2 in A549 cells at low and high multiplicities of infection (MOI).
They found that the rate of SARS-CoV-2 replication after low MOI was comparable to the replication rate after high MOI in ACE2-expressing A549 cells. 4 The original authors also observed that low MOI SARS-CoV-2 infection stimulated a relative muted proinflammatory response, which was ablated in high MOI SARS-CoV-2 infection in ACE2-expression A549 cells.
We specifically contrasted infection conditions by using low MOI SARS-CoV-2 in order to (1) stay consistent with previously published results and limit confounding effects from stoichiometry disruption of high SARS-CoV-2 components. The number of differentially expressed genes (DEGs) and NEM DEGs per biological source is listed in Table 1  downregulated every mtDNA-encoded protein gene and several tRNAs ( Figure 1A). Such strong regulatory effects of IAV on mtDNA-encoded genes was not consistent in cancerous cell lines. In A549 cells, while IAV downregulated mt-ND6 and mt-ATP6, we did not observe the same global downregulation as we did in primary NHBE cells, and we found 16S rRNA was upregulated only in these cancerous cells ( Figure 1B). HPV did not regulate expression of any mtDNA-encoded protein gene ( Figure 1B). RSV induced dramatic upregulation of every mtDNA-encoded protein along mt-rRNA and a few tRNAs ( Figure 1B). SARS-CoV-2 did not upregulate any mtDNA-encoded proteins but did upregulate 16S rRNA in Calu-3 and ACE2expressing A549 cells ( Figure 1B). In BALF, SARS-CoV-2 surprisingly downregulated nearly every mtDNA-encoded gene along with several mt-tRNAs ( Figure 1C). Overall, following SARS-CoV-2 infection, we observed downregulation of mtDNA-encoded genes in BALF, which is opposite to the minimal upregulation we observed in primary and cell lines. The complete list of significant mtDNA-encoded genes and fold changes are included in Supplementary Table: mtDNA Differentially Expressed Genes.

NEMS Sufficiently Classifies SARS-CoV-2
Given that MAVS and mtDNA-encoded gene expression differs among SARS-CoV-2, HPIV, RSV, and IAV, we hypothesized that the global NEM signature would sufficiently classify SARS-CoV-2. Therefore, we conducted a principal component analysis (PCA) exclusively on an NEMextracted gene set. As expected, the first two principal components sufficiently reduced NEM expression variance in a manner that classified SARS-CoV-2 in primary cells, cell lines, and clinical samples ( Figure 3). The amount of variance that the first two NEM-specific two principal components explain total 81%, 60%, and 56% for primary cells, cell lines, and clinical samples, respectively.

SARS-CoV-2 Specific NEM-Enriched Pathways
Since we showed that NEMs are sufficient to classify SARS-CoV-2, HPIV, RSV, and IAV, we attempted to unravel the biological processes that these NEMs modify. All NEMs were extracted from the complete list of statistically significant DEGs. Hierarchical clustering of NEMs showed distinct signatures by viral infection ( Figure 4A). We then then conducted gene enrichment analyses by inputting this set of significant NEMs against a universe background of all total significant DEGs. Four separate gene enrichment analyses were conducted (i.e., primary cells, cell lines, BALF, and lung).
The GO enriched terms of small molecule metabolism, phosphorus metabolism, oxidationreduction, and cellular amide metabolism were all shared between SARS-CoV-2 and IAV in NHBE primary cells ( Figure 4B; filtered by >20% NEM within gene set). SARS-CoV-2 particularly induced greater enrichment for mitochondrion organization and catabolism compared to IAV and IAVdNS1 ( Figure 4B). Furthermore, the top 10 most significant enriched for SARS-CoV-2 mapped back to mitochondrial translation, mitochondrial organization, and cellular respiration ( Figure 4C and 4D). SARS-CoV-2 not only induced global downregulation of the metabolic pathways shown in Figure 4C (Log2FC), but the degree of downregulation was greatest in SARS-CoV-2, as illustrated in the right panel of Figure 4C with hierarchical clustering scores colored (SARS-CoV-2 Score). Several mitochondrial ribosome protein genes were clear tissue and cell-specific differences across all analyses related to oxidation-reduction, small molecule metabolism, and carboxylic metabolism, suggesting SARS-CoV-2 may affect metabolism differently per cell type.

Discussion
We concentrated on the mitochondrial and nuclear-encoded mitochondrial gene host response to SARS-CoV-2 and other human respiratory viruses in multiple cell models and clinical samples. Our analyses showed that the metabolic transcriptome signature of SARS-CoV-2 infection included both shared and independent biological processes to IAV, HPIV, and RSV. We observed three notable metabolic transcriptomic signature differences including mitochondrial-gene expression, MAVS expression, and NEM biological processes.
First, we found that SARS-CoV-2 did not induce a dramatic mitochondrial-gene expression signature in NHBE, A549, and Calu-3 cells. In NHBE cells, we observed that IAV and IAVdNS1 induced near global downregulation of mitochondrial-encoded genes, including all 13 protein-encoding genes. Since IAVdNS1 (i.e., IAV without its interferon antagonizing protein NS1) also induced similar downregulation of mitochondrial-encoded genes as IAV, it is plausible that these mitochondrial-gene regulatory effects are independent from an interferon response. Furthermore, since SARS-CoV-2 did not drastically in uence mitochondrial-gene expression, this suggests (1) cells are not experiencing extreme energetic stress, (2) capacity for SARS-CoV-2 to evade mechanisms regulating mitochondrial-gene expression, or (3) compensation via the NEM network.
We did, however, observe dramatic downregulation of mitochondrial-gene expression in human BALF after SARS-CoV-2 infection, although we did not observe the same degree of downregulation in human lung. The cellular pro le of BALF after infection differs signi cantly compared to resting state conditions.
Prior reports have showed nearly a 900-fold increase in neutrophil in ltration after LPS stimulation. 13 Therefore, it is possible that decrease in mitochondrial-gene expression that we observed in SARS-CoV-2 patient BALF is a byproduct of a different cellular pro le, and it is also possible that SARS-CoV-2 affects mitochondria in immune cells differently than in lung cells. Furthermore, we are able to quantify expression mitochondrial-encoded genes in this traditional RNASeq work ow, but we cannot rule out that the effects we observed are partially in uenced by RNA with non-protein-encoding functions. Perhaps the differences (or lack thereof) in mitochondrial-gene expression are a byproduct of mitochondrial small open reading frame post-transcriptional regulation or the immune-signaling capacity of mitochondrial RNA and DNA released into the cytoplasm. [14][15][16] Second, we found that SARS-CoV-2 did not alter MAVS expression, whereas IAV, HPIV, and RSV all decreased MAVS expression. Although SARS-CoV-2 differentially regulates the host cell's interferon response, MAVS is not directly regulated by interferons. 17 MAVS is instead regulated by a negative feedback loop by reactive oxygen species (ROS). 18 MAVS expression typically decreases after RNA viral infection due to the rise of ROS. 19 MAVS is also regulated by calcium, ATP, cholesterol, ceramides, and phospholipids. 20 That SARS-CoV-2 does not affect MAVS expression suggests (1) SARS-CoV-2 may be capable of "hijacking" MAVS; (2) in the presence of SARS-CoV-2, MAVS is regulated differentially at the post-transcriptional or post-translational level; or (3) SARS-CoV-2 can propagate e ciently in the presence of unaltered MAVS expression, which could partially explain the muted antiviral response published previously.
Third, we observed that, across multiple cell and tissue types, SARS-CoV-2 reduced NEM expression related to cellular respiration and Complex 1 assembly. Reports on respiratory viruses such as RSV have indicated that Complex 1 inhibition could promote e cient viral replication. 21 We found that many NDUF Complex 1-associated proteins were downregulated after SARS-CoV-2 infection in NHBE cells, cell lines, BALF, and lung. Primary cells and BALF showed the greatest degree of downregulation of genes related to cellular respiration. SARS-CoV-2 not only decreased expression to many NDUF family of genes, but also several mitochondrial ribosome protein related genes in primary cells and BALF, observations of which did not carry over to cell lines. It is possible that the cancerous nature of A549 and Calu-3 cells limits our interpretation on metabolism due to naturally different cellular metabolic features. 22 In addition, in primary cells and cell lines, we found that SARS-CoV-2 uniquely affected a greater number of genes related to cellular catabolism, but it affected fewer genes related to cellular amide processing compared to other viruses, which could be due to the altered MAVS expression across viral infections. We also saw discordant shares of NEMs related to lipid processing and carboxylic acid metabolism across cell and tissues after SARS-CoV-2 infection perhaps due to the underlying biology of those cells. For example, we observed a greater share of NEMs related to carboxylic acid in primary cells after SARS-CoV-2 infection relative to other viruses, while, in contrast, we observed greater shares after other respiratory viral infections in cell lines.
In addition to the biological implications of SARS-CoV-2 infection presented here, our analyses suggest future research should carefully consider the in vitro model, especially if the desired outcome of interest is closely related to metabolic features. We highlighted more potent metabolic transcriptomic effects of SARS-CoV-2 and IAV in NHBE cells compared to A549 and Calu-3 cells. We observed several mitochondrial ribosome protein genes were downregulated after SARS-CoV-2 in NHBE cells, but we did not observe a similar phenomenon in cancerous cell lines and human BALF and lung. We found that IAV and IAVdNS1 potently downregulation mitochondrial-gene expression in primary NHBE cells but not to the same degree in A549 cells. We also identi ed signi cantly fewer overall genes and NEM genes that were differentially expressed in human lungs, which suggests tissues with heterogenous cellular pro les could be confounding and therefore hiding true metabolic effects.
Despite cell and tissue-speci c effects, our analyses showed consistent downregulation of MAVS and cellular respiration was a key, shared feature of SARS-CoV-2 infection in multiple cell and tissue sources.
Recent reports have predicted that SARS-CoV-2 RNA localizes to the mitochondrion. 23 An additional prereview report hinted that it is plausible that SARS-CoV-2 RNA can anneal with mitochondrial deubuiquitanse USP30, a subunit of ubuquiting protein ligase complex FBX021 (https://www.biorxiv.org/content/10.1101/2020.04.08.031856v3.full). SARS-CoV-2 RNA acting as an RNAi might explain many of the downregulatory expression effects that we observed. Downregulation of genes involved in cellular respiration could limit ROS production and, as a result, promote viral replication. Follow-up studies should consider the mechanisms by which SARS-CoV-2 affects mitochondrial biology. Analyses presented here have further characterized the pathogenesis of SARS-CoV-2 infection in context of mitochondrial biology and highlighted additional therapeutic directions.

Methods
The molecular preparation work ow for the RNASeq data used here has been reported previously. We downloaded sample FASTQ les from GEO or BIG Data Center. These FASTQ reads were mapped to rRNA sequences in order to remove cytoplasmic rRNA reads using STAR (v2.7.2b) with default parameters. Bioinformatically-ltered cytoplasmic rRNA FASTQ les were then mapped to the human genome (hg38) with GENCODE gene annotation (v33) using the following STAR parameters: sjdbScore 1, outFilterMultimapNmax 20, outFilterMismatchNmax 999, outFilterMismatchNoverReadLmax 0.04, alignIntronMin 20, alignIntronMax 1000000, alignMatesGapMax 1000000, and alignSJoverhangMin, alignSJDBoverhangMin 1. Aligned sorted BAM les were inputted into RR (v3.5.1) for counting. A count matrix and corresponding metadata were sorted in an S4 class derived from the SummarizedExperiment class of the GenomicRanges package in R. Count matrices were generated using the summarizeeOverlaps function of the GenomicAlignemnts package in R. Overlapping genomic features were resolved using the "Union" mode in the summarizeOverlaps function. Differential expression analyses were conducted using the DESeq R package. Signi cant differentially expressed genes were ltered by a padjusted value of 0.2. Genes under the "mitochondrion" Gene Ontology (GO) Term (GO:0005739) were downloaded and used for NEM enrichment analyses. There are 1842 mitochondrion GO terms, which we referred to as NEMs. This NEM gene set was used to extract such genes from the dds object generated from the DESEq2 R package. The dds object was transformed by variance stabilization using the varianceStabilizingTrasnformation function, which was used during the principal component analysis (plotPCA function in R) and euclidean heirarchial clustering (pheatmap function in R) scaled to each condition. Mitochondrial DNA expression heat maps were built using custom scripts in R. NEM enrichment was conducted using the clusterPro ler package in R. For GO analyses, the NEM-extracted gene set from statistically signi cant DEGs were tested against a universe background gene set of all statistically signi cant DEGs. This approach was implemented in order to identify biological processes that NEMS enriched compared to the background of all DEGs. The cut-off criteria for GO analysis were p < 0.05 and q < 0.20. Signi cant enriched terms were visualized using the two clusterPro ler package functions heatplot and cnetplot.

Con ict Of Interest
Pinchas Cohen is a consultant and stockholder of CohBar Inc. All other authors have no competing interests to declare.