Transcriptome profiling of different types of human respiratory tract cells infected by SARS-CoV-2 highlight an unique role for inflammatory and interferon response

ABSTRACT

The emergence of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) disease (COVID-19) at the end of 2019 has caused a large global outbreak and now become a major public health issue. However, there is currently a lack of data underlying how the human host interacts with SARS-CoV-2 virus. In the current study, We performed Venn-analysis, Gene ontology (GO), KEGG pathway analysis and Protein-protein interaction analysis of whole transcriptome studies with the aim of clarifying the genes and pathways potentially altered during human respiratory tract cell infection with SARS-CoV-2. We found 36 overlapping upregulated genes among different types of cells after viral infection. Further functional enrichment analysis revealed these Differential Expressed Genes (DEGs) are most likely involved in biological processes related to inflammatory response and response to cytokine, cell component related to extracellular space and I-kappa B/NF-kappa B complex, molecular function related to protein binding and cytokine activity. KEGG pathways analysis highlighted altered conical and casual pathways related to TNF, NF-kappa B, Cytokine-cytokine receptor interaction and IL-17 signaling pathways during SARS CoV-2 infection with CXCL1, CXCL2, CXCL3, CXCL8, CXCL10, IL32, CX3CL1, CCL20, IRF1, NFKB2 and NFKB1A up-regulated which may explain the inflammatory cytokine storms associated with severe cases of COVID-19.

KEYWORDS:

• SARS-CoV-2
• transcriptomics
• gene ontology
• pathway analysis
• inflammatory response

Introduction

The outbreak of Coronavirus Disease 2019 (COVID-19) caused by the novel coronavirus named severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2), was first reported at the end of 2019 (Tan et al. 2020; Zhu et al. 2020). The clinical characteristics of COVID-19 include respiratory symptoms, fever, cough, dyspnea, pneumonia (Chen et al. 2020; Huang et al. 2020; Li et al. 2020; Wang et al. 2020), acute respiratory distress syndrome (ARDS), and even multiple organ dysfunction in severe cases (Singhal 2020). Clinical characteristics varied in different COVID-19 patients, and a number of patients showed asymptomatic infection (Chan et al. 2020; Guan et al. 2020; Xu et al. 2020) Recent studies have revealed that asymptomatic cases could also transmit the virus to contacts around them, which makes it more difficult to control the spread of the virus (Rothe et al. 2020).

Several anti-viral drugs have been used to target the SARS-CoV-2 infection, including lopinavir/ritonavir, penciclovir, chloroquine, hydroxychloroquine, remdesivir, and ribavirin. However, most of them have limited efficacy to reduce the severity of COVID-19 patients. In addition, some also have serious side effects (Cao et al. 2020; Wang et al. 2020). The shortcomings of current tested drugs have resulted in scientists and physicians continuously seeking alternative, effective treatments for COVID-19 infection. Transcriptional profiling has provided remarkable opportunities for understanding the relationship between cellular function and metabolic pathways, as well as to define the possible implications of genetic variability and environmental conditions in many tissues and organisms after viral infection (Nowakowski et al. 2016; Li et al. 2019), which could contribute to the screening of drug targets against viruses. In bioinformatics, RNA-Seq analyses are challenging due to its complexity and voluminosity. However, many computational methods and techniques have been developed for the use of the analysis of transcriptomics data. In addition, network-based approach is also extensively used to analyze and identify the potential candidate genes after viral infection (Gratton et al. 2020).

The clinical picture of COVID-19 patients are characterized by an over-exuberant immune response with lung lymphomononuclear cells infiltration and proliferation that may account for tissue damage more than the direct effect of viral replication (Azkur et al. 2020; Blanco-Melo et al. 2020; Yang et al. 2020). Evidence also showed that SARS-CoV-2 infects human respiratory tract cells, prompting inflammation and tissue damage/loss as well as in other different cell types (Wong et al. 2020). Despite recent advances, little is known about the pathways involved in viral pathogenesis (Ziegler et al. 2020). In the present study, we performed an analysis of several transcriptome studies, aiming to clarify which genes and cellular networks were up- or down-regulated during SARS-CoV-2 infection in human lung cells. Moreover, we assessed a comprehensive GO enrichment and KEGG pathway analysis to predict how the modulation of these genes could affect the outcome of the disease.

Materials and methods

