Immune phenotypes that predict COVID-19 severity

Abstract

Severe COVID-19 causes profound immune perturbations, but pre-infection immune signatures contributing to severe COVID-19 remain unknown. Genome-wide association studies (GWAS) identified strong associations between severe disease and several chemokine receptors and molecules from the type I interferon pathway. Here, we define immune signatures associated with severe COVID-19 using high-dimensional flow cytometry. We measured the peripheral immune system from individuals who recovered from mild, moderate, severe or critical COVID-19 and focused only on those immune signatures returning to steady-state. Individuals that suffered from severe COVID-19 showed reduced frequencies of T cell, MAIT cell and dendritic cell (DCs) subsets and altered chemokine receptor expression on several subsets, such as reduced levels of CCR1 and CCR2 on monocyte subsets. Furthermore, we found reduced frequencies of type I interferon-producing plasmacytoid DCs and altered IFNAR2 expression on several myeloid cells in individuals recovered from severe COVID-19. Thus, these data identify potential immune mechanisms contributing to severe COVID-19.

Extended Data Figure 1:. Cohorts and timing of sample collection

a) Distribution of days between symptom onset and sample collection is shown as histograms and boxplots for individuals recovered from mild, moderate, severe and critical COVID-19 cases. Individuals are highlighted as dots within boxplot. Red dashed line indicates 60 days cutoff which was used for analysis shown in Figures 1b and 2a and Extended Data Figure 5a. Information about time between symptom onset and sample collection was unavailable for two samples from the mild COVID-19 group. b) Linear regression between length of hospitalization in days and age is shown for severe and critical COVID-19 cases. c) Boxplot shows length of hospitalization in days for severe and critical COVID-19 cases. Wilcoxon test was performed to determine significant difference between severe and critical COVID-19 cases. d) Length of hospitalization in days (x-axis) is shown as triangle and circles highlight length in days between symptom onset and sample collection. Donors are depicted in rows (y-axis). Symbols from hospitalized individuals are connected by colored bar. Blue or red bars highlight if length of hospitalization is shorter or longer, respectively. COVID-19 study groups based on severity are shown separately. e) Length of hospitalization (grey bar), admission to ICU (red symbol) and sample collection (blue symbol) based on days from symptom onset is shown for severe and critical COVID-19 cases. Death is indicated by cross.

Extended Data Figure 2:

Demographics summary * Gender information not available for one individual.

Extended Data Figure 3:. Comparison of individuals recovered from non-severe and severe COVID-19

a) Volcano plot shows comparison of individuals recovered from non-severe (mild/moderate) and severe (severe/critical) COVID-19. P-values were obtained from logistic regression, included correction for age and experiment and were corrected for multiple testing using Benjamini-Hochberg false discovery rate. Log2 fold change was calculated based on the mean of immune traits within non-severe and severe COVID-19 cases. P-values are shown as −log10. b) Bar graph shows FDR-adjusted −log10 P-values for significant immune traits with P < 0.01 derived from Extended Data Figure 3a. Bars are colored based on log2 fold change and split based on decrease (top) or increase (bottom) in individuals recovered from severe COVID-19. Bar on the left indicates trait type.

Extended Data Figure 4:. Comparison of analysis between all and non-hospitalized individuals at time of sample collection

a) Volcano plot shows comparison of individuals recovered from non-severe (mild/moderate) and severe (severe/critical) COVID-19. Only individuals not hospitalized or discharged at day of sample collection are included. P-values were obtained from logistic regression, included correction for age and experiment and were corrected for multiple testing using Benjamini-Hochberg false discovery rate. Log2 fold change was calculated based on the mean of immune traits within non-severe and severe COVID-19 cases. P-values are shown as −log10. b) Bar graph shows FDR-adjusted −log10 P-values for immune traits significantly different between non-severe and severe COVID-19 cases (cut-off for P-value < 0.012). Plot is similar to Extended Data Figure 3b but depicts P-values obtained with only individuals not hospitalized or released at day of sample collection. Bar on the left indicates the immune trait type. Color of bars indicate log2 fold change between non-severe and severe COVID-19 cases calculated as the ratio between the mean of immune traits between the two severity groups. c) Comparison of stable immune traits between non-severe and severe COVID-19 cases including either all individuals (x-axis) or only individuals not hospitalized or released at day of sample collection (y-axis) is shown. Plot shows FDR-adjusted −log10 P-values for manually gated immune traits (N = 801). P-values were obtained by logistic regression and corrected for age and experiment batch. Size of symbols is based on −log10 P-values from analysis including all individuals. Color depicts the type of trait. d) Venn graph depicts overlap of immune traits which differed between non-severe and severe COVID-19 group obtained from analysis including all (red circle, traits from Extended Data Fig. 3b) or only non-hospitalized individuals (blue circle, traits from Extended Data Fig. 4b). e) Example flow cytometry data and gating of MAIT cells is shown (left) for one donor recovered from mild and severe COVID-19. Boxplot (right) shows frequency of MAIT cells per group. More detailed gating information is shown in Supplementary Data 3. f) Boxplots show frequencies of CD4 (top) and CD8 (bottom) central memory (CM) cells of conventional CD4 and CD8 T cells, respectively, from all study groups.

