Gene Expression Profiling Reveals the Shared and Distinct Transcriptional Signatures in Human Lung Epithelial Cells Infected With SARS-CoV-2, MERS-CoV, or SARS-CoV: Potential Implications in Cardiovascular Complications of COVID-19

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is the causative virus for the current global pandemic known as coronavirus disease 2019 (COVID-19). SARS-CoV-2 belongs to the family of single-stranded RNA viruses known as coronaviruses, including the MERS-CoV and SARS-CoV that cause Middle East respiratory syndrome (MERS) and severe acute respiratory syndrome (SARS), respectively. These coronaviruses are associated in the way that they cause mild to severe upper respiratory tract illness. This study has used an unbiased analysis of publicly available gene expression datasets from Gene Expression Omnibus to understand the shared and unique transcriptional signatures of human lung epithelial cells infected with SARS-CoV-2 relative to MERS-CoV or SARS-CoV. A major goal was to discover unique cellular responses to SARS-CoV-2 among these three coronaviruses. Analyzing differentially expressed genes (DEGs) shared by the three datasets led to a set of 17 genes, suggesting the lower expression of genes related to acute inflammatory response (TNF, IL32, IL1A, CXCL1, and CXCL3) in SARS-CoV-2. This subdued transcriptional response to SARS-CoV-2 may cause prolonged viral replication, leading to severe lung damage. Downstream analysis of unique DEGs of SARS-CoV-2 infection revealed changes in genes related to apoptosis (NRP1, FOXO1, TP53INP1, CSF2, and NLRP1), coagulation (F3, PROS1, ITGB3, and TFPI2), and vascular function (VAV3, TYMP, TCF4, and NR2F2), which may contribute to more systemic cardiovascular complications of COVID-19 than MERS and SARS. The study has uncovered a novel set of transcriptomic signatures unique to SARS-CoV-2 infection and shared by three coronaviruses, which may guide the initial efforts in the development of prognostic or therapeutic tools for COVID-19.

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is the causative virus for the current global pandemic known as coronavirus disease 2019 (COVID- 19). SARS-CoV-2 belongs to the family of single-stranded RNA viruses known as coronaviruses, including the MERS-CoV and SARS-CoV that cause Middle East respiratory syndrome (MERS) and severe acute respiratory syndrome (SARS), respectively. These coronaviruses are associated in the way that they cause mild to severe upper respiratory tract illness. This study has used an unbiased analysis of publicly available gene expression datasets from Gene Expression Omnibus to understand the shared and unique transcriptional signatures of human lung epithelial cells infected with SARS-CoV-2 relative to MERS-CoV or SARS-CoV. A major goal was to discover unique cellular responses to SARS-CoV-2 among these three coronaviruses. Analyzing differentially expressed genes (DEGs) shared by the three datasets led to a set of 17 genes, suggesting the lower expression of genes related to acute inflammatory response (TNF, IL32, IL1A, CXCL1, and CXCL3) in SARS-CoV-2. This subdued transcriptional response to SARS-CoV-2 may cause prolonged viral replication, leading to severe lung damage. Downstream analysis of unique DEGs of SARS-CoV-2 infection revealed changes in genes related to apoptosis (NRP1, FOXO1, TP53INP1, CSF2, and NLRP1), coagulation (F3, PROS1, ITGB3, and TFPI2), and vascular function (VAV3, TYMP, TCF4, and NR2F2), which may contribute to more systemic cardiovascular complications of COVID-19 than MERS and SARS. The study has uncovered a novel

