We employed a search of the Gene Expression Omnibus (GEO) database to identify transcriptome datasets related to the seven HCoV infections in human respiratory cell lines. Of these, we selected four datasets related to SARS-CoV, SARS-CoV-2, MERS-CoV, and HCoV-229E for further analysis and listed them in Supplementary Table 1. Specifically, we focused only on the transcriptome data of the Calu3 and MRC-5 cell lines, as well as their corresponding mock controls, at the 24-h time point after infection with the four viruses in each dataset. The transcriptome data of PBMC samples from COVID-19 patients and corresponding healthy controls was obtained from GSE152418 dataset. All RNA-Seq data from our article were accessible through the NCBI GEO with the accession numbers GSE147507 (27), GSE56189, GSE148729 (28), GSE155986 (29), and GSE152418. The SARS-CoV-2 strains used in GSE147507 (USA-WA1/2020) and GSE148729 (Patient isolate, BetaCoV/Munich/BavPat1/2020| EPI_ISL_406862) were early original strains.

The FASTQ data from all these datasets was downloaded with the SRA toolkit and quality-controlled with fastp (30). The “fastq-dump” command was then employed to extract the fastq sequence information. Subsequently, we employed the “split-files” function to split the sequencing data, followed by applying a standardized quality control using the “fastp” tool with the parameter “-q 25 -u 30 thread = 5” to remove reads with substandard base quality. The human reference genome GRCh38.84 was used for alignment through HISAT2 (v2.2.0), and the mitochondrial genes were removed. Reads with alignment confidence scores exceeding 20 were extracted using “-F 4 -q 20,” resulting in the generation of the gene expression profiles required for downstream analyses. Gene expression levels of the sequencing data were quantified using the DESeq2 R package (version 1.18.1) and visualized through boxplots (Supplementary Figure 1A). Only genes with a count per million (CPM) greater than 1 in at least half of the samples were retained. The VennDiagram package in R (version 1.6.20) was implemented to plot Venn diagrams to visualize the relationship between genes in each infection group (Supplementary Figure 1B). The removal of the batch effects resulting from different links, such as methods, cell lines, and experimental designs, was achieved by the “removeBatchEffect” function of the limma R package (version 3.48.3). Finally, the results were checked through the fviz_pca_ind function and the plotted principal component analysis (PCA) diagrams (Supplementary Figure 1C).

The edgeR package (version 3.34.0) was executed for analyzing the differences in mRNA expression among the four HCoVs groups and the mock controls. The adjusted p-values of 0.05 together with absolute log2 (fold change) >1.5 were used as thresholds for identifying the differentially expressed genes (DEGs), and the up- and downregulated DEGs in each group were identified (Figure 1A). The count of the unique and common DEGs in each group was completed using the R package VennDiagram (Figure 1B). Volcano plots of the DEGs in each group were plotted using the ggplot2 R package (version 3.3.5), with the top 10 upregulated DEGs marked.

Function and pathway enrichment analysis of the DEGs were implemented with the software Ingenuity Pathway Analysis (IPA, Ingenuity Systems, Inc.) (31). For the IPA analysis, −log10 (p-value) >1.3 together with absolute Z-score >2 was set as the thresholds to identify the significantly up- or downregulated pathways and functional terms. The enrichment analysis results were visualized using positive and negative bar graphs (Figure 2A). The pathway and function terms were grouped into different categories based on functional characteristics (Figure 2B). Bubble plots of the enriched terms were generated with the ggplot2 R package. The heatmap of cytokine related DEGs was drawn using the R function pheatmap. The Cytoscape software was employed to exhibit the relationship between the genes and enriched functions. Tox function analysis was also conducted using IPA software.

In this study, we analyzed publicly available transcriptome data from cells infected with different coronaviruses (SARS-CoV-2, SARS-CoV, MERS-CoV, and HCoV-229E) (Supplementary Table 1). The results of PCA revealed distinct mRNA expression profiles among the four groups of infected and mock-infected cells (Supplementary Figure 1C). In particular, we observed that the SARS-CoV-2 cluster exhibited a closer distance to the SARS-CoV and MERS-CoV clusters than to the HCoV-229E cluster. This finding indicated marked variations in gene expression between highly pathogenic and lowly pathogenic HCoVs.

