Introduction

Infection of severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) is the cause of human coronavirus disease 2019 (COVID-19). The recent pandemic has caused devastation due to rapid spread of this viral infection. As a respiratory illness, the disease is readily transmitted. It also has a long incubation and can be carried asymptomatically, thus spreading through communities1. The COVID-19 pandemic has affected almost every country, regardless of their medical infrastructure and economic status. It has caused a healthcare crisis and created a devastating economic burden, including high unemployment, which exacerbates the effect of the disease2. At present, more than 70 million people from 213 countries have been infected with the virus (https://www.worldometers.info/coronavirus/). The rapid spreading of this respiratory infection has forced millions to shelter in their homes and has led to death of more than 1,580,000 individuals. Additionally, in the U.S., COVID-19 has disproportionately affected patients, particularly minorities and those with chronic problems such as hypertension, lung disease, diabetes, and immunocompromised conditions.

In the last three decades, the world has witnessed zoonotic transmission of various viruses (from animals to humans) leading to severe respiratory complications. These include H1N1, avian influenza, severe acute respiratory syndrome (SARS), and Middle East respiratory syndrome coronavirus (MERS‐CoV)3,4. Although infections of these viruses are often fatal, their effect is generally restricted to geographic locations such as Africa, Asia, and South America. As the world awaits a vaccine for SARS-CoV-2, efforts are being made to understand the molecular mechanisms of these infections5,6,7,8,9.

High-throughput technologies such as RNA sequencing and microarrays are useful in the detection of respiratory virus infections and in understanding their molecular effect on human lung epithelial cells10. Extensive data on corona virus sequences have been deposited in public repositories such as NCBI Gene Expression Omnibus (GEO) and EMBL Array Express11,12. Meta-analysis and mining of such data can aid in a) understanding the molecular impact of COVID-19, b) elucidating differences and similarities between SARS-CoV-2 and other respiratory virus infections, and c) identifying targets for drug development. In the current study, we performed comparative analysis of publicly available gene expression data related to human lung epithelial cells infected with a respiratory virus. The analyses identified genes specifically expressed by SARS-CoV-2 infections and those that are commonly altered due to infection of coronovirus-2 and/or other respiratory viruses. In particular, expression of CSF2 (colony-stimulating factor 2) appears to be involved in COVID-19 disease. Several COVID-19 clinical trials are currently focusing on inhibition of CSF2 [GM-CSF]13,14,15.

Results

RNA sequencing identified genes altered on SARS-CoV-2 infection of normal human lung epithelial cells

Differential expression analysis was conducted separately using RNA-seq data from GSE147507 and GSE153970. In total, 164 genes were up-regulated, and 76 genes were down-regulated on SARS-CoV-2 infection of normal human bronchial epithelial cells [GSE147507] (Fig. 1a). Similarly, 405 genes were up-regulated and 544 genes were down-regulated in SARS-CoV-2-infected human airway epithelium [GSE153970] (Fig. 1b). Gene list comparison identified 27 commonly down-regulated and 73 commonly up-regulated genes (Fig. 1c). Gene Ontology enrichment analysis of common differentially expressed genes showed “Inflammatory response”, “Neutrophil chemotaxis”, “Immune response”, “Cell chemotaxis”, and “Keratinization” as the top 5 enriched biological processes (Fig. 1d).

Figure 1
figure 1

High throughput gene expression analysis of SARS-CoV-2 infection in human lung epithelial cells. (a,b) Volcano plots showing the proportion of differentially expressed genes on SARS-CoV-2 infection in human lung epithelial cells [from GSE147507 and GSE153970, respectively]. (c) Venn diagram depicting genes commonly up-/down-regulated as per both studies. (d) Top 10 biological processes associated with common differentially expressed genes on SARS-CoV-2 infection of human lung epithelial cells.

Full size image

The impact of other viral infections on the human lung epithelial cell transcriptome was explored using publicly available microarray data

