Time-resolved systems immunology reveals a late juncture linked to fatal COVID-19

Abstract

COVID-19 exhibits extensive patient-to-patient heterogeneity. To link immune response variation to disease severity and outcome over time, we longitudinally assessed circulating proteins as well as 188 surface protein markers, transcriptome, and T cell receptor sequence simultaneously in single peripheral immune cells from COVID-19 patients. Conditional-independence network analysis revealed primary correlates of disease severity, including gene expression signatures of apoptosis in plasmacytoid dendritic cells and attenuated inflammation but increased fatty acid metabolism in CD56^dimCD16^hi NK cells linked positively to circulating interleukin (IL)-15. CD8⁺ T cell activation was apparent without signs of exhaustion. Although cellular inflammation was depressed in severe patients early after hospitalization, it became elevated by days 17-23 post symptom onset, suggestive of a late wave of inflammatory responses. Furthermore, circulating protein trajectories at this time were divergent between and predictive of recovery versus fatal outcomes. Our findings stress the importance of timing in the analysis, clinical monitoring, and therapeutic intervention of COVID-19.

Keywords: CITE-seq; COVID-19; exhaustion; immune juncture; mixed-effect statistical modeling; single cell multi-omics; time-resolved.

Graphical abstract

Figure 1

Study design and a multi-parameter disease severity metric (DSM) (A) COVID-19 patient cohort overview and sample collection timeline. The distribution of sample collection timing since symptom onset (TSO) is summarized on the top while the sampling for each of the 33 donors (age and gender indicated in parentheses) is plotted below indicating individual T0 (1^st) up to T3 (4^th) time points. One patient (P098) with unknown TSO is depicted here by using his hospital admission date as the reference. Some samples were excluded from analysis after initial quality control as they did not have a sufficient number of cells (see STAR methods). (B) Experimental design, analysis questions, and approach. Frozen PBMCs from the collected samples indicated in (A) and age-matched healthy controls were thawed and pooled, followed by combined surface protein/mRNA single-cell expression analysis (CITE-seq); data from individual samples were demultiplexed by a combination of SNP-based and “hashtag” antibody staining (see STAR methods). After automated clustering (plus additional targeted cell populations identified by manual gating) and cell subset identification based on surface protein profiles, pooled RNA data (“pseudobulk”) libraries were generated in silico for each sample per cell cluster/population and then used for downstream analysis using mixed effect statistical modeling to assess the effect of disease severity and TSO, while controlling for other variables such as TSO (when characterizing the effect of disease severity), age, and experimental batch. (C) Heatmap of 18 clinical and serum protein measurements of patients after correction for days since hospital admission. Measurement values are standardized across subjects. Samples are divided into four groups based on unsupervised hierarchical clustering. Patients admitted to the intensive care unit (ICU) during their hospitalization are denoted with an asterisk (^∗) next to their labels. None of the T0 samples were collected at the ICU. (D) Parameter importance of the fitted coefficient values from partial least squares (PLS) ordinal regression models of disease severity. Parameters are ordered by their absolute median coefficient values from independent training iterations in leave-one-out cross validation. Error bars indicate SD of coefficient distribution across all iterations. CXCL9/MIG, monokine induced by γ interferon; IP-10, interferon-γ-inducible protein 10. (E) Distribution of patient disease severity metric (DSM) groups based on the disease severity-outcome classification for the 30 patients with CITE-seq data (left panel, see STAR methods) and an independent set of 64 patients from Brescia (right panel). The DSM formula was trained using the CITE-seq cohort only. Significance of group difference is determined by two-tailed Wilcoxon test. ^∗p ≤ 0.05, ^∗∗p ≤ 0.01, and ^∗∗∗p ≤ 0.001. See also Figure S1 and Table S1.

Figure S1

