A single-cell atlas of the peripheral immune response in patients with severe COVID-19

Nature Medicine Published: 08 June 2020

Abstract

There is an urgent need to better understand the pathophysiology of  Coronavirus disease 2019 (COVID-19), the global pandemic caused by  SARS-CoV-2, which has infected more than three million people worldwide1. Approximately 20% of patients with COVID-19 develop severe disease and 5% of patients require intensive care2. Severe disease has been associated with changes in peripheral immune  activity, including increased levels of pro-inflammatory cytokines3,4 that may be produced by a subset of inflammatory monocytes (炎症单核细胞产生的炎症细胞因子升高)5,6, lymphopenia（淋巴细胞减少）7,8 and T cell exhaustion（T 细胞耗尽）9,10. To elucidate pathways in peripheral immune cells that might lead to  immunopathology or protective immunity in severe COVID-19, we applied  single-cell RNA sequencing (scRNA-seq) to profile peripheral blood  mononuclear cells (PBMCs) from seven patients hospitalized for COVID-19, four of whom had acute respiratory distress syndrome, and six healthy  controls. We identify reconfiguration of peripheral immune cell  phenotype in COVID-19, including a heterogeneous interferon-stimulated  gene signature, HLA class II downregulation and a developing neutrophil  population that appears closely related to plasmablasts appearing in  patients with acute respiratory failure requiring mechanical  ventilation正在发展的嗜中性粒细胞群体，该群体与发生在需要机械通气的急性呼吸衰竭患者出现的浆母细胞密切相关. Importantly, we found that peripheral monocytes and  lymphocytes do not express substantial amounts of pro-inflammatory  cytokines 外周单核细胞和淋巴细胞不表达大量促炎性细胞因子 . Collectively, we provide a cell atlas of the peripheral  immune response to severe COVID-19.

Fig. 1:  Expansion of plasmablasts and depletion of multiple innate immune cell  subsets in the periphery of patients with COVID-19.

a, Demographics, sample characteristics and disease course of patients with COVID-19. b, UMAP dimensionality reduction embedding of peripheral blood mononuclear cells (PBMCs) from all profiled samples (n = 44,721 cells) colored by donor of origin. IDs of patients with COVID-19 (n = 7) begin with ‘C’ and are colored in shades of orange (patients who were  not ventilated at the time of draw) or red (patients with ARDS who were  ventilated at the time of draw) and those of healthy donors begin with  ‘H’ (n = 6) and are colored in blues. c, UMAP embedding of the entire dataset colored by orthogonally generated clusters labeled by manual cell type annotation. d, Proportions of each cell type in each sample colored by donor of origin. The x axes correspond to the ventilation or ARDS status of each patient. Shown are exact two-sided P values by the Wilcoxon rank-sum test. n = 6, n = 4 and n = 4 biologically independent samples for Healthy, NonVent and ARDS,  respectively. Boxplot features: minimum whisker, 25th  percentile − 1.5 × inter-quartile range (IQR) or the lowest value  within; minimum box, 25th percentile; center, median; maximum box, 75th  percentile; maximum whisker, 75th percentile + 1.5 × IQR or greatest  value within.

Fig. 2: Robust HLA class II downregulation and type I interferon-driven  inflammatory signatures in monocytes are characteristics of SARS-CoV-2  infection.

a, UMAP embedding of all monocytes colored by sample of origin. n = 10,339 cells are plotted from n = 14 biologically independent samples. b, UMAP embedding of monocytes colored by CD14 and FCGR3A (encoding CD16a, to distinguish between CD14+ and CD16+ monocytes), HLA-DPB1 and HLA-DMA (illustrating HLA class II downregulation in patients with COVID-19) and S100A9 and IFI27 (demonstrating canonical inflammatory signatures in patients with COVID-19). c, UMAP embedding of monocytes colored by genes encoding pro-inflammatory  cytokines previously reported to be produced by circulating monocytes in severe COVID-196, namely TNF, IL6, IL1B, CCL3, CCL4 and CXCL2. d,g,h, Heatmaps of DE genes (d), differentially regulated canonical pathways (g) and differentially regulated predicted upstream regulators (h) between CD14+ monocytes of each COVID-19 sample compared to CD14+ monocytes of all healthy controls. The heatmap in d is colored by average log(fold-change), while heatmaps in g and h are colored by z-score. All displayed genes, pathways and regulators are statistically significant at the P < 0.05 confidence level by Seurat’s implementation of the Wilcoxon rank-sum  test (two-sided, adjusted for multiple comparisons using Bonferroni’s  correction, in d) or Ingenuity Pathway Analysis (IPA) implementation of the Fisher exact test (right-tailed, in g and h). The 50 genes (d), 25 pathways (g) or 50 regulators (h) with the highest absolute average log(fold-change) or z-score across all donors are labeled. Genes with a net positive average log(fold-change) or z-score are labeled in red; genes with a net negative average log(fold-change) or z-score are labeled in blue. DPS, days post-symptom onset; DTF, days from first reported or measured fever. e, Boxplot showing the mean HLA class II module score of CD14+ monocytes from each sample, colored by healthy donors (blue),  non-ventilated patients with COVID-19 (orange) or ventilated patients  with COVID-19 (red). Shown are exact P values by two-sided Wilcoxon rank-sum test. n = 6, n = 4 and n = 4 biologically independent samples for Healthy, NonVent and ARDS, respectively. f, Dot plot depicting percent expression and average expression of all detected HLA genes in CD14+ monocytes by donor. i, Boxplot showing the IFNA module score of each cell, colored by healthy donors (blue),  non-ventilated patients with COVID-19 (orange) or ventilated patients  with COVID-19 (red). j, Scatter plots depicting the correlation between the mean ISG module score of CD14+ monocytes in each sample and the patient age (top) and time–distance  from first measured or reported fever (bottom). Shown are Pearson’s r, exact two-sided P values and the 95% confidence interval. n = 8 (top) and n = 6 (bottom) independent biological samples. Number of cells for d,f–i: C1 A, 1,561; C1 B, 1,858; C2, 217; C3, 1,102; C4, 713; C5, 462; C6,  277; C7, 2,095; H1, 680; H2, 325; H3, 215; H4, 166; H5, 444; H6, 224.  For d,g–h, cells from all healthy controls (n = 2,054 cells) are used to generate comparisons with each COVID-19 sample. For e,i, boxplot features: minimum whisker, 25th percentile − 1.5 × IQR or the  lowest value within; minimum box, 25th percentile; center, median;  maximum box, 75th percentile; maximum whisker, 75th  percentile + 1.5 × IQR or greatest value within.

