Transcriptomic clustering of critically ill COVID-19 patients.

Infections caused by SARS-CoV-2 may cause a severe disease, termed COVID-19, with significant mortality. Host responses to this infection, mainly in terms of systemic inflammation, have emerged as key pathogenetic mechanisms, and their modulation has shown a mortality benefit.In a cohort of 56 critically-ill COVID-19 patients, peripheral blood transcriptomes were obtained at admission in an Intensive Care Unit (ICU) and clustered using an unsupervised algorithm. Differences in gene expression, circulating microRNAs (c-miRNA) and clinical data between clusters were assessed, and circulating cell populations estimated from sequencing data. A transcriptomic signature was defined and applied to an external cohort to validate the findings.We identified two transcriptomic clusters characterised by expression of either interferon-related or immune checkpoint genes, respectively. Steroids have cluster-specific effects, decreasing lymphocyte activation in the former but promoting B-cell activation in the latter. These profiles have different ICU outcome, in spite of no major clinical differences at ICU admission. A transcriptomic signature was used to identify these clusters in two external validation cohorts (with 50 and 60 patients), yielding similar results.These results reveal different underlying pathogenetic mechanisms and illustrate the potential of transcriptomics to identify patient endotypes in severe COVID-19, aimed to ultimately personalise their therapies.

Infections caused by SARS-CoV-2 have a wide range of severity, from asymptomatic to life-threatening cases. The most severe forms of Coronavirus-induced disease (termed COVID-19) [1] lead to respiratory failure fulfilling the acute respiratory distress syndrome (ARDS) criteria [2]. These critically ill patients often require mechanical ventilation and supportive therapy in an intensive care unit (ICU) and show mortality rates that range from 12 to 91% depending on patient and hospital factors [3].

Local and systemic inflammation are key pathogenetic mechanisms in severe COVID-19 [4]. Viral infection triggers a host response that involves not only anti-viral mechanisms, such as release of interferons, but may also activate a systemic, non-specific inflammatory response that has been related to multiple organ failure and death [5]. In addition to standard supportive care, the only treatments that have shown a survival benefit in critically-ill COVID-19 patients aim to modulate this inflammatory response [6]. However, it has been suggested that these treatments do not benefit patients with less severe forms of the disease or with only a mild activation of inflammation [7, 8].

There is increasing evidence that ARDS patients show different clinical features or systemic responses to severe diseases (phenotypes and endotypes respectively) [9]. Although the underlying causes responsible for this heterogeneity are not fully understood, clinical data showing different outcomes in response to a given treatment suggest that pathogenetic mechanisms may be different [10]. Therefore, identification of patient pheno/endotypes may be relevant not only for risk stratification, but also to design specific, personalised therapies in the ICU. Interestingly, whereas clustering of severe COVID-19 patients using respiratory data at ICU admission did not identify different phenotypes [11], addition of circulating biomarkers allowed the translation of the previously identified ARDS phenotypes to COVID-19 showed two groups of patients with different responses to steroid therapy [12], highlighting the relevance of the systemic response in this setting.

Transcriptomic profiling after sequencing of whole blood RNA may be useful to identify groups of critically-ill patients with different underlying pathogenetic mechanisms [13-15]. In addition, microRNAs have been proposed to confer robustness to biological processes by reinforcing transcriptional programs, with important pathophysiological consequences [16]. Preliminary results suggest that circulating micro-RNA (c-miRNA) expression could also play a role in this setting [17]. We hypothesized that clustering of COVID-19 patients using transcriptomics at ICU admission could help to identify subgroups with different pathogenesis. To test this hypothesis, we prospectively sequenced peripheral blood RNA and serum c-miRNA at ICU admission in a cohort of COVID-19 patients, applied an unbiased clustering algorithm and compared gene expression, clinical data and outcomes in the identified subgroups. Finally, we validated our findings in an external cohort.

Methods

Study design