Timeline of treatments relative to hospital admission date and DSM validation in an independent set of patients, related to Figure 1 (A) Timeline of treatments relative to hospital admission date. Dashed lines represent treatments with unknown start and/or end dates. Age and gender are listed next to each patient label. (B) Distribution of patient disease severity metric (DSM) grouped based on the WHO ordinal disease severity score of the CITE-seq cohort at the earliest time of PBMC sampling. The p value was computed using the Jonckheere trend test to assess the positive correlation between the WHO score and DSM. (C) Distribution of validation cohort disease severity metric (DSM) grouped by whether they were ever admitted to the ICU during the course of hospitalization for all patients (left) and only patients classified as Critical-Alive (right). Significance of difference between groups is determined by two-tailed Wilcoxon test.

Figure 2

A multimodal single peripheral immune cell atlas of COVID-19 (A) Average DSB-normalized protein expression (see STAR methods) in each coarse-level cell cluster. Only proteins with a DSB value >3 in at least 1% of the cells are shown. (B) UMAP visualization of single cells based on protein expression profiles for innate and adaptive groupings of cells labeled by the name of the corresponding coarse-level cluster. (C) Average expression of selected surface protein markers in example finer-resolution CD4⁺ T cell clusters (CM, central memory; TM, transitional memory; EM.TE, terminal/effector memory). (D–I) Transcript-based UMAP visualization of classical monocytes defined by surface proteins. Cells are colored according to donor class (D), disease severity group (only COVID-19 samples are shown) (E), days since symptom onset (only COVID-19 samples are shown) (F), surface marker CD163 expression (G), surface HLA-DR expression (H), and donor (I). See also Figure S2 and Tables S2 and S3.

Figure S2

Single immune cell atlas of COVID-19 reveals cell populations associated with COVID-19 disease status and severity, related to Figure 2 (A) Correlations of cell frequencies gated from CITE-seq and independently obtained 27-color flow cytometry data of the same samples. Shown are Pearson correlations and p values. Note that the cTfhs shown here were gated using the same strategy as flow data, i.e., based on CXCR5^hiCD45RA^-CD4⁺ T cells. Given the less-than-ideal staining of CXCR5 in both CITE-seq and flow cytometry, the cTfh being analyzed elsewhere in the paper were gated based on ICOS and PD1 in ICOS^hiPD1^hi CD4⁺ T cell RNA-based cluster (see STAR methods). (B) Frequencies of immune cell clusters/subsets in HC, DSM-low (less severe disease; DSM at or below median of DSM) and DSM-high (more severe disease; DSM above median) groups at T0 (near hospitalization). Shown are the populations with the most significant changes in either COVID-19 versus HCs or correlation with DSM. P values shown are unadjusted p values based on linear models controlling for age, experimental batch and TSO (TSO only controlled for in DSM model, see STAR methods). P values are reported for assessing differences between COVID-19 versus HCs (green) and for correlation with DSM (red). Additional cell clusters meeting the p value cutoff are shown in (B). (CD4_Mem_CM: central memory CD4⁺ T cell cluster as a fraction of total CD4⁺ T cells; CD8_Mem_CM.TM: central/transitional memory CD8⁺ T cell cluster as a fraction of total CD8⁺ T cells; the frequency of γδT, mucosal associated invariant T cells (MAIT), cDC, pDC, Non-Classical Monocytes and Platelets are expressed as fractions of total PBMCs.) 104 populations were tested in total. (C) Heatmap showing cell frequencies of major cell clusters/subsets in individual subjects (columns), grouped by HC and DSM. Plasmablast and B cell subsets are expressed as fractions of CD19⁺ B cells; CD4⁺ T cell subsets are expressed as fractions of CD4⁺ T cells; CD8⁺ T cell subsets are expressed as fractions of CD8⁺ T cells; the CD163^hi classical monocyte cluster is expressed as a fraction of total classical monocytes; others are shown as fractions of total PBMCs. Rows in the heatmap were scaled to have mean = 0 and variance = 1.

Figure 3