Extended Data Figure 5:. Dynamics of FlowSOM clusters in COVID-19

a) FlowSOM clusters affected by long-term perturbations were identified in individuals recovered from moderate (top) or severe COVID-19 (bottom) either by linear regression of cluster frequency and days between symptom onset and sample collection or Wilcoxon analysis of cluster frequency between early and late timepoints (cut-off >60 days between symptom onset and sample collection). −log10 P-values from both analyses are shown for All 388 FlowSOM clusters. P-value cutoff of 0.05 is shown by red line. Symbols are colored based on lineage and shaped based on CR1 (circle) or CR2 (triangle) panel. Symbol size is according to −log10 P-value from Wilcoxon analysis. b) Graph shows FDR-adjusted −log10 P-values derived from comparison of stable FlowSOM cluster (N = 291) frequencies between individuals recovered from non-severe and severe COVID-19. Analyses included either all individuals (x-axis) and only individuals not hospitalized or released at day of sample collection (y-axis). P-values were obtained by logistic regression correcting for age and experiment batch. Symbols are colored based on lineage and shape corresponds to CR1 or CR2 panel. Symbol size is based on FDR-adjusted −log10 P-value derived from analysis with all individuals. c) Venn graph shows overlap of significant FlowSOM clusters between individuals recovered from non-severe and severe COVID-19 from analysis including either all individuals (red circle) or only individuals not hospitalized or released at day of sample collection (blue circle).

Extended Data Figure 6:. Expression pattern of significant innate-like T cell clusters between non-severe and severe COVID-19

Expression (logicle-transformed fluorescence signal) of markers from CR1 (A) and CR2 (B) panel for innate-like T cell clusters are shown as overlaid histograms. T cell clusters are described in Figure 5c and d and are significantly different between individuals recovered from non-severe and severe COVID-19. All remaining clusters within innate-like T cells are depicted in grey and labeled as “Rest” as a reference population.

Extended Data Figure 7:. Expression pattern of significant myeloid cell clusters between non-severe and severe COVID-19

Expression (logicle-transformed fluorescence signal) of markers from CR1 (A) and CR2 (B) panel for myeloid cell clusters are shown as overlaid histograms. Myeloid cell clusters are described in Figure 6 and are significantly different between individuals recovered from non-severe and severe COVID-19. All remaining clusters within myeloid cells are depicted in grey and labeled as “Rest” as a reference population.

Figure 1:. Expression of functional markers and temporal dynamics of immune traits in COVID-19

a) Expression of chemokine receptors, CD40, CD86, IFNAR2, CD39 and TIGIT (rows) on immune cell populations (columns) is depicted. Median of mean fluorescence intensities (MFI) derived from 173 healthy individuals is visualized by min-max normalized color gradient. Dot size corresponds to median percentage of cells expressing these markers. Missing dots indicate that marker was not measured. b) Immune traits (N = 1779) at baseline or affected by long-term perturbations were distinguished in individuals recovered from moderate (left) and severe (right) COVID-19 cases. A combination of I) linear regression analysis between immune traits and days between symptom onset and sample collection and II) comparison of samples collected before and after 60 days of symptom onset using a Wilcoxon test was used as described in Online methods. Plot shows unadjusted −log10 P-values from both analyses. Dot size increases with significance from Wilcoxon test. Trait types are colored (Frequency of immune subset in blue, Frequency of expressing functional marker in orange and MFI values in green). Red line highlights threshold for unadjusted significance (P = 0.05).

Figure 2:. Long-term perturbations of immune traits in COVID-19