A search of microarray datasets in the NCBI GEO database led to identification of five studies involving human lung epithelial cells subjected to viral infection. GEO2R analyses of each study were performed separately to identify differentially expressed genes. Figure 2a shows differentially expressed probes identified from each comparative analysis. From GSE71766, there were 574 differentially expressed probes (281 genes) in H1N1-infected BEAS-2B cells compared to control; 166 probes (72 genes) were up-regulated by RV16 infection. The combined infection of RV16 and H1N1 altered expression of 589 probes (288 genes). In the case of GSE49840, there were 19,959 probes (14,892 genes) differentially expressed in H7N7-infected Calu-3 cells compared to mock infected cells; 20,155 (14,948 genes), 14,859 (11,362 genes), and 17,496 probes (13,012 genes) were differentially expressed by H5N1, H3N2, or H7N9 infections, respectively. For GSE17400, DOHV or SARS-CoV infection of Calu-3 cells led to altered expression of 447 (344 genes) and 221 probes (182 genes), respectively. Analysis of GSE48575 led to identification of 130 (99 genes) and 9 (7 genes) differentially expressed probes on seasonal (sH1N1) or pandemic (H1N1pdm) influenza virus infection, respectively, of NHBECs. Lastly, the processing of GSE47962 resulted in discovery of 13,414 (7938 genes), 1415 (997 genes), 7 (5 genes), and 8 (6 genes) differentially expressed probes after H1N1, SARS-CoV-dORF6, SARS-CoV-BAT, or SARS-CoV viral infections, respectively, of HAE cells.

Figure 2
figure 2

Comparative transcriptome analysis of collected gene expression data. (a) Histogram showing the differentially expressed genes as determined by GEO2R analysis of respiratory viral infection microarray data. (b) Heatmap showing 7 genes exclusively altered on SARS-CoV-2 infection.

Full size image

Comparative analysis of DEGs resulted in identification of SARS-CoV-2 infection-specific genes and those commonly affected by most of the viral infections

In order to identify genes that are exclusively affected on SARS-CoV-2 infection, we compared common DEGs (98 protein-coding genes of 100) from RNA-seq studies with GEO2R analysis results of GSE47962, GSE17400, GSE48575, GSE49840, and GSE71766 (Supplementary Table S1). GEO2R results from each microarray study include lists of differentially expressed probes, gene symbols, fold change, and adjusted p-values. If multiple probes related to the same gene were differentially expressed, we considered the probe with highest fold change value. Comparative analysis showed 7 genes exclusively altered on SARS-CoV-2 infection, including 5 up-regulated genes (CSF2, S100A8, MRGPRX3, S100A9 and MAB21L4) and two down-regulated genes (CXCL14 and PCDH7) (Fig. 2b). However, microarray platforms showed the absence of probes related to MAB21L4 in all microarray studies (Table 1).

Table 1 List of probes associated with COVID-19 exclusive genes in microarray platforms related to selected gene expression studies.
Full size table

Inflammation-related genes (such as IL6, IL1A, IL1B, CXCL2, CXCL6, CCL20, TNIP1, VNN1, TNFAIP3 and NFKBIZ) were commonly affected due to infection of SARS-CoV-2 and other SARS or avian/human influenza viruses (Supplementary Figure S1c).

Validation of SARS-CoV-2 exclusive genes in human bronchial organoid RNA-seq data

Processing and analysis of RNA-seq data related to SARS-CoV-2-infected human bronchial organoids (HBO) led to identification of 1532 differential expressed genes (861 up-regulated and 671 down-regulated genes) (Supplementary Figure S1a). With the exception of MAB21L4, all SARS-CoV-2 exclusive genes from comparative transcriptome analysis showed the same expression pattern in HBO cells (Supplementary Figure S1b).

Protein–protein interaction analysis of altered genes on SARS-CoV-2 infection revealed hub genes

The common 98 protein-coding, differentially expressed genes on SARS-CoV-2 infection were first queried in the STRING database to identify known/predicted interactions among them. The database returned a PPI network of 333 interactions (edges) between 72 genes (nodes) (Fig. 3). The PPI network was downloaded as a simple interaction format (SIF) file, visualized with Cytoscape, and analyzed with Cytohubba plugin to identify hub genes. The top 50 genes were obtained based on three network parameters: degree, closeness, and betweenness separately. The 43 genes featured in all three lists were considered as hub genes (Table 2). The hub genes were IL6, IL1B, CXCL8, MMP9, CXCL1, CSF2, CCL20, ICAM1, IL1A, CXCL2, CSF3, SAA1, NFKB2, PI3, TNFAIP3, CXCL3, CXCL5, CXCL6, EDN1, HBEGF, BCL2A1, NFKBIZ, S100A12, PLAUR, BIRC3, IL36G, LIF, SERPINB2, SPRR1A, SPRR1B, IVL, DDX58, SPRR2A, ZC3H12A, S100A9, PDGFB, IL7R, MAP3K8, IKBKE, MAF, ADAM8 and GBP5. Among hub genes, most were affected by infection of human or avian influenza viruses such as H7N7, H1N1, H7N9, H3N2, and H5N1 (Fig. 4); CSF2 and S100A9 were exclusive to SARS-CoV-2 infection.

Figure 3
figure 3

Protein–protein interaction network of genes altered by SARS-CoV-2 infection. The network with 333 edges among 72 nodes was obtained from the STRING database. Interactions with scores of ≥ 0.4 were considered.