Cell-type-specific gene expression signatures of COVID-19 disease status and severity (A) Gene set enrichment analysis (GSEA) result of COVID-19 versus HCs at T0. Selected gene sets (rows, see STAR methods) are grouped into functional/pathway categories. Differential expression moderated t statistics were computed from a linear model accounting for age, and experimental batch (TSO not accounted for as TSO does not exist for healthy individuals). Dot color denotes normalized gene set enrichment score (NES) and size denotes -log10(adjusted p value). p values were adjusted using the Benjamini-Hochberg method. The sign of the NES corresponds to increase/decrease of change in COVID-19 compared to HCs. Cell populations (columns) are grouped by lineage/type. See Table S4A for detailed results of all gene sets which passed the adjusted p value cutoff of 0.2. (B) Similar to (A) but the enrichment analysis was performed based on association with DSM at T0, and the model controlled for TSO. The sign of the NES score denotes the sign of the association with DSM. See also Table S4B. (C) Heatmap of type I IFN response in classical monocytes. Heatmap showing the scaled average mRNA expression (within the indicated cell cluster/population for a given subject) of shared leading-edge (LE) genes from the GSEA analysis of COVID-19 versus HCs (see A) and association with DSM (see B). Selected top LE genes are labeled by gene symbol and only the sample from the first time point for each patient (T0) is shown. For gene expression heatmaps, subjects are grouped by HC and DSM classes; also shown are the DSM value and severity-outcome classification of each patient (top of the heatmaps). Subjects with a small number of cells not included in the GSEA test are marked by a # (cell number <8). (D) Per-sample gene set signature scores of the GO_Response to type I IFN gene set versus TSO (in days) in DSM-low and DS-high groups. Classical monocyte is shown as an example. Individual samples are shown as dots and longitudinal samples from the same individual are connected by gray lines. Trajectories (using LOESS smoothing, see STAR methods) were fitted to the groups separately with the shaded areas representing standard error. The median score of age- and gender-matched HCs is shown as a green dotted line. Gene set scores were calculated using the union of LE genes from both the timing (TSO) and disease severity (DSM) association models using gene set variation analysis (see STAR methods). (E) Scatterplot comparing the effect size of association between TSO and the GO_Response to type I IFN gene set in the DSM-high (y axis) and DSM-low (x axis) patients. Each dot corresponds to a cell type/cluster. The effect size reflects the change of type I IFN signature enrichment with time (effect size <0 corresponds to the enrichment decreasing with time). Cell types are colored by their statistical significance (adjusted p value <0.05) of whether the association with time is significantly different between the two DSM groups (based on a model with a DSM-TSO interaction term, see STAR methods). The area framed by blue borders corresponds to cell types with more precipitous drop of IFN signature over time in the DSM-low than DSM-high groups. (F) Heatmap of apoptosis/cell death signature in pDCs. Similar to (C), but instead of the shared LE genes, all LE genes from COVID-19 versus HCs and association with DSM are shown. Shared LE genes and genes in the REACTOME_Oxidative stress-induced senescence gene set are annotated on the right. (G) Heatmap showing the sample-level pairwise Pearson correlations among serum IFN-α2a level, pDC frequency, the apoptosis signature score in pDCs, as well as the IFN-I and protein translation signature scores in classical monocyte, CD56^dimCD16^hi NK, and CD8 memory T cell clusters (^∗p value <0.05). Selected scatterplots are shown, and the corresponding entry in the heatmap is indicated using bold borders. Samples were filtered to remove those with fewer than 8 cells or less than 40,000 unique molecular identifiers (UMI) in the pseudobulk library. Note that not all subjects had IFN-α2a measurements. Each dot indicates a subject and is colored by the severity-outcome category. See also Figure S3 and Table S4.

Figure S3

Cell-type-specific gene expression signatures association with time since symptom onset and disease severity, related to Figure 3 (A and B) Similar to Figure 3B, but here showing GSEA results (of select gene sets) based on association with time since symptom onset (TSO) in DSM-low (A) and DSM-high (B) groups using all time points. See Tables S4C and S4D for detailed results of these selected gene sets and all sets that passed the adjusted P value cutoff of 0.2. Dot color denotes NES and size denotes -log10(adjusted P value). The sign of NES corresponds to the sign of correlation with TSO. Cell clusters/populations are grouped by lineage/type. (C) Similar to Figure 3C. Heatmap of translation/ribosomal gene signature in classical monocytes at T0. (D) Similar to Figure 3D. Time course of gene set signature scores of REACTOME_Translation and KEGG_Ribosome gene sets in DSM-low and DSM-high groups, respectively. Classical monocyte is shown as an example. (E and F) Similar to Figure 3F. Heatmap of apoptosis/cell death gene signature in pDCs of validation cohorts (Schulte-Schrepping et al., 2020) cohort 1 (E) and cohort 2 (F). The LE genes from the GSEA analysis of COVID-19 versus HCs and association with DSM in our Brescia CITE-seq cohort are shown (i.e., same genes as Figure 3F). Only the first time point (T0 sample) of each subject is shown. (G) GSEA results of Schulte-Schrepping et al. (2020) cohorts for pDC apoptosis/cell death signature identified in Brescia cohort. Given the limited pDC cell numbers in most samples, particularly in cohort 1, the enrichment was calculated without filtering based on cell number or library size cutoffs.

