Comparisons of the immunological landscape between COVID-19, influenza, and respiratory syncytial virus patients by clustering analysis.

BACKGROUND: COVID-19 has stronger infectivity and a higher risk for severity than most other contagious respiratory illnesses. The mechanisms underlying this difference remain unclear.
METHODS: We compared the immunological landscape between COVID-19 and two other contagious respiratory illnesses (influenza and respiratory syncytial virus (RSV)) by clustering analysis of the three diseases based on 27 immune signatures' scores.
RESULTS: We identified three immune subtypes: Immunity-H, Immunity-M, and Immunity-L, which displayed high, medium, and low immune signatures, respectively. We found 20%, 35.5%, and 44.5% of COVID-19 cases included in Immunity-H, Immunity-M, and Immunity-L, respectively; all influenza cases were included in Immunity-H; 66.7% and 33.3% of RSV cases belonged to Immunity-H and Immunity-L, respectively. These data indicate that most COVID-19 patients have weaker immune signatures than influenza and RSV patients, as evidenced by 22 of the 27 immune signatures having lower enrichment scores in COVID-19 than in influenza and/or RSV. The Immunity-M COVID-19 patients had the highest expression levels of ACE2 and IL-6 and lowest viral loads and were the youngest. In contrast, the Immunity-H COVID-19 patients had the lowest expression levels of ACE2 and IL-6 and highest viral loads and were the oldest. Most immune signatures had lower enrichment levels in the intensive care unit (ICU) than in non-ICU patients. Gene ontology analysis showed that the innate and adaptive immune responses were significantly downregulated in COVID-19 versus healthy individuals.
CONCLUSIONS: Compared to influenza and RSV, COVID-19 displayed significantly different immunological profiles. Elevated immune signatures are associated with better prognosis in COVID-19 patients.

Background

The coronavirus disease 2019 (COVID-19) caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has been reported to exceed 126 million cases and resulted in more than 2.7 million deaths globally as of March 26, 2021 [1]. COVID-19 shares similar clinical presentations with influenza and respiratory syncytial virus (RSV), such as cough, fever, severe lung infections, and occasional respiratory failure causing deaths [2], nevertheless, COVID-19 has also demonstrated distinct clinical characteristics, such as gastrointestinal discomfort, nerve injuries, and hypogeusia [3], [4]. Compared to other contagious respiratory illnesses, such as influenza and RSV, COVID-19 has stronger infectivity and a higher risk for severity, such as respiratory distress and systematic inflammation [3], [5]. More seriously, COVID-19 often has a long incubation period that significantly increases its transmissibility [3], [5], [6]. Unlike many other respiratory viruses, SARS-CoV-2 infection may result in host hyperinflammatory response, characterized by a high level of circulating inflammatory cytokines and interleukine-6 (IL-6) [7]. Several studies have shown that a large amount of type I interferon promoted inflammatory cytokines production and resulted in lethal pneumonia [8], [9]. Immune dysfunction is also a characteristic of COVID-19 [10]; as shown in severe COVID-19 cases, the total numbers of T and B cells are reduced, resulting from overproduced pro-inflammatory cytokines [11].

To uncover potential mechanisms underlying the significantly different clinical characteristics of COVID-19, we compared the immunological landscape between COVID-19, influenza, and RSV based on their gene expression profiles. We performed a clustering analysis to identify immune subtypes of COVID-19, influenza, and RSV based on the enrichment levels of 27 immune signatures. We identified three immune subtypes with high, medium, and low immune signatures, respectively, and dissected these subtypes’ molecular features. We also explored correlations between these immune signatures and various clinical characteristics in COVID-19 patients. This study aimed to provide new insights into the immunological landscape of COVID-19.

Methods

Datasets

