Skip to main content

ORIGINAL RESEARCH article

Front. Genet., 16 March 2021
Sec. Computational Genomics
This article is part of the Research Topic Coronavirus Disease (COVID-19): Molecular Mechanisms, Translational Approaches and Therapeutics View all 118 articles

Network Analysis and Transcriptome Profiling Identify Autophagic and Mitochondrial Dysfunctions in SARS-CoV-2 Infection

  • 1Bioinformatics and Computational Biology Laboratory, National Heart, Lung, and Blood Institute, National Institutes of Health, Bethesda, MD, United States
  • 2Laboratory of Mitochondrial Biology and Metabolism, National Heart, Lung, and Blood Institute, National Institutes of Health, Bethesda, MD, United States

Analyzing host cells' transcriptional response to severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection will help delineate biological processes underlying viral pathogenesis. First, analysis of expression profiles of lung cell lines A549 and Calu3 revealed upregulation of antiviral interferon signaling genes in response to all three SARS-CoV-2, MERS-CoV, or influenza A virus (IAV) infections. However, perturbations in expression of genes involved in inflammatory, mitochondrial, and autophagy processes were specifically observed in SARS-CoV-2-infected cells. Next, a validation study in infected human nasopharyngeal samples also revealed perturbations in autophagy and mitochondrial processes. Specifically, mTOR expression, mitochondrial ribosomal, mitochondrial complex I, lysosome acidification, and mitochondrial fission promoting genes were concurrently downregulated in both infected cell lines and human samples. SARS-CoV-2 infection impeded autophagic flux either by upregulating GSK3B in lung cell lines or by downregulating autophagy genes, SNAP29, and lysosome acidification genes in human samples, contributing to increased viral replication. Therefore, drugs targeting lysosome acidification or autophagic flux could be tested as intervention strategies. Finally, age-stratified SARS-CoV-2-positive human data revealed impaired upregulation of chemokines, interferon-stimulated genes, and tripartite motif genes that are critical for antiviral signaling. Together, this analysis has revealed specific aspects of autophagic and mitochondrial function that are uniquely perturbed in SARS-CoV-2-infected host cells.

Introduction

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a beta-coronavirus and the cause of the Coronavirus Disease 2019 (COVID-19) pandemic. Lack of effective treatment strategy and vaccine makes SARS-CoV-2 infection a big threat to human health and well-being. Understanding the cellular processes impacted in the host cells responding to a virus infection will be necessary to understand the virus pathogenesis and to develop drug intervention strategies.

COVID-19 presents as a wide range of clinical manifestations, ranging from asymptomatic to respiratory failure or multiorgan and systemic manifestations (Cascella et al., 2020; Ludwig and Zarbock, 2020; Wang C. et al., 2020; Zhu et al., 2020). This viral pneumonia outbreak caused by SARS-CoV-2 was first identified in Wuhan, China, in December 2019 (Chan et al., 2020). Since then, the virus has continued to spread globally, with a current transmissibility estimate (R0) between 3 and 4 (Fung et al., 2020; Yuen et al., 2020). Several drugs are currently under various phases of clinical trials, and management strategies include supportive medical care for existing cases and social distancing for prevention. Understanding this novel pathogen and the host response it elicits is crucial to combatting the emerging threat to public health.

SARS-CoV-2 is the 7th and most recent addition to human coronaviruses (hCoVs), which include four globally endemic hCoVs that cause a substantial portion of upper respiratory infections (229E, OC43, HKU1, and NL63), as well as two other highly pathogenic strains that have also caused recent pandemics [SARS-CoV and MERS-CoV (Fung et al., 2020; Raoult et al., 2020) in 2002–2003 and 2012, respectively (De Wit et al., 2016)]. All seven hCoVs are single-stranded, positive-sense RNA viruses (Ludwig and Zarbock, 2020). They all have zoonotic origins, with bats as the evolutionary reservoir host of five viruses (229E, NL63, SARS-CoV, MERS-CoV, and SARS-CoV-2). Although SARS-CoV-2 is phylogenetically similar to both MERS-CoV, and SARS-CoV (Wu et al., 2020), there are biological differences. Notably, although SARS-CoV-2 has a lower, but yet undetermined mortality rate, it is distinctly more contagious than these other highly pathogenic hCoVs, causing vastly different epidemiological dynamics.

The hCoVs differentially infect the human respiratory tract. The low pathogenic hCoVs infect the upper respiratory tract, and the highly pathogenic hCoVs infect the lower respiratory tract (Channappanavar and Perlman, 2017). Consistent with this, SARS-CoV, SARS-CoV-2, and MERS-CoV were shown to differentially infect the lung alveolar cell subtypes in cynomolgus macaques (Rockx et al., 2020) and SARS-CoV elicited distinct immune response in different tissues (To et al., 2004). Furthermore, cell tropism study of the SARS-CoV and SARS-CoV-2 in different cell type cultures could partially explain the symptomatic differences of these two virus infections (Chu et al., 2020). Single cell (sc) transcriptomic data of the COVID-19 lung tissue have been analyzed to identify the subset of cells most prone to the SARS-CoV-2 infection and the marker genes associated with the infected cells. One such study intriguingly identified upregulation of the receptor-angiotensin-converting enzyme 2 (ACE2) in the SARS-CoV-2-infected type II pneumocyte population of the lung cells as a potential mechanism facilitating virus infection (Ziegler et al., 2020). Several studies have linked the expression of ACE2 and TMPRSS2 with increased susceptibility to viral entry (Hoffmann et al., 2020; Walls et al., 2020; Yan et al., 2020). In addition to ACE2 and TMPRSS2, another study has shown that the neuropilin-1 (NRP1) could act as a host cofactor and facilitate viral entry (Daly et al., 2020).

One of the intriguing aspects of COVID-19 is that the infection leads to a wide range of the symptoms. The infected patients may be asymptomatic or may present mild to severe symptoms, which likely arise from altered immune response (Garcia, 2020). A dysregulated immune response is caused by rapid viral replication, cytokine storms, delayed interferon response, and macrophage infiltration and excessive proinflammatory cytokines (Channappanavar and Perlman, 2017; Garcia, 2020). This immunopathogenesis mechanism is supported by the observation of decreased viral loads occurring with increased disease severity (De Wit et al., 2016). Severity of illness for SARS-CoV-2 infections is likely impacted by both the direct cytotoxic effects of the virus, and the effectiveness of the complex host response (Mar et al., 2018; Astuti and Ysrafil, 2020; Chen G. et al., 2020). However, efforts to understand the molecular mechanisms driving different clinical outcomes require further study to help develop appropriate drug intervention strategies.

To delineate the host cell transcriptional response to the viral infection and potentially identify genes and biological processes (BPs) specifically impacted by SARS-CoV-2 infection, we have utilized gene expression information from several datasets from cell lines (Blanco-Melo et al., 2020), and from human nasopharyngeal samples (Lieberman et al., 2020) classified into young and old groups that were positive for high or low viral loads (see Figure 1, study schema). To facilitate SARS-CoV-2 viral entry, A549 cells transduced with human ACE2 (hACE2) infected with SARS-CoV-2 were utilized (see Methods). The transformed A549 cells and Calu3 cells both revealed viral reads when infected with SARS-CoV-2 (Blanco-Melo et al., 2020). Viral infection in the transduced A549 cells and Calu3 cells presented in this dataset was confirmed by evaluating the percent reads that aligned with the viral genome for each of the infected samples and has been published (Blanco-Melo et al., 2020). Analysis of gene expression profiles of cells infected with either SARS-CoV-2, MERS-CoV, or influenza A virus (IAV) comparisons revealed upregulation of interferon signaling genes. However, the SARS-CoV-2 infection uniquely elicited differential expression of genes involved in inflammation, autophagy, and mitochondrial processes. To validate the findings from cell lines, expression profile of the SARS-CoV-2-positive human nasopharyngeal samples was analyzed. Consistent with the cell line data, the differentially expressed (DE) genes from control vs. infected human samples also annotated to inflammation, autophagy, and mitochondrial processes. Notably, mTOR, mitochondrial ribosomal, mitochondrial complex I, lysosome acidification genes, and mitochondrial fission promoting genes were downregulated in both cell line and human datasets. Perturbation of the mitochondrial function and autophagy could negatively impact the host cells' immune response against the viral infection leading to systemic inflammation (Won et al., 2015; Jang et al., 2019). Decreased expression of mitochondrial fission promoting genes may contribute to hyperfused mitochondria and impaired interferon response (Barbier et al., 2017; Das and Chakrabarti, 2020). In the cell lines, the autophagy flux impeding GSK3B was upregulated. In the infected human samples, several autophagy genes, p62 and SNAP29, were downregulated. Together, these gene expression changes support the idea that the autophagic flux is likely decreased in SARS-CoV-2-infected cells, which may contribute to viral propagation. Therefore, drugs increasing autophagic flux, or lysosome acidification, could be tested as treatment strategies. Furthermore, the gene expression profile of A549 cell line strongly correlated with the lung epithelial lineage basal and ionocyte cell types from the lung single-cell (sc) RNA-seq data. This correlation suggests that the gene expression profile of A549 cells likely reflects lung cells' response to the SARS-CoV-2 infection. Using the age-stratified human nasopharyngeal expression data, we have also delineated some age-specific changes in antiviral signaling, which may provide more insight into the age-dependent differences in viral pathogenicity. Therefore, from this analysis, we have identified some key aspects of autophagy and mitochondrial processes that are uniquely impacted in SARS-CoV-2 infection and are likely representative of host cells' response to the infection.

FIGURE 1
www.frontiersin.org

Figure 1. Study schema. A schema highlighting the various datasets used in the study and the downstream analysis performed for each dataset is shown as a flow chart. The datasets utilized in this study are divided into two groups: Blue boxes highlight the cell lines used in the study. The viruses used for infection are indicated between gray horizontal lines. The green box represents the human nasopharyngeal datasets and the age and viral load criteria used to classify these samples. For each pair of infected and control samples, DE analysis was performed to identify DE genes and pathway enrichment analysis was performed to reveal the biological processes to which the DE genes annotated to (shown in the flow chart). From the cell line data, the DE genes and pathway enrichment results were compared to identify BPs that were either uniquely perturbed in SARS-CoV-2 infection or were commonly perturbed in all viral infection (blue lines with arrowhead leading to the blue box). For validating the findings from the cell line data, the DE results from infected human nasopharyngeal samples were analyzed. Specifically, the DE results from the A549 dataset (highlighted in orange box) was compared with the DE results from old age control vs. high viral load positive human nasopharyngeal samples (highlighted in orange box). This comparison will help identify concordant gene expression changes and BPs impacted in both datasets (gray line with arrowhead leading to orange box). Finally, using the age-stratified human nasopharyngeal datasets, age-specific gene expression changes were delineated (green boxes). DE, differential expression; BPs, biological processes; MOI, multiplicity of infection.

Results

To perform an in-depth analysis of the host cells' transcriptional response to SARS-CoV-2 infection, several gene expression datasets from different cell lines and human samples were used. An overview of the datasets analyzed and compared are depicted in Figure 1 (study schema).

Interferon Autophagy and Mitochondrial Processes Are Impacted in A549 Cells Infected With SARS-CoV-2

SARS-CoV-2 (High Viral Titer) vs. Mock

Differential gene expression analysis of hACE2 receptor-transduced A549 lung epithelial cell line that was either mock-infected or infected with SARS-CoV-2 at higher viral titer (2MOI) [see Methods (Blanco-Melo et al., 2020)] identified ~8,000 DE genes. The volcano plot profiles both upregulated and downregulated genes in the SARS-CoV-2-infected cells (Supplementary Figure 1A, Supplementary Table 3: Sheet 1). Pathway enrichment analysis of the DE genes showed enrichment in a wide range of BPs (Supplementary Figure 1B, Supplementary Table 4: Sheet 1). These DE genes were classified into upregulated or downregulated following SARS-CoV-2 infection and analyzed by pathway enrichment analysis. Upregulated DE genes annotated to a wide range of pathways, notably including the interferon signaling, NFkB/cytokine signaling processes, and proteasomal degradation (Figure 2A, Supplementary Table 4: Sheet 2). Heatmaps highlight the upregulation of genes in interferon and cytokine processes, and perturbation of genes in the autophagy pathways (Figures 2B–D, respectively). DE genes downregulated in the SARS-CoV-2-infected cells annotated to pathways primarily involving cell cycle and mitochondrial processes (Figure 2E, Supplementary Table 4: Sheet 3). A heatmap shows that the expression of the genes in mitochondria-related processes, electron transport chain, and respiration was mostly downregulated in SARS-CoV-2-infected cells (Figure 2F).

FIGURE 2
www.frontiersin.org