This prospective observational study was reviewed and approved by the regional ethics committee (Comité de Ética de la Investigación Clínica del Principado de Asturias, ref 2020.188). Informed consent was obtained from each patient's next of kin. Fifty-six consecutive patients admitted to one of the participant ICUs at Hospital Universitario Central de Asturias (Oviedo, Spain) from April to December 2020 were included in the study. Inclusion criteria were ICU admission and PCR-confirmed COVID-19. Exclusion criteria were age <18, any condition that could explain the respiratory failure other than COVID-19, do-not-resuscitate orders or terminal status, refusal to participate or severe comorbidities that may alter the systemic response (immunosuppression, history of organ transplantation, disseminated neoplasms). All patients were managed following a standardised written clinical protocol.

Sample acquisition and processing

After inclusion, two samples of peripheral blood were drawn in the first 72 h after ICU admission. One sample was collected in Tempus Blood RNA tubes (Thermo Fisher) to facilitate cell lysis, precipitate RNA and prevent its degradation. The other sample was immediately centrifuged to obtain serum and mixed with TRI reagent for serum RNA precipitation. These tubes were stored at −80°C until processing. Whole blood RNA was extracted by isopropanol precipitation and sequenced in an Ion S5 GeneStudio sequencer using AmpliSeq Transcriptome Human Gene Expression kits that amplify canonical human transcripts (18,574 coding genes and 2,228 non-coding genes with a complete annotation in RefSeq [https://www.ncbi.nlm.nih.gov/refseq/]). Details on RNA extraction and sequencing have been provided elsewhere [8]. FASTQ files containing RNA sequences were pseudoaligned using a reference transcriptome (http://refgenomes.databio.org) and salmon software [18] to obtain transcript counts.

Total serum RNA was extracted using miRNEasy kit (Qiagen), following manufacturer's instructions, and c-miRNA isolated and sequenced at BGI Genomics (Wuhan, China). c-miRNA readouts were mapped using bowtie2 [19], with an index built using the hg38 human reference genome. Quantification of sequenced miRNAs was performed using miRDeep2 [20] with reference human mature and hairpin miRNA sequences downloaded from miRBase (release 22, https://www.mirbase.org).

Clustering

Clustering of RNA samples was performed following a previously described protocol [21]. Briefly, log2-transformed gene expression data (expressed as transcripts per million reads) were filtered to keep the 5% of features with the largest variance. Clusters were built based on Euclidean distances following the Ward clustering algorithm [22] and represented after dimensionality reduction using a uniform manifold approximation and projection (UMAP) algorithm [23]. Cluster p-values, indicating how strong the cluster is supported by the data (this is, the p-value with the alternative hypothesis that the cluster does not exist), were calculated by multiscale bootstrap resampling using the pvclust package [24] for R.

Analysis of differentially expressed genes and circulating microRNA

Gene raw counts obtained after pseudoalignment were compared between clusters using DESeq2 [25]. Log2 fold change for each gene between clusters and the corresponding adjusted p-value (corrected using a false discovery rate of 0.05) were calculated. Genes with an absolute log2 fold change above 2 and an adjusted p-value lower than 0.01 were used for Gene Set Enrichment Analysis (GSEA) using the clusterProfiler R package [26].

A correlation analysis was performed in genes annotated to a Gene Ontology category involved in interferon pathway. Correlation coefficients between each gene pair were transformed to z-scores and the p-values for each comparison calculated using the DGCA package for R [27]. Genes with opposite correlations in each cluster were selected and the networks defined by their significant correlations traced.

Differentially expressed genes between clusters were also matched with the c-miRNAs expressed for each group using the MicroRNA Target Filter tool from Ingenuity Pathway Analysis (Qiagen Digital Insights), to identify predicted interactions. Intersected mRNA and miRNA datasets were filtered to explicitly pair opposed and reciprocal expression changes. Only experimentally observed predictions were considered. Key mRNA-miRNA relationships identified were overlayed onto the networks of interest to explore the predicted functionality in our datasets. Pathways related to humoral, and T and B cellular immune responses were selected as relevant. miRNAs with <3 targeted mRNAs were filtered out from the network.

Changes in gene expression after steroid therapy

To study the effects of steroids in each cluster, peripheral blood gene expression after 4 days in the ICU was assessed in 27 patients (18 assigned to CTP1 and 9 to CTP2), comparing those receiving treatment with dexamethasone (6 mg/12 h) to those who were not treated with this steroid. RNA extraction, sequencing and analysis were performed as described. Pathways with a differential response to steroids were identified after GSEA as those with significant enrichment scores with opposite signs.

Clinical data

Demographics and comorbidities were collected at ICU admission (day 1). Data on gas exchange, respiratory support, hemodynamics, received treatments and results from routine laboratory analyses were prospectively collected at days 1 and 7 after ICU admission. Patients were followed up to ICU discharge. During this period, duration of ventilatory support and vital status were collected for outcome analysis.

Circulating cell populations

Proportions of transcriptionally active circulating cells in each sample were estimated using Immunostates [28], a previously published deconvolution algorithm. From the original reference matrix, cell populations not commonly identified in peripheral blood (Mast cells and macrophages) were removed. Using this modified reference matrix containing expression of 318 genes for 16 different blood cell types, the percentage of each one of these types was estimated from the bulk RNAseq.

Validation

To validate our results in an external cohort, we used two publicly available dataset of 50 and 60 transcriptomes from severe and critically-ill COVID-19 patients [29, 30]. Sample acquisition was performed at enrolment. Clinical data and gene counts were downloaded from Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/, accession number GSE157103) [29] or Zenodo (https://zenodo.org/record/6120249) [30]. First, we identified differentially expressed genes that best discriminate between clusters in our data, as those with an Area Under the Receiver Operating Characteristic curve (AUROC) above 0.95. A transcriptomic score was calculated as the geometric mean of these genes, and the AUROC for this score determined and a threshold between clusters was defined. Finally, raw gene expression data from validation cohorts were normalised using DESeq2, transcriptomic scores calculated, and each sample assigned to one cluster using the previously established threshold, scaled to the range of obtained values (to account for the variability in sequencing techniques). Clinical data, outcomes and estimated cell populations (by bulk RNA-seq deconvolution as previously described) were compared between clusters.

Statistical analysis

Given the observational nature of the study and the lack of previous results, no formal sample size calculations were done. Data are expressed as median and interquartile range. Missing data were not imputed. Differences between clusters were assessed using two-tailed Wilcoxon or chi-square tests (for quantitative and qualitative data respectively). For survival analysis, patients were followed up to ICU discharge, with ICU discharge alive and spontaneously breathing being the main outcome measurement. Differences in this outcome between clusters were assessed using a competing risk model as previously described [8], and hazard ratio for the main outcome, with the corresponding 95% confidence interval, was calculated. All the analyses were performed using R v4.1.1 [31] and packages ggplot2 [32], pROC [33] and survival [34], in addition to those previously cited. All the code and raw data can be found at https://github.com/Crit-Lab/COVID_clustering.

Results

Patient clustering

Peripheral gene expression was sequenced in 56 consecutive critically-ill patients (20% female, age 68 [61-75] year) admitted to one of the participant ICUs. Amongst 16 903 genes counted, 1727 were used for hierarchical clustering (figure 1a). The two main branches of the obtained clustering tree showed the highest p values for an alternative hypothesis that the clusters do not exist (figure 1b and supplementary Figure 1). Therefore, the sample was divided in two mutually exclusive groups, termed COVID-19 transcriptomic profiles (CTP) 1 and 2. Bidimensional representation of the study population using a UMAP algorithm confirmed the separation of the two clusters (figure 1c). Supplementary figure 2 shows a heatmap with the expression of the genes used for clustering.

FIGURE 1

Patient clustering. a) Clustering strategy based on peripheral blood RNAseq, using the 5% genes with the highest variance among samples. b) Hierarchical clustering tree, showing the p-values (corresponding to the alternative hypothesis that the cluster does not exist) of the two main clusters. c) Uniform manifold approximation and projection (UMAP) showing a bidimensional representation of all the samples and clusters. TPM: Transcripts per million reads.