We collected eight gene expression profiling datasets, including four COVID-19 RNA-seq datasets [12], [13], [14], [15], one RSV RNA-seq dataset from the Gene Expression Omnibus database (GEO, https://www.ncbi.nlm.nih.gov/geo/), one RSV microarray dataset [16], one influenza microarray dataset from GEO, and one microarray dataset for both RSV and influenza [17] (Table 1). First, we collected raw data of these datasets, composed of 1,036 patients and 273 healthy controls. The 1,036 patients included 842 COVID-19, 44 influenzas, and 150 RSV patients. The median age of these patients was 53 years (interquartile range: 38–67 years). Among these patients, 404 (38.96%) were females, 348 (33.61%) were males, and 284 (27.43%) were unidentified. There were 100 COVID-19 cases in the dataset GSE157103, including 50 intensive care unit (ICU) cases and 50 non-ICU cases. The four COVID-19 datasets included patients from the USA. The two influenza datasets included patients from Canada and the UK, respectively, and the two RSV datasets included patients from the UK and the USA. The patients’ tissues in these datasets originated from the upper airway, peripheral blood mononuclear cell, leukocytes, and blood. Next, we performed quality control analysis to detect outliers in each dataset using WGCNA [7]. Finally, we merged these datasets using the “merge” function in the R package “base” and performed the adjustment of batch effects and normalization of the combined data using the “normalizeBetweenArrays” function in the R package “limma.” Because this study aimed to compare immune signatures between COVID-19, influenza, and RSV, the bacterial and other anonymous viruses-associated samples were excluded.

Table 1

A summary of the datasets analyzed.

Dataset	Disease	Platform	Number of patients	Number of controls	Excluded outliers	Source of samples
GSE152075	COVID-19	RNA-Seq	430	54	3	upper airway
GSE152418	COVID-19	RNA-Seq	76	69	3	Peripheral blood mononuclear cell
GSE156063	COVID-19	RNA-Seq	234	41	1	upper airway
GSE157103	COVID-19	RNA-Seq	102	26	0	leukocytes
GSE155237	RSV	RNA-Seq	46	32	0	upper airway
GSE100161	RSV	Microarray	70	14	1	blood
GSE21802	Influenza	Microarray	19	4	0	blood
GSE42026	Influenza, RSV	Microarray	59	33	0	blood

Clustering

We performed an extensive literature search for immune signatures associated with respiratory viruses’ infection and collected 27 immune signatures involving 136 marker genes [18], [19], [20], [21]. We first quantified the enrichment levels of the 27 immune signatures in each sample by the single-sample gene-set enrichment analysis (ssGSEA) [22]. Based on the ssGSEA scores of the 27 immune signatures, we performed hierarchical clustering of the samples in the merging dataset based on the Euclidean metrics.

Gene-set enrichment analysis

We identified gene ontology (GO) and KEGG pathways upregulated and downregulated in COVID-19, influenza, RSV, and their immune subtypes by GSEA [23] (R implementation) with a threshold of adjusted P value < 0.05.

Identification of gene modules enriched in immune subtypes

We used WGCNA [7] to identify gene modules significantly enriched in immune subtypes based on gene co-expression network analysis. Based on the expression correlations between the gene modules’ hub genes, we identified the GO terms significantly correlated to specific traits.

Statistical analysis

In the comparison of different groups of data, we used Student’s t-test or ANOVA test if they were normally distributed; otherwise, we used the Mann–Whitney U test or Kruskal–Wallis test. We used the Pearson correlation coefficient (R) to evaluate the expression correlation between two genes and the Spearman correlation coefficient (ρ) to assess correlations between other variables. We adjusted P values in multiple tests using the false discovery rate (FDR), which were calculated by the Benjamini-Hochberg method. We performed all statistical analyses using R programming (version 4.0.2).

Results

The immune subtypes of COVID-19, influenza, and RSV

Based on the 27 immune signatures’ enrichment scores, we hierarchically clustered the 809 virus-infected patient samples based on the Euclidean distance metric. We obtained three clear clusters, termed Immunity-H, Immunity-M, and Immunity-L, which displayed the highest, medium, and lowest immune signatures levels, respectively (Fig. 1A). In fact, among the 27 immune signatures, 21 displayed significantly higher enrichment scores in Immunity-H than in Immunity-M and Immunity-L (one-way ANOVA test, P < 0.05) (Fig. 1B). We found 20%, 35.5%, and 44.5% of COVID-19 cases included in Immunity-H, Immunity-M, and Immunity-L, respectively; all influenza cases were included in Immunity-H; 66.7% and 33.3% of RSV cases belonged to Immunity-H and Immunity-L, respectively. These results indicate that most COVID-19 patients have weaker immune signatures than the patients infected with other viruses. This was evidenced by that 22 of the 27 immune signatures displayed lower enrichment scores in COVID-19 than in influenza and/or RSV patients (Fig. 1C).

Fig. 1

Identification of immune subtypes of COVID-19, influenza, and RSV. (A) Hierarchical clustering yields three subtypes: Immunity-H, Immunity-M, and Immunity-L, based on the enrichment scores of 27 immune signatures. (B) Comparisons of the enrichment scores of 27 immune signatures between the three immune subtypes. (C) Comparisons of the enrichment scores of 27 immune signatures between COVID-19, influenza, and RSV patients. (D) Comparisons of the expression levels of ACE2 and IL-6 between the three immune subtypes and between the three immune subtypes containing COVID-19 only. (E) Comparisons of COVID-19 patients’ ages between the three immune subtypes. (F) Comparisons of viral loads in COVID-19 patients between the three immune subtypes. (G) Spearman correlation between the expression levels of ACE2 and IL-6 and viral loads in COVID-19 patients. (H) Pearson correlation between the expression levels of ACE2 and IL-6 in COVID-19 patients. * P < 0.05, ** P < 0.01, *** P < 0.001.

Angiotensin-converting enzyme 2 (ACE2) is the host cell receptor of SARS-CoV and SARS-CoV-2 [3]. We found that ACE2 expression levels were the highest and lowest in Immunity-M and Immunity-H, respectively (two-tailed Student’s t-test, P < 0.05) (Fig. 1D). IL-6 has been indicated as an adverse prognosticator of COVID-19. Notably, the expression levels of IL-6 were the highest and lowest in Immunity-M and Immunity-H, respectively (P < 0.05) (Fig. 1D). When the cases of influenza and RSV were excluded, the expression levels of ACE2 and IL-6 were still the highest and lowest in Immunity-M and Immunity-H, respectively (P < 0.05) (Fig. 1D). It indicated that the COVID-19 cases in Immunity-M and Immunity-H may have different prognoses. Indeed, the COVID-19 patients were the youngest and oldest in Immunity-M and Immunity-H, respectively (Kruskal–Wallis test, P = 0.003) (Fig. 1E). Moreover, the COVID-19 patients had the lowest and highest viral loads in Immunity-M and Immunity-H, respectively (P = 0.004) (Fig. 1F). The expression levels of ACE2 and IL-6 were inversely correlated with viral loads in COVID-19 patients (Spearman correlation, ρ = − 0.29, P < 0.001, ρ = − 0.12, P < 0.015, respectively) (Fig. 1G). Meanwhile, we observed a strong expression correlation between ACE2 and IL-6 in COVID-19 patients (Pearson correlation R = 0.65, P < 0.001) (Fig. 1H).

Overall, these results indicate that the host immune response to SARS-CoV-2 infection is likely to be more heavily dysregulated than the response to the other viruses’ infection.

Differentially expressed genes and pathways between COVID-19 patients and influenza and RSV patients

ACE2 and TMPRSS2 are two crucial molecules for SARS-CoV-2 invading host cells [24]. We found that ACE2 expression levels were significantly higher in COVID-19 patients than in healthy controls, influenza, and RSV patients (P < 0.05) (Fig. 2A). TMPRSS2 expression levels were significantly lower in COVID-19 patients than in healthy controls, while it was significantly higher in COVID-19 patients than in influenza and RSV patients (P < 0.05) (Fig. 2A). Furthermore, we identified 41 genes showing significantly different expression levels between COVID-19 patients and healthy controls (FDR < 0.1, fold change greater than 2). Among them, 12 genes encoding transcriptional factors (TFs) were upregulated in COVID-19 patients. The 12 TFs included five viral replication and cellular response inhibitors (IFIT3, DDX58, MX1, OAS1, and TRIM22), one chemokine (CXCL10), two interferon-gamma signaling mediators (XAF1 and GBP1), one macrophage activation modulator (EPSTI1), two immune response modulators (STAT1 and EIF2AK2), and one interferon’s antiviral activity modulator (IFITM3). In influenza and RSV patients, 14 and 40 TF genes were upregulated relative to their healthy controls, respectively. IFITM3 was the only TF gene commonly upregulated in COVID-19, influenza, and RSV patients. However, IFITM3 had significantly lower expression levels in COVID-19 than in influenza and RSV patients and showed strong positive expression correlations with immune signatures (Fig. 2B). Again, these results indicate the differences in host immune response in COVID-19 versus influenza and RSV patients.

Fig. 2

Identification of differentially expressed genes and gene ontology (GO) between COVID-19 patients and influenza and RSV patients. (A) Comparisons of the expression levels of ACE2 and TMPRSS2 between virus-infected patients and healthy controls. (B) IFITM3 shows significantly lower expression levels in COVID-19 than in influenza and RSV patients and strong positive expression correlations with immune-promoting signatures. GO terms upregulated and downregulated in COVID-19 (C), influenza (D), and RSV (E) patients versus healthy controls. (F) GO terms highly enriched in COVID-19 versus influenza. * P < 0.05, ** P < 0.01, *** P < 0.001.

GSEA identified 81 GO terms enriched in COVID-19 patients, which were mainly involved in the regulation of metabolism, extracellular matrix, Toll-like receptors (TLRs), Hippo signaling, and immune response. In contrast, two GO terms involved in regulating ubiquitin and innate immunity were downregulated in COVID-19 patients (Fig. 2C). In addition, GSEA identified 272 GO terms enriched in influenza patients, most of which were involved in the regulation of inflammatory and immune responses (Fig. 2D). 49 GO terms were enriched in RSV patients, most of which were associated with immune regulation (Fig. 2E). Surprisingly, none of these GO terms were commonly upregulated or downregulated among the three groups of virus-infected patients, indicating unique molecular alterations in each group of virus-infected patients.

We also conducted a comparison of gene expression profiles between the COVID-19 patients and influenza and RSV patients. We found 126 GO terms highly enriched in COVID-19, most of which were involved in the regulation of metabolism, chemokines, and cytokines. In addition, 32 GO terms were downregulated in COVID-19, which were mainly associated with genomic stability regulation (Fig. 2F and Supplementary Fig. S1A). Besides, we compared gene expression profiles between Immunity-H, Immunity-M, and Immunity-L. We found 107 GO terms significantly upregulated in Immunity-H versus Immunity-L, and 136 GO terms significantly upregulated in Immunity-M versus Immunity-H (Supplementary Fig. S1B, C, D). As expected, these GO terms were mainly involved in the regulation of the immune system. Collectively, these results indicate that the host immune response to SARS-CoV-2 is significantly different from that of the other viruses and that SARS-CoV-2 infection may affect the host metabolic process and genomic stability.

Gene modules enriched in immune subtypes and diseases

We used WGCNA [7] to identify similar co-expression gene modules that differentiated patients by subtype (Immunity-H, Immunity-M, and Immunity-L) or disease (COVID-19, influenza, and RSV). We identified 12 gene modules highlighted in different colors: turquoise, green, purple, tan, magenta, yellow, blue, brown, black, pink, green-yellow, and red (Fig. 3A). As expected, the immune response, a representative GO term for the yellow module, was highly enriched in Immunity-H versus Immunity-M and Immunity-L (r = 0.26, P = 4.0 × 10−14). The viral transcription was also highly enriched in Immunity-H (r = 0.77, P = 3.0 × 10−160), while it was impoverished in Immunity-M (r = −0.28, P = 2.0 × 10−16) and Immunity-L (r = −0.47, P = 2.0 × 10−46). It is reasonable because high levels of viral replication may cause a stronger host immune response. The extracellular matrix was highly enriched in Immunity-M (r = 0.47, P = 2.0 × 10−46), while it was downregulated in Immunity-H (r = −0.52, P = 1.0 × 10−57). The defense response to virus was also upregulated in Immunity-M (r = 0.12, P = 5.0 × 10−4), while it was downregulated in Immunity-L (r = −0.15, P = 1.0 × 10−5) and showed no significant correlation with Immunity-H (r = 0.045, P = 0.2). It indicated that Immunity-M and Immunity-L had the strongest and weakest defense response to virus infection among the three subtypes. The epithelial cell proliferation was highly enriched in Immunity-L and Immunity-M, while it was downregulated in Immunity-H.

Fig. 3

Gene modules enriched in immune subtypes and diseases. 12 gene modules (indicated in turquoise, green, purple, tan, magenta, yellow, blue, brown, black, pink, green-yellow, and red colors) significantly differentiated these patients by immune subtype (A) and disease (B) identified by WGCNA. The representative gene ontology (GO) for gene modules, correlation coefficients, and P values (in parenthesis) are shown. Gray module contains unassigned genes and has no GO associated with it. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Likewise, WGCNA identified 12 gene modules that significantly differentiated these patients by disease (COVID-19, influenza, and RSV) (Fig. 3B). Notably, ten gene modules displayed higher enrichment in COVID-19 relative to influenza and RSV. Interestingly, the immune response was highly enriched in COVID-19 versus influenza and RSV (r = 0.46, P = 8.0 × 10−44). The defense response to virus was also enriched in COVID-19 compared to influenza and RSV (r = 0.1, P = 0.004). It suggests that COVID-19 patients have a stronger immune response to virus infection than influenza and RSV patients, which may be described as an inflammatory response to antigenic stimulus. In addition, two gene modules for which representative GO terms were nucleus and chromosome, respectively, were significantly upregulated in COVID-19. In contrast, the viral transcription was significantly downregulated in COVID-19 versus influenza and RSV (r = −0.72, P = 6.0 × 10−127).

Associations between immune signatures and clinical features in COVID-19 patients

We further analyzed associations between the 27 immune signatures and clinical features in COVID-19 patients in GSE157103. The clinical features included ICU, the Acute Physiology and Chronic Health Evaluation (APACHE II), the Sequential Organ Failure Assessment (SOFA), Charlson score, C-Reactive Protein (CRP), D-dimer, ferritin, fibrinogen, hospital free days (45 days follow-up), lactate, procalcitonin, and ventilator-free days. Both APACHE II and SOFA scores measure the severity of ICU patients [25]. The Charlson score predicts the ten-year mortality for a patient with a range of comorbid conditions [26]. The CRP, D-dimer, ferritin, fibrinogen, lactate, and procalcitonin are laboratory measurements whose high levels were correlated with the severity of COVID-19 [27]. Notably, we found that 20 immune signatures had significantly lower enrichment levels (ssGSEA scores) in ICU than in non-ICU patients (P < 0.05) (Fig. 4A). Only two immune signatures (type II IFN reponse and inflammatory cytokines) had significantly higher ssGSEA scores in ICU than in non-ICU pateints. Moreover, most of these immune signatures’ ssGSEA scores were negatively correlated with the levels of CRP, D-dimer, ferritin, fibrinogen, lactate, and procalcitonin (Fig. 4B). In constrast, a majority of these immune signatures’ ssGSEA scores were positively correlated with hospital free days and ventilator-free days. Collectively, these results suggest that elevated immune signatures are likely to be associated with better prognosis in COVID-19 patients.

Fig. 4

Associations between immune signatures and clinical features in COVID-19 patients. (A) Comparisons of the enrichment levels of 27 immune signatures between ICU and non-ICU cases. Two-tailed Student’s t-test P values are indicated. (B) Heatmap showing Spearman correlations between the enrichment levels of 27 immune signatures and clinical parameters in COVID-19. * P < 0.05, ** P < 0.01, *** P < 0.001. The results in (A, B) were obtained by analyzing the dataset GSE157103.

Discussion and conclusions

Respiratory tract infection is liable for undue disease burden [28]. The innate and adaptive immune responses to respiratory virus infection are critical for sustaining a healthy respiratory system and preventing pulmonary disease [29]. Although COVID-19 shares certain primary symptoms with influenza and RSV [30], [31], its transmissibility and destructivity far exceed influenza and RSV [1]. The unique host immune response in COVID-19 patients could greatly contribute to this disease’s high transmissibility and destructivity [32]. For example, SARS-CoV-2 infection may not always employ interferon production but delay or suppress interferon II or III production in lung tissues [33].

Upon infection, the influenza virus is recognized by host cells through pathogen-associated molecular patterns (PAMPs) [34]. PAMPs are recognized by TLRs, RIG-I, and NLRP3 to generate pro-inflammatory cytokines and type I interferon production, which has strong antiviral activity by limiting virus replication in host cells. Type I interferons can induce STAT activation and transcription of the interferon-stimulated gene to influence innate and adaptive immune responses. Furthermore, host immune cells recruit alveolar macrophages to activate apoptotic influenza virus-infected cells, while blood-derived macrophages infected with influenza viruses produce large amounts of pro-inflammatory cytokines [34]. RSV infection elicits a strong innate immune response, especially respiratory tract neutrophil response [35]. The peak of neutrophil response often concurs with maximum clinical severity and viral load. However, because RSV infection can reduce natural killer (NK) cell count to regulate the inflammatory response, a lower NK cell count is correlated with more serious infection [36]. Upon acute infection, the adaptive immune response among RSV patients is activated via recruiting CD8+ T cells to clear viruses. Meanwhile, B cells play an important role in the production of protective antibodies. RSV infection also triggers interferon II production, which plays a strong protective role against virus infection [37]. Our results demonstrated that compared to COVID-19, influenza and RSV patients elicited a stronger innate immune response, e.g., more interferons I production, macrophages, NK cells, and dendritic cells. Furthermore, influenza and RSV patients elicited a stronger adaptive immune response described by more B and T cell production.

Compared to influenza and RSV, COVID-19 is likely to display lower immune signatures and higher IFN II production, leading to poor prognosis and worse outcomes. The innate and adaptive immune responses are significantly lower in COVID-19 than in the other viruses-infected patients, while chemoattractant protein-1 (IL-6) is upregulated in COVID-19 patients. Meanwhile, we found that pro-inflammatory cytokines were remarkably higher in COVID-19 than in influenza and RSV patients. In contrast, anti-inflammatory cytokines were significantly lower in COVID-19 patients, which may describe the cytokines storm among COVID-19 patients. Interestingly, COVID-19 patients displayed a higher expression of cytolytic perforin than the patients infected with the other viruses. However, since perforin’s important role lies in the cytotoxic activity against pathogens, abnormal production of this protein could be responsible for resistance to the cytotoxic activity of NK and CD8+ T cells. It has been shown that target pathogen cells could employ various escape mechanisms against perforin activity [38].

The clustering analysis showed that Immunity-M and Immunity-L harbored 80% of COVID-19 patients, suggesting that most COVID-19 patients have an inadequate immune response to virus infection. Furthermore, we found that molecular and clinical features were significantly different between COVID-19 patients in Immunity-M and Immunity-H. For example, Immunity-M patients had significantly higher expression levels of ACE2 and IL-6 while lower viral loads than Immunity-H patients. Immunity-M patients were younger than Immunity-H patients. The lower viral loads in Immunity-M than in Immunity-H patients could explain why the former had a weaker immune response than the latter. In addition, we found that the expression levels of ACE2 and IL-6 had a negative correlation with viral loads in COVID-19 patients. The negative association between ACE2 and viral loads could be attributed to that ACE2 can reduce viral loads [39].

Interestingly, younger COVID-19 patients have a lower immune response and fewer viral loads in Immunity-M than older patients in Immunity-H. A potential explanation could be that the immune system is often dysfunctional and heterogeneous in older people [40]. Through aging, the immune system often changes in two ways [40]. One is immunosenescence; that is, dysfunctional cells become epigenetically locked into a pro-inflammatory state, secreting cytokines and chemokines. The other one is inflamm-aging arising from the overactive and ineffective alert system [40]. In addition, when elderly alveolar TLRs detect the pathogen, they respond by producing type I interferons to attract immune cells to the site of infection. Furthermore, although the number of macrophages increases with age, their plasticity to convert between pro- and anti-inflammatory states is reduced [33]. Therefore, the higher immune signatures in Immunity-H are likely to be associated with dysfunctional immune responses in older COVID-19 patients, leading to worse clinical outcomes, as observed in the previous results. Besides, we found that IFITM3, a gene playing a critical role in promoting the immune system’s defense against viruses, was significantly downregulated in COVID-19 versus influenza and RSV patients. Again, it suggests an inadequate immune response in COVID-19 versus influenza and RSV patients, although the levels of pro-inflammatory cytokines are remarkably higher in COVID-19 than in influenza and RSV patients.

In addition to the pathways of immune responses, several metabolic pathways were dysregulated in COVID-19 patients. These metabolic pathways are involved in the abnormal activity of mitochondrial respiratory chain, negative regulation of signaling receptor activity, dysregulation of extracellular matrix binding, dysregulation of nucleoside phosphate catabolic process, regulation of cytokinesis and chemotaxis, regulation of TOR signaling, and decreased pulmonary function. These data may explain why many COVID-19 patients exhibited severe cardiovascular and pulmonary injury [41]. Moreover, the extracellular matrix regulation is enriched in several COVID-19 patients, which plays a vital role in establishing and maintaining immune homeostasis [42]. Also, we found a strong correlation between metabolic dysregulation and immunogenesis in COVID-19 patients. For example, the enrichment levels of MHC class 1 were positively correlated with decreased pulmonary function, extracellular matrix binding, abnormal activity of mitochondrial respiratory chain, abnormality of the diaphragm, positive regulation of DNA metabolic process, and negative regulation of ToR signaling. In contrast, they were negatively correlated with the regulation of positive chemotaxis. Regulation of positive chemotaxis and positive regulation of cytokinesis were negatively associated with most of the 27 immune signatures, while the abnormal activity of the mitochondrial respiratory chain was positively associated with a majority of these immune signatures (Supplementary Fig. S2).

Our results also showed that some COVID-19 patients had a stronger immune response than other COVID-19 patients, while these patients often had more severe clinical outcomes. The main reason could be that their hyperinflammatory reaction to SARS-CoV-2 infection, known as cytokine storm [43], [44], resulted in serious organ damage. Besides, our results showed that ACE2 expression levels were significantly lower in immunity-H than in immunity-M and immunity-L. It indicates that ACE2 upregulation could lead to antiviral immunosuppression to cause a worse prognosis in COVID-19 patients [45], [46].

Although our study provides a strong base for further investigation into different immunological profiles of COVID-19, influenza, and RSV patients, it should be noted that this study has several limitations. First, there are insufficient clinical data (e.g., disease severity) associated with our analyses, especially for influenza and RSV patients. Second, this study was conducted without considering the patients’ underlying disease history due to limited detailed data available. Finally, this study’s underlying bias in study design, patient populations, source tissues, and technical errors cannot be ruled out as residual confounding. However, our data still provide a clear insight into the difference in immune responses between COVID-19, influenza, and RSV patients. The significantly different immune responses in different COVID-19 patients could be responsible for significantly different clinical outcomes between COVID-19 cases.

Ethics approval and consent to participate

Ethical approval and consent to participate were waived since we used only publicly available data and materials in this study.

Consent for publication

Not applicable.

Funding

This work was supported by the Science and Technology Project of Traditional Chinese Medicine in Zhejiang Province (grant number: 2017ZKL003).

CRediT authorship contribution statement

Zeinab Abdelrahman: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Data curation, Visualization, Writing - original draft. Zuobing Chen: Methodology, Software, Validation, Formal analysis, Investigation, Data curation, Visualization, Funding acquisition. Haoyu Lyu: Software, Validation, Formal analysis, Investigation, Data curation, Visualization. Xiaosheng Wang: Conceptualization, Methodology, Resources, Investigation, Writing - review & editing, Supervision, Project administration, Funding acquisition.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.