Figure 2. SARS-CoV-2 infection of hACE2 transduced lung epithelial A549 cells impacts expression of genes in interferon, cytokine, and autophagic processes. (A) Top 25 pathways from the pathway enrichment analysis of genes upregulated in the SARS-CoV-2-infected A549 cell is presented as a horizontal bar plot, where the x axis represents the –log10-transformed q-value and the color of the horizontal bar is scaled blue to red representing low to high q-values, respectively. (B) Heatmap highlighting the expression of genes in the interferon processes in mock and infected cells. The red and green color bands represent up- and downregulated genes, respectively. This heatmap shows that cytokine-related genes were predominantly upregulated in infected cells. (C) Heatmap highlighting the expression of genes in the cytokine processes in mock and infected cells. The red and green color bands represent up- and downregulated genes, respectively. This heatmap shows that interferon-related genes were predominantly upregulated in infected cells. (D) Heatmap highlighting the expression of genes in the autophagy-related processes in mock and infected cells. The red and green color bands represent up- and downregulated genes, respectively. This heatmap shows that autophagy-related genes were perturbed in infected cells. (E) Top 25 pathways from the pathway enrichment analysis of gene upregulated in SARS-CoV-2 infection is presented as a horizontal bar plot, where the x axis represents the –log10-transformed q-value and the color of the horizontal bar is scaled blue to red representing low to high q-values, respectively. (F) Heatmap highlighting the expression of genes in the mitochondrial organization and translation in mock and infected cells. The red and green color bands represent up- and downregulated genes, respectively. This heatmap shows that mitochondrial process-related genes were predominantly downregulated in infected cells. (G) Pathway enrichment summary map for mock vs. SARS-CoV-2 at high MOI (blue nodes) and low MOI (red nodes) comparisons. Each node represents a pathway/biological process (BP). The node size is proportional to the number of DE genes overlapping with the BP. The nodes that share genes are connected with edges. Single color nodes are pathways that are distinctly enriched by DE genes from one comparison. Two colored nodes are pathways enriched by DE genes from both comparisons. The label above each black circle summarizes the gene ontology (GO) terms of similar BPs present inside the circle. Notable groups of BPs associated with antigen processing, autophagy and mitochondria that were predominantly enriched by DE genes from Mock vs. SARS-CoV-2 (high MOI) comparison and are highlighted in red. BP associated interferon and UPR processes were also predominantly enriched by DE genes from Mock vs. SARS-CoV-2 (high MOI) comparison. MOI, multiplicity of infection; DE, differentially expressed; UPR, unfolded protein response.

SARS-CoV-2: Low Viral Titer vs. High Viral Titer

Differential gene expression analysis of hACE2-transduced A549 cells infected with mock and a 10-fold lower viral titer of SARS-CoV-2 [see Methods (Blanco-Melo et al., 2020)] was also performed. The resulting DE genes could be compared to the DE genes from mock vs. SARS-CoV-2 infection at higher viral titer (Figure 2). Given the exposure of cells to a lower viral titer, the number of DE genes from this comparison was smaller (4,494 genes) vs. the comparison of high titer SARS-CoV-2 against mock (~8,000 genes) (Supplementary Figure 1C, Supplementary Table 3: Sheet 2). Analysis of the 4,494 DE genes showed significant enrichment in inflammation, autophagy and mitochondrial processes (Supplementary Figure 1D, Supplementary Table 4: Sheet 4). To further assess the extent of overlap BPs between low MOI and high MOI A549 cell datasets, the pathway enrichment results were graphically summarized and presented in a single map. This pathway summary map overlays the pathway enrichment results of mock vs. high titer SARS-CoV-2-infected cells on top of the mock vs. low virus titer infected cells to show processes exclusively (single-color nodes) or commonly (double-color nodes) enriched between the two datasets. This analysis confirmed that perturbation in autophagy, inflammation, and mitochondrial processes were enriched by DE genes from both datasets (i.e., high and low MOI infected A549 cells) (Figure 2G).

Correlation of Expression Profiles Between Cell Lines and Lung Cell Types

The results presented in the section above and in subsequent sections used gene expression profiles of A549 (that were transduced with hACE2) and Calu3 lung cell lines. How close these cell lines are to lung cells was assessed by evaluating the correlation of the gene expression profiles of cell lines with lung cells using the scRNA-seq information from lung. First, using the scRNA-seq data, cell-type expression profile was computed as the mean expression across cells within each cell type. The top 1,000 genes with the highest variance among the 57 cell-type expression profiles were selected as highly variable genes, which were presumably informative for differentiating the 57 cell types. Next, the expression profiles of lung cell lines were compared with the expression profile of (hACE2-transduced) A549 and Calu3 cell lines. Correlation between the highly variable genes from lung scRNA-seq data and either A549 or Calu3 cells was calculated and plotted (Supplementary Figure 1E). This analysis revealed that the hACE2-transduced A549 cell gene expression strongly correlated with the basal and ionocyte lung cell subpopulations, which both represent lung epithelial cell lineage (Morrisey, 2018; Schiller et al., 2019). These data suggest that the A549 cells' response to SARS-CoV-2 infection likely reflects lung cells' response to the same infection. The Calu3 cells showed a similar pattern but a lower correlation with the lung cell types analyzed (Supplementary Figure 1E).

The correlation analysis (presented above) suggests that the BPs impacted in SARS-CoV-2-infected A549 cells is likely impacted in the SARS-CoV-2-infected lung epithelial cells too. However, given the limitations of analyzing lung cell line data, gene expression analysis of lung samples from patients with severe or mild COVID-19 will help test if these processes are differently impacted depending on the severity of the disease. Together, these results support that SARS-CoV-2 infection impacts the expression of genes involved in the cytokine signaling, autophagy, and mitochondria/respiration.

Comparing SARS-CoV-2 Infection in hACE2-Transduced A549 and Calu3 Cell Lines

A number of studies have been published that focused on identifying receptors used by SARS-CoV-2 to delineate viral entry mechanisms. Several of these studies have identified Angiotensin-Converting Enzyme 2 (ACE2) as the receptor that interacts with the SARS-CoV-2 spike protein to mediate viral entry (Li et al., 2007; Shang et al., 2020). Furthermore, TMPRSS2 and TMPRSS4, which are two membrane-bound serine proteases, were found to facilitate viral entry into the cells (Iwata-Yoshikawa et al., 2019; Hoffmann et al., 2020; Zang et al., 2020). Analysis of the gene expression datasets of A549 and Calu3 cells revealed that ACE2 and TMPRSS2 genes are highly expressed in the latter and showed robust gene expression changes in response to SARS-CoV-2 infection (Supplementary Table 3: Sheet 3). To facilitate SARS-CoV-2 infection, A549 cells were transduced with hACE2 vector (Blanco-Melo et al., 2020). The gene expression profile of SARS-CoV-2-infected A549 and Calu3 cells was compared and the correlation of gene expression between infected A549 and Calu3 cells was determined. We found significant correlation (R = 0.68, p < 2.2e−16) between gene expression of SARS-CoV-2-infected A549 and Calu3 cells (Supplementary Figure 2A). Furthermore, 65% of DE genes from mock vs. SARS-CoV-2 in Calu3 comparison overlapped with DE genes from the respective A549 comparison (Supplementary Figure 2B). Finally, we plotted a pathway enrichment summary map by using the pathway enrichment results from mock vs. SARS-CoV-2 comparison in A549 and Calu3 cells (Supplementary Figure 2C). Overlaying the pathway analysis results from A549 over Calu3 revealed an overlap of a wide range of BPs including the interferon, neutrophils, mitochondrial, and autophagy processes between the two datasets (Supplementary Figure 2C). Together, these results show that, upon ACE2 expression, the gene expression changes in infected A549 cells is highly correlated with infected Calu3 cells.

Network Analysis Identified Protein–Protein Interaction Subnetworks of Genes Involved in Interferon, Inflammation, and Mitochondrial Translation

SARS-CoV-2 vs. Mock: Network Analysis

To further understand the potential BPs in play during SARS-CoV-2 infection, we performed a consensus weighted gene coexpression network analysis (WGCNA) (Langfelder and Horvath, 2008) on combined, batch-corrected (see Methods), gene expression values of mock and SARS-CoV-2 infected at low and high titer cells, to identify clusters/modules of correlated gene. WGCNA identified more than 47 coexpression modules. The overlap of genes in each of these modules with significant DE genes from mock vs. SARS-CoV-2 infected at high titer is presented in the cluster dendrogram where each correlated module is represented by a color, and their overlap with DE genes is shown in horizontal bars (Figure 3A, Supplementary Table 1).

FIGURE 3
www.frontiersin.org

Figure 3. Consensus network analysis of hACE2-transduced mock and SARS-CoV-2-infected (high and low MOI) A549 cells. (A) Cluster dendrogram showing correlated genes grouped into clusters marked by different colors on the horizontal block labeled “Consensus Module Colors.” The DE genes in each cluster is marked as a black vertical line in the horizontal block labeled “DE genes.” The up- and downregulated genes are shown as red and blue vertical lines in a block labeled “Upregulated/downregulated genes,” respectively. (B) Top pathways from the pathway enrichment analysis of correlated DE genes in the blue module is presented as a horizontal bar plot, where the x axis represents the –log10-transformed q-value and the color of the horizontal bar is scaled blue to red representing low to high q-values, respectively. (C) Protein–protein interaction (PPI) subnetworks in the blue module is presented where each node represents a gene and the border color of the nodes indicate upregulation (red) and downregulation (blue) in SARS-CoV-2-infected A549 cells (high MOI) compared to mock-infected cells. The edge between the nodes indicate interaction based on the GeneMANIA database information. The network shows a highly connected interactome of interferon-stimulated genes (ISGs) that are coordinately upregulated in infected cells. (D) Another PPI subnetwork identified in the blue module shows several highly interconnected mitochondrial ribosomal (MRP) and complex I (NDUF) genes. Each node represents a gene and the border color of the nodes indicate upregulation (red color) and downregulation (blue) in SARS-CoV-2-infected A549 cells (high MOI) compared to mock-infected cells. Additional interactome of MYD88 TRAF3 and TLR3 genes that were coordinately upregulated in infected cells can be seen (top left). (E) NDUFV1 transcript level reported as fold change in beta-coronavirus-infected HCT-8 cells compared to control (no infection) HCT-8 cells. Student's t-test p-value = 2E-07. (F) Top pathways from the pathway enrichment analysis of correlated DE genes in the turquoise module is presented as a horizontal bar plot, where the x axis represents the –log10-transformed q-value and the color of the horizontal bar is scaled blue to red representing low to high q-values, respectively. (G) MRPL20 and MRPL43 transcript levels reported as fold change in beta-coronavirus-infected HCT-8 cells compared to control (no infection) HCT-8 cells. Student's t-test p-value = 8E-08 and 0.007. (H) A PPI subnetwork of correlated DE genes in the turquoise module shows a well-connected interactome of genes encoding mitochondrial ribosomal proteins, mitochondrial coiled-coil-helix-coiled-coil-helix domain proteins, and cytochrome oxidase. Each node represents a gene and the border color of the nodes indicates upregulation (red color) and downregulation (blue) in SARS-CoV-2-infected cells (high MOI) compared to mock-infected cells. DE, differentially expressed; MOI, multiplicity of infection.

First, pathway enrichment analysis of the correlated DE genes in the blue module showed significant annotation to the mitochondria, immunity, and mRNA/transcription processes (Figure 3B, Supplementary Table 4: Sheet 5). Using the GeneMANIA (Franz et al., 2018) database, protein–protein interaction (PPI) subnetworks for the DE genes in this module/cluster were identified. This analysis identified two PPI subnetworks of genes involved in interferon signaling (Figure 3C), TLR3 and MYD88 (Figure 3D), and mitochondrial translation and complex I genes (Figure 3D). Consistent with the expression level of complex I genes in the RNA-seq data, the mRNA level of NDUFV1 was down in beta-coronavirus-infected HCT-8 cells when measured by quantitative real-time polymerase chain reaction (qRT-PCR) (Figure 3E). Together, from these data, we concluded that the interferon signaling and inflammation genes were upregulated, and mitochondrial genes were downregulated in SARS-CoV-2-infected cells.

Next, pathway enrichment analysis of DE genes from the turquoise module revealed significant annotation to viral gene expression and apoptosis processes (Figure 3F, Supplementary Table 4: Sheet 6). Using GeneMANIA database, a PPI subnetwork of genes encoding the mitochondrial ribosomal proteins that were mostly downregulated in SARS-CoV-2-infected cells was also identified (Figure 3G). Consistent with the RNA-seq data, mRNA levels of the mitochondrial ribosomal genes MRPL20 and MRPL43 were downregulated in beta-coronavirus-infected HCT-8 cells when measured by qRT-PCR (Figure 3H). Together, these data suggest that SARS-CoV-2 infection results in a coordinated change in the interferon signaling, inflammation, and mitochondrial processes.

Gene Expression Changes Associated With SARS-CoV-2 Infection Is Distinct From IAV Infection With Minor Overlaps

SARS-CoV-2 vs. IAV

To compare the expression profile of SARS-CoV-2-infected cells with other virus-infected cells, DE analysis of mock- vs. IAV-infected cells was performed and the up- and downregulated genes are presented in a volcano plot (Supplementary Figure 3A, Supplementary Table 3: Sheet 4). The pathway analysis of the DE genes from this comparison showed enrichment in protein translation, localization, and anti-viral responses (Supplementary Figure 3B, Supplementary Table 4: Sheet 7). Additionally, genes upregulated in the IAV-infected cells annotated to pathways for virus response, protein trafficking, and unfolded protein response (UPR) (Figure 4A, Supplementary Table 4: Sheet 8). Genes that were downregulated in IAV-infected cells enriched in vacuole- and lysosome-related processes (Figure 4B, Supplementary Table 4: Sheet 9). Interestingly, few DE genes from the mock vs. SARS-CoV-2 overlapped with the DE genes from mock vs. IAV comparison (Supplementary Figure 3C).

FIGURE 4
www.frontiersin.org