INTRODUCTION
The novel severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) belongs to the Coronaviridae family of viruses (coronaviruses) and is responsible for the coronavirus disease 2019 (COVID- 19) pandemic (1). Along with its other accomplices, Middle East respiratory syndrome coronavirus (MERS-CoV) and severe acute respiratory syndrome coronavirus (SARS-CoV), SARS-CoV-2 can jump species barrier followed by human-to-human transmission via droplet infection. In late December 2019, initial reports suggested the origin of SARS-CoV-2 in a seafood and wild animal trading market in Wuhan, China (2). To date, the pandemic has caused more than 83 million infections and more than 1.8 million deaths worldwide (https://www.worldometers.info/coronavirus/). SARS-CoV-2 leads to more cardiovascular complications than do MERS-CoV and SARS-CoV; however, what causes these major differences remains obscure (3,4).
The initial genome identification of SARS-CoV-2 suggested that it has a ∼80% similarity with SARS-CoV and 96% identical to a bat coronavirus; however, there are differences in its pathogenicity and host response (2). The virus nucleic acid shedding patterns in both symptomatic and asymptomatic patients of SARS-CoV-2 are similar, which explains the transmission potential of otherwise asymptomatic carriers (5). In contrast, the viral burden in the upper respiratory tract in SARS-CoV infection peaks at around 10 days after the initial exposure (6). On the contrary to SARS-CoV-2, viral load in MERS-CoVinfected individuals peak at week 2 of the onset of infection (7). This suggests the difference in the virulence and host response of these three strains.
Upon entry, next steps are viral replication, amplification, and spread in the host, which largely depend on similarities and/or uniquenesses in transcriptional signature of these viruses. Patients with SARS-CoV-2 manifest a few different but aggravated symptoms, particularly major cardiovascular complications, from SARS-CoV and MERS-CoV, which may be attributable to the difference in their transcriptional signatures. Reports of aggravated blood coagulation in COVID-19 patients suggested the mechanism of prominent elevation of D-dimer and fibrin/fibrinogen degradation products (8). Higher mortality rate is reported in COVID-19 patients with thromboembolic events (9), and treatment with anticoagulant-heparin has improved prognosis (10).
The present comparative analysis has determined key differences in transcriptional changes in lung epithelial cells induced by these virus strains. To further examine transcriptional responses of SARS-CoV-2 and other two coronaviruses, we analyzed a comprehensive map of lung epithelial cells infected with these three coronaviruses and explored pathological host responses unique to SARS-CoV-2. Our findings may help to understand potential mechanisms by which SARS-CoV-2 causes more cardiovascular complications than do two other coronaviruses and to establish molecular bases for the development of therapies against COVID-19. Figure 1A depicts the workflow of gene expression analysis. For differential gene expression analysis of SARS-CoV-2 infection, raw expression counts were downloaded from Gene Expression Omnibus (GEO) accession number GSE147507 (11). RNA sequencing (RNAseq) dataset was generated on Illumina Nextseq 500 platform. The raw read counts were normalized by log2 transformation, before and after normalization box plot; principal component analysis (PCA) and density plot are shown in Supplementary Figure 1. Using INMEX tool that employs DESeq (12,13), differential expression analysis was performed and differentially expressed genes (DEGs) were characterized for each sample with adjusted p < 0.05 [false discovery rate (FDR) corrected by Benjamini-Hochberg method]. GSE81909, the dataset we used for analysis of MERS-CoV, was generated on Agilent-Whole Human Genome Microarray 4x44K G4112. After downloading the raw read counts from GEO, we normalized the dataset using variance-stabilizing normalization followed by quantile normalization (14). Before and after normalization box plot, PCA and density plot are shown in Supplementary Figure 2. Similarly, we downloaded raw read counts of GSE17400 (15) for analysis of SARS-CoV infection. This dataset was generated on Affymetrix Human Genome U133 plus 2.0 Array. After normalization of dataset using variancestabilizing normalization followed by quantile normalization (Supplementary Figure 3), both microarray datasets, GSE81909 and GSE17400, were subjected to DEGs analysis using LIMMA algorithm (16). DEGs were characterized for each sample with adjusted p < 0.05 (FDR corrected by Benjamini-Hochberg method). Heatmap visualization of a subset of 25 overexpressed and underexpressed genes was constructed using heatmap.2 from the gplot package in R. Volcano plots were constructed using custom scripts in R, and PCA was performed on log2 fold-change values using PMA package in R (17) (Supplementary Figure 4). It is worthwhile mentioning that all three datasets used in this study are collected from different laboratories and using different cell lines, as well as experimental techniques (e.g., microarrays, RNAseq). We selected RNAseq dataset for SARS-CoV-2 as there were especially no microarray datasets on SARS-CoV-2 in humans. Therefore, we did take the present analysis strategy of comparing each dataset with its internal control to call the DEGs for each coronavirus infection. Further, we compared the three sets of DEGs to find the shared and unique genes between coronavirus infections; we applied this strategy to minimize data variabilities (e.g., operator and platform biases). GSE147507 was generated in primary human lung epithelium (NHBE); GSE81909 was generated in human airway epithelial cells, whereas GSE17400 was generated in human bronchial epithelial cells. Supplementary Table 1 provides detailed information of each dataset and sequencing/microarray platform used.

Functional Gene Set Enrichment Analysis of DEGs
To discern the implication of DEGs called from transcriptome analysis of coronavirus infection in lung epithelial cells, we performed a functional analysis using the EnrichR platform (18). This web-based software product evaluates significantly enriched pathways/terms in an input gene list with the help of its extensive gene set libraries, which includes Gene Ontology (GO) (19) and various pathway analysis libraries such as Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, Reactome pathway, wikipathway, Panther, and Biocarta. We retrieved tables of enriched pathways from each database and prepared a comprehensive table of most significant pathways for each coronavirus infection based on the adjusted p value (ranking derived from Fisher exact test for gene sets) significance.

Common DEGs Analysis Between Coronaviruses
We created a coronavirus-gene network for better visualization of the shared genes between the coronaviruses using Cytoscape software (20). The network was generated by utilizing the list of DEGs from three coronaviruses studies in which coronaviruses are the source nodes; genes are the target nodes, and the connections between them are the edges in the network. The network core represents the coronaviruses, whereas the innercircle genes in the network are the shared ones, and outercircle genes are unique to each coronavirus ( Figure 2C). A Venn diagram representing the shared and unique DEG portion between three coronaviruses was generated using VENNY 2.1 tool (https://bioinfogp.cnb.csic.es/tools/venny/index.html). A heatmap represents the expression profiles for common DEGs between coronaviruses. Clustering of selected genes on the heatmap was performed by hierarchical clustering algorithm utilizing Euclidean distance measure.

Pathway Clustering and Network-Based Hub Gene Analysis
For visualization and interpretation of the biological relevance of unique DEGs to SARS-CoV-2 DEGs, Cytoscape v3.1 plugin was used for analysis. Biological pathway clustering analysis was done using BinGO (21). BinGO analyzes GO terms and functional groups association within the biological networks. We performed biological pathway clustering analysis to see collective function of these genes. The size of a node is proportional to the number of targets in the biological process category. The color represents enrichment significance-the deeper the color on a color scale, the higher the enrichment significance. Hub gene network analysis was performed using NetworkAnalyst (13), which created a protein-protein interaction (PPI) network by integrating the InnateDB interactome with the original seed of 221 DEGs. This tool supports integrative analysis of gene expression data through statistical, visual, and networkbased analysis approaches by taking the advantage of common functions for network topology and module analysis approaches. Briefly, the complete list of unique DEGs from SARS-CoV-2 was uploaded into the web-based server of NetworkAnalyst. Network construction was restricted to contain all the original seed proteins in order to visualize the connections. To help identify highly interconnected hub nodes, topological measures (e.g., degree and betweenness centrality) were used. Expression of the genes was considered as the network feature, where redcolored nodes are genes with increased expression, green nodes are genes with decreased expression, and gray nodes are genes not expressed in our data.

Expression2Kinases Analysis of Regulatory Gene Networks
ChEA is a comprehensive databases of kinases and transcription factors (22), and it is used in background of Expression2Kinases (X2K) (23), the tool we used to understand the upstream regulatory molecules of DEGs in SARS-CoV-2 infection. The 10 most significant transcription factors and kinases were extracted based on Fisher exact test p value enrichment scoring. We downloaded the ".graphml" file generated from the analysis to create and visualize regulatory network on Cytoscape environment. This ensures that the protein network obtained during network expansion is properly connected by automatically increasing the path length, so that there are more intermediate proteins used to connect the transcription factors. In the network, a yellow node represents intermediate proteins in the PPI regulatory network. Node size represents the significance of protein based on adjusted p value; the bigger the nodes size, the higher the significance value.

Statistical Analyses
For differential expression analysis of SARS-CoV-2 dataset GSE147507, read counts were subjected to differential expression analysis using INMEX, which utilizes the Rpackage DESeq (13). Genes with adjusted p < 0.05 were considered significant. The p value adjustment for multiple comparisons was done by the Benjamini-Hochberg method. For MERS-CoV-GSE81909 and SARS-CoV-GSE17400, differential expression analysis was performed with LIMMA algorithm for each dataset, independently using adjusted p < 0.05, based on the FDR using the Benjamini-Hochberg method and moderated t test. Significantly enriched GO terms were identified using hypergeometric tests, and p ≤ 0.05 was applied as a cutoff for statistical significance.

Selection of Eligible Gene Expression Datasets for Coronavirus Infection in Human Lung Epithelial Cells
We selected three studies from the GEO accession numbers: GSE147507 for SARS-CoV-2, GSE81909 for MERS-CoV, and GSE17400 for SARS-CoV. The search was limited to transcriptome data generated in human lung epithelial cells. Figure 1A depicts the overall workflow of the analysis in this study. A total of 3/3, 20/20, and 9/9 control/infected cell culture replicates SARS-CoV-2, MERS-CoV, and SARS-CoV, respectively, were used in this analysis. GSE147507 dataset was RNAseq data and generated on Illumina Nextseq 500, and we utilized only six samples (GSM4432378, GSM4432379, GSM4432380, GSM4432381, GSM4432382, and GSM4432383); these were independent biological triplicates of primary human lung epithelium (NHBE), which were mock treated or infected with SARS-CoV-2 (USA-WA1/2020). Of note, the other two datasets were generated by microarray using Affymetrix Human Genome U133A series (GSE17400-SARS-CoV) and Agilent-014850 Whole Human Genome Microarray 4x44K G4112F (GSE81909-MERS-CoV). Sample sources of all three datasets were of human lung epithelial cells and primary lung cells infected with coronaviruses. Supplementary Table 1 provides detailed information of each dataset and sequencing/microarray platform used.

Analysis of Differentially Expressed Genes (DEGs) in the SARS-CoV-2 Dataset Led to Perturbation of Inflammatory, Coagulation, and Apoptotic Pathways
In SARS-CoV-2-infected dataset (GSE147507), we identified a total of 338 DEGs with adjusted p < 0.05 (Supplementary Dataset 1). Among these 338 DEGs, 92 genes increased, and 246 decreased. Figure 1B depicts the heatmap of expression of top significant DEGs among the samples. Table 1 lists the top 20 increased and decreased DEGs from our analysis of SARS-CoV-2 infection. Interferon (IFN)induced transmembrane protein 10 (IFITM10), C-X-C motif chemokine ligand 14 (CXCL14), and myosin light chain kinase (MYLK) were among the most significantly increased genes, while small proline-rich protein 2D (SPRR2D), interleukin 36 gamma (IL36G), and serum amyloid A2 (SAA2) were the most decreased genes in our analysis of SARS-CoV-2-infected lung epithelial cells compared to mock controls. When these DEGs were subjected to the analysis of overrepresented biological pathways and enriched terms, several pathways related to inflammation, apoptosis, blood coagulation, and lung fibrosis were enriched ( Table 2). Enriched terms and biological pathways were significantly overrepresented in the gene list if they showed an adjusted p < 0.05. DEGs from SARS-CoV-2 infection were associated with the KEGG pathways such as IL-17 signaling pathway (hsa04657) with database overlap of 21/93 (which means of 93 genes associated with this pathway reported in KEGG, 21 are present among our DEGs) and adjusted p =  Table 2, SSX family member 2 (SSX2), fos protooncogene, AP-1 transcription factor subunit (FOS), and early growth response 1 (EGR1) were among the most significantly increased genes, whereas transcription elongation regulator 1 (TCERG1), G protein-coupled receptor 55 (GPR55), and casein kappa (CSN3) were the most decreased genes in our analysis of MERS-CoV-infected lung epithelial cells compared to controls. Similarly, Supplementary Table 3 shows that pentraxin 3 (PTX3), early growth response 1 (EGR1), and EGR2 were among the most significantly increased genes, whereas aldo-keto reductase family 1 member B10 (AKR1B10), IFN-α-inducible protein 6 (IFI6), and matrix metallopeptidase 1 (MMP1) are the most decreased genes in our analysis of SARS-CoV-infected lung epithelial cells compared to mock controls. When these DEGs from both MERS-CoV and SARS-CoV were subjected to the analysis of overrepresented biological pathways and enriched terms, several pathways related to inflammation and apoptosis were commonly enriched (Supplementary Tables 4, 5).

Muted Expression of Acute Inflammatory Genes Was Observed in the SARS-CoV-2 When Compared to MERS-CoV-2 and SARS-CoV
Previous studies revealed that lung epithelial cells, dendritic cells, and macrophages all express cytokines to some extent during major viral infection causing cytokine storm. However, little is known about the situation in COVID-19. Earlier studies showed IFN-γ-related cytokine storm in SARS-CoV infection, whereas MERS-CoV infection had delayed induction of proinflammatory cytokines and suppression of innate antiviral response. It is crucial to identify the primary source of the cytokine storm in response to SARS-CoV-2 infection and the underlying virological mechanisms. Our analysis of shared DEGs between three coronaviruses resulted in 17 shared DEGs (Figure 2A), among which most genes are related to acute inflammation. Figure 2B indicates that the expression levels of these genes were lower in SARS-CoV-2 when compared with the other two coronaviruses. Our identification of a muted transcriptional response to SARS-CoV-2 supports a model in which initial failure to rapidly respond to infection results in prolonged viral replication and subsequent recruitment of proinflammatory cells as the infection progresses to induce alveolar damage in COVID-19. Table 3 depicts the expression value of DEGs shared in all three coronaviruses. It is evident that the expression of critical acute inflammatory genes including TNF-α-induced protein 3 (TNFAIP3), C-X-C motif chemokine ligand 1 (CXCL1), and TNF was lower in the SARS-CoV-2 dataset compared to two other coronaviruses. These results may suggest that lung epithelial cells do not directly contribute to the cytokine storm during COVID-19 and that other immune cells appear to participate in this process.

SARS-CoV-2 Elicits Suppressed Type I IFN Response and Activation of Apoptotic Gene Signature
Dissecting the DEGs involved in IFN response to coronavirus infection in primary human lung epithelial cells revealed that SARS-CoV-2 elicits a muted response that lacks robust induction of a subset of cytokines including the type I IFN compared to the response to MERS-CoV and SARS-CoV ( Figure 3A). Furthermore, our analysis revealed that in desperate bid to control the viral propagation, SARS-CoV-2 infection induced several apoptosis-related genes in human lung epithelial cells compared to responses to MERS-CoV and SARS-CoV ( Figure 3B).

Downstream Analysis of DEGs Unique to SARS-CoV-2 and Network-Based Meta-Analysis Led to Pathways and Hub Genes Related to Inflammation and Vascular Dysfunction
A set of 221 DEGs unique to SARS-CoV-2 underwent downstream analysis to understand the enriched pathways and associated hub genes. Using BinGO enrichment clusters of biological process (GO terms) associated with unique DEGs of SARS-CoV-2 was generated, which revealed enriched pathway clusters associated with immune responses/chemotaxis, blood coagulation, apoptosis, vascular remodeling, and vascular cell proliferation ( Figure 4A). Supplementary Dataset 4 compiles the complete list of GO: terms enriched in DEGs associated with SARS-CoV-2 infection. We generated a PPI network by integrating the InnateDB interactome with the original seed of 221 DEGs. An expanded PPI network was generated with 2,542 nodes representing the proteins and 4,457 edges representing the interaction between these proteins. Network-based hub DEG analysis ( Figure 4B) identified ribosomal protein L9 (RPL9) and SMAD family member 3 (SMAD3) to be the most highly ranked hub genes that increased and decreased among the DEGs, respectively, based on betweenness centrality and degree score. The list of top 15 hub genes based on network topology scores is shown in Table 4.

Identification of the Transcription Factors and Regulatory Kinases Network Upstream to the Unique DEGs Obtained From SARS-CoV-2
To understand what lies to the upstream of the unique DEGs identified from the SARS-CoV-2 infection, we used X2K bioinformatics tool. The regulatory gene network analysis resulted in identification of transcription factors and kinases related to the DEGs. Network in Figure 4C shows the top kinases and transcription factors related to our DEGs. The list of top 10 ranked transcription factors and protein kinases is shown in Supplementary Table 4. This analysis revealed the most important regulatory gene candidates that may be involved in the formation of regulatory complexes. Mitogen-activated protein kinase 1 (MAPK1) and MAPK3 are among the top kinases, whereas SMAD3 and SMAD2 are among the top transcription factors associated with the unique DEGs from lung epithelial cells infected with SARS-CoV-2.

Distinct Pathways and Gene Signatures Associated With SARS-CoV-2
After we removed the DEGs shared by three coronaviruses (17 DEGs) and those shared between SARS-CoV-2 and MERS-CoV (91 DEGs); and SARS-CoV-2 and SARS-CoV (9 DEGs), a total of 221 DEGs remained that was specific to SARS-CoV-2 (Figure 2A and Supplementary Dataset 1). Our interest was to identify the distinct pathways that may participate in the pathogenesis of COVID-19. Using the list of SARS-CoV-2specific DEGs, we conducted biological process (GO) analysis on the unique set of DEGs using a Cytoscape plugin, BinGO tool. Table 5 depicts the most distinct pathways and their associated representative genes with its known function and expression fold change that may have implication on the disease pathogenesis.

DISCUSSION
In the present study, we focus on defining transcriptional responses to SARS-CoV-2 relative to MERS-CoV and SARS-CoV. A major goal was to discover unique cellular responses to SARS-CoV-2 among these three coronaviruses. We specifically selected datasets generated in cultured human lung epithelial cells infected with each of these three coronaviruses as this cell type is the major interface between the environment and the host and defends the lung against foreign substances and pathogens. In general, our data show that overall transcriptional footprints to SARS-CoV-2 infection were distinct from those to the other two coronaviruses. Despite the decreased expression of acute inflammatory and type I IFN genes in response to SARS-CoV-2, we observed increased expression of several genes associated with interleukin signaling, complement pathways, and chemokines. This finding echoes with the previously published study, which conducted RNAseq analysis to understand host transcriptional response to influenza A virus and SARS-CoV-2 in primary human bronchial epithelial cells (11). We used a publicly available subset of RNAseq data (GSE147507) from this study to compare it with independent datasets for SARS-CoV and MERS-CoV. It is worth mentioning that the list of DEGs unique to SARS-CoV-2 infection generated in our analysis was associated with coagulation and vascular function, which may explain why COVID-19 causes more systemic cardiovascular complications than do MERS and SARS (4,8).
By analyzing RNAseq dataset of lung epithelial cells infected with SARS-CoV-2, we defined transcriptional signatures of 338 DEGs, including 92 increased and 246 decreased genes across the datasets. Among the top 10 increased DEGs, IFITM10, CXCL14, and MYLK are the most significantly increased genes. CXCL14 is a cytokine involved in immunoregulatory and inflammatory  processes by mediating the chemotactic activity for monocytes and therefore can be implicated in the immune cell infiltration in the lung during SARS-CoV-2 infection (24). IFITM proteins family inhibit the entry of a large number of viruses; however, the exact role of IFITM10 as an antiviral agent remains unknown (25). SAA2 is the most significantly decreased gene in our analysis. SAA2 is a useful inflammatory marker in acute viral infections such as influenza, but its decreased expression in our analysis is consistent with the aberrant inflammatory response of SARS-CoV-2 infection (26). Viral infection is marked by the activation of immune system, which is evident from the enrichment of several pathways, including IL-17, TNF, and apoptosis signaling pathways among others (27). Patients who are infected with COVID-19 may develop pneumonia and progress to severe respiratory failure termed acute respiratory distress syndrome, which may result in the development of lung fibrosis   (28). Consistent with this report, several genes such as colonystimulating factor 3 (CSF3), endothelin 1 (EDN1), plasminogen activator, urokinase (PLAU), and MMP9 were reported to be differentially expressed in SARS-CoV-2 infection.
Previous studies revealed that lung epithelial cells, macrophages, and dendritic cells express cytokines to some extent during major viral infection causing cytokine storm. Evidence for molecular mechanisms of cytokine storm in COVID-19 remains limited. Earlier studies have shown IFNγ-related cytokine storm in SARS patients (29), while delayed induction of proinflammatory cytokines and suppression of innate antiviral response by the MERS-CoV (30). It is crucial to identify the primary source of the cytokine storm in response to SARS-CoV-2 infection and the virological mechanisms behind the cytokine storm. Consistent to this model, our analysis of shared DEGs between three coronaviruses resulted in 17 DEGs, among which most molecules are related to acute inflammation. It is evident that important acute inflammatory genes (e.g., TNFAIP3, CXCL1, and TNF) are decreased in SARS-CoV-2infected cells compared to other coronaviruses. These results may suggest that the major source of cytokine storm in COVID-19 is not lung epithelial cells, but possibly immune cell types.
These aberrant transcriptional responses to SARS-CoV-2 may indicate low responses to infection, resulting in prolonged viral replication and serious lung damage in COVID-19 (31).
Our study linked DEGs unique to SARS-CoV-2 infection with pathway clusters related to immune responses, blood coagulation, apoptosis, and vascular remodeling. Apoptosis, which is a defense mechanism of hosts against the viral infection, depends on the rapid programmed cell death to curtail viral spread. Previous studies reported that SARS-CoV has evolved sophisticated molecular strategies to trigger host cell apoptotic defenses (32,33). In our data, SARS-CoV-2 infection increased the expression of NRP1 and FOXO1. NRP1 has a regulatory role of apoptotic pathways, whereas FOXO1 is an important regulator of cell death acting downstream of several signaling pathways, including CDK1, PKB/AKT1, and STK4/MST1 (34).
Evidence suggests elevation of D-dimer and fibrin/fibrinogen degradation products in patients with COVID-19, highlighting aggravated blood coagulation (8). A recent clinical study examined seven lungs obtained during autopsy from patients with COVID-19 (35). The study observed vascular endothelial injury and widespread thrombus formation in pulmonary vessels. Immunohistochemical staining of pulmonary vasculature of COVID-19 showed alveolar capillary microthrombi were nine times as prevalent in patients with COVID-19 compared with influenza patients. In consistent with these reports, we identified coagulation-related genes in SARS-CoV-2 infection, including tissue factor that initiates the external coagulation cascades. In contrast, SARS-CoV-2 suppressed the expression of antithrombotic gene TFPI2, which inhibits a variety of serine proteases including factor VIIa/tissue factor complex.
Despite the clinical impact, the information on the mechanisms of COVID-19 and its cardiovascular complications remains limited. A systems approach, involving unbiased bioinformatics and network analysis, may help to identify causative genes and integrated pathways as drug targets for the improvement in disease management (36). Network-based analysis of hub genes in the DEGs dataset unique to SARS-CoV-2 infection resulted in prioritization of RPL9 as the most highly ranked DEG that had increased expression, based on betweenness centrality and degree score. The increased expression of RPL9, a ribosomal protein, can be attributed to the fact that virus hijacks the translational machinery of the host for its survival by the mechanisms such as ribosome shunting and phosphorylation of ribosomal proteins (37,38).
Regulatory gene network analysis helps to understand what lies upstream of the DEGs in cells infected with SARS-CoV-2. It is important to find out the regulatory kinases and transcription factors as they participate in the pathogenesis and the progression of the virus infection (39). Among several kinases regulating the expression of DEGs expression in our analysis, genes involved in MAPK cascades (MAPK1, MAPK2, MAPK8, and MAP3K7) have roles in host response to viral infection (40). SMAD3, an effector molecule in the transforming growth factor-β signaling pathway, is also an interesting candidate in our analysis as it was the top hub gene in our network analysis and also the most significant transcription factor in our regulatory gene analysis. In the present study, SARS-CoV-2 reduced SMAD3 expression, which is consistent with previous findings on the decreased expression of SMAD3 during viral infection to overtake the host innate antiviral mechanism (41,42).
In conclusion, our study provides the snapshot of transcriptional host responses to SARS-CoV-2 infection, in which expression of various inflammatory genes is decreased. These results may explain why many infected individuals show no symptoms. Furthermore, our study revealed that SARS-CoV-2 elicits muted antiviral type I IFN response, which may result in prolonged viral replication. These findings may also explain why SARS-CoV-2 infection has a longer incubation period than other coronavirus infections. Our analysis revealed expression of several genes related to apoptosis, coagulation, and vascular function, which may contribute to cardiovascular complications. Furthermore, our study has identified a novel set of candidate transcriptomic signatures unique to SARS-CoV-2 infection, which may guide the initial efforts in the development of diagnostic or therapeutic tools for COVID-19.

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.