Study search and RNA-Seq data collection

We searched for RNA-Seq experiments deposited in the Gene Expression Omnibus (GEO) database and NCBI PubMed related to SARS-CoV-2 infection in different types of human respiratory tract cells: 1. Calu-3 was infected with SARS-CoV-2 (USA-WA1/2020/NR-52281) at an MOI of 2 for 24 h, sequencing libraries were sequenced on an Illumina NextSeq 500 platform. ID: doi/10.3390/biology9090260 (Shaath and Alajez 2020); 2. A549^ACE2 was infected with SARS-CoV-2 (USA-WA1/2020/ NR-52281) at an MOI of 0.2 for 24 h, sequencing libraries were sequenced on an Illumina NextSeq 500 platform. ID: doi/10.1101/2020.07.12.199687 (GSE154613) (Daisy et al. 2020); 3. A lung organoid model driving by human pluripotent stem cells was infected with SARS-CoV-2 (USA-WA1/2020/ NR-52281) at an MOI of 0.01 for 24 h, sequencing libraries were sequenced on an Illumina NovaSeq 6000 platform. ID: doi/10.1038/s41586-020-2901-9 (GSE155241) (Tatsuya et al. 2020); 4. Human bronchial organoids generated by primary human bronchial epithelial cells was infected with SARS-CoV-2 (Hu/DP/Kng/19-020) of 5.0 × 10⁴ PFU for 5 days, sequencing libraries were sequenced on an Illumina NextSeq 500 platform. ID: doi/10.1101/2020.05.25.115600v2 (GSE150819) (Han et al. 2020) which carried out in human cells (ex vivo) or human cell lines (in vitro). Genes with fold change >1.3 or <0.77 and p-values <0.05 were considered to be statistically significant.

Functional analysis of differential expressed genes