Differences between transcriptomic profiles

Then we assessed the overall differences in gene expression. Using an adjusted p-value cut-off point of 0.01, there were 9700 differentially expressed genes (Supplementary file 1), with 3640 having an absolute log2 fold change above 2 (figure 2a). Interestingly, most of these genes were downregulated in CTP2. Then, GSEA was used to identify the molecular pathways involving these differentially expressed genes. One hundred and ten biological processes with significant differences between clusters were identified (Supplementary Figure 3). Among these, several categories related to the interferon-mediated response and lymphocyte activation were identified (figure 2b), and participating genes were plotted (figure 2c-e). Patients included in CTP1 showed an enrichment of several interferon genes, linked to the activation of a number of immune populations related to innate and adaptative responses (figure 2c), whereas CTP2 was enriched in genes involved in B-cell receptor signaling (figure 2d) and regulatory T-cell differentiation (figure 2e).

FIGURE 2

Differentially expressed genes between COVID-19 transcriptomic clusters (CTP). a) Volcano plot showing fold-change for each gene and their significance level. Genes with an adjusted p-value lower than 0.01 and an absolute log2 fold change above 2 are colored in orange. Differentially expressed genes included in interferon-dependent pathway are labelled. b) Enrichment of Gene Ontology categories related to Interferon signaling in COVID-19 Transcriptomic Profile 2 (CTP2, n=14 compared to CTP1, n=42). C-E: Networks combining pathways and genes with differential expression between clusters; involving Interferon-dependent lymphoid activation (c) upregulated in CTP1, and B-cell receptor signaling (d) and regulatory T-cell differentiation (e), upregulated in CTP2.