Figure 4

Conditional independence network analysis points to IL-15-associated fatty acid metabolism and attenuated inflammation in CD56^dimCD16^hi NK cells as primary correlates of disease severity (A) Disease severity network showing cell-type-specific gene set signatures directly connected with DSM (see STAR methods). Direct connections between the nodes to each other are also shown in a lighter shade. Edge value and width indicate Spearman correlation and its statistical significance, respectively, between DSM and the gene set signature score. In the legend: % leading-edge (LE) genes denotes the proportion of genes from the gene set that are in the LE based on gene set enrichment analysis. The top three LE genes based on effect size (association with DSM; see STAR methods) and other selected genes representative of the biology are shown for each gene set. (B, C, G, and I) Scatterplots showing the correlation of circulating IL-15 level versus REACTOME_Fatty acid metabolism signature score (B) and DSM (C); REACTOME_Fatty acid metabolism score versus HALLMARK_TNFα signaling via NF-κB score and GO_Cellular response to IL-1 score (G); REACTOME_Fatty acid metabolism score versus HALLMARK_mTORC1 signaling score and normalized IFNG mRNA expression (I). Pearson correlation (R) and associated p value were computed using pairwise complete observations (see STAR methods). Each dot indicates a subject and is colored by the severity-outcome category. (D) Heatmap of REACTOME_Fatty acid metabolism LE genes from the GSEA analysis of DSM association in CD56^dimCD16^hi NK cells at T0 (see Figure 3B). LE genes are grouped based on annotations obtained from Gene Ontology. Genes with an asterisk (^∗) are those with both biosynthetic and catabolic annotations. Genes in red are those in the fatty acid oxidation set. (E and F) Similar to (D), heatmaps of inflammation related gene sets in CD56^dimCD16^hi NK cells at T0: HALLMARK_TNFα signaling via NF-κB (E), GO_Cellular response to IL1 and KEGG_Chemokine signaling pathway (see A) (F). Heatmaps showing the scaled average mRNA expression (within the indicated cell cluster/population for a given subject) of LE genes from the GSEA analysis of DSM association in CD56^dimCD16^hi NK cells (see Figure 3B). Selected top LE genes are labeled by gene symbol. For all gene expression heatmaps (D–F), subjects are grouped by HC and DSM classes; also shown are the DSM value and severity-outcome classification of each patient (top of the heatmaps in D). Subjects with a small number of cells not included in the GSEA test are marked by a # (cell number <8). (H) Average IFNG mRNA expression of CD56^dimCD16^hi NK cells; (D), (E), (F), and (H) are aligned column wise. (J and K) Per-sample gene set signature scores of REACTOME_Fatty acid metabolism (J), HALLMARK_TNF-α signaling via NF-κB and GO_Cellular response to IL-1 (K) versus TSO (in days) in DSM-high (red) and DSM-low (blue) groups of CD56^dimCD16^hi NK cells. Individual samples are shown as dots and longitudinal samples from the same individual are connected by gray lines. Trajectories (using LOESS smoothing) were fitted to the groups separately with the shaded areas representing standard error. The median score of age- and gender-matched HCs is shown as a green dotted line. Gene set scores were calculated using the union of LE genes from both TSO association and DSM association models reflecting the change over time and severity differences at T0 (see STAR methods). (L) Circulating IL-15 levels in DSM-low and DSM-high groups versus TSO. See also Figures S4 and S5 and Table S5.

Figure S4