a) Top immune traits affected by long-term perturbations are depicted. Traits are derived from analysis in Figure 1b and selected for P-value <0.001 in one of both analysis (linear regression and/or Wilcoxon test). Bars pointing to the left and right are derived from linear regression and Wilcoxon test, respectively, and are colored based on trait type (Frequency of immune subset in blue, Frequency of expressing functional marker in orange and MFI values in green). Colored bar on the left depicts severity group from which the significant trait is derived. R² and slope from linear regression are shown as colored bars on the right. Values in the right bar are slope values from linear regression. Red and black dashed lines show P-value cut-off of 0.001 and 0.05, respectively. b) Frequencies of switched (top row) and naïve (bottom row) B cells of total B cells are shown. Plot on the left shows frequencies as boxplots for healthy subjects (grey) and individuals recovered from mild (purple), moderate (burgundy), severe (orange) or critical COVID-19 (yellow). The two plots on the right show the frequency of cells as a function of time between symptom onset and sample collection for individuals recovered from moderate and severe COVID-19. Far right plot shows the distribution of the traits in 173 healthy individuals. Similar to Fig. 2b, dynamics of c) CD38⁺HLA-DR⁻ (left) and CD38⁻ HLA-DR⁻ of CD4 naïve T cells (right), d) frequencies of cDC1s of total DCs, e) CCR3 MFI of basophils and f) CD95 MFI of early NK and NK2 cells are shown. Age-corrected residuals from linear regression were used for statistical analysis. For comparison between groups, one-way ANOVA was used on residuals to test for overall significant difference prior to Wilcoxon test with Bonferroni correction. Second and third plot show dot plots with linear regression (red line) and 95% confidence interval for individuals recovered from moderate and severe COVID-19, respectively. * P < 0.05, ** P < 0.01, *** P < 0.001, **** P < 0.0001

Figure 3:. Potential immune features at baseline predicting COVID-19 severity

a) Flow cytometry data (left) depicts CX3CR1 expression on early NK cells from an individual recovered from mild (top) and severe (bottom) COVID-19. Quantification of mean fluorescence intensity (MFI) of CX3CR1 on early NK cells is shown as boxplot for all study groups. b) Flow cytometry dot plots on top row (left and middle plot) depicts expression of CCR4 and CCR9 for a mild and severe COVID-19 case. Histogram overlay shows CXCR3 expression for the same cell subset and donors. MFI values for the same receptors are shown as boxplots for all groups (bottom row). c) Flow cytometry plot depicts expression of CCR4 on CD4⁺ naïve (top row), CD4⁺ transitional memory (TM, second row), CD8⁺ naïve (third row) and CD8⁺ TM (bottom row) T cells from an individual recovered from mild (left column) or severe (right column) COVID-19. Quantification of these subsets in all study groups are shown as boxplots (right). d) Flow cytometry plot depicts TIGIT expression on CD8⁺ naïve (top row), CD8⁺ stem-cell like memory (TSCM, second row), CD8⁺ central memory (CM, third row), CD8⁺ terminal effector* (TE*, fourth row) T cells and MAIT cells (bottom row) from an individual recovered from mild (left column) and severe (right column) COVID-19. Quantification of these subsets in all study groups are shown as boxplots (right). e) Flow cytometry data (left) depicts CD38 and HLA-DR expression on CD8 effector (EM; top) and terminal (TM; bottom) memory T cells from an individual recovered from mild (left) and critical (right) COVID-19. The gate defines CD38⁺HLA-DR⁺ activated T cells. Quantification of these subsets in all study groups are shown as boxplots (right). f) Flow cytometry example data (left) for gating of marginal zone (MZ) B cells from total B cells in an individual recovered from mild and severe COVID-19 is shown. Boxplot (right) shows frequencies of MZ B cells of total B cells in all study groups. Residuals from linear regression between immune trait and age were used to calculate statistics on age-corrected data. ANOVA with subsequent Wilcoxon test and Bonferroni correction on residuals was performed for statistics highlighted in boxplots. * P < 0.05, ** P < 0.01, *** P < 0.001, **** P < 0.0001

Figure 4:. Innate immune signatures predict COVID-19 severity

a) Example flow cytometry data for frequencies of pDCs from myeloid cells and inflammatory CD14⁺ DC3s of total DC3s is shown from an individual recovered from mild (top) and critical (bottom) COVID-19. Corresponding enumeration for all subjects based on study group are shown as boxplots (right). Precise delineation of pDCs is shown in Supplementary Data 1. b) Mean fluorescence intensity (MFI) of CX3CR1 on cDC1s (left) and pDCs (right) is shown for all study groups as boxplot (top row). Example flow cytometry data for CX3CR1 signal (red peak) on cDC1s (first column) and pDCs (second column) is shown as histogram for an individual recovered from mild (top row) and severe (bottom row) COVID-19. B cells (grey) and Monocytes (blue) are overlaid as reference populations known to lack and express CX3CR1, respectively. Numbers in histogram plots highlight MFI. c) Flow cytometry data (left) depicts CCR1 and CCR2 expression on classical (top), intermediate (middle) and non-classical (bottom) monocytes from a patient recovered from mild (left column) and critical (right column) COVID-19. Boxplots (right) show MFI values of CCR1 (first column) and CCR2 (second column) on the same monocyte populations for all study groups. d) Expression of IFNAR2 from an individual recovered from mild (red) and severe (blue) COVID-19 is shown as overlaid histogram (left) for classical monocytes (top), CD14⁺ DC3s (second row), pDCs (third row) and naïve B cells (bottom). Plot on the right depicts fold change of median IFNAR2 expression of each disease severity group compared to median IFNAR2 expression of healthy individuals on all defined myeloid (top) and non-myeloid (bottom) subsets. e) Boxplots show IFNAR2 MFI for all study groups for classical monocytes, CD14+ DC3s, pDCs and naïve B cells. Residuals from linear regression between immune trait and age were used to calculate statistics on age-corrected data. ANOVA with subsequent Wilcoxon test and Bonferroni correction on residuals was performed for statistics highlighted in boxplots. * P < 0.05, ** P < 0.01, *** P < 0.001, **** P < 0.0001