In addition to these quantitative changes in expression of interferon-related genes, we explored the existence of qualitative differences between clusters. We calculated the linear correlation coefficients among the 145 genes included in the Gene Ontology categories involving interferon signaling in each cluster. There was a significant difference between the two correlation matrices (figure 3, p<0.001 calculated using a Chi-square test), thus demonstrating differences in the orchestration/structure of IFN responses between groups. In addition, pairwise differences in correlation coefficients for each gene pair were assessed. Gene pairs with correlation coefficients with an adjusted p-value for their difference below 0.05 and opposite signs in each cluster were selected, and networks including these genes traced (figure 3 and Supplementary Figure 4). These results suggest that both clusters have a qualitatively different activation of the interferon pathway, with some genes such as HSP90AB1 and JAK1 acting as hubs with opposite correlations. Of note, CTP1 was hallmarked by strong, positive correlations among effector IFN proteins, whereas this was not the case for CTP2.

FIGURE 3

Correlation between genes included in Interferon-dependent pathways. Correlograms (bottom) and gene networks (top) showing correlations with opposite sign between genes in each COVID-19 transcriptomic cluster (CTP, n=42 and 14 for CTP1 and 2 respectively). Only Pearson correlation coefficients with a p-value lower than 0.05 are shown.

Differences in circulating cell populations

The previous results suggest that the identified clusters may have a different circulating lymphocyte profile. To further explore this finding, cell populations were estimated by deconvolution of RNAseq data. This analysis revealed a higher granulocyte proportion in patients assigned to CTP1, a lower proportion of lymphocytes and no differences in monocytes or NK cells (figure 4a-d). Although no differences in absolute lymphocyte counts were found (645 [483—948] versus 730 [580-908]/mm3, p=0.71, table 1), deconvolution and adjustment by total lymphocyte fraction revealed a higher proportion of CD4+ T cells (figure 4e) and a lower proportion of CD8+ T cells (figure 4f) and naïve B-cells (figure 4g), with no differences in memory B-cells (figure 4h) in this group. Detailed data on other cell populations can be found in the online supplement (supplementary figure 5).

FIGURE 4

Estimated circulating cell populations. a-d) Proportions of blood cells were estimated from RNA-seq using a deconvolution algorithm. e-h) Lymphocyte subpopulations expressed as percentage of the absolute number of lymphocytes. Points represent individual patient data. In boxplots, bold line represents the median, lower and upper hinges correspond to the first and third quartiles (the 25th and 75th percentiles) and upper and lower whiskers extend from the hinge to the largest or smallest value no further than 1.5 times the interquartile range. p-values were calculated using a two-tailed Wilcoxon test.

TABLE 1

Clinical differences between COVID-19 transcriptomic profiles (CTP). BMI: Body mass index. COPD: Chronic Obstructive Pulmonary Disease. APACHE-II: Acute Physiology and Chronic Health disease Classification System II. PBW: Predicted body weight (according to height). PEEP: Positive end-expiratory pressure. IL-6: Interleukin-6. Data are expressed as median (interquartile range) or count (percentage). p-values were calculated using a Wilcoxon test (quantitative data) or Chi-square test (proportions)

	CTP1 (n=42)	CTP2 (n=14)	p_value