Figure 4. SARS-CoV-2 infection of A549 lung epithelial cells results in distinct gene expression changes that are not seen in IAV infection. (A) Top 25 pathways from the pathway enrichment analysis of DE genes upregulated in IAV-infected cells is presented as a horizontal bar plot, where the x axis represents the –log10-transformed q-value and the color of the horizontal bar is scaled blue to red representing low to high q-values, respectively. (B) Top 25 pathways from the pathway enrichment analysis of DE genes downregulated in IAV-infected cells is presented as a horizontal bar plot, where the x axis represents the –log10-transformed q-value and the color of the horizontal bar is scaled blue to red representing low to high q-values, respectively. (C) Pathway enrichment summary map for mock vs. SARS-CoV-2 (blue nodes) and mock vs. IAV (red nodes) comparisons. Each node represents a pathway/biological process (BP). The node size is proportional to the number of DE genes overlapping with the BP. The nodes that share genes are connected with edges. The label above each black circle summarizes the gene ontology (GO) terms of similar BPs present inside the circle. Single-color nodes are pathways that are distinctly enriched by DE genes from one comparison. Two-colored nodes are pathways enriched by DE genes from both comparisons. Notable groups of BPs associated with NFkB, ROS, autophagy and mitochondria are highlighted in red. DE genes from Mock vs. SARS-CoV-2 comparison exclusively enriched in highlighted BPs.

A pathway enrichment summary map was created by overlaying the pathway enrichment results of the mock vs. IAV comparison on top of the mock vs. SARS-CoV-2 comparison. Consistent with the DE genes comparison (Supplementary Figure 3C), the enrichment map also highlighted little overlap of pathways between the two comparisons. DE genes from both comparisons commonly enriched in a subset of pathways associated with protein trafficking (Figure 4C). Furthermore, only a subset of the interferon pathway genes and a few chemokine genes that were upregulated in SARS-CoV-2-infected cells were also upregulated in IAV-infected cells, while the autophagy and inflammation genes remained mostly unchanged in the latter. Therefore, upregulation of cytokine/inflammation, changes in autophagy, and downregulation of the mitochondrial processes were uniquely observed in SARS-CoV-2-infected cells. Upregulation of DE genes involved in the cytokine/inflammation processes is consistent with cytokine storm observed in severe cases of COVID-19 patients. Since these observations were made by analyzing the gene expression changes in a lung cell line, future studies profiling gene expression changes in severe COVID-19 patients' lung samples will be needed to confirm these findings.

SARS-CoV-2-Infected Cells Share Some Gene Expression Signature With MERS-CoV-Infected Cells With Few Exceptions

SARS-CoV-2 vs. MERS-CoV

Comparison of gene expression profiles revealed that SARS-CoV-2-infected cells are distinct from those of IAV-infected cells (Figure 4). Although these are both viruses, they are not phylogenetically close. Therefore, we next compared the gene expression profiles of SARS-CoV-2- and MERS-CoV-infected cells, since both are hCoVs. DE analysis of the mock- vs. SARS-CoV-2-infected Calu3 lung carcinoma cells identified several up- and downregulated genes (Supplementary Figure 4A, Supplementary Table 3: Sheet 3). Pathway enrichment analysis showed annotation of the DE genes to cell cycle, inflammation, and apoptosis processes (Supplementary Figure 4B, Supplementary Table 4: Sheet 10). A pathway enrichment summary map for mock vs. SARS-CoV-2 and mock vs. MERS-CoV comparisons was generated to assess the extent of overlap of pathways between the two datasets (Figure 5A). Notably, the DE genes from both comparisons enriched in the mitochondria, autophagy, cell cycle, and UPR processes. However, DE genes from mock vs. SARS-CoV-2 comparison predominantly enriched in inflammation, cytokine signaling, and immunity-related processes (Figure 5A). Consistently, genes upregulated in the SARS-CoV-2-infected Calu3 cells enriched in inflammation and nuclear factor kappaB (NFkB) processes (Figure 5B, Supplementary Table 4: Sheet 11), while upregulated genes from both hCoV-infected cells annotated to protein trafficking and small GTPase signaling (Figures 5B,C, Supplementary Table 4: Sheets 11 and 12). On the other hand, genes downregulated in both comparisons commonly annotated to mitochondrial processes (Figures 5D,E, Supplementary Table 4: Sheets 13 and 14). These findings suggest that perturbation of autophagy and mitochondrial genes are common gene expression signatures associated with hCoV infection, but the SARS-CoV-2 virus almost exclusively impacts the cytokine/inflammatory processes in the lung cells. It is likely that perturbation of mitochondrial processes and autophagy may lead to a dysfunctional immune response (Won et al., 2015; Jang et al., 2019). Further studies will be required to understand if and how these processes may contribute to inflammation during SARS-CoV-2 infection.

FIGURE 5
www.frontiersin.org

Figure 5. SARS-CoV-2 and MERS-CoV infection have some common and some distinct gene expression signatures. (A) Pathway enrichment summary map for mock vs. SARS-CoV-2 (blue nodes) and mock vs. MERS-CoV (red nodes) comparisons in Calu 3 cells. Each node represents a pathway/biological process (BP). The node size is proportional to the number of DE genes overlapping with the BP. The nodes that share genes are connected with edges. The label above each black circle summarizes the gene ontology (GO) terms of similar BPs present inside the circle. Single-color nodes are pathways that are distinctly enriched by DE genes from one comparison. Two-colored nodes are pathways enriched by DE genes from both comparisons. Notable groups of BPs associated with immunity, ROS, autophagy, and mitochondria are highlighted in red. These notable groups of BPs were commonly enriched by DE genes from both SARS-CoV-2 and MERS-CoV comparisons. The DE genes from mock vs. SARS-CoV-2 comparison predominantly enriched in inflammation and immunity-related processes. (B) Top 25 pathways from the pathway enrichment analysis of DE genes upregulated in SARS-CoV-2-infected Calu3 cells is presented as a horizontal bar plot, where the x axis represents the –log10-transformed q-value and the color of the horizontal bar is scaled blue to red representing low to high q-values, respectively. (C) Top 25 pathways from the pathway enrichment analysis of DE genes upregulated in MERS-CoV-infected Calu3 cells is presented as a horizontal bar plot, where the x axis represents the –log10-transformed q-value and the color of the horizontal bar is scaled blue to red representing low to high q-values, respectively. (D) Top 25 pathways from the pathway enrichment analysis of DE genes downregulated in SARS-CoV-2-infected Calu3 cells are presented as a horizontal bar plot, where the x axis represents the –log10-transformed q-value and the color of the horizontal bar is scaled blue to red representing low to high q-values, respectively. (E) Top 25 pathways from the pathway enrichment analysis of DE genes downregulated in MERS-CoV-infected Calu3 cells is presented as a horizontal bar plot, where the x axis represents the –log10-transformed q-value and the color of the horizontal bar is scaled blue to red representing low to high q-values, respectively. DE: differentially expressed.

Validation of Findings From Cell Line Data in SARS-CoV-2-Positive Human Samples

Comparing SARS-CoV-2-Positive vs. -Negative Human Nasopharyngeal Expression Profile With A549 Dataset

The gene expression profiling of the A549 and Calu3 cell line data revealed that perturbations in inflammatory, autophagy, and mitochondrial processes were unique to coronavirus infections. To validate the findings from SARS-CoV-2-infected lung cell lines, we analyzed the gene expression profiles of SARS-CoV-2-positive human nasopharyngeal samples. DE analysis of SARS-CoV-2-positive vs. -negative samples revealed up- and downregulated genes (Supplementary Table 3: Sheet 5). Concurrent with the cell line data, the pathway enrichment analysis of the DE genes revealed significant annotation to inflammation, autophagy, and mitochondrial processes (Supplementary Figure 5A, Supplementary Table 4: Sheet 15). Approximately 60% of the DE genes from the SARS-CoV-2-infected human samples comparison overlapped with the infected A549 cell data (Supplementary Figure 5B). Of this, ~50% of the DE genes were concurrently downregulated in both datasets (i.e., human samples and A549 cells) and significantly annotated to the autophagy, immunity, and mitochondrial processes (Supplementary Figure 5C). Since age (Wang D. et al., 2020) and viral load (Liu et al., 2020) can determine the severity of the COVID-19 outcome, the human samples were further classified into young (<40 years) and old (>60 years) with low or high viral load positive samples. Consistent with the SAR-CoV-2-positive vs. -negative human data, the DE genes from the control vs. high viral load in both old and young subjects significantly annotated to inflammation, autophagy, and mitochondrial processes [Figure 6A (old subject), Supplementary Table 3: Sheets 7 and Sheet 9, and Supplementary Table 4: Sheets 16 and 17]. When the pathway enrichment result from the old human subjects (that were high viral load positive) was overlaid on the pathway enrichment results from A549 cells, transduced hACE2 infected with high MOI of SARS-CoV-2, autophagy, NFkB, oxidative stress, and mitochondrial processes were commonly perturbed in both datasets (Figure 6B). Together, these observations suggest that across cell types, SARS-CoV-2 infection alters the cells' autophagy and mitochondrial processes and that perturbations in these processes may impede an effective immune response leading to severe outcomes.

FIGURE 6
www.frontiersin.org

Figure 6. Gene expression profiling of human nasopharyngeal samples. (A) Top 25 pathways from the pathway enrichment analysis of DE genes from control vs. high viral load positive old age human samples is presented as a horizontal bar plot, where the x axis represents the –log10-transformed q-value and the color of the horizontal bar is scaled blue to red representing low to high q-values, respectively. (B) Pathway enrichment summary map for control vs. high viral load positive human samples (blue nodes) and mock vs. SASR-CoV-2 (red nodes) comparisons in hACE2-transduced A549 cells. Each node represents a pathway/biological process (BP). The node size is proportional to the number of DE genes overlapping with the BP. The nodes that share genes are connected with edges. The label above each black circle summarizes the gene ontology (GO) terms of similar BPs present inside the circle. Single-color nodes are pathways that are distinctly enriched by DE genes from one comparison. Two-colored nodes are pathways enriched by DE genes from both comparisons. Notable groups of BPs associated with immunity, autophagy, and mitochondria are highlighted in red. These notable groups of BPs were commonly enriched by DE genes from both comparisons. (C) Heatmap of the mean expression values of the indicated genes in young and old human samples that were negative (control) or positive with either high or low viral loads of SARS-CoV-2 virus. The indicated genes are broadly grouped into four different processes using a vertical bar present on the right side of the heatmap. The orange and blue color bands represent upregulated and downregulated genes, respectively. (D) CCL4 and CXCR5 chemokine and receptor transcript levels reported as fold change in beta-coronavirus-infected HCT-8 cells compared to control (no infection) HCT-8 cells. Student's t-test p-value = 3E-04 and 1.7E-05.

Delineating the Age and Viral Load Impact on SARS-CoV-2 Response in Human Nasopharyngeal Samples

The human nasopharyngeal samples were classified into young and old samples with low or high viral load to delineate either age specific or viral load specific gene expression changes in response to infection. The analysis in control vs. infected (with either low or high viral load) young or old samples helped identify DE genes that significantly annotated to autophagy, neutrophils, and mitochondrial processes (Supplementary Table 3: Sheets 6–9), which is consistent with the pathway enrichment results from SARS-CoV-2-infected A549 cell line data (Figure 2). However, in-depth gene expression analysis of the human data revealed some concurrent gene expression changes between cell line and human datasets, and some sample specific changes. Genes encoding the mTOR (Figure 6C), mitochondrial ribosomal genes, mitochondrial complex I genes, and lysosome acidification genes were downregulated in all of the SARS-CoV-2-positive human samples as well as in SARS-CoV-2-infected A549 cells (Supplementary Table 3: Sheet 1). Additionally, autophagy initiation, nucleation genes, p62 (SQSTM1), SNAP29, and MITF were specifically downregulated in SARS-CoV-2-positive human samples, which indicates decreased autophagic flux in infected samples (Figure 6C and Supplementary Figure 6D). Furthermore, antioxidant encoding SOD1 (Fukai and Ushio-Fukai, 2011) gene, the oxidative stress sensor PARK7 (DJ-1) (Wang et al., 2016), and mitochondria fission promoting SOCS6 (Lin et al., 2013) were also downregulated in infected human samples (Figure 6C). Finally, several cytokine/inflammation process genes were downregulated in infected human samples, which was opposite to the upregulation of these genes in SARS-CoV-2-infected A549 cells (Supplementary Table 3: Sheet 1).

There were several age- and viral load-dependent gene expression changes observed in the human nasopharyngeal samples, which are described below. Several interferon pathway genes were upregulated in high viral infected young and old samples. However, the upregulation of these genes in high viral load positive old samples was muted compared to the corresponding young samples (Supplementary Figure 6A). The difference in interferon genes induction was most prominent in low viral load positive samples, as unlike the young samples where most of these genes were upregulated, in old samples, these interferon genes were significantly downregulated (Supplementary Figure 6A). Similarly, ADAP2 [which is an interferon-stimulated gene (ISG)], antiviral tripartite motif family E3 ligase TRIM38, chemokine CCL4 and its receptor CCR5, and WARS1 encoding the Tryptophanyl-tRNA Synthetase were all robustly down in old, low viral dose positive samples compared to the age-matched control (Figure 6C). In response to the high viral titer infection, both CCL4 and CXCR5 were upregulated more robustly in young subjects compared to old subjects. Consistent with the expression in the nasopharyngeal samples, robust upregulation of CCL4 and CXCR5 mRNAs was also observed in HCT-8 cells infected with the beta-coronavirus compared to the control cells by qRT-PCR (Figure 6D). Upregulation of chemokine signaling may be a common host cell response to coronaviruses. Together, these data suggest age-specific differences in host transcriptional response to SARS-CoV-2 infection. However, further analysis of the infected lung samples from old and young patients will be required to test if these differences may be causing more severe outcomes in older patients.

Gene Expression Analysis of a Severe COVID-19 Lung Sample Shows Exaggerated Immune/Inflammation Response

Healthy vs. COVID-19 Lung Biopsy Samples

