We included patients who consented to donate blood to the University of Virginia (UVA) Intensive Care Unit (ICU) Biorepository. We then selected patients with confirmed COVID-19 respiratory failure, viral, non-COVID-19 respiratory failure, and patients with non-viral causes of respiratory failure (presumed bacterial infections) in the UVA ICU Biorepository to serve as a comparison cohort. Control patients were patients admitted on mechanical ventilation without respiratory failure (usually intubated for airway protection). We followed the Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) guidelines (37) and our study complied with all principles outlined in the Declaration of Helsinki. All study protocols were approved by the UVA Institutional Review Board for Health Sciences Research (Protocol #21101). Respiratory failure was defined as patients with acute respiratory distress syndrome (ARDS) using Berlin criteria (38) who were on mechanical ventilation in the ICU. COVID-19 diagnoses were confirmed by RealTime SARS-CoV-2 assay performed on the m2000 system (Abbott Molecular Inc.; Des Plaines, IL).

Serum cytokine and chemokine levels were measured by Merck Millipore MILLIPLEX Human Cytokine/Chemokine/Growth Factor Panel A (HCYTA-60K) and Panel II (HCYP2MAG) assays.

Blood was collected into PAXgene® Blood RNA tubes upon admission to the ICU, and after 24 and 72 hours. RNA was prepared and sequenced by Genewiz (Azenta Life Sciences). RNA was isolated using a Qiagen total RNA isolation kit, followed by RNA-seq library preparation using standard Illumina protocols with rRNA and globin depletion. RNA was sequenced on an Illumina HiSeq 4000.

Three independent whole blood datasets were analyzed. In the first dataset of patients from the UVA ICU Biorepository, RNA-seq data was obtained from 13 COVID-19 AHRF patients, 8 viral AHRF, 5 non-viral AHRF, and 5 control ICU patients. COVID-19 patients were subdivided into 2 groups (COVID Group 1 and COVID Group 2) based on separation by PCA of the top 500 variable genes. Additional publicly available RNA-seq datasets (GSE161731; GSE172114), were obtained and analyzed as confirmation studies. Additional dataset details can be found in Supplementary Table 1 .

For RNA-seq analysis, the quality of raw FASTQ reads was analyzed using FASTQC (39) to identify the poor-quality reads and the adaptor contamination. Adaptors and low-quality sequencing reads were trimmed using Trimmomatic (40) and reads before 14bp were discarded. The clean raw sequencing reads were aligned to human reference genome (Gencode hg38) using STAR(v2) (41). The SAM files were converted into BAM files using sambamba (42). The aligned BAM files were fed to read summarization program featureCounts (43), to assign the sequencing reads to genomic features. The differential gene expression between COVID-19 and normal patients was carried out using the R/Bioconductor package DESeq2. The raw counts were normalized within the DESeq2 analysis pipeline using the median of ratios method and genes with low expression were filtered using HTSFilter (44). The p-values and adjusted p-values (FDR) were calculated for each gene. Genes with FDR < 0.05 were determined as significant differentially expressed genes. PCA plots of log₂ transformed gene expression values were generated using the plotPCA function from the DESeq2 package.

The R/Bioconductor package GSVA (45) (v1.25.0) was used as a non-parametric, unsupervised method to estimate the variation in enrichment of pre-defined gene sets in RNA-seq dataset samples as previously described (36) (www.bioconductor.org/packages/release/bioc/html/GSVA.html). In brief, a matrix of log₂ transformed gene expression values for each sample and pre-defined gene sets were used as inputs for the GSVA algorithm. Then enrichment scores (GSVA scores) for each gene set were calculated using a Kolmogorov Smirnoff (KS)-like random walk statistic. GSVA scores for each patient and control were calculated and normalized to scores between -1 (no enrichment) and +1 (enriched). Significance of gene set enrichment between cohorts was calculated using a Welch’s t-test and p-value < 0.05 was considered significant.

Input gene sets used for GSVA analysis were previously used for the analysis of COVID-19 patient datasets (36) and can be found in Supplementary Table 4 .

Multivariable linear regression analysis was performed with MaAsLin 2 (46), a Bioconductor R package that helps to determine the association between gene expression and complex metadata features. MaAsLin 2 is a statistical method that relies on general linear regression models which can test for the association between various functional and cell specific modules versus individual discrete and categorical clinical variables. Computed GSVA scores and patient metadata were used as input for the MaAsLin 2 function in R with normalization method and transformation method applied “NONE”, analysis method “LM”, and correction method “BH”. The significant associations with clinical variables were visualized using scatterplots and box plots.

Additional linear regression analyses for individual patient cohorts and between PC GSVA scores and log₂ expression of Ig heavy chain transcripts were performed in GraphPad Prism (v 9.1.0; www.graphpad.com). For each analysis, the r² value indicating the Goodness of Fit and the p-value testing the significance of the slope are displayed.

Patient demographic data from COVID Group 1 and Group 2 were compared using an unpaired t-test with Welch’s correction for continuous variables. Pearson’s chi-squared tests were used to analyze count data between groups. Cytokine data for COVID Group 1, COVID Group 2, and control ICU patients was compared using Brown-Forsythe and Welch ANOVA tests with Dunnett T3 test for multiple comparisons. Significant differences in gene set enrichment between COVID-19 patients and non-COVID-19 patient cohorts were calculated using Brown-Forsythe and Welch ANOVA tests with Dunnett T3 test for multiple comparisons.

All statistical analyses were performed in GraphPad Prism (v 9.1.0; www.graphpad.com) and a two-tailed p-value of 0.05 was used for statistical significance.

The RNA-seq data generated in the current study are publicly available through the National Center for Biotechnology Information (NCBI) under SRA project number PRJNA777938. The confirmation study datasets are publicly available under accessions GSE161731 and GSE172114.

To identify differences in immune pathology among severe COVID-19 patients, we analyzed whole blood transcriptomes from 13 individuals with COVID-19-induced AHRF, 6 with other viral and 7 with non-viral-induced AHRF, and 5 controls, who were patients admitted to the ICU on mechanical ventilation but without evidence of AHRF ( Supplementary Table 1 , Supplementary Table 2 ). Out of 18,130 total transcripts, we found 344 differentially expressed genes (DEGs) between control and COVID-19 ICU patients, 248 DEGs between viral AHRF and COVID-19 patients, and 650 DEGs between non-viral AHRF and COVID-19 patients ( Supplementary Table 3 and Figure 1A ). Notably, expression of none of the genes previously associated with severe disease or increased mortality and reported as therapeutic targets for treatment of COVID-19 was significantly changed between controls and COVID-19 patients in the ICU ( Supplementary Figure 1B ) (47–49).

To examine differences in inflammatory pathways between COVID-19 patients and control ICU patients, we carried out Gene Set Variation Analysis (GSVA) (45) using a set of immune cell and pathway gene signatures ( Figure 1B ). Gene expression from COVID-19 patients was enriched for signatures of granulocytes, including inflammatory neutrophils and low-density granulocytes (LDGs) as well as plasma cells (PCs) and CD40 activated B cells. In addition, as compared to controls, COVID-19 patients exhibited decreased enrichment for signatures of dendritic cells (DCs), activated T cells, and Tregs.

Principal component analysis (PCA) using the log₂ transformed gene expression values differentiating COVID-19 and control ICU patients indicated a division among COVID-19 patients into two groups ( Figure 2A ). Notably, the two COVID-19 groups differed in expression of specific COVID-19 associated genes ( Supplementary Figure 1C ). COVID Group 1 patients tended to show an increase in the innate immune checkpoint molecule CD24, whereas COVID Group 2 patients had increased expression of the anti-viral response genes OAS1, OAS2, and OAS3.

We then utilized GSVA to examine inflammatory pathways in the two gene expression-derived COVID-19 patient groups in greater detail ( Figure 2B ). Enrichment of PCs and de-enrichment of DCs was conserved between both COVID-19 groups compared to controls. However, the majority of signatures were differentially enriched in the two groups, revealing distinct immune profiles. Specific granulocyte population signatures were enriched in the COVID-19 patient groups with increased LDGs in COVID Group 1 and increased inflammatory and suppressive neutrophils in COVID Group 2. In addition, COVID Group 1 was uniquely enriched for signatures of CD40 activated B cells, the alternative complement pathway, the cell cycle, glycolysis, and the NFkB complex and de-enriched for activated T cell signatures. In COVID Group 2, natural killer (NK) cell, general interferon (IFN), IFNA2, and IFNB1, but not IFNG signatures were significantly increased, whereas no inflammatory signatures were decreased compared to controls.

We also compared gene signature enrichment from ICU patients with COVID-19-induced AHRF to patients admitted to the ICU with respiratory failure from other non-SARS-CoV-2 viral infections or non-viral causes ( Figure 3A ). The PC gene signature was consistently enriched in both COVID Group 1 and 2 compared with the non-viral and viral AHRF cohorts. Numerous differences in gene set enrichment patterns were noted between the COVID-19 groups and those with other causes of AHRF. Many of the differences in immune cell and pathway enrichment between COVID Group 1 and 2 and non-viral AHRF patients were consistent with the differences from control ICU patients, whereas, in general, viral and COVID-19 AHRF were more similar. In COVID Group 1, CD40 activated B cells and the cell cycle were increased over the non-viral AHRF group. In COVID Group 2, suppressive neutrophils, NK cells, T cells, IFN, IFNA2, and IFNB1 were increased, whereas granulocytes and glycolysis were decreased. as compared to non-viral AHRF. The most consistent difference between COVID Group 1 or COVID Group 2 and viral AHRF patients was the increased PC signature in the COVID patients.

To correlate GSVA gene signature enrichment with AHRF patient cohort and clinical features, we performed multivariable linear regression analysis using MaAsLin 2 and plotted relationships with significant correlations (46) ( Figure 3B ; Supplementary Table 2 ). As a result, we found the TNF, IFNG, inflammatory neutrophil, and suppressive neutrophil signatures had significant positive correlations with length of stay. In addition, the LDG signature had a significant positive correlation and the activated T cell signature a significant negative correlation with length of intubation. However, when linear regression analysis was carried out individually for each cohort, we found that the gene signature to clinical feature correlations were not significant for all patient cohorts ( Figure 3B ). None of these gene signature to clinical feature correlations were significant when considering the control ICU patients or non-viral AHRF patients alone, indicating that these relationships were specific for patients with virus-induced respiratory failure. Among viral and COVID AHRF cohorts, the positive correlation between inflammatory neutrophils and length of stay was significant for both COVID and viral AHRF patients while the correlations between TNF, IFNG, or suppressive neutrophils and length of stay were only significant for COVID but not viral AHRF cohorts. Furthermore, the negative correlation between activated T cells and length of intubation was only significant for the viral AHRF cohort, but not for the COVID groups. Therefore, subsets of ICU patients exhibit differences in immune signatures indicative of a worse clinical prognosis, but this is not unique for COVID-19 patients.

COVID-19 patients, whether Group 1 or Group 2, had a significant increase in PCs over all other ICU patients. This result was further probed by multivariable linear regression analysis, in which the most significant correlation between GSVA enrichment score and patient cohort was for the PC signature, which was uniquely associated with COVID-19 patient groups ( Figure 4A ). To investigate the immunoglobulin (Ig) heavy chain(s) expressed by AHRF COVID-19 patient PCs, we carried out linear regression using PC GSVA scores and Ig heavy chain gene expression ( Figures 4B, C ). As a whole, COVID-19 patient PC GSVA scores were significantly correlated with IGHG3 and IGHA1 Ig heavy chain isotypes ( Figure 4B ). However, when the COVID groups were analyzed separately, only COVID Group 1 showed a significant correlation with expression of IgHA1 (r^{2 =} 0.8, p=0.004). This was not found with COVID group 2 ( Figure 4C ).

To determine whether gene expression-derived groups of hospitalized COVID-19 patients also differed in the level of disease severity, we compared clinical feature and cytokine data of COVID Group 1 and Group 2 patients ( Figure 5 ; Supplementary Table 3 ). We found that baseline demographic data as well as length of stay in ICU and length of intubation were similar between COVID Group 1 and 2 ( Figures 5A, B ). Notably, however, these cohorts varied widely in a number of clinical features indicating that COVID Group 2 had more severe disease ( Figure 5B ). On average, COVID Group 2 had two fewer days of symptoms before admission to the ICU and thus had accelerated disease onset. Upon admission, ferritin and AST levels were over 2X and 1.5X higher, respectively, in Group 2 patients whereas their lung function, as measured by mean PF ratio, was lower. Furthermore, maximum ferritin and aspartate aminotransferase (AST) levels were even more elevated in COVID Group 2 than at admission, indicative of rapid disease progression in these patients. In contrast to clinical features, pro-inflammatory cytokines were only modestly elevated in COVID Group 1 and 2 over controls and in COVID Group 2 over Group 1 ( Figure 5C ). COVID Group 1 and 2 exhibited modest increases in IL6, IL8, and TNF, although these differences did not reach statistical significance. In addition, COVID Group 1 had slightly elevated CD40L and VEGF and COVID Group 2 had significantly elevated levels of the myeloid chemokines CCL2 and CXCL10 as well as IFNA2 and IFNG. In many cases, severe COVID-19 patients are thought to have had greater viral exposure and thus greater viral load in relation to mild cases (50, 51). However, we found no significant differences in viral loads between the COVID-19 patient groups despite their clear differences in clinical features and gene expression profiles ( Figure 5D ).

To examine the persistence of gene signature enrichment over time in the ICU, several COVID-19 AHRF, Viral AHRF, and Non-viral AHRF patients were also sampled at 24- and 72-hours post-admission and individualized trajectories of gene expression were assessed ( Figure 6 and Supplementary Figure 2 ). Gene signatures for all AHRF cohorts remained largely stable over time, but among individuals with COVID-19, Group 1 patients appeared to have greater variation than Group 2 patients. In particular, the IFN and neutrophil signatures that were uniquely enriched in COVID Group 2 were decreased in Group 1 patients over time, whereas the LDG signature was increased. Importantly, Group 2 patient 142, who succumbed, displayed very little change in gene expression over the 72 hours.

Our initial dataset of COVID-19 patients consisted entirely of severe AHRF cases admitted to the ICU. Therefore, we wanted to characterize the immune profiles of COVID-19 patients at different stages of diseases and severity (non-hospitalized vs hospitalized) as compared to healthy controls. To do this, we analyzed a second publicly available COVID-19 transcriptomic dataset (GSE161731, Supplementary Table 1 ), which sampled 16 COVID-19 patients at early-stage (< 10 days), 30 at mid-stage (11-21 days), and 17 at late-stage (> 21 days) disease (52). In total, there were 82 common DEGs between individuals with COVID-19 and healthy controls in this study and our analysis of COVID-19 ICU patients compared to controls in the ICU ( Supplementary Figure 3 ). GSVA analysis revealed that many gene signatures enriched in AHRF COVID-19 patients were selectively enriched in the early and mid-stage, but not late-stage disease cohorts ( Figure 7A ). Furthermore, early-stage patients most resembled the COVID Group 2 cohort, whereas mid-stage disease patients resembled COVID Group 1. Early stage COVID-19 patients were enriched for suppressive neutrophil, monocyte, PC, IFN, CD40 activated B cell, cell cycle, and NFkB gene signatures. Mid-stage patients were enriched for PC, CD40 activated B cell, alternative complement pathway, and cell cycle gene signatures. Late-stage patients were de-enriched for all of these signatures as compared to the early and mid-stage disease cohorts and had no significant differences from healthy controls. Notably, while early and mid-stage cohorts included both mild, non-hospitalized and severe, hospitalized patients, none of the patients with late-stage disease required hospitalization.

In addition to differences in disease stage, patients with early and mid-stage disease were further differentiated by disease severity based on whether they were hospitalized. Comparing 35 non-hospitalized and 11 hospitalized COVID-19 patients ( Figure 7B ) revealed that a few signatures were commonly enriched in all COVID-19 patients regardless of disease severity, including PCs, CD40 activated B cells, the alternative complement pathway, and the cell cycle. Interestingly, linear regression analysis of PC GSVA scores with IgH chain gene expression revealed that both non-hospitalized and hospitalized COVID-19 patients had significant correlations with IgG and IgA chain genes, although the non-hospitalized patient correlations were stronger than those of hospitalized patients ( Supplementary Figure 4 ). The NFkB complex signature was the only one uniquely enriched in non-hospitalized COVID-19 patients as compared to healthy controls. In contrast, many inflammatory signatures were specific to hospitalized COVID-19 patients, including increased granulocyte, inflammatory neutrophil, suppressive neutrophil, LDG, monocyte, IFN, classic complement pathway, and anti-inflammation signatures. In addition, only hospitalized patients had decreased DC and T cell gene signatures. To provide further support for these results, we analyzed a third publicly available dataset (GSE172114) of 23 non-critical and 46 critical COVID-19 patients. Critical COVID-19 patients from this study shared 103 and 2573 DEGs with our previous analyses of COVID-19 ICU patients and COVID-19 patients with varying disease severity respectively ( Supplementary Figure 3 ). Notably, the enriched gene signatures previously noted in hospitalized compared to non-hospitalized COVID-19 patients were also enriched in this group of critical compared with non-critical COVID-19 patients ( Supplementary Figure 5 ). Therefore, severe cases of COVID-19, which require hospitalization, have conserved immune profiles as measured by inflammatory gene signatures, but upon further dissection reveal patient heterogeneity indicative of risk for more severe disease.