To further explore the differences in gene expression patterns among HCoVs, we identified the DEGs among the four groups of infected cells. In total, 949, 171, 278, and 1,198 DEGs were observed in SARS-CoV-2, SARS-CoV, MERS-CoV, as well as HCoV-229E, respectively (Figure 1A). Notably, we observed that a large proportion of DEGs (70.7 and 88.1%) were specifically identified in the SARS-CoV-2 and lowly pathogenic HCoV-229E group, respectively (Figure 1B). These findings indicated significant differences between highly and lowly pathogenic HCoVs. Besides these, the results also revealed that the majority of DEGs were significantly upregulated in the highly pathogenic HCoV-infected groups, particularly in the SARS-CoV-2 group, where 90.3% of the 949 DEGs showed an upregulation (Figure 1A). These observations suggested that infected cells may mount a resistance response to highly pathogenic HCoVs by upregulating the expression of certain genes.

To identify the key genes involved in the host response to each virus, we analyzed the top 10 upregulated DEGs in each group (Figure 1C). In the highly pathogenic SARS-CoV-2 group, 8 out of the top 10 upregulated DEGs were interferon cytokine genes (IFNL2, IFNL3, IFNB1, IFNL1, and IFNL4) and IFN-stimulated genes (GBP5, TNF, and ZBP1), indicating a robust interferon response in SARS-CoV-2 infection. Similarly, in the SARS-CoV group, 6 out of the top 10 upregulated DEGs were interferon cytokine genes (IFNs) and ISGs, including IFNL2, IFNL3, IFNB1, GBP5, TNF, and ZBP1. In the MERS-CoV group, 4 out of the top 10 upregulated DEGs were IFNs and ISGs, including IFNL2, IFNL3, IFNB1, and TNF. Additionally, MRC1, which is related to antigen recognition (32), was also among the top upregulated DEGs in the MERS-CoV group. In the HCoV-229E group, the top upregulated DEGs, such as FAM83B, AATBC, Linc02257, and SPATA12, were mainly related to cell proliferation and migration (33–36) and DNA oxidative damage (37).

To gain insights into the molecular mechanisms underlying the pathogenicity of different HCoVs and to identify enriched pathways and functional categories, we conducted an IPA analysis. The analysis revealed a total of 66, 11, 10, and 15 pathways, as well as 167, 72, 64, and 56 functional terms, in the highly pathogenic SARS-CoV-2, SARS-CoV, MERS-CoV, and lowly pathogenic HCoV-229E groups, respectively (Figure 2A). Notably, the majority of enriched pathways and functions were significantly activated across all four groups [Z-score ≥2 and −log (p-value) >1.3], with the SARS-CoV-2 group displaying the highest number of activated terms (Figure 2A).

To provide a comprehensive overview of the enriched pathways and functional terms, we have grouped them based on their functional characteristics (Figure 2B). The results showed that both the top-ranked pathway or function term category in all three highly pathogenic HCoV groups was “Immune response.” Specifically, 27, 7, and 4 pathways and 100, 43, and 46 functions associated with “Immune response,” were significantly activated across all three highly pathogenic HCoV groups (Figure 2B). In contrast, only one significantly inhibited immune response-related pathway (“Interferon Signaling”) was discovered in the lowly pathogenic HCoV-229E. Moreover, none of the immune-related function terms were significantly activated or inhibited [Z-score ≥2 and −log (p-value) >1.3]. Our findings suggested that highly pathogenic HCoVs induced a robust immune response, while lowly pathogenic HCoVs might not activate extra immune responses or elicit a weaker immune response.

To determine the precise immunological response of cells to the highly pathogenic HCoVs, we intensively investigated the innate and adaptive immune responses of infected cells (Figure 3). For innate immunity, the SARS-CoV-2 infection significantly activated 20 out of the 21 enriched pathways and all of the 45 enriched functions (Figures 3A, B). Of these, 14 pathways and 28 functions were unique to the SARS-CoV-2 infection, such as “Toll-like Receptor Signaling” (Z = 2.333), “Production of Nitric Oxide and Reactive Oxygen Species in Macrophages” (Z = 3.606), “iNOS Signaling” (Z = 2.449), “Natural Killer Cell Signaling” (Z = 2.887), “Interferon Signaling” (Z = 3.5), and “Role of PKR in Interferon Induction and Antiviral Response” (Z = 2.673). In contrast, the SARS-CoV and MERS-CoV groups exhibited fewer enriched innate immune pathways and functions (Figures 3A, B). Only three common innate immune-associated pathways, including “Phagosome Formation,” “Role of Pattern Recognition Receptors in Recognition of Bacteria and Viruses,” and “Necroptosis Signaling Pathway,” as well as 10 innate immune-associated functions, were significantly activated in all three highly pathogenic HCoV groups (Supplementary Figure 2A). Our results indicated that the innate immune response triggered by SARS-CoV-2 infection was distinct from that of SARS-CoV and MERS-CoV and may be potentially capable of inducing stronger innate immune responses than the other two.