Sex			0.174
Male	36 (86%)	9 (64%)
Female	6 (14%)	5 (36%)
Age (years)	69 (63–75)	63.5 (59–69)	0.147
BMI (Kg·m−2)	29 (25–33)	29 (27–31)	0.781
Race			0.582
Caucasian	38 (90%)	14 (100%)
Black	2 (5%)	0
Latino	2 (5%)	0
Chronic kidney disease	4 (10%)	0	0.549
COPD	5 (12%)	1 (7%)	1
Liver cirrhosis	1 (2%)	0	1
Arterial hypertension	26 (62%)	6 (43%)	0.35
Diabetes	9 (21%)	3 (21%)	1
Dyslipemia	18 (43%)	6 (43%)	1
Day 1
APACHE-II score	18 (14–21)	16 (13–17)	0.120
FiO2	0.5 (0.4–0.6)	0.45 (0.3–0.5)	0.438
PaO2/FiO2	197 (157–245)	188 (151–278)	0.863
PaCO2 (mmHg)	43 (39–47)	41 (39–42)	0.091
Respiratory rate (/min)	18 (16–21)	18 (17–22)	0.859
pH	7.37 (7.32–7.41)	7.42 (7.36–7.43)	0.099
Lactate (mEq·L−1)	1.3 (1.08–1.8)	1.1 (0.9–1.2)	0.040
Tidal volume (mL)	479 (455–504)	500 (475–514)	0.499
Tidal volume/PBW (mL·Kg−1)	7.5 (6.9–8.3)	8 (7.5–8.7)	0.239
Plateau pressure (cmH2O)	27 (24–29.75)	25 (22–29)	0.776
PEEP (cmH2O)	14 (12–15)	12 (10–12)	0.088
Driving pressure (cmH2O)	14 (11–15)	15 (12–15)	0.568
Compliance (mL/cmH2O)	36 (31–43)	31 (29–44)	0.697
Creatinin (mg·dL−1)	0.92 (0.68–1.23)	0.71 (0.59–0.97)	0.130
Creatine kinase	96 (60–279)	87 (62–177)	0.446
Lactate dehydrogenase	440 (396–521)	459 (390–493)	0.773
Aspartate aminotransferase	47 (37–73)	45 (31–55)	0.458
Alanine aminotransferase	35 (21–55)	29 (20–58)	0.893
Procalcitonin (ng·mL−1)	0.23 (0.14–0.6)	0.14 (0.13–0.27)	0.250
C Reactive Protein	19 (9–24)	16 (2–25)	0.699
IL-6 (pg·mL−1)	113 (54–276)	164 (36–250)	0.784
Ferritin (ng·mL−1)	1329 (968–1606)	1673 (856–2182)	0.576
D-dimer (ng·mL−1)	1495 (842–3304)	1084 (750–2126)	0.501
Leukocytes (/μL)	9010 (6750–11 825)	5440 (4418–6453)	0.002
Neutrophils (/μL)	7170 (4990–10 190)	4180 (3340–6200)	0.007
Monocytes (/μL)	330 (180–470)	260 (160–480)	0.656
Lymphocytes (/μL)	645 (482.5–948)	730 (580–908)	0.705
Neutrophil-to-Lymphocyte ratio	10.5 (7.8–17.3)	7.1 (3.5–10.7)	0.010
Days from hospital to ICU admission	2 (0–3)	2 (1–4)	0.5
Treatments during ICU stay
Mechanical ventilation	38 (90%)	11 (79%)	0.484
Prone ventilation	23 (61%)	8 (73%)	0.981
Neuromuscular blockade	23 (61%)	6 (55%)	0.643
Extracorporeal membrane oxygenation	1 (3%)	0	1
Vasoactive drugs			0.204
None	17 (40%)	6 (43%)
One	25 (60%)	7 (50%)
Two or more	0	1 (7%)
Steroid therapy	19 (45%)	5 (36%)	0.755
ICU evolution
IL-6 at day 7	54 (11–171)	42 (16–130)	0.713
Ferritin at day 7	1100 (698–1504)	1544 (805–1908)	0.745
D-dimer at day 7	2068 (1249–4586)	1541 (988–3370)	0.422
Ventilator-free days at day 28	12 (0–19)	19 (9–23)	0.050