The gene expression analysis of the SAR-CoV-2-infected A549 and Calu3 cell lines revealed upregulation of the cytokine/inflammatory processes. However, this observation was inconsistent with the human nasopharyngeal expression profile where inflammatory process genes were downregulated. This difference could arise due to disparate cell type being compared, or that severe inflammatory symptoms may arise in severe case of COVID-19. Since the human nasopharyngeal data did not include clinical symptoms indicating the severity of COVID-19 in the samples whose sequence were analyzed, we could not stratify the samples based on disease severity. However, to test if the cytokine/inflammatory processes were also impacted in COVID-19 lungs, RNA-seq data from the healthy and COVID-19 lung biopsy were analyzed. It should be noted that COVID-19 lung biopsies were technical replicates and therefore statistical significance of this analysis is limiting. Future studies involving bigger sample size will be required to confirm these observations. Nevertheless, the gene expression profile of a SARS-CoV-2-infected lung was distinct from healthy lungs with up- and downregulated genes highlighted in the volcano plot (Supplementary Figure 7A, Supplementary Table 3: Sheet 10). The pathway enrichment summary map and the plot showed that the DE genes predominantly annotated to the inflammation, ROS, and leukocyte/monocyte-related pathways (Supplementary Figure 7B). Furthermore, the DE genes upregulated in the COVID-19 lungs enriched in the anti-viral response processes, cytokine secretion, immune cell proliferation/migration, and inflammation (Supplementary Figure 7C, Supplementary Table 4: Sheet 18). The downregulated genes were significantly enriched in protein trafficking, RNA metabolism, and oxygen-sensing processes (Supplementary Figure 7D, Supplementary Table 4: Sheet 19). It is likely that perturbations in the oxygen-sensing processes are reflective of the severe respiratory distress often seen in severe COVID-19 patients due to reduced oxygenation ability of the failing lungs.

Discussion

Highly pathogenic hCoVs are known to infect the lower respiratory airways and cause severe acute respiratory syndrome (SARS) (Channappanavar and Perlman, 2017). The recently discovered SARS-CoV-2 virus is the cause of COVID-19 (Tammaro et al., 2020). The clinical manifestations of this virus infection include fever, cough, fatigue, respiratory distress, and cardiac injury (Chen N. et al., 2020; Guan et al., 2020; Huang et al., 2020). While some patients with COVID-19 suffered from mild symptoms, other patients had increasingly life-threatening symptoms (Guan et al., 2020). Age and underlying medical conditions such as diabetes and hypertension are likely to determine the severity of the symptoms (Wang D. et al., 2020). Analyzing the gene expression profiles of host cells infected with SARS-CoV-2 will be necessary to decipher the subcellular functions perturbed by this virus and to inform drug development strategies.

Here, we present an in-depth differential expression analysis of A549 and Calu3 cell lines, comparing mock to infection with either SARS-CoV-2, IAV, or MERS-CoV. Since A549 cells lacked expression of ACE2, TMPRSS2, or TMPRSS4 that is required for SARS-CoV-2 viral entry into the cells (Iwata-Yoshikawa et al., 2019; Hoffmann et al., 2020; Zang et al., 2020), A549 cells transduced with hACE2 were used. Upon SARS-CoV-2 infection at low and high MOI, viral transcripts were detected in these cells indicating infection (Blanco-Melo et al., 2020). Concurrently, we also observed strong correlation between SARS-CoV-2-infected Calu3 and hACE2-transduced A549 cells, suggesting that SARS-CoV-2 infection of the lung cell lines results in similar gene expression profile. Furthermore, we report that (i) SARS-CoV-2 infection impacted the expression of genes in inflammation, cell cycle, reactive oxygen species (ROS), autophagy, and mitochondrial processes, which were absent in IAV-infected cells; (ii) while perturbation in autophagy and mitochondrial processes is common in hCoV infections (SARS-CoV-2 and MERS-CoV), we found that increased expression of the inflammatory/cytokine signaling genes was exclusively observed in SARS-CoV-2-infected lung cells; (iii) coexpression network analysis helped identify a cluster of genes involved in inflammation and mitochondrial translation process that were either coordinately up- or downregulated in SARS-CoV-2-infected cells, respectively. Together, these data suggest that perturbation in the autophagy, mitochondrial processes in SARS-CoV-2-infected lung cells could hinder an effective immune response (Won et al., 2015; Jang et al., 2019) and increase inflammation, which is often seen in severe COVID-19 patients suffering from cytokine storm (Channappanavar and Perlman, 2017; Mehta et al., 2020). Since these conclusions were made using the data from virus-infected lung cell lines, the correlation between these cells' expression profiles and marker gene expression from different lung cell types were determined. While the A549 cells showed robust correlation with lung epithelial lineage basal and ionocyte cells, Calu3 cells showed a similar pattern but lower correlation with these cell types. Therefore, the processes delineated in SARS-CoV-2 A549 cells likely represent the lung epithelial cells response to SARS-CoV-2 infection. To further substantiate the cell line findings, the gene expression profiles of SARS-CoV-2-positive human nasopharyngeal samples were used for validation. Analysis of control vs. SARS-CoV-2-positive samples helped identify DE genes that significantly annotated to the autophagy, NFkB, oxidative stress, and mitochondrial processes. Using the patient information available from this dataset, the samples were grouped into young (<40 years) and old (>60 years) that were either positive with low or high viral load. Comparing the gene expression and pathway enrichment results of old age samples with that of A549 high MOI data revealed a wide range of BPs that were commonly perturbed in both datasets. Therefore, this analysis has delineated several BPs, discussed in more detail below, that are impacted in the SARS-CoV-2-infected host cells.

Some complement genes (C1S and C1R) were specifically upregulated in high viral titer SARS-CoV-2-infected cell lines. Consistently, the C1q/TNF-related protein 6, a glycoprotein that regulates complement activation, was downregulated in both SARS-CoV-2-infected cells and human samples. This gene is implicated in arthritis, and intra-articular injection of the recombinant C1qTNF6 protein was shown as an effective strategy in improving arthritis and inflammation in C1qtnf6–/– mice (Murayama et al., 2015). An elevated complement response could likely lead to excessive inflammation, which was also observed in MERS-CoV infection of the hDPP4-transgenic mouse model (Jiang et al., 2018). Additionally, several past studies have highlighted the interplay between the complement and coagulation systems (Skoglund et al., 2010; Oikonomopoulou et al., 2012). It is likely that the increased thrombosis in COVID-19 patients (Bikdeli et al., 2020) is a result of excessive complement activation. Further assessment of complement activation in COVID-19 patients will be required to confirm this. Together, these observations suggest that inhibition of the complement system as potential treatment strategies could be tested.

Infection of A549 cells with SARS-CoV-2 at higher viral titer perturbed autophagy; upregulated genes in the interferon, cytokine, nuclear factor kappaB (NFkB), and reactive oxygen species (ROS) processes; and downregulated genes in the mitochondrial and electron transport chain processes. Consistently, analysis of DE genes in one of the correlated clusters from WGCNA showed significant enrichment in the interferon signaling processes. Additionally, GeneMANIA analysis of the correlated DE genes in two modules revealed PPI subnetworks of genes involved in ISGs and NFkB, which were both mostly upregulated in the infected cells. In addition to the ISGs, the JAK-STAT signal transduction genes, which play critical role in type I cytokine (such as IL6) signaling and inflammation (Leonard and O'shea, 1998; O'brown et al., 2015; Banerjee et al., 2017), were also upregulated in the SARS-CoV-2-infected A549 cells (Figure 4C). IL6, a pleotropic cytokine, was shown to be elevated in critically ill COVID-19 patients (Chen X. et al., 2020). Consistently, IL6 was upregulated in the SARS-CoV-2-infected cells. IL6 acts via the JAK-STAT signaling through SOCS3 protein kinase (also upregulated in SARS-CoV-2-infected cells) to activate the immune response (Brocke-Heidrich et al., 2004). Excessive IL6 causes excessive inflammation as seen in arthritis (Srirangan and Choy, 2010). However, the inflammatory/cytokine gene expression profile in the virus-infected human nasopharyngeal samples was distinct from the cell line data. In the infected nasopharyngeal samples, most of the cytokine/inflammatory process genes were significantly downregulated compared to the control. These discordant inflammatory gene expression profiles between two datasets may be due to different cell types being compared (A549 is a lung cell line, and human nasopharyngeal samples are predominantly squamous epithelial cells), or the upregulation of cytokine/inflammatory processes genes in SARS-CoV-2-infected A549 cells may represent a severe COVID-19 infection state. Due to lack of clinical information describing the disease state of subjects whose nasopharyngeal samples were analyzed, we could not test this possibility in the nasopharyngeal dataset. However, analysis of the COVID-19 lung biopsy samples revealed significant upregulation of genes enriched in cytokine/inflammatory processes. Therefore, upregulation of IL6 and NFkB genes may contribute to the inflammatory symptoms observed in severe COVID-19 patients (Channappanavar and Perlman, 2017; Mehta et al., 2020). These data support a central role for cytokine signaling in COVID-19 pathogenesis. Treatment strategies aimed at mitigating the cytokine effects or complement system could be tested in treatment of COVID-19. One such clinical trial aimed at mitigating the IL6 effects is already underway (NCT04322773). Furthermore, another study showed decreased mortality in patients treated with tocilizumab, which blocks IL6 (Somers et al., 2020).

Analysis of young (<40 years) and old (>60 years) nasopharyngeal samples also revealed some age-specific changes in the gene expression profile. Notably, the ISGs (IFIT1, IFIT2, and IFIT3) were upregulated in both low and high viral load SARS-CoV-2-positive young samples. Upregulation of these genes in high viral load infected old age samples were, however, less robust. Furthermore, in the low viral load positive old samples, most of the interferon genes (except IFIT3) were downregulated. In addition to the interferon genes, ADAP2, which is an ISG (Shu et al., 2015), TRIM5, which is a retroviral restriction factor (De Silva and Wu, 2011), tripartite motif TRIM22, and TRIM38 were significantly downregulated in low viral load positive old samples compared to age-matched control. The latter two genes are involved in innate immunity and in restricting viral infections (Barr et al., 2008; Lian and Sun, 2017). Moreover, the Tryptophanyl-tRNA Synthetase encoding WARS1, which stimulates immunity against viral infection (Lee et al., 2019), chemokines CXCL11, CCL4, and CCL4 receptor CCR5 were also significantly downregulated in low viral load positive old samples compared to age-matched control. Chemokine CXCL10 expression was unchanged in the low viral load positive old samples. All these chemokines and the receptor were upregulated in the low viral load positive young samples. Consistent with this, qRT-PCR analysis of the CCL4 and CXCR5 transcript levels in beta-coronavirus-infected HCT-8 cells also showed a robust upregulation of these genes in infected cells. It should be noted that upregulation of CXCL10, CXCL11, and IFIT2 in the nasopharyngeal samples has been proposed to accurately predict the presence of respiratory virus infection (Landry and Foxman, 2018). Lack of induction of these genes in old age infected subjects may indicate an inability of the host cells to detect viral entry. This, in combination with downregulation of genes involved in antiviral immunity, likely contributes to severe disease outcomes, which is consistent with severe COVID-19 manifestations in older patients. However, more studies using the lung biopsy samples of SARS-CoV-2-infected young and old age patients will be required to confirm these findings.

What processes may be causing/contributing to defective immune response in the SARS-CoV-2-infected cells? Autophagy- and mitochondria-related processes were two other prominent categories of the BPs that were exclusively impacted in the hCoV-infected cells.

Gene expression analysis of SARS-CoV-2-infected A549 and Calu3 cells revealed differential expression of genes involved in autophagic processes. In contrast, most of the DE genes from human nasopharyngeal sample comparison that annotated to autophagy processes were downregulated. SARS-CoV-2-infected A549 and Calu3 cell lines and human nasopharyngeal samples expressed significantly low levels of mTOR and LAMTOR genes. The regulatory associated protein of mTOR (RAPTOR) expression was also significantly decreased in all SARS-CoV-2-infected samples. mTOR inhibits autophagy (Kapuy et al., 2014; Rabanal-Ruiz and Korolchuk, 2018). Consistent with the mTOR function, the autophagy inducing microphthalamia-associated transcription factors (MiTF/TFE) and unc-51-like autophagy activating kinase 1 (ULK1) were upregulated in the infected A549 cells. Concurrent with these data, inhibition of mTORC1 in SARS-CoV-2-infected human bronchial epithelial cells NCI-H1299 and monkey kidney cells (VeroFM) has also been reported (Gassen et al., 2020). These observations suggest that downregulation of mTOR may result in autophagy induction upon SARS-CoV-2 infection in A549 cells. However, contrary to the expression profiles of infected A549 cells, the autophagy initiation and nucleation genes were downregulated in infected human nasopharyngeal cells compared to control samples. It is worth noting that a past study in mouse embryonic cells showed that coronavirus infection induced autophagy and that coronavirus mouse hepatitis virus (MHV) replication was impaired in atg5–/– cells (Prentice et al., 2004). While the gene expression data from SARS-CoV-2-infected cell lines and human samples indicate changes in autophagy gene expression, very little can be inferred about the changes to the autophagic flux, which is the capacity of the cells to degrade the autophagosome by fusing with the lysosomes. Studies have shown that dengue and enteroviral infections inhibit autophagic flux in host cells by decreasing p62 or cleaving SNAP29, respectively, to facilitate infection (Metz et al., 2015; Mohamud et al., 2018). Consistent with this observation, we found that the autophagic flux inducing p62 (SQSTM1) and SNAP29 gene expression level were down in SARS-CoV-2-positive human samples. Moreover, inhibition of S-phase kinase-associated protein 2 (SKP2) as E3 ligase decreased Beclin1 (BECN1) degradation and increased autophagic flux, which, in turn, decreased MERS-CoV (Gassen et al., 2019; Yang and Shen, 2020). However, SKP2 expression in both SARS-CoV-2-infected cell lines and human samples was significantly low. Nevertheless, drugs known to increase autophagic flux has been shown to impede SARS-CoV-2 infection (Gassen et al., 2020; Gorshkov et al., 2020). It is likely that autophagic flux is decreased in human samples by downregulation of p62 and SNAP29 and in SARS-CoV-2-infected A549 cells potentially through upregulation of glycogen synthase kinase 3 beta (GSK3B), which has been shown to impair lysosome acidification (Weikel et al., 2016). Consistently, GSK3B inhibition has been shown to increase autophagy flux in mice liver and human pancreatic cancer cells (Marchand et al., 2015; Ren et al., 2016). Additionally, several lysosome acidification genes were concordantly downregulated in both SARS-CoV-2-infected A549 cells and human samples. Since impaired lysosome acidification has been associated with impaired autophagic flux (Yim and Mizushima, 2020), drugs targeting lysosome reacidification or increasing autophagic flux could potentially be tested as a therapeutic intervention to SARS-CoV-2 infection. Another independent study also showed that beta-coronavirus deacidifies lysosomes as a mechanism to facilitate viral infection and egress from the host cells (Ghosh et al., 2020). Together, these findings suggest that autophagy and regulation of autophagic flux may be central to SARS-CoV-2 infection and that their differential regulation may be key to the control of the infection.