Supporting data for dissecting primary gene expression signature correlates inferred by conditional independence network analysis, related to Figure 4 (A) Scatterplot of REACTOME_Oxidative stress-induced senescence signature score and GO_Apoptotic signaling signature score in pDCs (Pearson correlation and p value were computed using pairwise complete observations [samples filtered to those with at least 8 cells and 40,000 UMI detected in the pseudobulk library]). Each dot denotes a subject and is colored by the severity-outcome category. (B) Similar to (A), but between circulating IL-15 level and fatty acid metabolism signature score in CD56^dimCD16^hi NK cells after regressing out their associations with DSM. (C and D) Similar to Figure 4D. Heatmaps of REACTOME_Fatty acid metabolism in NK cells of two validation cohorts (Schulte-Schrepping et al., 2020) cohort 1 (C) and cohort 2 (D). The LE genes from the GSEA analysis of association with DSM in each cohort are shown. Only the first time point (T0 sample) of each subject is shown. (E) GSEA results of Schulte-Schrepping et al. (2020) cohorts for NK cell REACTOME_Fatty acid metabolism. (F) Similar to Figure 4G. Scatterplot of REACTOME_Fatty acid metabolism score and HALLMARK_TNFα signaling via NF-κB score in the validation cohorts (Schulte-Schrepping et al., 2020). (G and H) Similar to Figure 4E. Heatmaps of inflammation related gene sets in classical monocytes: HALLMARK_TNFα signaling via NF-κB (G) and HALLMARK_Inflammatory response (H).

Figure S5

Exogeneous corticosteroid treatment is not a major driver of cell-type-specific gene expression signatures associated with disease severity, related to Figure 4 (A) GSEA results for glucocorticoid response signature (see STAR methods on how the signature was derived) in DSM association model. Positive enrichment corresponds to higher level in the DSM-high samples. (B and C) Scatterplot showing the correlations between the indicated signature scores (computed using GSVA) and the glucocorticoid response signature score (B) or the TSC22D3 mRNA expression level (C) in CD56^dimCD16^hi NK cells. Pearson correlations and p values are reported for patients with corticosteroid use (orange; steroid use TRUE) and those without (cyan; steroid use FALSE) as well as for all samples (black). Each dot indicates a subject, shaped by steroid use status and colored by the severity-outcome category. (D) TSC22D3 mRNA expression levels of CD56^dimCD16^hi NK cells and classical monocytes in HC, no steroid-use and steroid-use COVID-19 patient groups. P values from testing for differences between the no steroid-use and steroid-use groups are shown. P values were obtained from an ANOVA test accounting for DSM, TSO, age and experimental batch.

Figure S6