Fig. 3: Heterogeneous patterns of NK cell exhaustion and IFN response in COVID-19.

a, UMAP embedding of CD4+ T cells, CD8+ T cells and NK cells colored by sample of origin. b, UMAP embedding colored by lineage genes (CD3D, CD3G, CD4, CD8A, FCGR3A and NCAM1) and selected functional/phenotypic markers (GZMB and MKI67). For a,b, n = 22,016 cells are plotted from n = 14 biologically independent samples. c, Boxplots depicting proportions of CD56dim NK cells, CD56bright NK cells and proliferating lymphocytes among total T and NK cells by  sample of origin. The cells used to calculate each proportion are  highlighted in bold black in the adjacent UMAP embeddings and were  identified by manually labeling clusters generated by clustering CD4+ T cells, CD8+ T cells and NK cells alone. Shown are exact two-sided P values from the Wilcoxon rank-sum test. n = 386 (top), n = 4,899 (middle), n = 781 (bottom) total from n = 6, n = 4 and n = 4 biologically independent samples for Healthy, NonVent and ARDS, respectively. d, Boxplot showing the mean expression score by only NK cells of three canonical markers of NK cell exhaustion: LAG3, PDCD1 (encoding PD-1) and HAVCR2 (encoding TIM-3). Shown are exact two-sided P values by Wilcoxon rank-sum test. e, Boxplot showing the mean expression score by only NK cells of four canonical NK cell cytokine genes (CCL3, CCL4, IFNG and TNF). Shown are exact P values by Wilcoxon rank-sum test. For d,e, n = 6, n = 4 and n = 4 biologically independent samples for Healthy, NonVent and ARDS, respectively. In c–e, boxplot features: minimum whisker, 25th percentile − 1.5 × IQR or the  lowest value within; minimum box, 25th percentile; center, median;  maximum box, 75th percentile; maximum whisker, 75th  percentile + 1.5 × IQR or greatest value within. f,g, Heatmaps of DE genes (f) and differentially regulated predicted upstream regulators (g) between NK cells of each COVID-19 sample compared to NK cells of all healthy controls. As in Fig. 2, f is colored by average log(fold-change), while g is colored by z-score. All displayed genes and regulators are statistically significant at the P < 0.05 confidence level by Seurat’s implementation of the Wilcoxon rank-sum  test (two-sided, adjusted for multiple comparisons using Bonferroni’s  correction, f) or IPA’s implementation of Fisher exact test (right-tailed, g). The 50 genes or regulators with the highest absolute average log(fold-change) or z-score across all donors are labeled. Genes with a net positive average log(fold-change) or z-score are labeled in red; genes with a net negative average log(fold-change) or z-score are labeled in blue. DPS, days post-symptom onset; DTF, days from first reported or measured fever. Number of cells for f,g: C1 A, 354; C1 B, 387; C2, 271; C3, 328; C4, 104; C5, 518; C6, 58; C7, 130; cells from all healthy controls (n = 4,707 cells) were used to generate comparisons with each COVID-19 sample. h,i, Heatmaps of differentially upregulated ISGs (h; Supplementary Table 25) and cytokines (i; Supplementary Table 25) in donors with COVID-19, colored by the number of COVID-19 samples in  which the gene was differentially expressed relative to all healthy  controls. DE genes used to construct these heatmaps are provided in  Supplementary Tables 4–10. An ISG or cytokine was counted as differentially expressed if it had an average log(fold-change) > 0.25 and an adjusted two-sided P value < 0.05 by Seurat’s implementation of the Wilcoxon rank-sum test. n = 8 biologically independent COVID-19 samples compared to n = 6 biologically independent healthy controls.