In addition to changes in expression of genes involved in autophagic processes, several genes involved in the mitochondrial processes were downregulated in SARS-CoV-2-infected A549 cells and human samples. This is consistent with our current understanding that viruses either induce or inhibit various mitochondrial processes as part of their replication and dissemination efforts (Anand and Tikoo, 2013). Infection of cells with SARS-CoV-2 at higher viral titer downregulated the genes in the mitochondrial processes. Consistently, WGCNA of SARS-CoV-2-infected A549 datasets also identified a DE gene cluster that annotated to the mitochondrial organization and translation processes. Subsequent GeneMANIA analysis identified a PPI subnetwork of genes involved in mitochondrial translation, which were coordinately downregulated in SARS-CoV-2-infected cells. Since mitochondrial import and translation are interlinked (Mokranjac and Neupert, 2010; Sanchez-Caballero et al., 2016), we found that several mitochondrial complex I and translocase genes were downregulated in the SARS-CoV-2-infected cells. Given the extensive crosstalk between autophagy and mitochondrial function (Graef and Nunnari, 2011; Rambold and Lippincott-Schwartz, 2011), it is likely that perturbations in autophagy and mitochondrial processes observed in SARS-CoV-2-infected cells are interlinked. It is worth noting that the mTORC1 complex, which comprises the mTOR and RAPTOR, stimulates synthesis of mitochondrial ribosomal, complex I proteins (Morita et al., 2013, 2015) and mitochondrial fission process 1 (MTFP1) (Morita et al., 2017). Consistently, we found decreased expression of mitochondrial ribosomal and complex I genes, which is likely a result of decreased mTOR and RAPTOR expression in SARS-CoV-2-infected cells. Additionally, decreased expression of MTFP1 that may impede mitochondrial fission resulting in hyperfused mitochondria (Tondera et al., 2005) was specifically observed in SARS-CoV-2-infected A549 cells. While MTFP1 gene expression was undetectable in human nasopharyngeal samples, another mitochondrial fission promoting SOCS6 (Lin et al., 2013) was downregulated in these infected samples. Past studies have shown that hyperfused or elongated mitochondria in dengue- and SARS-infected cells can suppress interferon signaling and innate immunity (Shi et al., 2014; Barbier et al., 2017; Das and Chakrabarti, 2020). Furthermore, reduced complex I expression has been found in many cancer cells and is shown to affect the oxidative phosphorylation, which also impacts the immune cell function (Simonnet et al., 2003; Bonora et al., 2006; Baracca et al., 2010; Won et al., 2015; Angajala et al., 2018). Together, these data highlight that SARS-CoV-2-infected cells have decreased mTOR expression and perturbation in autophagy and mitochondrial processes, which, in turn, could properly impede the immune response to infection. Further studies will, however, be required to delineate which of these perturbations are the direct result of SARS-CoV-2 infection and/or contribute to pathogenesis/severe clinical manifestations.

Summary

In summary, we have presented a detailed DE and coexpression network analysis of the RNA-seq data from SARS-CoV-2-infected A549 cells. Using the gene expression profiles of A549 and Calu3 cells infected with IAV or MERS-CoV, we concluded that perturbations in cytokine signaling and inflammation processes, downregulation of genes in the mitochondrial processes, and perturbation of autophagy were uniquely observed in novel coronavirus-infected cells. To validate the findings from the cell line data, gene expression analysis of control and SARS-CoV-2-positive human nasopharyngeal samples was performed. Consistent with the cell line data, DE genes from human data also significantly annotated to inflammation, autophagy, and mitochondrial processes. It is likely that perturbation of autophagy and mitochondrial processes may impede an effective immune response leading to severe outcomes. Furthermore, age-stratified human nasopharyngeal was used to analyze gene expression changes in control vs. high viral load positive old age (>60 years) samples. This analysis also revealed several BPs that were concordantly impacted in both cell line and human datasets, with few differences. Specifically, genes encoding mTOR were downregulated in infected cells, which likely caused downregulation of mitochondrial ribosomal genes in both cell line and human datasets. Additionally, genes encoding mitochondrial complex I and lysosome acidification were also concurrently downregulated in infected cells from both datasets. Several inflammation process genes and autophagy genes were discordantly regulated in both datasets. It is likely that the autophagic flux is impeded in infected cells lines and human samples due to increased expression of GSK3B gene or downregulation of p62 and SNAP29, respectively, which may further promote viral propagation. These data suggest that drugs that enhance autophagic flux or increase lysosome acidification could be tested as intervention strategies. Increased expression of inflammation process genes in A549 cells is likely due to these cells representing a severe COVID-19 infection state. Consistently, we found upregulation of inflammatory process genes in expression profiles of the COVID-19 lung biopsy sample. These data support a central role for cytokine/inflammation processes in COVID-19 pathogenesis. Finally, using the age-stratified expression profiles of infected human nasopharyngeal samples, we identified muted or downregulation of several chemokines, ISGs, and tripartite motif genes that are critical for innate immunity and antiviral signaling. It is likely that defective antiviral response in old age patients in combination with perturbations in autophagy and mitochondrial processes could result in severe COVID-19 disease state often seen in older population. In summary, using gene expression data of SARS-CoV-2-infected cells, we show that viral infection of host cells results in perturbations in specific aspects of autophagy and mitochondrial processes. Future studies focusing on how these perturbations either contribute to viral propagation or impede an effective immune response will be required to gain more understanding of the viral pathogenesis.

Materials and Methods

Data Collection

Raw gene count matrix for bulk RNA-seq was obtained from GEO (accession number GSE147507) (Blanco-Melo et al., 2020). The data contained gene expression count matrix of two lung carcinoma cell lines A549 (Lieber et al., 1976; Holownia et al., 2016) and Calu3 (Shen et al., 1994; Martens et al., 2018). In this dataset, the A549 treatment conditions included mock, infection with IAV (N = 2 per group), and infection with SARS-CoV-2 at 2 (high titer, N = 3 per group) and 0.2 (low titer, N = 3 per group) multiplicity of infection (MOI), after transduction with a vector expressing human ACE2 (hACE2) (N = 3 per group). The adenovirus Ad-GFP-h-ACE2 from Vector Biolabs was used for transduction of A549 cells at 500 MOI (Blanco-Melo et al., 2020). Subsequently, cells were infected with SARS-CoV-2 [Isolate USA-WA1/2020 (NR-52281)] at 0.2 or 2 MOI as indicated (Blanco-Melo et al., 2020). From this dataset, the raw gene count matrix for two healthy human lung biopsy samples and one COVID-19 sample (two technical replicates) and Calu3 cells infected with SARS-CoV-2 at 2 MOI (N = 3 per group) were also analyzed. Since the A549 cells infected with SARS-CoV-2 at high and low MOI were used in the network analysis, and DE genes from A549 and Calu3 comparisons were used, a boxplot of the normalized counts showed that all samples were comparable with no outlier (data not shown). Additionally, gene expression data in FPKM were downloaded from GEO (accession number GSE139516) (Zhang et al., 2020) for the Calu3 cell line infected with MERS-CoV and mock (N = 3 per group). Differential expression analysis was performed on mock- and MERS-CoV-infected cells for 24 h. Human lung single-cell RNA-seq (scRNA-seq) data with 57 annotated cell types were downloaded from Synapse (accession syn21041850) (Travaglini et al., 2020). For validation study, human nasopharyngeal gene expression matrix was obtained from GEO (accession number GSE152075) (Lieberman et al., 2020). Using patient age and viral load information, the samples were further classified into young (<40 years) and old age (>60 years) that were either negative or positive with low or high viral loads.

Data Analysis

RNA-Seq Analysis and Network Analysis

Differential expression analysis was performed using limma-voom and limma trends (Law et al., 2014; Ritchie et al., 2015). Genes with adjusted p < 0.05 were considered DE. The DE genes were tested for pathway enrichment using clusterProfiler and pathways with q-values (i.e., p-values corrected for multiple comparison) < 0.05 were considered significant (Yu et al., 2012). To perform the consensus WGCNA (Langfelder and Horvath, 2008), first the pooled control and hACE2-transduced SARS-CoV-2-infected A549 cells at low and high titer expression data were batch corrected using combat (from the sva R package) (Leek et al., 2012), and the low expressing genes with count < 5 in four out of six samples were removed. The batch-corrected normalized count was analyzed using the WGCNA R package (Langfelder and Horvath, 2008) with default parameters. A total of 50 coexpressing modules were identified, of which DE genes in correlated modules > 50 genes in size were selected for downstream analysis. For subnetwork analysis, GeneMANIA database (Mostafavi et al., 2008) was used to identify potential PPI between the DE genes from the correlated modules. The PPI networks were then overlaid with the fold-change information using Cytoscape (Shannon et al., 2003). Prior to DE analysis of the human nasopharyngeal gene expression data (Lieberman et al., 2020), the gene count matrix was subset to only include the SARS-CoV-2-positive and -negative samples that were either <40 years (young) or >60 years (old) and exclude samples that were positive with medium load viral titer. Using the sva R package and batch information, the count matrix was batch corrected across three experimental groups: control (or negative), low viral load, and high viral load in both young and old groups. DE analysis of this dataset was performed as described above. DE analysis was also performed on batch-corrected negative (control) vs. SARS-CoV-2-positive human nasopharyngeal samples.

Code Availability

The differential gene expression and pathway enrichment results for all the comparisons are included in the supplement. All the R codes used to analyze the data and make figures can be found at https://github.com/NHLBI-BCB/COVID-19_Transcriptomics. There is no restriction on data access.

Generating Pathway Enrichment Summary Map

The pathway enrichment summary map was generated using the indicated pathway enrichment results presented in Supplementary Table 4. To compare the pathway enrichment results from two different comparisons, the “Description,” “q-value,” and “GeneID” information from each pathway enrichment table were used to compute similarity of the pathways and overlap of genes between the pathways. Using this information, the enrichment Map (Merico et al., 2010) app in Cytoscape creates a graphical network, where the nodes represent a pathway. If the genes annotated to pathway is shared with another pathway (which arises due to genes usually getting annotated to more than one pathway), an edge is drawn between the pathway nodes to show this information. Additionally, if DE genes from two different comparisons were significantly annotated to the same pathway/node, that node is highlighted with two colors. The common terms in the pathway description were then used to annotate/label a group of nodes using the Enrichment Map (Merico et al., 2010) and AutoAnnotate (Kucera et al., 2016) apps in Cytoscape. R (Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/) was used for data visualization. The datasets supporting the results of this article are also available from Figshare (https://doi.org/10.6084/m9.figshare.12272351.v8) (Mehdi, 2020).

Single-Cell RNA-Seq Analysis

Analysis of human lung single-cell RNA-seq (scRNA-seq) data with 57 annotated cell types was performed in R (v3.6) using Seurat (v3.1.1) (Stuart et al., 2019). The UMI (Unique Molecular Identifier) count matrix was filtered for genes expressed in < 3 cells and normalized using SCTransform implemented in Seurat. DE genes were computed for 57 cell types using FindAllMarkers implemented in Seurat with default parameters. The UMAP plot was plotted using the top 50 principal components computed from the expression of highly variable genes selected by SCTransform.

Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR) Analysis

HCT-8 cell line (CCL-244) and the RNA of HCT-8 cells (VR-1558D) infected with OC43, a beta-coronavirus 1 strain, were purchased from ATCC. The control HCT-8 RNA was isolated after culturing the cells in RPMI medium supplemented with 10% fetal bovine serum (Gibco). The total RNA was isolated using the miRNeasy isolation Kit (Qiagen). The cDNA was synthesized using SuperScriptIII (Invitrogen). The transcripts were measured using FastStart Essential DNA Green Master (Roche). The 18s transcript measurements were used as control. All the primers were ordered from Qiagen QuantiTect assay. Four technical replicates are reported. p-values from Student's t-test are reported.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