Full size image
Table 2 Hub genes involved in SARS-CoV-2 infection of normal human bronchial epithelial cells.
Full size table
Figure 4
figure 4

Association of COVID 19 hub genes with other respiratory viral infections. A network generated by use of Cytoscape depicts the expression of COVID-19 hub genes on infection by various respiratory viruses. The green arrows indicate down-regulation, and the red arrows show up-regulation of these genes.

Full size image

Discussion

Microarray meta-analysis and comparative transcriptome analysis have been useful bioinformatic approaches for maximum utilization of publicly available gene expression data16,17. Since the advent of high-throughput technologies such as microarrays and RNA-seq, researchers have performed in-depth transcriptome analyses of various biological conditions, leading to various discoveries18. Data from thousands of such experiments are being deposited in public repositories such as GEO and ArrayExpress. Selection and combinatorial analysis of such data can aid researchers in understanding the molecular mechanisms of a disease and in discovering biomarkers19,20,21. We observed the availability of various microarray studies related to human respiratory viral infection and sensed the opportunity to compare the effect of SARS-CoV-2 and other respiratory viral infections on the human lung transcriptome.

The comparative transcriptome analysis led to identification of genes that were altered exclusively after SARS-CoV-2 infection. Among these genes, S100 calcium-binding protein A9 (S100A9) and S100 calcium-binding protein A8 (S100A8) are calcium- and zinc-binding proteins that are elevated in inflammatory lung disorders22. S100A8 and S100A9 form a heterodimer complex called Calprotectin (CLP). Elevated levels of CLP is found in bronchoalveolar lavage fluid (BALF), serum and lung tissue of pneumonia patients23. Serum CLP could be potential biomarker for COVID-19 and further research is required to test this hypothesis24. Colony stimulating factor 2 (CSF2) is a cytokine-coding gene associated with respiratory diseases such as pulmonary alveolar proteinosis25. CSF2 (GM-CSF) is known to be pro-inflammatory cytokine produced by wide variety of cells such as macrophages, T-cells, fibroblast, tumor cells, endothelial cells, and with primary production at the inflammation site26,27. GM-CSF influences activation and proliferation of immune cells such as macrophages, monocytes, dendritic cells, neutrophils, eosinophils28. Although GM-CSF plays important role in maintaining immune homeostasis, its over-expression in lung could lead to fibrotic reactions and severe immune cell infiltrations29. Recent reports show that a) COVID-19 patients requiring intensive care unit (ICU) have increased level of GM-CSF in China and b) drugs targeting CSF2 (GM-CSF) or its receptor (such as lenzilumab, namilumab, gimsilumab, and otilimab) are being evaluated in clinical trials with COVID-19 patients13,14,15,30. MAS-related GPR family member X3 (MRGPRX3), a member of the mas-related/sensory neuron specific subfamily of G protein coupled receptors, is down-regulated in human airway epithelial cells exposed to smoke from electronic cigarettes31. In oral cancers, lung cancers, and head and neck cancers, C-X-C Motif chemokine ligand 14 (CXCL14) functions as a tumor suppressor; it also induces growth of prostate and breast cancers32,33,34,35,36. Protocadherin 7 (PCDH7) is involved in cell–cell recognition and adhesion37. Mab-21-like 4 (MAB21L4) has no known association with lung disorders or respiratory virus infections.

The PPI analysis of genes differentially expressed by SARS-CoV-2 infection identified 43 hub genes. CSF2 and S100A9 were the only hub genes to show a SARS-CoV-2-exclusive gene expression pattern. Almost all other hub genes were affected by infection of other respiratory viruses. In conclusion, both PPI analysis and comparative transcriptome analysis point to a role of CSF2 in the molecular mechanism of SARS-CoV-2 infections of human lung epithelium. The current study also highlights the exclusivity of known lung inflammation disorder genes such as S100A8 and S100A9 with respect to SARS-CoV-2 infection.

Our current study utilizes exclusively the transcriptomic data. Similar comparative analysis as well as integrative analysis with proteomics and other “omic” data will be useful as more data becomes public. The current bioinformatics approach makes use of available data to identify potential molecular targets for treatment of COVID-19. Future studies will focus on experimental validation to establish the exclusive association of CSF2, S100A8 and S100A9 with SARS-CoV-2 infections in lung epithelial cells.

Methods

RNA sequencing data analysis