The Go terms enrichment, KEGG pathways and protein–protein interaction analysis of the DEGs were performed using OmicsBean (http://www.omicsbean.cn) and Cytoscape v. 3.1.0 software (http://www.cytoscape.org/) (Saito et al. 2012) as a complementary and more comprehensive approach for identifying the central hub genes. Gene ID and Fold change are required to upload to complete the relevant analysis, the results were visualized using the ClueGO v. 2.1.1 application in Cytoscape, which creates clusters of functionally related genes using the following GO databases: KEGG, Wiki Pathways, REACTOME, and GO.

Results

Overlapping and functional enrichment analysis of differential expression genes among different lung cells infected with SARS-CoV-2

To highlight the changes in cellular host genes in response to SARS-CoV-2 infection and further explore the mechanism of viral infection and pathogenesis, we searched the electronic databases (GEO databases or PubMed) to identify the related studies. We included 4 studies in our qualitative study, they are all about transcriptome profiling of human respiratory tract cells (Calu-3, doi:10.3390/biology9090260; ACE2-A549, GSE154613; Lung organoid model driving by human pluripotent stem cells, GSE155241; Primary human bronchial epithelial cells, GSE150819).

Using Venn diagram analysis for the groups’ data, the following was found: there were 36 overlapping up-regulated and 0 overlapping down-regulated mRNAs among the four infection groups (Figure 1(A,B)). We further performed a functional enrichment analysis using the 36 up-regulated genes. Results revealed 2,506 terms associated with biological processes (BP), 164 terms associated with cell components (CC), 226 terms associated with molecular functions (MF), and 102 terms associated with Kyoto Encyclopedia of Genes and Genomes pathways (KEGG). Of these, a number of terms were significantly enriched including 1,907 for BP, 78 for CC, 98 for MF and 39 for KEGG (Figure 2(A)). The 10 most enriched terms are presented in detail for BP, CC and MF (Figure 2(B)).

Figure 1. Comparative analysis of differentially expressed mRNAs in SARS-CoV-2 infected lung cells. A,B Venn analysis of all up/down expressed mRNAs among the four SARS-CoV-2 infection groups. a: Calu-3 infected with SARS-CoV-2 (doi:10.3390/biology9090260), b: ACE2-A549 infected with SARS-CoV-2 (GSE154613), c: Lung organoid model using human pluripotent stem cells infected with SARS-CoV-2 (GSE155241), d: Primary human bronchial epithelial cells infected with SARS-CoV-2 (GSE150819).

Figure 2. GO term enrichment analysis of SARS-CoV-2 induced up-regulated genes in all four groups. A, Number of enriched GO terms and KEGG pathways. B, Overview of significantly enriched biological processes, cell components, and molecular functions.

Biological process analysis of overlapping all up-regulated genes

The 1,907 different biological process GO terms were further analyzed using OmicsBean software. According to the different levels of the GO terms (Figure 3(A,C)), most differential genes were related to positive regulation of immune system process, regulation of immune system process and positive regulation of response to stimulus (level 3 GO), and continued subdivision of the grades showed that the differential genes were involved in the defense response, response to cytokine and cytokine-mediated signaling pathways (level 3–7 GO). According to the number of enriched genes (Figure 3(B)), they mainly enriched in signal transduction (29 differential genes), defense response (24 differential genes), regulation of cell communication (24 differential genes), response to organic substance (23 differential genes), multicellular organism development (22 differential genes) and cellular response to chemical stimulus (20 differential genes). A sorting of the different GO terms based on the P-value and gene number (Figure 3(C–E)) revealed that these up-regulated genes are most likely involved in biological processes related to inflammatory response and response to cytokine, which indicates that SARS-CoV-2 infection could cause a strong host immune response.

Figure 3. Analysis of 1907 biological process significantly enriched by all-up-regulated genes. A, Levels of enriched biological processes. B, Percent of genes associated with enriched biological processes. C, Top 10 of significantly enriched biological processes. D, Expressed genes associated with enriched biological processes. E, Pie chart of expressed proteins of enriched biological processes.

Cell component analysis of all overlapping up-regulated genes

The analysis of the 78 different cell components-related GO terms showed that according to the classification of the different GO levels (Figure 4(A)), most up-regulated genes were related to continued subdivision identified extracellular region part, membrane-enclosed lumen (level 2 GO) and organelle lumen (level 3 GO), cytoplasm, cytoplasmic part and cytosol (level 4–6 GO). According to the number of enriched genes (Figure 4(B)), the most abundantly enriched up-regulated genes were associated with cytoplasm (27 genes), intracellular organelle (22 genes), intracellular organelle part (20 genes) and vesicle (9 genes). By distinguishing the difference in the P-value of the GO terms and gene number (Figure 4(C,D)), we found that these up-regulated genes were enriched in extracellular space and I-kappa B/NF-kappa B complex of cell components.

Figure 4. Analysis of 78 cell components significantly enriched by all up-regulated genes. A, Levels distribution of enriched cell components. B, Percent of genes associated with enriched cell components. C, Top 10 of significantly enriched cell components. D, Pie chart of expressed proteins of enriched cell components.

Molecular function analysis of overlapping all upregulated genes

The analysis of the 98 different molecular functions according to the number of the different GO classification levels (Figure 5(A)), the genes showing the highest differential expression are associated with protein binding (level 3 GO), and further subdivisions showed that the up-regulated genes are associated with receptor binding, cytokine receptor binding, G-protein coupled receptor binding, chemokine receptor binding and CXCR chemokine receptor binding (level 4–7 of GOs). Further analysis of enriched genes (Figure 5(B)) revealed that the most abundant up-regulated genes were enriched in receptor binding (12 genes), identical protein binding (10 genes) and enzyme binding (6 genes). Sorting the different GO terms based on ascending P-values and gene number (Figure 5(C,D)) revealed that the genes showing the greatest differential expression are likely enriched in the molecular functions of the protein binding and cytokine activity.

Figure 5. Analysis of 98 molecular function enriched by all up-regulated genes. A, Levels distribution of enriched molecular functions. B, Percent of genes associated with enriched molecular functions. C, Top 10 of significantly enriched molecular functions. D, Pie chart of expressed proteins of enriched molecular functions.

Based on the above results, it can be hypothesized that SARS-CoV-2 infection could activate the immune response in human respiratory tract cells. A statistical analysis of immune regulation GO terms and pathways (Supplemental Tables 1–4) showed that many up-regulated genes, such as, PTGS2, LIF, NFKBIA, ICAM1, CSF1, NFKB2, NFKBIZ, CX3CL1, CCL20, CXCL3, BIRC3 and so on, are involved in.

KEGG pathway and protein–protein interaction analysis of overlapping all up-regulated genes

Through a KEGG pathway cluster analysis, we further found that the up-regulated genes mainly regulated the TNF signaling pathway, cytokine-cytokine receptor interaction, viral protein interaction with cytokine and cytokine receptor and IL-17 signaling pathway (Figure 6(A)). Sorting the different KEGG pathways based on ascending P-values and gene number (Figure 6(B,C)) revealed that the up-regulated genes mainly enriched in TNF signaling pathway, NF-kappa B signaling pathway and IL-17 signaling pathway.

Figure 6. KEGG pathway cluster analysis and correlation analysis of all up-regulated genes. A, Percent of genes associated with enriched KEGG pathway. B, Distribution of enriched KEGG pathways. C, Pie chart of significant enriched KEGG pathway. D, Correlation analysis of up-regulated genes associated with key pathways. E, SARS-CoV-2 infection regulated the IL-17 signaling pathway. The key up-regulated factors are marked in red.

A protein–protein interaction analysis (Figure 6(D)) was also explored. The results showed that the key genes regulating these processes are CXCL1, CXCL2, CXCL3, CXCL8, CXCL10, IL32, CX3CL1, CCL20, IRF1, NFKB2 and NFKB1A, NFKB1E and their certain interaction relationship was also shown. Furthermore, an analysis of the most significant IL-17 signaling pathway in KEGG (Figure 6(E)) revealed that the TNFAIP3, NFKBIA (upstream) gene in the pathway and the CXCL1, CXCL2, CXCL8. CXCL10, CCL20 and PTGS2 (downstream) genes in the pathway play an important role in immunity of host cells to SARS-CoV-2 infection. In addition, the details of NF-kappa B and TNF signaling pathway activated by SARS-CoV-2 infection were also shown in Figure S1–2. This finding is consistent with the results of the GO analysis, which showed that SARS-CoV-2 infection mainly activated the ability of respiratory tract cells to respond to the invasion of extracellular pathogens.

Discussion

As of 13 Oct 2020, a total of 37,423,660 confirmed cases of COVID-19 with 1,074,817 deaths had been reported (WHO 2020), the case numbers are still increasing rapidly due to the second wave of epidemics worldwide. The underlying biological processes and mechanisms impacted by this viral infection of the host are still not clear. 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.

RNA-Seq has contributed to increase the knowledge relative to the biological and cellular processes involved during viral infection, especially in the context of an emergent pathogen (Li et al. 2019; Gratton et al. 2020), such as SARS-CoV-2. Venn analysis of several RNA-Seq assays increases the statistical power of samples and we applied this method for re-examining past evidence to better understand the pathogenesis of infections. In our current study, we explored recent SARS-CoV-2 infected cell-derived transcriptome datasets to dissect the cellular and molecular changes in response to SARS-CoV-2 infection. We finally included four available studies and analyzed the gene expression patterns using Venn analysis (Figure 1).

The RNA-Seq results of respiratory tract cells infected by SARS-CoV-2 showed that the up-regulated genes identified from the comparison of two groups such as CXCL1, CXCL2, CXCL3, CXCL8, CXCL10, IL32, CX3CL1, CCL20, IRF1, NFKB2 and NFKB1A, NFKB1E are involved in the regulation of TNF signaling pathway, NF-kappa B signaling pathway and IL-17 signaling pathway. Tumor necrosis factor (TNF) was first described as a glycoprotein induced in response to endotoxin and has the capacity to kill tumor cells (Carswell et al. 1975). Then, TNF was shown to be a T helper type 1 (Th1) cytokine produced by several cell types, and to be involved in pro-inflammatory responses like interleukin (IL)-1β and IL-6 (Beutler and Cerami 1986). In concert with chemokines, TNF signaling pathway is involved in the regulation of inflammatory processes of infectious diseases including HIV (Human immunodeficiency virus), JEV (Japanese encephalitis virus), WNV (West Nile virus) and USUV (Ustuuvirus) infection (Cheng et al. 2004; Bradley 2008; Pasquereau et al. 2017; Obdulio et al. 2019). The nuclear factor NF-kappa B pathway has also long been considered a prototypical pro-inflammatory signaling pathway, largely based on the role of NF-kappa B signaling pathway in the expression of pro-inflammatory genes including cytokines and chemokines (Lawrence et al. 2001). Therefore, NF-kappa B has been considered as a target for new anti-inflammatory drugs for a long time. In our study, we found SARS-CoV-2 infection could activate NF-kappa B signaling pathway and up-regulate some pro-inflammatory genes, we further speculate drugs targeting the NF-kappa B signaling pathway may be one of the candidate drugs for the treatment severe COVID-19 patients.

It is worth mentioning that, IRF1 in respiratory tract cells was up-regulated by SARS-CoV-2 infection. IRF1 is a member of the IFN regulatory factor (IRF) family that plays an important role in immunomodulation. IRF1 prominently participates in antiviral defense by regulating early expression of IFNs (Panda et al. 2019). In addition, as a transcription factor, IRF1 regulates constitutive expression of other antiviral genes by binding to their promoter directly. Studies revealed that IRF1 deficient mice are more susceptible to RNA virus infections (Kimura et al. 1994; Reis et al. 1994; Maloney et al. 2012; Mboko et al. 2014; Nair et al. 2017). Thus, we think IRF1 may play a critical role in regulating constitutive antiviral gene networks to confer resistance against SARS-CoV-2 infections in human respiratory epithelial cells.

Inflammation is the first line of defense against tissue damage caused viral infection, involving activating both the innate and adaptive immune responses (D’Elia et al. 2013). However, exuberant immune responses after viral infection, which also called cytokine storm, has been found to be associated with excessive levels of pro-inflammatory cytokines and also some tissue damage (Thompson et al. 2017). Cytokine storm has been proved to exist after several influenza viruses infection (Wang et al. 2010; Zhou et al. 2013), as well as coronavirus (Zhang et al. 2004; Jiang et al. 2005; Ayman et al. 2019), and contributes to acute lung injury or ARDS. Preliminary studies have revealed that SARS-CoV-2 infection could also trigger cytokine storm in severe cases, accompanied by an increase of various cytokines (Liu et al. 2020). In the current study, we found that several cytokine-related genes were upregulated after SARS-CoV-2 infection which could provide explanation about the cytokine storms associated with severe COVID-19 cases.

Conclusions

Our results highlighted activation of TNF, NF-kappa B, cytokine-cytokine interaction, IL-17 signaling pathway and other inflammatory response as the main hallmark associated with SARS-CoV-2 infection in different types of human respiratory epithelial cells in vitro. We further identified several differentially expressed genes during the course of infection such as IRF1, CXCL1-3, IL32, CXCL10, CX3CL1, CCL20, NFKB2 and others. Upregulation of these genes which may serve as useful biomarkers of COVID-19 infection. Interferon and inflammatory response to SARS-CoV-2 infection might provide explanation to cytokine storms associated with severe COVID-19 cases. However, precise mechanisms of the host response to SARS-CoV-2 require further detailed investigation.

Acknowledgements

This work was supported by Shenzhen Natural Science Foundation (JCYJ20190809152415652), China Postdoctoral Science Foundation (2019M660836), the National Science and Technology Major Project (2018ZX10711001, 2017ZX10103011, 2018ZX09711003, and 2020YFC0841700), Shenzhen Science and Technology Research and Development Project (202002073000001). Yingxia Liu, Yu Wang, Minghui Yang and Luping Lei contributed to the study design. Minghui Yang, Luping Lei, Yang Yang and Jun Wang contributed to the manuscript and figures. Qiumei Cao, Xiao Jiang, Jinzhi Lai, Ling Qing and Kun Huang contributed to the proof reading. All the authors have read and approved the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

We searched for RNA-Seq experiments deposited in the Gene Expression Omnibus (GEO) database and NCBI PubMed related to SARS-CoV-2 infection in different types of human respiratory tract cells (Calu-3, doi:10.3390/biology9090260; ACE2-A549, GSE154613; a lung organoid model driving by human pluripotent stem cells, GSE155241; Primary human bronchial epithelial cells, GSE150819) which carried out in human cells (ex vivo) or human cell lines (in vitro).

Additional information

Funding

This work was supported by Shenzhen Natural Science Foundation [grant number JCYJ20190809152415652], China Postdoctoral Science Foundation [grant number 2019M660836], Shenzhen Key Technology Research Project [grant number JSGG20200207161928126], the National Science and Technology Major Project [grant numbers 2018ZX10711001, 2017ZX10103011, 2018ZX09711003, and 2020YFC0841700], and Shenzhen Science and Technology Research and Development Project [grant number 202002073000001].