KS and MP conceived the project. KS, YC, and MP designed the workflow and performed the analyses. KS, FS, and MP collected the data. KS, SH, KH, and MNS validated the results. KS and JJ drafted the manuscript. All authors contributed to the analysis and evaluation of the results. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the Intramural Research Programs (IRPs) of the National Heart, Lung, and Blood Institute (NHLBI) grant number: 1ZICHL006228-03 to MP.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Dr. Mark Knepper, head of the Epithelial Systems Biology Laboratory at the NHLBI, for providing his input in data analysis. This work utilized the computational resources of the NIH HPC Biowulf cluster (http://hpc.nih.gov). A previous version of this manuscript has been released as a preprint at BioRxiv (Singh et al., 2020).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.599261/full#supplementary-material

Supplementary Figure 1. (A) Volcano plot showing DE genes that were up (red color dots) and down regulated (blue color dots) in hACE2 transduced A549 cells infected with SARS-CoV-2 (high MOI). (B) Top 25 pathways from the pathway enrichment analysis of the DE genes from the mock vs. SARSCoV-2 (high MOI) comparison is presented as a horizontal bar plot, where x axis represents the log10 transformed q-value and the color of the horizontal bar is scaled blue to red representing low to high q-values, respectively. (C) Volcano plot showing DE genes that were up (red color dots) and down regulated (blue color dots) in SARS-CoV-2 (low MOI) infected A549 cells that were transduced with hACE2. (D) Top 25 pathways from the pathway enrichment analysis of the DE genes from the mock vs. SARS-CoV-2 (low MOI) comparison is presented as a horizontal bar plot, where x axis represents the –log10 transformed q-value and the color of the horizontal bar is scaled blue to red representing low to high q-values, respectively. DE, differentially expressed; MOI, multiplicity of infection. (E) Plot showing correlation between marker genes from different lung subpopulations (on x-axis) and hACE2 transduced A549 and Calu3 cells lines (color coded independent samples with legend at the bottom of the plot).

Supplementary Figure 2. (A) Correlation plot between mean gene expression from SARS-CoV-2 infected hACE2 transduced A549 (x axis) and Calu3 (y axis) cells. (B) Venn diagram showing overlap between DE genes from mock vs. SARS-CoV-2 A549 and Calu3 cell comparisons. (C) Pathway enrichment summary map for mock vs. SARS-CoV-2 comparisons in Calu3 (blue nodes) and hACE2 transduced A549 (red nodes) cells. Each node represents a pathway/biological process (BP). The node size is proportional to the number of DE genes overlapping with the BP. The nodes that share genes are connected with edges. The black circle outlines group the gene ontology (GO) terms of similar BPs. Single color nodes are pathways that are distinctly enriched by DE genes from one comparison. Two colored nodes are pathways enriched by DE genes from both comparisons. The DE genes from both comparisons enriched in inflammation, ROS, mitochondria, and autophagy processes.

Supplementary Figure 3. (A) Volcano plot showing DE genes that were up (red color dots) and down regulated (blue color dots) in IAV infected A549 cells. (B) Top 25 pathways from the pathway enrichment analysis of the DE genes from the mock vs. IAV comparison is presented as a horizontal bar plot, where x axis represents the –log10 transformed q-value and the color of the horizontal bar is scaled blue to red representing low to high q-values, respectively. (C) Venn diagram showing DE genes overlap between mock vs. SARS-CoV-2 (High MOI) and mock vs. IAV comparisons. DE, differentially expressed; MOI, multiplicity of infection.

Supplementary Figure 4. (A) Volcano plot showing DE genes that were up (red color dots) and down regulated (blue color dots) in SARS-CoV-2 infected Calu3 cells. (B) Top 25 pathways from the pathway enrichment analysis of the DE genes from the mock vs. SARS-CoV-2 comparison in Calu3 is presented as a horizontal bar plot, where x axis represents the –log10 transformed q-value and the color of the horizontal bar is scaled blue to red representing low to high q-values, respectively.

Supplementary Figure 5. (A) Top 25 pathways from the pathway enrichment analysis of the DE genes from the positive (infected) vs. negative human nasopharyngeal samples comparison is presented as a horizontal bar plot, where x axis represents the –log10 transformed q-value and the color of the horizontal bar is scaled blue to red representing low to high q-values, respectively. (B) Venn diagram showing DE genes overlap between control vs. high viral load old age human samples and mock vs. SARS-CoV-2 infected A549 cells comparisons. (C) Pathway enrichment result of common DE genes indicates in figure (B), that were concordantly downregulated in both datasets is presented as a horizontal bar plot, where x axis represents the –log10 transformed q-value and the color of the horizontal bar is scaled blue to red representing low to high q-values, respectively. DE, differentially expressed.

Supplementary Figure 6. Heatmap of the mean expression values of the indicated genes in young and old human samples that were negative (control) or positive with either high or low viral loads of SARS-CoV-2 virus is presents. (A) Heatmap of interferon signaling genes. (B) Heatmap of mitochondrial ribosomal genes. (C) Heatmap of mitochondrial complex I genes. (D) Heatmap of lysosome acidification genes. The orange and blue color bands represent upregulated and downregulated genes, respectively.

Supplementary Figure 7. DE genes from COVID-19 lung compared to healthy lungs show robust upregulation of immunity, cytokines, and inflammatory processes. (A) Volcano plot showing DE genes that were up (red color dots) and down regulated (blue color dots) in COVID-19 lung biopsy samples compared to healthy samples. (B) Pathway enrichment summary map for healthy vs. COVID-19 lungs (technical replicates) (blue nodes). Each node represents a pathway/biological process (BP). The node size is proportional to the number of DE genes overlapping with the BP. The nodes that share genes are connected with edges. The black circle summarizes the gene ontology (GO) terms of similar BPs. The DE genes from healthy vs. COVID-19 lung comparison predominantly enriched in inflammation and immunity related processes. (C) Pathway enrichment result of DE genes upregulated in COVID-19 lung vs. healthy lung biopsy samples is presented as a horizontal bar plot, where x axis represents the –log10 transformed q-value and the color of the horizontal bar is scaled blue to red representing low to high q-values, respectively. (D) Pathway enrichment result of DE genes downregulated in COVID-19 lung vs. healthy lung biopsy samples is presented as a horizontal bar plot, where x axis represents the log10 transformed q-value and the color of the horizontal bar is scaled blue to red representing low to high q-values, respectively.

Supplementary Table 1. Consensus module name and size (gene numbers) and its overlap with significant genes.

Supplementary Table 2. List of marker genes of lung cell types from single cell RNA-seq data.

Supplementary Table 3. Gene Differential Expression (GDE) Analysis results in COVID-19 transcriptomics comparisons.

Supplementary Table 4. Pathway Enrichment analysis results of GDE genes in COVID-19 transcriptomics comparisons.

References

Anand, S. K., and Tikoo, S. K. (2013). Viruses as modulators of mitochondrial functions. Adv. Virol. 2013:738794. doi: 10.1155/2013/738794

CrossRef Full Text | Google Scholar

Angajala, A., Lim, S., Phillips, J. B., Kim, J. H., Yates, C., You, Z., et al. (2018). Diverse roles of mitochondria in immune responses: novel insights into immuno-metabolism. Front. Immunol. 9:1605. doi: 10.3389/fimmu.2018.01605

PubMed Abstract | CrossRef Full Text | Google Scholar

Astuti, I., and Ysrafil (2020). Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2): an overview of viral structure and host response. Diabetes Metab. Syndr. 14, 407–412. doi: 10.1016/j.dsx.2020.04.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Banerjee, S., Biehl, A., Gadina, M., Hasni, S., and Schwartz, D. M. (2017). JAK-STAT signaling as a target for inflammatory and autoimmune diseases: current and future prospects. Drugs 77, 521–546. doi: 10.1007/s40265-017-0701-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Baracca, A., Chiaradonna, F., Sgarbi, G., Solaini, G., Alberghina, L., and Lenaz, G. (2010). Mitochondrial complex I decrease is responsible for bioenergetic dysfunction in K-ras transformed cells. Biochim. Biophys. Acta 1797, 314–323. doi: 10.1016/j.bbabio.2009.11.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Barbier, V., Lang, D., Valois, S., Rothman, A. L., and Medin, C. L. (2017). Dengue virus induces mitochondrial elongation through impairment of Drp1-triggered mitochondrial fission. Virology 500, 149–160. doi: 10.1016/j.virol.2016.10.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Barr, S. D., Smiley, J. R., and Bushman, F. D. (2008). The interferon response inhibits HIV particle production by induction of TRIM22. PLoS Pathog. 4:e1000007. doi: 10.1371/journal.ppat.1000007

PubMed Abstract | CrossRef Full Text | Google Scholar

Bikdeli, B., Madhavan, M. V., Jimenez, D., Chuich, T., Dreyfus, I., Driggin, E., et al. (2020). COVID-19 and thrombotic or thromboembolic disease: implications for prevention, antithrombotic therapy, and follow-up. J. Am. Coll. Cardiol. 75, 2950–2973. doi: 10.1016/j.jacc.2020.04.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Blanco-Melo, D., Nilsson-Payant, B. E., Liu, W. C., Uhl, S., Hoagland, D., Moller, R., et al. (2020). Imbalanced host response to SARS-CoV-2 drives development of COVID-19. Cell 181, 1036–1045.e1039. doi: 10.1016/j.cell.2020.04.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Bonora, E., Porcelli, A. M., Gasparre, G., Biondi, A., Ghelli, A., Carelli, V., et al. (2006). Defective oxidative phosphorylation in thyroid oncocytic carcinoma is associated with pathogenic mitochondrial DNA mutations affecting complexes I and III. Cancer Res. 66, 6087–6096. doi: 10.1158/0008-5472.CAN-06-0171

PubMed Abstract | CrossRef Full Text | Google Scholar

Brocke-Heidrich, K., Kretzschmar, A. K., Pfeifer, G., Henze, C., Loffler, D., Koczan, D., et al. (2004). Interleukin-6-dependent gene expression profiles in multiple myeloma INA-6 cells reveal a Bcl-2 family-independent survival pathway closely associated with Stat3 activation. Blood 103, 242–251. doi: 10.1182/blood-2003-04-1048

PubMed Abstract | CrossRef Full Text | Google Scholar

Cascella, M., Rajnik, M., Cuomo, A., Dulebohn, S. C., and Di Napoli, R. (2020). “Features, Evaluation and Treatment Coronavirus (COVID-19),” in StatPearls. Treasure Island, FL: StatPearls Publishing LLC.

PubMed Abstract | Google Scholar

Chan, J. F., Yuan, S., Kok, K. H., To, K. K., Chu, H., Yang, J., et al. (2020). A familial cluster of pneumonia associated with the 2019 novel coronavirus indicating person-to-person transmission: a study of a family cluster. Lancet 395, 514–523. doi: 10.1016/S0140-6736(20)30154-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Channappanavar, R., and Perlman, S. (2017). Pathogenic human coronavirus infections: causes and consequences of cytokine storm and immunopathology. Semin. Immunopathol. 39, 529–539. doi: 10.1007/s00281-017-0629-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, G., Wu, D., Guo, W., Cao, Y., Huang, D., Wang, H., et al. (2020). Clinical and immunological features of severe and moderate coronavirus disease 2019. J. Clin. Invest. 130, 2620–2629. doi: 10.1172/JCI137244

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, N., Zhou, M., Dong, X., Qu, J., Gong, F., Han, Y., et al. (2020). Epidemiological and clinical characteristics of 99 cases of 2019 novel coronavirus pneumonia in Wuhan, China: a descriptive study. Lancet 395, 507–513. doi: 10.1016/S0140-6736(20)30211-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X., Zhao, B., Qu, Y., Chen, Y., Xiong, J., Feng, Y., et al. (2020). Detectable serum SARS-CoV-2 viral load (RNAaemia) is closely associated with drastically elevated interleukin 6 (IL-6) level in critically ill COVID-19 patients. medRxiv. doi: 10.1101/2020.02.29.20029520

CrossRef Full Text | Google Scholar

Chu, H., Chan, J. F. W., Yuen, T. T. T., Shuai, H., Yuan, S., Wang, Y., et al. (2020). Comparative tropism, replication kinetics, and cell damage profiling of SARS-CoV-2 and SARS-CoV with implications for clinical manifestations, transmissibility, and laboratory studies of COVID-19: an observational study. Lancet Microbe. 1, e14–e23. doi: 10.1016/S2666-5247(20)30004-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Daly, J. L., Simonetti, B., Antón-Plágaro, C., Kavanagh Williamson, M., Shoemark, D. K., Simón-Gracia, L., et al. (2020). Neuropilin-1 is a host factor for SARS-CoV-2 infection. bioRxiv 370, 861–865. doi: 10.1126/science.abd3072

PubMed Abstract | CrossRef Full Text | Google Scholar

Das, R., and Chakrabarti, O. (2020). Mitochondrial hyperfusion: a friend or a foe. Biochem. Soc. Trans. 48, 631–644. doi: 10.1042/BST20190987

PubMed Abstract | CrossRef Full Text | Google Scholar

De Silva, S., and Wu, L. (2011). TRIM5 acts as more than a retroviral restriction factor. Viruses 3, 1204–1209. doi: 10.3390/v3071204

PubMed Abstract | CrossRef Full Text | Google Scholar

De Wit, E., Van Doremalen, N., Falzarano, D., and Munster, V. J. (2016). SARS and MERS: recent insights into emerging coronaviruses. Nat. Rev. Microbiol. 14, 523–534. doi: 10.1038/nrmicro.2016.81

PubMed Abstract | CrossRef Full Text | Google Scholar

Franz, M., Rodriguez, H., Lopes, C., Zuberi, K., Montojo, J., Bader, G. D., et al. (2018). GeneMANIA update 2018. Nucleic Acids Res. 46, W60–W64. doi: 10.1093/nar/gky311