Potential regulatory miRNAs

To identify c-miRNA potentially related to the observed changes in RNA expression and immune cell populations, we analysed miRNA content using the MiRNA Target Filter tool included in Ingenuity Pathway Analysis. After filtering by experimentally confirmed miRNA-gene relationships, and only opposed changes in miRNA/gene expression levels, 83 miRNAs targeting 608 genes were identified in our dataset of differentially expressed genes. Given the observed differences in lymphocyte populations, we focused on miRNAs involved in humoral and cellular immune regulation (29 miRNAs and 151 genes). Paired miRNA-gene networks are depicted in Supplementary Figure 6 (104 downregulated genes/18 predicted upregulated miRNAs) and figure 5 (47 upregulated genes/11 predicted downregulated miRNAs), with an overlay including differentially expressed genes between CTP1 and CTP2. miRNAs predicted to regulate expression of these genes were identified and compared (figure 5b-h). Among these, counts of miR-145a-5p and miR-181-5p were significatively lower in CTP2 (figure 5c and d respectively).

FIGURE 5

Regulation of gene expression by micro-RNAs. a) Micro-RNAs potentially regulating genes with increased differential expression were identified and a network built. b-h) Counts of hub micro-RNAs (defined as those regulating 3 or more differentially expressed genes) in serum. Points represent individual patient data. In boxplots, bold line represents the median, lower and upper hinges correspond to the first and third quartiles (the 25th and 75th percentiles) and upper and lower whiskers extend from the hinge to the largest or smallest value no further than 1.5 times the interquartile range. p-values were calculated using a two-tailed Wilcoxon test.

Cluster-specific effects of steroids

To assess cluster-specific effects of steroids, we compared gene expression after 4 days of ICU stay in patients with and without steroids in each cluster. Although steroids modified the transcriptomic profile in both clusters, the overlap in differentially expressed genes between clusters was minimal (figure 6a-b). When pathways with divergent responses were assessed (figure 6c), we found that steroids downregulated T- and B- cell activation and IL production and activated JAK-STAT signaling only in patients from the CTP1 cluster. In opposite, steroid therapy was related to B-cell activation in patients assigned to CTP2.

FIGURE 6

Cluster-specific effects of steroids. Differences in peripheral blood gene expression after 4 days in the ICU between patients treated or not with dexamethasone, stratified by cluster. a-b) Euler diagrams showing the number of genes up- (a) and down-regulated in patients receiving steroids. c) Pathways with divergent activation/suppression response to steroids between clusters.

Clinical differences and outcome

Clinical differences between clusters at ICU admission were studied (table 1). There were no significant differences in demographic and clinical variables other than a higher neutrophil count in cluster CTP1, with no differences in lymphocyte counts. Patients assigned to CTP2 cluster showed more ventilator-free days during the first 28 days in ICU (table 1). In the survival analysis, after adjusting for age, sex, and need for intubation during the ICU stay, assignation to CTP2 increased the probability of ICU discharge alive and spontaneously breathing (HR 2.00 [1.08–3.70], p=0.028, figure 7). Other used biomarkers as neutrophil count, neutrophil-to-lymphocyte ratio or C-reactive protein showed only a moderate performance for cluster assignment (AUROCs of 0.74 [0.61–0.88], 0.73 [0.57–0.89] and 0.53 [0.31–0.77], respectively). Replacing cluster assignment with neutrophil-to-lymphocyte ratio or C-reactive protein did not yield a statistically significant HR in the survival analyses (HR 0.958 [0.908–1.011], p=0.122 for neutrophil-to-lymphocyte ratio; HR 0.986 [0.957–1.016], p=0.365 for C-reactive protein).

FIGURE 7