Single-cell and clonal expansion analysis in CD8⁺ T cells and exhaustion assessment in clonal CD8⁺ T cells, related to Results section “Extensive single-cell and clonal expansion heterogeneity without signs of exhaustion in CD8+ T cells (A) Heatmap showing the gene expression profile of CD8⁺ T cell clusters identified based on single-cell mRNA expression of the leading-edge genes of selected pathways presented in Figure 3 (see text and STAR methods); only COVID-19 T0 samples are shown. Average scaled mRNA expression per cluster of genes (columns) in each of the CD8⁺ T cell clusters (rows) is shown. Gene memberships in selected gene sets are indicated by the color bars at the bottom of the heatmap. Total number of cells in each cluster is indicated on the right side of the heatmap. Bar plot shows the fraction of cells within each cluster that are defined as expanded (> 1 cell detected per CDR3 α-β pair). Two small clusters (< 32 cells each) are not shown. (B) Average expression of selected surface proteins in the clusters from (A). (C) Results of linear model accounting for age, and experimental batch, comparing the frequency of CD8⁺ T cell clusters from (A) between COVID-19 and HC samples. Shown is the difference in mean frequency between COVID-19 and HC (x axis) versus the -log10 (p value); horizontal line denotes an unadjusted p value of 0.05. 18 CD8 clusters were tested in linear model. (D) Fraction of overlap between RNA based clustering (from A) and surface protein based CD8⁺ naive and memory T cell cluster annotations (based on high resolution clustering). (E) CD8⁺ T cell clonality (median rarefied Simpson index, see STAR methods) in HC, DSM-low and -high groups. p values shown are from linear model adjusted for experimental batch, age, and TSO (green for assessing differences between COVID-19 versus HCs and red for correlation with DSM; TSO only adjusted for in the model assessing DSM association). (F) Coefficients from linear models of mean surface protein expression of canonical exhaustion markers fitted to COVID-19 patients and HCs. Positive coefficients (red bars) correspond to higher expression in COVID-19 versus HCs (above) or higher expression in DSM-high versus -low (below) and vice versa for negative coefficients (green bars). p values are indicated with red indicating significance at the 0.05 level (unadjusted). Among them, only CTLA4 was mildly significant but the association was in the opposite direction expected for exhaustion: its expression was lower in patients, particularly in the more severe. Similarly, insignificant changes or inconsistent directions of change were detected for the mRNA of these protein markers (data not shown). (G) Association of proportion of CD39⁺PD1⁺ cells with COVID-19 versus HCs and DSM in clonal CD8⁺ memory T cells using different cutoffs for CD39 and PD1 surface protein expression DSB counts (0.5, 1). p values are from Wilcoxon tests. Similarly, we also tested whether the frequency of cells co-expressing multiple exhaustion markers in (F) was different. Independent of the surface marker combination or protein expression cutoff used, we saw no signs of exhaustion or association with DSM (data not shown). (H) GSEA results for assessing enrichment for known exhaustion signatures in DE genes for expanded CD8⁺ T cells in COVID-19 versus HCs and DSM-high versus DSM-low comparisons. Positive enrichment corresponds to higher level in the first group. Both exhaustion gene sets were not enriched in the DSM-high versus DSM-low comparison, indicating that exhaustion was not associated with disease severity. The enrichment in the COVID versus HC comparisons could reflect cellular activation (“up” signature – see (I) below) and depression of translation genes in COVID-19 patients negatively associated with IFN-I signatures (“down” signature) – see also Figure 3A. (I) Gene set enrichment of Wherry et al. (2007) up and down genes in KEGG, GO BP, REACTOME, and Li BTM’s. Only top hits (ranked by adjusted p values from Fisher’s exact test) are shown. The “exhaustion up” genes appeared to reflect CD8⁺ T cell activation as they are enriched for NF-κB, JAK-STAT signaling, and IL-2 signaling genes. (J) GSEA result of Schulte-Schrepping et al. (2020) cohorts for exhaustion signatures of COVID-19 versus HCs and severe versus mild comparisons at T0. Differential expression moderated t statistics were computed from a linear model accounting for age, experimental batch and TSO (TSO only accounted for in the severe versus mild comparison). Dot color denotes normalized gene set enrichment score (NES) and size denotes -log10(adjusted P value). Dot shape indicates the significance of adjusted P values (adjusted P value < 0.05). Both cohorts 1 and 2 and both CD8⁺ T cell clusters in the original data and memory CD8⁺ T cell populations using transferred cell type labels from our CITE-seq data are shown. (K) Similar to (H). GSEA results of Schulte-Schrepping et al. (2020) cohorts for Wherry et al. (2007) exhaustion signatures of COVID-19 versus HCs and severe versus mild comparisons at T0. Positive enrichment corresponds to higher level in the first group. Only memory CD8 population based on transferred labels from our CITE-seq data is shown as an example.

Figure 5