CrossRef Full Text | Google Scholar

Fukai, T., and Ushio-Fukai, M. (2011). Superoxide dismutases: role in redox signaling, vascular function, and diseases. Antioxid Redox Signal. 15, 1583–1606. doi: 10.1089/ars.2011.3999

PubMed Abstract | CrossRef Full Text | Google Scholar

Fung, S. Y., Yuen, K. S., Ye, Z. W., Chan, C. P., and Jin, D. Y. (2020). A tug-of-war between severe acute respiratory syndrome coronavirus 2 and host antiviral defence: lessons from other pathogenic viruses. Emerg. Microbes Infect. 9, 558–570. doi: 10.1080/22221751.2020.1736644

PubMed Abstract | CrossRef Full Text | Google Scholar

Garcia, L. F. (2020). Immune response, inflammation, and the clinical spectrum of COVID-19. Front. Immunol. 11:1441. doi: 10.3389/fimmu.2020.01441

PubMed Abstract | CrossRef Full Text | Google Scholar

Gassen, N. C., Niemeyer, D., Muth, D., Corman, V. M., Martinelli, S., Gassen, A., et al. (2019). SKP2 attenuates autophagy through Beclin1-ubiquitination and its inhibition reduces MERS-Coronavirus infection. Nat. Commun. 10:5770. doi: 10.1038/s41467-019-13659-4

PubMed Abstract | CrossRef Full Text

Gassen, N. C., Papies, J., Bajaj, T., Dethloff, F., Emanuel, J., Weckmann, K., et al. (2020). Analysis of SARS-CoV-2-controlled autophagy reveals spermidine, MK-2206, and niclosamide as putative antiviral therapeutics. bioRxiv [Preprint]. doi: 10.1101/2020.04.15.997254

CrossRef Full Text | Google Scholar

Ghosh, S., Dellibovi-Ragheb, T. A., Kerviel, A., Pak, E., Qiu, Q., Fisher, M., et al. (2020). β-coronaviruses use lysosomes for egress instead of the biosynthetic secretory pathway. Cell 183, 1520–1535.e1514. doi: 10.1016/j.cell.2020.10.039

PubMed Abstract | CrossRef Full Text | Google Scholar

Gorshkov, K., Chen, C. Z., Bostwick, R., Rasmussen, L., Xu, M., Pradhan, M., et al. (2020). The SARS-CoV-2 cytopathic effect is blocked with autophagy modulators. bioRxiv [Preprint]. doi: 10.1101/2020.05.16.091520

PubMed Abstract | CrossRef Full Text | Google Scholar

Graef, M., and Nunnari, J. (2011). Mitochondria regulate autophagy by conserved signalling pathways. EMBO J. 30, 2101–2114. doi: 10.1038/emboj.2011.104

PubMed Abstract | CrossRef Full Text | Google Scholar

Guan, W. J., Ni, Z. Y., Hu, Y., Liang, W. H., Ou, C. Q., He, J. X., et al. (2020). Clinical Characteristics of Coronavirus Disease 2019 in China. N. Engl. J. Med. doi: 10.1056/NEJMoa2002032

CrossRef Full Text | Google Scholar

Hoffmann, M., Kleine-Weber, H., Schroeder, S., Kruger, N., Herrler, T., Erichsen, S., et al. (2020). SARS-CoV-2 cell entry depends on ACE2 and TMPRSS2 and is blocked by a clinically proven protease inhibitor. Cell 181, 271–280.e278. doi: 10.1016/j.cell.2020.02.052

PubMed Abstract | CrossRef Full Text | Google Scholar

Holownia, A., Wielgat, P., Rysiak, E., and Braszko, J. J. (2016). Intracellular and extracellular cytokines in A549 cells and THP1 cells exposed to cigarette smoke. Adv. Exp. Med. Biol. 910, 39–45. doi: 10.1007/5584_2016_214

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, C., Wang, Y., Li, X., Ren, L., Zhao, J., Hu, Y., et al. (2020). Clinical features of patients infected with 2019 novel coronavirus in Wuhan, China. Lancet 395, 497–506. doi: 10.1016/S0140-6736(20)30183-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Iwata-Yoshikawa, N., Okamura, T., Shimizu, Y., Hasegawa, H., Takeda, M., and Nagata, N. (2019). TMPRSS2 contributes to virus spread and immunopathology in the airways of murine models after coronavirus infection. J. Virol. 93, e01815–18. doi: 10.1128/JVI.01815-18

PubMed Abstract | CrossRef Full Text | Google Scholar

Jang, Y. J., Kim, J. H., and Byun, S. (2019). Modulation of autophagy for controlling immunity. Cells 8:138. doi: 10.3390/cells8020138

CrossRef Full Text | Google Scholar

Jiang, Y., Zhao, G., Song, N., Li, P., Chen, Y., Guo, Y., et al. (2018). Blockade of the C5a-C5aR axis alleviates lung damage in hDPP4-transgenic mice infected with MERS-CoV. Emerg. Microbes Infect. 7:77. doi: 10.1038/s41426-018-0063-8

PubMed Abstract | CrossRef Full Text

Kapuy, O., Vinod, P. K., and Banhegyi, G. (2014). mTOR inhibition increases cell viability via autophagy induction during endoplasmic reticulum stress - An experimental and modeling study. FEBS Open Biol. 4, 704–713. doi: 10.1016/j.fob.2014.07.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Kucera, M., Isserlin, R., Arkhangorodsky, A., and Bader, G. D. (2016). AutoAnnotate: a cytoscape app for summarizing networks with semantic annotations. F1000 Res. 5:1717. doi: 10.12688/f1000research.9090.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Landry, M. L., and Foxman, E. F. (2018). Antiviral response in the nasopharynx identifies patients with respiratory virus infection. J. Infect. Dis. 217, 897–905. doi: 10.1093/infdis/jix648

CrossRef Full Text | Google Scholar

Langfelder, P., and Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9:559. doi: 10.1186/1471-2105-9-559

CrossRef Full Text | Google Scholar

Law, C. W., Chen, Y., Shi, W., and Smyth, G. K. (2014). voom: Precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biol. 15:R29. doi: 10.1186/gb-2014-15-2-r29

CrossRef Full Text

Lee, H. C., Lee, E. S., Uddin, M. B., Kim, T. H., Kim, J. H., Chathuranga, K., et al. (2019). Released tryptophanyl-tRNA synthetase stimulates innate immune responses against viral infection. J. Virol. 93, e01291–18. doi: 10.1128/JVI.01291-18

PubMed Abstract | CrossRef Full Text | Google Scholar

Leek, J. T., Johnson, W. E., Parker, H. S., Jaffe, A. E., and Storey, J. D. (2012). The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics 28, 882–883. doi: 10.1093/bioinformatics/bts034

PubMed Abstract | CrossRef Full Text | Google Scholar

Leonard, W. J., and O'shea, J. J. (1998). Jaks and STATs: biological implications. Annu. Rev. Immunol. 16, 293–322. doi: 10.1146/annurev.immunol.16.1.293

CrossRef Full Text | Google Scholar

Li, W., Sui, J., Huang, I. C., Kuhn, J. H., Radoshitzky, S. R., Marasco, W. A., et al. (2007). The S proteins of human coronavirus NL63 and severe acute respiratory syndrome coronavirus bind overlapping regions of ACE2. Virology 367, 367–374. doi: 10.1016/j.virol.2007.04.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Lian, Q., and Sun, B. (2017). Interferons command Trim22 to fight against viruses. Cell. Mol. Immunol. 14, 794–796. doi: 10.1038/cmi.2017.76

PubMed Abstract | CrossRef Full Text | Google Scholar

Lieber, M., Smith, B., Szakal, A., Nelson-Rees, W., and Todaro, G. (1976). A continuous tumor-cell line from a human lung carcinoma with properties of type II alveolar epithelial cells. Int. J. Cancer 17, 62–70. doi: 10.1002/ijc.2910170110

PubMed Abstract | CrossRef Full Text | Google Scholar

Lieberman, N. A. P., Peddu, V., Xie, H., Shrestha, L., Huang, M. L., Mears, M. C., et al. (2020). In vivo antiviral host response to SARS-CoV-2 by viral load, sex, and age. bioRxiv [Preprint]. doi: 10.1101/2020.06.22.165225

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, H. Y., Lai, R. H., Lin, S. T., Lin, R. C., Wang, M. J., Lin, C. C., et al. (2013). Suppressor of cytokine signaling 6 (SOCS6) promotes mitochondrial fission via regulating DRP1 translocation. Cell Death Differ. 20, 139–153. doi: 10.1038/cdd.2012.106

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y., Yan, L. M., Wan, L., Xiang, T. X., Le, A., Liu, J. M., et al. (2020). Viral dynamics in mild and severe cases of COVID-19. Lancet Infect. Dis. 20, 656–657. doi: 10.1016/S1473-3099(20)30232-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Ludwig, S., and Zarbock, A. (2020). Coronaviruses and SARS-CoV-2: a brief overview. Anesth. Analg. 131, 93–96. doi: 10.1213/ANE.0000000000004845

CrossRef Full Text | Google Scholar

Mar, K. B., Rinkenberger, N. R., Boys, I. N., Eitson, J. L., Mcdougal, M. B., Richardson, R. B., et al. (2018). LY6E mediates an evolutionarily conserved enhancement of virus infection by targeting a late entry step. Nat. Commun. 9:3603. doi: 10.1038/s41467-018-06000-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Marchand, B., Arsenault, D., Raymond-Fleury, A., Boisvert, F. M., and Boucher, M. J. (2015). Glycogen synthase kinase-3 (GSK3) inhibition induces prosurvival autophagic signals in human pancreatic cancer cells. J. Biol. Chem. 290, 5592–5605. doi: 10.1074/jbc.M114.616714

CrossRef Full Text | Google Scholar

Martens, K., Hellings, P. W., and Steelant, B. (2018). Calu-3 epithelial cells exhibit different immune and epithelial barrier responses from freshly isolated primary nasal epithelial cells in vitro. Clin. Transl. Allergy 8:40. doi: 10.1186/s13601-018-0225-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Mehdi, P. (2020). Supplemental material for Singh et al. 2020. doi: 10.6084/m9.figshare.12272351.v8

CrossRef Full Text

Mehta, P., Mcauley, D. F., Brown, M., Sanchez, E., Tattersall, R. S., Manson, J. J., et al. (2020). COVID-19: consider cytokine storm syndromes and immunosuppression. Lancet 395, 1033–1034. doi: 10.1016/S0140-6736(20)30628-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Merico, D., Isserlin, R., Stueker, O., Emili, A., and Bader, G. D. (2010). Enrichment map: a network-based method for gene-set enrichment visualization and interpretation. PLoS ONE 5:e13984. doi: 10.1371/journal.pone.0013984

PubMed Abstract | CrossRef Full Text | Google Scholar

Metz, P., Chiramel, A., Chatel-Chaix, L., Alvisi, G., Bankhead, P., Mora-Rodriguez, R., et al. (2015). Dengue virus inhibition of autophagic flux and dependency of viral replication on proteasomal degradation of the autophagy receptor p62. J. Virol. 89, 8026–8041. doi: 10.1128/JVI.00787-15

PubMed Abstract | CrossRef Full Text | Google Scholar

Mohamud, Y., Shi, J., Qu, J., Poon, T., Xue, Y. C., Deng, H., et al. (2018). Enteroviral infection inhibits autophagic flux via disruption of the SNARE complex to enhance viral replication. Cell Rep. 22, 3292–3303. doi: 10.1016/j.celrep.2018.02.090

PubMed Abstract | CrossRef Full Text | Google Scholar

Mokranjac, D., and Neupert, W. (2010). The many faces of the mitochondrial TIM23 complex. Biochim. Biophys. Acta 1797, 1045–1054. doi: 10.1016/j.bbabio.2010.01.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Morita, M., Gravel, S. P., Chenard, V., Sikstrom, K., Zheng, L., Alain, T., et al. (2013). mTORC1 controls mitochondrial activity and biogenesis through 4E-BP-dependent translational regulation. Cell Metab. 18, 698–711. doi: 10.1016/j.cmet.2013.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Morita, M., Gravel, S. P., Hulea, L., Larsson, O., Pollak, M., St.-Pierre, J., et al. (2015). mTOR coordinates protein synthesis, mitochondrial activity and proliferation. Cell Cycle 14, 473–480. doi: 10.4161/15384101.2014.991572

PubMed Abstract | CrossRef Full Text | Google Scholar

Morita, M., Prudent, J., Basu, K., Goyon, V., Katsumura, S., Hulea, L., et al. (2017). mTOR controls mitochondrial dynamics and cell survival via MTFP1. Mol. Cell 67, 922–935.e925. doi: 10.1016/j.molcel.2017.08.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Morrisey, E. E. (2018). Basal cells in lung development and repair. Dev. Cell 44, 653–654. doi: 10.1016/j.devcel.2018.03.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Mostafavi, S., Ray, D., Warde-Farley, D., Grouios, C., and Morris, Q. (2008). GeneMANIA: a real-time multiple association network integration algorithm for predicting gene function. Genome Biol. 9(Suppl. 1):S4. doi: 10.1186/gb-2008-9-s1-s4

PubMed Abstract | CrossRef Full Text

Murayama, M. A., Kakuta, S., Inoue, A., Umeda, N., Yonezawa, T., Maruhashi, T., et al. (2015). CTRP6 is an endogenous complement regulator that can effectively treat induced arthritis. Nat. Commun. 6:8483. doi: 10.1038/ncomms9483