Intensive Care Unit (ICU) stay. Cumulative incidence of the main outcome (ICU discharge alive and spontaneously breathing), modelled using a competing risk model (with death as a competitive risk) and adjusted by age, sex and need for mechanical ventilation during the ICU stay. The hazard ratio (HR) for COVID-19 transcriptomic profile 2 (CTP-2) with its 95% confidence interval is shown.

Definition of a transcriptomic signature and external validation

To apply our findings to an external cohort, we first developed a characteristic gene signature that allows assignation to one cluster using gene expression data. We focused on genes upregulated in CTP2, as they constitute a relatively small group, given the massive gene downregulation in this group. Among these 117 upregulated genes, 15 (BCL2, CARD11, CD247, CD7, CD81, CLSTN1, E2F6, MCM5, PARP1, PNPO, RASGRP1, RCC2, RPTOR, RUNX3 and ZAP70) had an AUROC to identify CTP2 higher than 0.95. Expression of these genes was synthesized into a transcriptomic score. As expected, the score was higher in CTP2 (Supplementary figure 7A), with an AUROC of 0.99 (95% CI 0.97–1) (Supplementary figure 7B). In a Cox-regression analysis including this transcriptomic score, age, sex and need for mechanical ventilation, the score was correlated to ICU discharge (HR 1.202 [1.041–1.387] per 100 points increase in transcriptomic score, p=0.012, supplemental figure 8). Based on these results, a cut-off point of 250 in this score, aimed to include all CTP2 cases, was chosen.

Then, this transcriptomic score was calculated in two external cohorts. Regarding the first cohort (n=50), 13 patients were classified as CTP1 and 37 as CTP2. Comparisons between these clusters are shown in table 2. In spite of no significant differences in age, sex, APACHE-II or SOFA scores, patients assigned to CTP2 showed more ventilator-free days at day 28 of ICU stay, and the percentage of patients with zero ventilator-free days at day 28 was lower in CTP2. In the second validation cohort (n=60), 22 patients were classified as CTP1 and 38 as CTP2, after rescaling the cut-off point to account for differences in sequencing technology and depth. Resembling the previous results, there were no differences in age and sex, but mortality by day 28 was higher in patients assigned to CTP1 (table 2). Deconvolution of peripheral blood transcriptomes in both validation cohorts recapitulated some of the differences observed in the discovery cohort, including higher neutrophil counts and lower proportions of CD8+ T-cells in CTP1 (Supplementary Figures 9 and 10).

TABLE 2

Clinical data and outcomes in the validation cohort. APACHE-II: Acute Physiology and Chronic Health disease Classification System II. SOFA: Sequential Organ Failure Assessment. VFD: Ventilator-free days. p-values were calculated using a Wilcoxon test (quantitative data) or Chi-square test (proportions)

	CTP1	CTP2	p_value