This was also true for adaptive immunity in response to highly pathogenic HCoV infections (Figures 3A, C). The SARS-CoV-2 group exhibited 7 significantly activated adaptive immunity pathways and 29 functions, including 6 pathways and 23 functions that were specific to the SARS-CoV-2 infection. Interestingly, three out of the seven adaptive immunity pathways were related to IL-17 (Figure 3A). In contrast, the SARS-CoV and MERS-CoV groups displayed less enrichment in adaptive immune pathways and functions, with only “IL-17 signaling” (Z = 2) and “IL-15 production” (Z = 2.236) being significantly activated in the SARS-CoV and MERS-CoV groups, respectively. Furthermore, there were no common adaptive immunity-associated pathways among the three groups, and only three adaptive immunity-associated functions were significantly activated in all three groups, which were “Cell movement of lymphocytes,” “Lymphocyte migration,” and “T cell migration” (Supplementary Figure 2B). Collectively, these results uncovered that SARS-CoV-2 infection evoked more specific and robust innate and adaptive immune responses than SARS-CoV and MERS-CoV.

The COVID-19 severity is linked to an excessive inflammatory response characterized by the generation of a substantial amount of pro-inflammatory cytokines, a phenomenon known as the cytokine storm (38). To explore the differences in cytokine responses among the highly pathogenic SARS-CoV-2, SARS-CoV, MERS-CoV, and the lowly pathogenic HCoV-229E, we investigated the cytokines and cytokine-related pathways in each group (Figure 4).

Our results unveiled that SARS-CoV-2 infection significantly activated 14 cytokine-related pathways, whereas only two cytokine-related pathways were activated in the SARS-CoV and MERS-CoV group, respectively (Figure 4A). The “Role of Hypercytokinemia/hyperchemokinemia in the Pathogenesis of Influenza” pathway displayed the highest enrichment among all three highly pathogenic HCoV groups. The Z-scores for this pathway in the SARS-CoV-2, SARS-CoV, and MERS-CoV groups were 5.568, 3.162, and 2.828, respectively. Notably, no activated cytokine-related pathways and only one inhibited pathway, “Interferon signaling” (Z = −2.449), were identified in the lowly pathogenic HCoV-229E group.

To further explore the cytokine profiles during infection with these HCoVs, we analyzed the expression levels of the 41 cytokines and 13 cytokine receptor genes that were significantly regulated during infections (Figure 4B). Our analysis signified that the SARS-CoV-2 group exhibited significantly higher cytokine and cytokine receptor gene expression levels than the other three HCoV groups. Among these, the upregulated chemokines and interleukins were the most common cytokine genes in the SARS-CoV-2 group. Specifically, three chemokine genes (CXCL10, CXCL8, and CX3CL1) and five interleukin genes (IL-6, IL-12A, IL-1B, IL-23A, and TSLP) were found to be specifically upregulated. In contrast, no interleukin-related genes were observed to be upregulated in the SARS-CoV or MERS-CoV group. Furthermore, three interleukin genes were found to be downregulated in the HCoV-229E group, with CSF2 (log2FC = −15.644) exhibiting the most significant downregulation among the interleukin genes. These findings suggested that chemokines and interleukins might play a pivotal role in the pathogenesis of SARS-CoV-2 infection.

The IFNs and ISGs, which are mainly involved in innate immunity, have been identified as being highly upregulated in response to infection with highly pathogenic HCoVs (Figure 1C). To further investigate this phenomenon, we analyzed the expression levels of IFNs and ISGs in the four HCoV groups. Our results demonstrated that IFNs were significantly upregulated in cells infected with the highly pathogenic SARS-CoV-2, SARS-CoV, and MERS-CoV, whereas no significant upregulation was observed in cells upon infection with the lowly pathogenic HCoV-229E (Figure 4B). Of particular interest, the top three upregulated IFNs in highly pathogenic HCoVs were type III IFNL2, IFNL3, and type I IFNB1. Notably, our analysis revealed that type III IFNL1 was commonly upregulated in both SARS-CoV-2 and SARS-CoV groups, while type III IFNL4 was uniquely identified in the SARS-CoV-2 group. Furthermore, we noted that the “Interferon signaling” pathway was significantly activated (Z = 3.5) in the SARS-CoV-2 group, whereas this pathway was significantly inhibited (Z = −2.449) in the lowly pathogenic HCoV-229E (Figure 4A). These results highlighted the critical role of interferon in the host’s immune response to highly pathogenic HCoV infections.