PubMed Abstract | CrossRef Full Text | Google Scholar

O'brown, Z. K., Van Nostrand, E. L., Higgins, J. P., and Kim, S. K. (2015). The inflammatory transcription factors NFkappaB, STAT1 and STAT3 drive age-associated transcriptional changes in the human kidney. PLoS Genet. 11:e1005734. doi: 10.1371/journal.pgen.1005734

PubMed Abstract | CrossRef Full Text | Google Scholar

Oikonomopoulou, K., Ricklin, D., Ward, P. A., and Lambris, J. D. (2012). Interactions between coagulation and complement–their role in inflammation. Semin. Immunopathol. 34, 151–165. doi: 10.1007/s00281-011-0280-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Prentice, E., Jerome, W. G., Yoshimori, T., Mizushima, N., and Denison, M. R. (2004). Coronavirus replication complex formation utilizes components of cellular autophagy. J. Biol. Chem. 279, 10136–10141. doi: 10.1074/jbc.M306124200

PubMed Abstract | CrossRef Full Text | Google Scholar

Rabanal-Ruiz, Y., and Korolchuk, V. I. (2018). mTORC1 and nutrient homeostasis: the central role of the lysosome. Int. J. Mol. Sci. 19:818. doi: 10.3390/ijms19030818

PubMed Abstract | CrossRef Full Text | Google Scholar

Rambold, A. S., and Lippincott-Schwartz, J. (2011). Mechanisms of mitochondria and autophagy crosstalk. Cell Cycle 10, 4032–4038. doi: 10.4161/cc.10.23.18384

CrossRef Full Text | Google Scholar

Raoult, D., Zumla, A., Locatelli, F., Ippolito, G., and Kroemer, G. (2020). Coronavirus infections: Epidemiological, clinical and immunological features and hypotheses. Cell Stress 4, 66–75. doi: 10.15698/cst2020.04.216

PubMed Abstract | CrossRef Full Text | Google Scholar

Ren, F., Zhang, L., Zhang, X., Shi, H., Wen, T., Bai, L., et al. (2016). Inhibition of glycogen synthase kinase 3beta promotes autophagy to protect mice from acute liver failure mediated by peroxisome proliferator-activated receptor alpha. Cell Death Dis. 7:e2151. doi: 10.1038/cddis.2016.56

CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Rockx, B., Kuiken, T., Herfst, S., Bestebroer, T., Lamers, M. M., De Meulder, D., et al. (2020). Comparative pathogenesis Of COVID-19, MERS And SARS in a non-human primate model. bioRxiv [Preprint]. doi: 10.1101/2020.03.17.995639

CrossRef Full Text | Google Scholar

Sanchez-Caballero, L., Guerrero-Castillo, S., and Nijtmans, L. (2016). Unraveling the complexity of mitochondrial complex I assembly: a dynamic process. Biochim. Biophys. Acta 1857, 980–990. doi: 10.1016/j.bbabio.2016.03.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Schiller, H. B., Montoro, D. T., Simon, L. M., Rawlins, E. L., Meyer, K. B., Strunz, M., et al. (2019). The human lung cell atlas: a high-resolution reference map of the human lung in health and disease. Am. J. Respir. Cell Mol. Biol. 61, 31–41. doi: 10.1165/rcmb.2018-0416TR

PubMed Abstract | CrossRef Full Text | Google Scholar

Shang, J., Ye, G., Shi, K., Wan, Y., Luo, C., Aihara, H., et al. (2020). Structural basis of receptor recognition by SARS-CoV-2. Nature 581, 221–224. doi: 10.1038/s41586-020-2179-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, B. Q., Finkbeiner, W. E., Wine, J. J., Mrsny, R. J., and Widdicombe, J. H. (1994). Calu-3: a human airway epithelial cell line that shows cAMP-dependent Cl- secretion. Am. J. Physiol. 266, L493–501. doi: 10.1152/ajplung.1994.266.5.L493

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, C. S., Qi, H. Y., Boularan, C., Huang, N. N., Abu-Asab, M., Shelhamer, J. H., et al. (2014). SARS-coronavirus open reading frame-9b suppresses innate immunity by targeting mitochondria and the MAVS/TRAF3/TRAF6 signalosome. J. Immunol. 193, 3080–3089. doi: 10.4049/jimmunol.1303196

PubMed Abstract | CrossRef Full Text | Google Scholar

Shu, Q., Lennemann, N. J., Sarkar, S. N., Sadovsky, Y., and Coyne, C. B. (2015). ADAP2 is an interferon stimulated gene that restricts RNA virus entry. PLoS Pathog. 11:e1005150. doi: 10.1371/journal.ppat.1005150

PubMed Abstract | CrossRef Full Text | Google Scholar

Simonnet, H., Demont, J., Pfeiffer, K., Guenaneche, L., Bouvier, R., Brandt, U., et al. (2003). Mitochondrial complex I is deficient in renal oncocytomas. Carcinogenesis 24, 1461–1466. doi: 10.1093/carcin/bgg109

PubMed Abstract | CrossRef Full Text | Google Scholar

Singh, K., Chen, Y. C., Judy, J. T., Seifuddin, F., Tunc, I., and Pirooznia, M. (2020). Network analysis and transcriptome profiling identify autophagic and mitochondrial dysfunctions in SARS-CoV-2 infection. bioRxiv [Preprint]. doi: 10.1101/2020.05.13.092536

PubMed Abstract | CrossRef Full Text | Google Scholar

Skoglund, C., Wettero, J., Tengvall, P., and Bengtsson, T. (2010). C1q induces a rapid up-regulation of P-selectin and modulates collagen- and collagen-related peptide-triggered activation in human platelets. Immunobiology 215, 987–995. doi: 10.1016/j.imbio.2009.11.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Somers, E. C., Eschenauer, G. A., Troost, J. P., Golob, J. L., Gandhi, T. N., Wang, L., et al. (2020). Tocilizumab for treatment of mechanically ventilated patients with COVID-19. medRxiv. doi: 10.1101/2020.05.29.20117358

PubMed Abstract | CrossRef Full Text | Google Scholar

Srirangan, S., and Choy, E. H. (2010). The role of interleukin 6 in the pathophysiology of rheumatoid arthritis. Ther. Adv. Musculoskelet. Dis. 2, 247–256. doi: 10.1177/1759720X10378372

PubMed Abstract | CrossRef Full Text | Google Scholar

Stuart, T., Butler, A., Hoffman, P., Hafemeister, C., Papalexi, E., Mauck, W. M. III., Hao, Y., et al. (2019). Comprehensive integration of single-cell data. Cell 177, 1888–1902.e1821. doi: 10.1016/j.cell.2019.05.031

CrossRef Full Text | Google Scholar

Tammaro, A., Adebanjo, G. a,.R, Parisella, F. R., Pezzuto, A., and Rello, J. (2020). Cutaneous manifestations in COVID-19: the experiences of Barcelona and Rome. J. Eur. Acad. Dermatol. Venereol. 34, e210–e211. doi: 10.1111/jdv.16530

PubMed Abstract | CrossRef Full Text | Google Scholar

To, K. F., Tong, J. H., Chan, P. K., Au, F. W., Chim, S. S., Chan, K. C., et al. (2004). Tissue and cellular tropism of the coronavirus associated with severe acute respiratory syndrome: an in-situ hybridization study of fatal cases. J. Pathol. 202, 157–163. doi: 10.1002/path.1510

PubMed Abstract | CrossRef Full Text | Google Scholar

Tondera, D., Czauderna, F., Paulick, K., Schwarzer, R., Kaufmann, J., and Santel, A. (2005). The mitochondrial protein MTP18 contributes to mitochondrial fission in mammalian cells. J. Cell Sci. 118, 3049–3059. doi: 10.1242/jcs.02415

PubMed Abstract | CrossRef Full Text | Google Scholar

Travaglini, K. J., Nabhan, A. N., Penland, L., Sinha, R., Gillich, A., Sit, R. V., et al. (2020). A molecular cell atlas of the human lung from single-cell RNA sequencing. Nature 587, 619–625. doi: 10.1101/742320

PubMed Abstract | CrossRef Full Text | Google Scholar

Walls, A. C., Park, Y. J., Tortorici, M. A., Wall, A., Mcguire, A. T., and Veesler, D. (2020). Structure, function, and antigenicity of the SARS-CoV-2 spike glycoprotein. Cell 181, 281–292.e286. doi: 10.1016/j.cell.2020.02.058

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, B., Cai, Z., Tao, K., Zeng, W., Lu, F., Yang, R., et al. (2016). Essential control of mitochondrial morphology and function by chaperone-mediated autophagy through degradation of PARK7. Autophagy 12, 1215–1228. doi: 10.1080/15548627.2016.1179401

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, C., Horby, P. W., Hayden, F. G., and Gao, G. F. (2020). A novel coronavirus outbreak of global health concern. Lancet 395, 470–473. doi: 10.1016/S0140-6736(20)30185-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, D., Hu, B., Hu, C., Zhu, F., Liu, X., Zhang, J., et al. (2020). Clinical characteristics of 138 hospitalized patients with 2019 novel coronavirus-infected pneumonia in Wuhan, China. JAMA. 323, 1061–1069. doi: 10.1001/jama.2020.1585

PubMed Abstract | CrossRef Full Text | Google Scholar

Weikel, K. A., Cacicedo, J. M., Ruderman, N. B., and Ido, Y. (2016). Knockdown of GSK3beta increases basal autophagy and AMPK signalling in nutrient-laden human aortic endothelial cells. Biosci. Rep. 36:e00382. doi: 10.1042/BSR20160174

PubMed Abstract | CrossRef Full Text | Google Scholar

Won, J. H., Park, S., Hong, S., Son, S., and Yu, J. W. (2015). Rotenone-induced impairment of mitochondrial electron transport chain confers a selective priming signal for NLRP3 inflammasome activation. J. Biol. Chem. 290, 27425–27437. doi: 10.1074/jbc.M115.667063

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, A., Peng, Y., Huang, B., Ding, X., Wang, X., Niu, P., et al. (2020). Genome composition and divergence of the novel coronavirus (2019-nCoV) originating in China. Cell Host Microbe 27, 325–328. doi: 10.1016/j.chom.2020.02.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, R., Zhang, Y., Li, Y., Xia, L., Guo, Y., and Zhou, Q. (2020). Structural basis for the recognition of SARS-CoV-2 by full-length human ACE2. Science 367, 1444–1448. doi: 10.1126/science.abb2762

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, N., and Shen, H. M. (2020). Targeting the endocytic pathway and autophagy process as a novel therapeutic strategy in COVID-19. Int. J. Biol. Sci. 16, 1724–1731. doi: 10.7150/ijbs.45498

PubMed Abstract | CrossRef Full Text | Google Scholar

Yim, W. W., and Mizushima, N. (2020). Lysosome biology in autophagy. Cell Discov. 6:6. doi: 10.1038/s41421-020-0141-7

CrossRef Full Text | Google Scholar

Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284–287. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuen, K. S., Ye, Z. W., Fung, S. Y., Chan, C. P., and Jin, D. Y. (2020). SARS-CoV-2 and COVID-19: The most important research questions. Cell Biosci. 10:40. doi: 10.1186/s13578-020-00404-4

PubMed Abstract | CrossRef Full Text

Zang, R., Gomez Castro, M. F., Mccune, B. T., Zeng, Q., Rothlauf, P. W., Sonnek, N. M., et al. (2020). TMPRSS2 and TMPRSS4 promote SARS-CoV-2 infection of human small intestinal enterocytes. Sci Immunol. 5:eabc3582. doi: 10.1126/sciimmunol.abc3582

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X., Chu, H., Wen, L., Shuai, H., Yang, D., Wang, Y., et al. (2020). Competing endogenous RNA network profiling reveals novel host dependency factors required for MERS-CoV propagation. Emerg. Microbes Infect. 9, 733–746. doi: 10.1080/22221751.2020.1738277

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, N., Zhang, D., Wang, W., Li, X., Yang, B., Song, J., et al. (2020). A novel coronavirus from patients with pneumonia in China, 2019. N. Engl. J. Med. 382, 727–733. doi: 10.1056/NEJMoa2001017

PubMed Abstract | CrossRef Full Text | Google Scholar

Ziegler, C. G. K., Allon, S. J., Nyquist, S. K., Mbano, I. M., Miao, V. N., Tzouanas, C. N., et al. (2020). SARS-CoV-2 receptor ACE2 is an interferon-stimulated gene in human airway epithelial cells and is detected in specific cell subsets across tissues. Cell. 181, 1016–1035.e19. doi: 10.1016/j.cell.2020.04.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: COVID-19, SARS-CoV-2, transcriptomics, co-expression network analysis, inflammation

Citation: Singh K, Chen Y-C, Hassanzadeh S, Han K, Judy JT, Seifuddin F, Tunc I, Sack MN and Pirooznia M (2021) Network Analysis and Transcriptome Profiling Identify Autophagic and Mitochondrial Dysfunctions in SARS-CoV-2 Infection. Front. Genet. 12:599261. doi: 10.3389/fgene.2021.599261

Received: 26 August 2020; Accepted: 15 February 2021;
Published: 16 March 2021.

Edited by:

Pier Paolo Piccaluga, University of Bologna, Italy

Reviewed by:

Biju Issac, Leidos Biomedical Research, Inc., United States
Yuhai Chen, Chinese Academy of Sciences, China

Copyright © 2021 Singh, Chen, Hassanzadeh, Han, Judy, Seifuddin, Tunc, Sack and Pirooznia. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Mehdi Pirooznia, mehdi.pirooznia@nih.gov

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.