The NCBI GEO database was searched for microarray or RNA-sequencing data related to SARS-CoV-2 infections of human lung epithelial cells. We evaluated two studies, GSE147507 and GSE1539706,38. GSE147507 included normal human bronchial epithelial cells subjected to mock treatment (n = 3) or to SARS-CoV-2 infection (n = 3); GSE153970 included primary human airway epithelial cultures infected with mock (n = 3) or SARS-CoV-2 (n = 3). For the purpose of validation, we also selected RNA-seq data related to SARS-CoV-2-infected human bronchial organoids [GSE150819].

Raw sequencing data related to selected samples were downloaded from Sequence Read Archive (SRA) using fastq-dump of sratoolkit v2.9.6 (http://ncbi.github.io/sra-tools/). First, raw sequencing reads were trimmed to remove adapter sequences and low-quality regions using Trim Galore! (v0.4.1) (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). Trimmed reads were subjected to quality control analysis using FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Tophat v2.1 was used to map trimmed raw reads to the human reference genome (hg38)39. All bam files from multiple runs related to the same samples were merged and sorted using SAMtools (Version: 1.3.1)40. Finally, raw read counts were enumerated for each gene in each sample using GTF (gene transfer file) from Ensembl [Homo.sapiens_GRCh38.82.gtf] and HTSeq-count41.

Analysis of differential expression was performed using DESeq2 according to a standard protocol (https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html)42. Genes with adj.P-value < 0.05 and absolute fold change ≥ 1.5 were considered as significantly differentially expressed. Common up-regulated and down-regulated genes from GSE147507 and GSE153970 were obtained using the Venny online tool (https://bioinfogp.cnb.csic.es/tools/venny/).

Gene ontology enrichment analyses of the common Differentially Expressed Genes (DEGs) were accomplished by use of the Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8 online tool43. Gene ontology (GO) biological processes with P-values < 0.05 and gene counts > 2 were considered as significantly enriched.

Microarray data collection and analysis

The NCBI GEO database was queried for microarray data related to SARS-CoV infections of human lung epithelial cells. A query (SARS-CoV) AND "Homo sapiens"[porgn] AND ("gse"[Filter] AND ("Expression profiling by array"[Filter])) led to 15 search results. After screening, two studies (GSE47962, GSE17400) were selected. To find microarray data related to SARS-CoV infections of human lung epithelial cells, the GEO database was queried using (((Human lung epithelium) OR (Human bronchial epithelial) AND "Homo sapiens"[porgn] AND ("gse"[Filter] AND "Expression profiling by array"[Filter]))) AND (viral infection AND ("gse"[Filter] AND "Expression profiling by array"[Filter])) AND ("gse"[Filter] AND "Expression profiling by array"[Filter]) AND ("Expression profiling by array"[Filter]). This led to 38 search results, three of which (GSE49840, GSE71766, and GSE48575) were selected for analysis. Table 3 provides sample, platform, and cell line details for all five studies.

Table 3 Summary of publicly available, high-throughput experimental data considered in the current study.
Full size table

Since SARS-CoV-2 RNA-seq data included transcriptome profiling after 24 h of infection, in all microarray studies, we considered only samples after 24 h of viral infection.

GSE47962 included samples from human airway epithelium (HAE) cells infected with SARS-CoV, influenza virus (H1N1), or variants of SARS-CoV (SARS-dORF6 and SARS-BatSRBD)44. GSE71766 comprised human bronchial epithelial cells (BEAS-2B) infected with rhino virus (RV), influenza virus (H1N1), or both (RV + H1N1)45. Bronchial epithelial cell line 2B4 (a clonal derivative of Calu-3 cells) infected with SARS-CoV or Dhori virus (DOHV) was part of GSE1740046. GSE49840 included polarized calu-3 (cultured human airway epithelial cells) infected with human influenza virus (H3N2) or avian influenza viruses (H7N9, H5N1, and H7N7)47. GSE48575 consisted of normal human bronchial epithelial cells (NHBEC) infected with seasonal H1N1 influenza A (sH1N1) or pandemic H1N1 influenza A (H1N1pdm)48. GEO2R was used to identify differentially expressed genes for each of these studies independently12. Probes with adj. P-value < 0.05 and absolute fold change ≥ 1.5 were considered as statistically significant and were compared with DEGs of SARS-CoV-2 infection from RNA-seq data.

Protein–protein interaction analysis

STRING, a database of known or predicted protein–protein interactions (PPIs) was used to obtain interactions between genes altered on SARS-CoV-2 infection49. Output from the STRING database was uploaded to Cytoscape v3.7.2 in simple interaction format, and the Cytohubba app was employed to identify hub genes50,51,52. The top 50 genes were obtained separately from the PPI network based on three network parameters (closeness, degree, and betweenness) then common genes among these were selected as hub genes. We also checked the DrugBank database to determine if a drug is available to target them53.