Further analysis of the expression levels of 628 known ISGs (39) in the four HCoV groups revealed that SARS-CoV-2 exhibited significantly higher levels of ISG expression than the other three HCoVs (Figure 5 and Supplementary Table 2). Specifically, we observed that 115 ISGs were distinctly upregulated in cells with the SARS-CoV-2 infection, implying a robust interferon response (Figure 5A). Notably, 95 of these ISGs (82.6%) were uniquely upregulated in the SARS-CoV-2 group (Supplementary Table 2). Our investigation further revealed an upregulation of several ISGs with known antiviral functions, including IFIT, IFITM, RSAD2, ZNFX1, TRIM21, ISG15, and ISG20 (40–45). Moreover, we observed an upregulation of GBP genes (GBP3, GBP4, GBP5, and GBP6), which have been demonstrated to modulate extensive innate immune responses against various pathogens (46). In addition, the upregulation of genes that potentiated IFN signaling was also observed, including TAP1, STAT1, and XAF1 (47–50). In contrast to SARS-CoV-2, infections with SARS-CoV and MERS-CoV induced significantly fewer ISGs (17 and 12, respectively), and infection with the lowly pathogenic HCoV-229E resulted in the downregulation of 17 out of 28 ISGs (Figure 5A).

Inhibition of virus replication has been demonstrated in many ISGs. To investigate the effectiveness of ISGs in inhibiting virus replication, we conducted a functional analysis that identified 52 ISGs that were significantly enriched in ten functional terms associated with virus replication (Figures 5B, C). These functional terms displayed significant inhibition in the highly pathogenic SARS-CoV-2 group, including “Viral life cycle” (Z = −4.491), “Replication of virus” (Z = −4.392), and “Replication of coronavirus” (Z = −3.517). Furthermore, it was found that more than half of the DEGs (52–64%) enriched in each functional term related to virus replication were ISGs (Figure 5B). This indicated a strong IFN/ISG response that contributed to inhibiting SARS-CoV-2 replication. In contrast, only three functional terms related to virus replication were inhibited in the SARS-CoV group (Figure 5C), including “Replication of Herpesviridae” (Z = −2.214), “Replication of vesicular stomatitis virus” (Z = −2.425), and “Replication of viral replicon” (Z = −2.236), which included three ISGs (TNF, ZBP1, and CCL5). Interestingly, no such virus replication-related functions were found in the MERS-CoV group.

To further investigate the putative role of ISGs in shaping the clinical course of COVID-19, we conducted an analysis using publicly accessible PBMC transcriptome data sourced from the GSE152418 dataset. This dataset comprises PBMC samples collected from a cohort of 15 COVID-19 patients exhibiting various degrees of symptom severity, with categorizations ranging from moderate, severe, to critical, alongside a control group consisting of 17 healthy individuals. We identified a total of 107 ISGs, demonstrating distinctive expression profiles across the spectrum of COVID-19 symptom severity, with marked upregulation in moderate cases and conspicuous downregulation in severe patients, particularly those requiring ICU admission (Supplementary Figure 3 and Supplementary Table 3). Of notable significance, among these ISGs, we found that 13 ISGs exhibited a significant upregulation in patients with moderate COVID-19 symptoms, while 17 ISGs displayed a significant downregulation in patients necessitating ICU admission. This observed pattern in ISG expression harmoniously reinforces the conclusions drawn from our prior investigation of lung cell lines. Notably, these findings suggested the potential role of insufficient ISG expression and its subsequent suppression as plausible contributory factors to the manifestation of more severe clinical symptoms in COVID-19 patients.

Overall, our findings demonstrated that SARS-CoV-2 infection induced a distinctive and strong interferon response characterized by the upregulation of multiple ISGs with antiviral functions. These results provided valuable insights into the pathological characteristics of SARS-CoV-2, which showed a relatively milder pathogenesis and a lower fatality rate in contrast to SARS-CoV and MERS-CoV.

We evaluated the potential multi-tissue injury using IPA tox-function analyses. Our results revealed that SARS-CoV-2 infection significantly upregulated multi-organ damage related functions, particularly in the liver, kidney, and vasculature (Figure 6). In contrast, fewer relevant terms were significantly enriched in the other three HCoV groups.

In the SARS-CoV-2 group, four functions associated with the angiocarpy system were significantly upregulated (Figure 6). These activated functions were all associated with elevated blood pressure. In contrast, the other four downregulated functions in this group were associated with congenital heart disease, including atrial and ventricular septal defects (Figure 6). These findings suggested that SARS-CoV-2 infection can increase blood pressure and reduce the risk of congenital heart malformations, which may be realized through the effect of ACE2 on the renin-angiotensin system (RAS) (51). In general, our findings unveiled that the SARS-CoV-2 infection had a more pronounced effect on multi-organ toxicity than other HCoV infections.