Analyses of timing effects suggest a late immune response juncture (A) Time course of monocyte subset frequencies in DSM-low and DSM-high groups. Classical monocyte is expressed as fraction of total PBMCs; the CD163^hi classical monocyte subset is expressed as a fraction of classical monocytes. The median cell frequency of age- and gender-matched HCs is shown as a green dotted line. Individual samples are shown as dots and longitudinal samples from the same individual are connected by gray lines. P values shown are from mixed effect linear models indicating the statistical significance of the timing effect (i.e., TSO) accounting for age and experimental batch in DSM-low and DSM-high groups, respectively. Trajectories (using LOESS smoothing) were fitted to DSM-low and DSM-high groups separately. TSO = days 17–23 period is highlighted with purple. (B) Similar to (A) but showing the absolute blood neutrophil and monocyte counts. The two green dotted lines mark the approximate reference range of cell counts in healthy adults. The shaded areas around trajectories denote standard error. (C) Effect size (normalized enrichment score from GSEA) comparison of the period before day 17 (TSO dimCD16^hi NK cells and classical monocytes. The gene set signature scores were calculated using the LE genes identified from the enrichment analysis shown in (C) above to highlight differences between the DSM-high versus DSM-low groups during the days 17–23 period. The median score of age- and gender-matched HCs is shown as a green dotted line. The shaded areas around trajectories denote standard error. (F) Time course of serum protein levels from DSM-low and DSM-high patients, respectively. Similar to (E). Top 4 serum proteins significantly different between DSM-high versus DSM-low groups during the days 17–23 period are shown. See Table S6 for the full list of differentially expressed proteins. See also Figure S7 and Tables S4 and S6.

Figure S7

Supporting data for critical juncture analysis, related to Figures 5 and 6 (A) Similar to Figure 3A, but here showing GSEA results for assessing the differences of delta between disease severity groups (DSM-high versus DSM-low) between the days 17-23 time window and the period before (TSO < day 17). See Table S4F for detailed results of these selected gene sets and the entire set that passed the adjusted p value cutoff of 0.2. (B) Time course of blood neutrophil and monocyte counts in recovered and deceased groups. Corresponding range of HCs are shown with green dotted lines. Longitudinal samples from the same individual are connected by gray lines. Trajectories were fitted to recovered and deceased groups separately with the shaded areas representing standard error. (C) Effect size comparison of DSM-high versus DSM-low (CITE-seq cohort) and deceased versus recovered (critical patients with distinct outcome subcohorts – see Figure 6A) at the days 17-23 period. Effect sizes were derived from the coefficients of the group comparison term in mixed-effect models. Proteins are colored by whether they were significantly different (with an unadjusted p value of less than or equal to 0.05) in the deceased versus recovered (“outcome”) and the DSM-high versus -low comparisons (“DSM”). See also Table S7 for the full list of DE proteins. (D) Similar to (C). Effect size comparison of Brescia deceased versus recovered and an independent US cohort (Yale cohort) (Lucas et al., 2020) deceased versus recovered patients (See STAR methods) for 38 overlapping circulating proteins/cytokines at the juncture period (TSO days 17-23). Proteins in red are those significantly different between the deceased and recovered patient groups in the Brescia cohort (unadjusted ANOVA p value < 0.05).

Figure 6

Divergence of deceased and recovered patients at the late juncture (A) Approach for assessing and validating the late immune response juncture hypothesis by using serum protein profiles of critical ill patients with either recovery or deceased outcomes. (B) Effect size plots of circulating serum proteins comparing the difference between critical deceased versus recovered patients before (days 7–16), during (days 17–23), and after (days 24–30) the juncture period. Mixed effect models were fitted to assess whether this difference between deceased and recovered groups changed significantly between (1) pre-juncture and juncture (top panel showing effect sizes in each period and its right bar plots showing the p value assessing the temporal difference), or (2) pre-juncture and after the juncture (bottom panel). The size of the circle denotes the statistical significance of the difference between deceased and recovered groups in the indicated period (p < 0.05 is marked by solid outlines). See also Table S7. (C) Outcome prediction performance (recovered versus fatal) at (17–23 days; purple) or post (24–30 days; blue) juncture using leave-one-out cross-validation. Feature selection was performed using the same procedure that identified proteins shown in (B). Area under the curve (AUC) and permutation p values are shown. (D) Similar to Figure 6F but showing serum protein levels of critical ill patients with recovered or deceased outcomes (see A). Trajectories (using LOESS smoothing) were fitted to the recovered versus deceased patient groups separately. Top proteins showing the largest temporal shifts in their differences between the deceased versus recovered patients in (B) are shown. (E) Similar to (D) but for antibody measurements against SARS-CoV-2 spike and nucleocapsid proteins in critically ill patients with recovered or deceased outcomes. LU, light unit. See also Figure S7 and Tables S4, S6, and S7.