Validation cohort 1
Sample size	13	37
Transcriptomic score	216 (197–228)	365 (311–470)
Age (years)	63 (55–73)	64 (55–72)	0.842
Sex (male/female)	8/5	25/12	0.741
APACHE-II score	23 (20–34)	21 (14–25)	0.097
SOFA score	7 (6–13)	8 (6–10)	0.35
VFDs at day 28	0 (0–20)	18 (2–28)	0.016
Zero VFDs at day 28	8 (62%)	8 (22%)	0.014
Validation cohort 2
Sample size	22	38
Transcriptomic score	1430 (1215–1506)	2194 (1782–2503)
Age≥50 years	18 (	31	1
Sex (male/female)	13/9	21/17	0.986
Death at day 28	7	3	0.042

Discussion

Our results show that unsupervised transcriptomic clustering of critically ill COVID-19 patients at ICU admission, results in two groups with different immune profiles, response to steroids and outcome. Application of a cluster-specific score to two independent cohorts confirmed this result. These findings suggest there are specific COVID-19 endotypes with different underlying immunopathogenesis and outcomes.

Clustering strategies have been proposed to identify different subgroups of critically ill patients with respiratory failure that may help to personalise treatments. In ARDS, a hyperinflammatory/reactive phenotype [9, 35], characterised by markers of acute inflammation and tissue hypoxia, has been linked to higher mortality rates and could specifically benefit from fluid restriction, higher PEEP levels or protective ventilation [36], in contraposition to the uninflamed phenotype. Of note, causes of ARDS were different between phenotypes, with a higher incidence of sepsis in the hyperinflamed/reactive group. Clustering of COVID-19 patients using respiratory data failed to identify phenotypes at ICU admission [11]. Addition of clinically available biomarkers allowed a direct translation of the inflammatory/reactive framework in two cohorts [12, 37]. Focusing on a single disease (COVID-19) rather than a syndrome (ARDS) could result in a reduced phenotypic variability, thus increasing the informative value of the systemic response evaluated by circulating biomarkers. Clinically available markers as neutrophil-to-lymphocyte ratio or C-reactive protein are usually elevated in severe forms of COVID-19 [38], but their role to stratify patients is yet to be determined and, in our study, failed to predict outcome or cluster assignment.

In this setting, transcriptomic clustering may offer several advantages by including a large number of features for classification, reduced intervention times and absence of imputed or not available data, although the superiority of this approach remains to be demonstrated. Increasing evidence points at c-miRNA as biomarkers with pathogenetic implications given their role as modulators of gene expression. Point-of-care devices under development would allow the quantification of our 15-gene signature or validated cluster-specific c-miRNA at the bedside, to rapidly identify these patient endotypes and predict outcomes [39].

Bulk peripheral blood RNAseq has been used to study COVID-19 pathogenesis, by comparing cases with different severity or against healthy controls [40-42]. Our approach included only severe cases, revealing two different clusters that include quantitative and qualitative differences in the regulation of the immune response to SARS-CoV-2 infection and different responses to steroids. Of note, these biological disparities occur despite no differences in clinical variables, suggesting that clusters reflect endotypes with specific pathogenetic mechanisms and may outperform clinical diagnostic instruments.

CTP1 is characterised by an interferon-driven response and CD4+ T lymphocyte activation, that have been linked to a worse outcome [43, 44]. miR-145a-5p and miR-181-5p, which play key roles promoting granulopoiesis [45] and CD4+ T-cell maturation [16, 46] respectively, were upregulated in this cluster. Steroids downregulated genes involved in lymphocyte activation, but upregulated the JAK/STAT pathway in CTP1, which can promote further overexpression of miR-181 family members [47, 48]. JAK/STAT pathway can be activated by IL-6 and has been related to mortality in COVID-19 patients [49, 50].

CTP2 cluster, with better outcome, is characterised by B-cell and Treg activation and upregulation of immune checkpoints such as BCL2 and immunoglobulin and TNF superfamilies. Of note, BCL2 and TNF superfamilies are targeted by miR-181, which is decreased in this cluster and has been described to induce immunoparalysis and block immune checkpoints [51]. Dysregulation of other immune checkpoints has been also linked to mortality in COVID-19 patients [52]. In this group, steroids further promoted B-cell activation. Collectively, these results raise the hypothesis that steroids may help to further regulate the inflammatory response in CTP2, but activate JAK/STAT- dependent immune response in CTP1, which in turn could partly explain our differences in ICU outcome and the proposed synergic effects of steroids and IL-6 / JAK blockade in COVID-19 [53].

Our results have several limitations. First, the sample size is reduced, so we cannot exclude the existence of additional clusters with other underlying pathogenetic mechanisms, or that different clustering parameters or strategies may yield different results. However, unbiased p-values associated to the identified clusters were high, and the results confirmed in two independent validation cohorts. It must be noted that time of sampling differs among cohorts (first 72 h after ICU admission in our study and in the COMBAT cohort [30], between days 1 and 6 in the study by Overmyer et al. [29]). We do not have data to define specific time windows. However, the consistency of the results among studies reinforces the external validity of our clustering strategy. Third, cell populations were estimated by deconvolution of the bulk transcriptome and should be confirmed using single cell RNA seq or flow cytometry. Finally, although our data show different effects of steroids in each cluster, it is unclear if therapeutic immunomodulation may impact outcomes in a cluster-specific manner.

In summary, our results show that transcriptomic clustering using peripheral blood RNA at ICU admission allows the identification of two groups of critically-ill COVID-19 with different immune profile and outcome. These findings could be useful for risk stratification of these patients and help to identify specific profiles that could benefit from personalised treatments aimed to modulate the inflammatory response or its consequences.