Figure 5:. Unsupervised analysis of immune system in individuals recovered from non-severe and severe COVID-19

a) FDR-adjusted −log10 P-values of FlowSOM clusters (N = 55) which differ significantly (P < 0.05) between individuals recovered from non-severe (mild/moderate) and severe (severe/critical) COVID-19 are shown. Bars are colored based on lineage (B cells, purple; CD4 T cells, orange; CD8 T cells, red; innate-like T cells, green; myeloid cells, blue; NK cells, pink). Bar on the left indicates whether traits originate from chemokine receptor panel 1 (CR1, grey) or 2 (CR2, black). b) Volcano plots show FDR-adjusted −log10 P-values and log2 fold change derived from comparison of FlowSOM clusters between individuals recovered from non-severe and severe COVID-19 cases. Main lineages are depicted in separated plots and contain FlowSOM clusters from both panels CR1 (circle) and CR2 (triangle). Data point size corresponds to −log10 P-values and color indicates log2 fold change. c) Heatmap depicts normalized median fluorescence intensity (MFI) values for lineage, differentiation and functional markers from top significant innate-like T cell clusters (Fig. 5a, P < 0.01). Values derived from CR1 (top) and CR2 (bottom) panels are separated. Heatmaps on the right highlight expression of markers specific for CR1 and CR2 panels including chemokine receptors, co-stimulatory markers and IFNAR2. Values are normalized based on trimmed 1–99% percentile values. Complete heatmaps for all innate-like T cell clusters are shown in Supplementary Data 12. d) Frequencies for same clusters described in Figure 5c are shown as boxplots based on study group. Values are log10(+ 1) transformed and plotted on linear scale. Logistic regression with correction for age and experiment batch was used to identify significant clusters between non-severe and severe COVID-19. Only FlowSOM clusters (N = 291) which did not show temporal changes within moderate and severe COVID-19 cases are shown as described in the Online methods section and results (Extended Data Fig. 5a). Residuals from linear regression between immune trait and age were used to calculate statistics on age-corrected data. ANOVA with subsequent Wilcoxon test and Bonferroni correction on residuals was performed for statistics highlighted in boxplots. * P < 0.05, ** P < 0.01, *** P < 0.001, **** P < 0.0001

Figure 6:. Myeloid cell populations from FlowSOM analysis as potential predictor for disease outcome

a) Frequencies (log10 +1) of myeloid cell clusters among top hits (P < 0.01) described in Figure 5a are shown for CR1 (top row) and CR2 panel (bottom row). Values are log10(+ 1) transformed and plotted on linear scale. Residuals from linear regression between immune trait and age were used to calculate statistics on age-corrected data. ANOVA with subsequent Wilcoxon test and Bonferroni correction on residuals was performed for statistics highlighted in boxplots. * P < 0.05, ** P < 0.01, *** P < 0.001, **** P < 0.0001 b) Heatmaps showing normalized median fluorescence intensity (MFI) values for clusters described in Figure 6a are shown. Heatmaps on the left show markers used to delineate immune cell subsets. On the right, heatmaps depict CR panel-specific markers. Values are normalized based on trimmed 1–99% percentile values. c) tSNE plots with myeloid cells from individuals recovered from mild (left column) or severe (right column) COVID-19 are shown. Data from panels CR1 and CR2 are shown in the top and bottom row, respectively. Each plot contains 50’000 subsampled myeloid cells (gating shown in Supplementary Data 1). Dots are colored based on FlowSOM cluster annotation and full data is shown in Supplementary Data 9. Clusters described in Figs. 6a and b are annotated and highlighted. d) Spearman analysis of normalized MFI values between clusters described in Figure 6a is shown in order to estimate the phenotypic overlap between CR1 and CR2 panel. Heatmap depicts Spearman correlation coefficient.