COVID-19 immune features revealed by a large-scale single-cell transcriptome atlas

Abstract

A dysfunctional immune response in coronavirus disease 2019 (COVID-19) patients is a recurrent theme impacting symptoms and mortality, yet a detailed understanding of pertinent immune cells is not complete. We applied single-cell RNA sequencing to 284 samples from 196 COVID-19 patients and controls and created a comprehensive immune landscape with 1.46 million cells. The large dataset enabled us to identify that different peripheral immune subtype changes are associated with distinct clinical features, including age, sex, severity, and disease stages of COVID-19. Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) RNA was found in diverse epithelial and immune cell types, accompanied by dramatic transcriptomic changes within virus-positive cells. Systemic upregulation of S100A8/A9, mainly by megakaryocytes and monocytes in the peripheral blood, may contribute to the cytokine storms frequently observed in severe patients. Our data provide a rich resource for understanding the pathogenesis of and developing effective therapeutic strategies for COVID-19.

Keywords: B cell receptor sequencing; COVID-19; SARS-CoV-2; T cell receptor sequencing; cell-cell interaction; cytokine storm; host cell range; ligand-receptor interaction; single-cell RNA-seq; single-cell transcriptomics.

Graphical abstract

Figure 1

Multi-tissue and multi-stage single-cell atlas of COVID-19 patients and healthy controls (A) Flowchart depicting the overall experimental design of this study. Cells circled with dash lines were enriched in samples from the disease progression stage. (B) Overview of the cell clusters in the integrated single-cell transcriptomes of 1,462,702 cells derived from COVID-19 patients and healthy controls. Clusters were named based on the cluster-specific gene expression patterns, in which we used “high” or “low” labels to indicate the relative expression levels in the corresponding clusters. Genes without high or low labels were specifically expressed in the corresponding clusters. (C) Tissue preference of each cluster measured by the ratio of observed to randomly expected cell numbers (R_O/E) (Zhang et al., 2018). (D) Patient group preference of each cluster measured by R_O/E. See also Figure S1 and Tables S1 and S2.

Figure S1

Basic characteristics of the integrated dataset and selected markers of cell subsets in different major cell lineages, related to Figure 1 (A) The age distribution of the dataset (color-coded by disease conditions). (B) Distribution of sex. Chi-square test. (C-E) Distribution of unique molecular identifier (UMI) counts per cell (C), gene counts per cell (D), and percentage of mitochondrial transcripts per cell (E) detected for cells in various tissue types. PBMC, peripheral blood mononuclear cells; BALF, bronchoalveolar lavage fluid; PFMC/Sputum, pleural effusion/sputum. (F-J) Violin plots of selected marker genes (rows) for cell subsets (columns) within each cell lineage, including 6 B/plasma B cell clusters (F), 23 Myeloid cell clusters (G), 3 NK cell clusters (H), 4 Epithelial cell clusters (I) and 28 T cell clusters (J).

Figure S2

Comparison of different immune cell types among patient groups, related to Figure 2 (A) Comparison on the major cell type level based on 159 unsorted PBMC samples with at least one thousand cells available in the scRNA-seq data. NK, natural killer cells; Mono, monocytes; DC, dendritic cells; Mega, megakaryocytes. (B) State transition between B_c05-MZB1-XBP1 and other B cell sub clusters quantified by the STARTRAC algorithm based on TCR clonotypes (Zhang et al., 2018). Clonotypes with more than 5 cells were shown in the right panel. (C) Percentage of B_c06-MKI67 in PBMC across disease conditions based on the same cohort with (A). (D) RNA velocity analysis shows the transition potential from B_c03-CD27-AIM2 to B_c05-MZB1-XBP1. Cell pairs transiting from B_c03-CD27-AIM2 to B_c05-MZB1-XBP1 or vice versa were quantified in the bar plot. (E) Percentage of B_c03-CD27-AIM2 in PBMC across disease conditions based on the same cohort with (A). (F) State transition quantified by STARTRAC (Zhang et al., 2018) between T_CD4_c13-MKI67-CCL5^low proliferating cells and other CD4 cell sub-clusters (left) and clones containing T_CD4_c13-MKI67-CCL5^low cells with more than 5 cells (right). (G) Percentage of T_CD4_c04−ANXA2 across disease conditions based on the same cohort with (A). (H) Sex differences of T_CD4_c04-ANXA2. Single-side Wilcoxon test. (I and J) Percentage of T_gdT_c14-TRDV2 and T_CD8_c09-SLC4A10 across disease conditions based on the same cohort with (A). Adjusted P-values smaller than 0.05 are indicated (two-sided unpaired Wilcoxon test).

Figure 2

Associations of patient age, sex, COVID-19 severity, and stage with cellular compositions in PBMCs (A) Heatmap for q values of ANOVA. Sample type, fresh or frozen; sample time, days after symptom onset. (B) Composition comparison for plasma B cells (B_c05-MZB1-XBP1) based on 159 unsorted PBMC samples with at least 1,000 cells available in the scRNA-seq data. (C) Classes of heavy chains for B_c05-MZB1-XBP1. (D–G) Composition comparison for DC_c4−LILRA4, Neu_c3−CST7, T_CD4_c13-MKI67-CCL5 low, and T_CD8_c10-MKI67-GZMK. (H) Associations between age and T_CD8_c01−LEF1 (Spearman’s correlation). (I) Sex differences of T_CD4_c08−GZMK−FOS^high. Adjusted p values < 0.05 are indicated (two-sided unpaired Wilcoxon tests). See also Figure S2 and Table S3.

Figure S3

Effects of sampling time and sample processing methods (fresh or frozen) on immune cell composition and the BCR/TCR diversity, related to Figures 2 and 3 (A) Gross relationship between B_c05−MZB1−XBP1 frequency in PBMC and sampling days. ANOVA rejected the association between B_c05−MZB1−XBP1 frequency and sampling days after incorporating age, sex, COVID-19 severity and stage (Figure 2A). (B) Gross relationship between DC_c4−LILRA4 frequency in PBMC and sampling days. (C) Gross relationship between Neu_c3−CST7 frequency in PBMC and sampling days. (D-G) Comparison among patient groups for T_CD4_c02−AQP3, T_CD8_c01−LEF1, T_CD8_c02−GPR183, and T_CD4_c08−GZMK−FOS^high via separating fresh and frozen PBMC samples. (H-K) Gross relationship of sampling time with frequencies of T_CD4_c02−AQP3, T_CD4_c08−GZMK−FOS^high, T_CD8_c01−LEF1, and T_CD8_c02−GPR183.

Figure 3

Associations of patient age, sex, COVID-19 severity, and stage with the diversity of B and T cell repertoires in PBMCs (A) Heatmap for q values of ANOVA. Sample type, fresh or frozen; sample time, days after symptom onset. (B–E) Comparison for T_CD4_c02−AQP3, T_CD4_c08−GZMK−FOS^high, T_CD8_c01−LEF1, and T_CD8_c02−GPR183. (F) Sex differences for T_CD4_c08−GZMK−FOS^high. (G–I) Age associations of the TCR diversity of T_CD8_c01−LEF1, T_CD8_c05-ZNF683, and T_CD8_c09-SLC4A10 (Spearman’s correlation). (J) V gene usage of published SARS-CoV-2 neutralized antibodies and their relationship with those differentially used IGHV genes in our dataset. Gini-index was used to quantify the skewness of the V gene usage of the published SARS-CoV-2 neutralized antibodies. IGHV genes differentially used by moderate or severe COVID-19 patients compared with healthy controls and their intersections are shown with different colors. Venn diagram is used to show their overlaps with those published SARS-CoV-2 antibodies. Adjusted p values < 0.05 are indicated (two-sided unpaired Wilcoxon test). See also Figure S3 and Table S3.

Figure 4

Cell types with SARS-Cov-2 RNA detected (A) 3,085 cells with SARS-CoV-2 RNA detected (UMI > 0) from BALF (6/12) and sputum (2/22) samples. No cells from PBMCs or PFMCs were detected as SARS-CoV-2 positive. (B) Markers used to determine cell types. Goblet and basal cells were merged as secretory epithelial cells for convenience in the subsequent analyses. (C) Viral load in each cell quantified by log(CPM). (D) Expression levels of host factors reported to associate with SARS-CoV-2 infection in literature. (E) Pearson’s correlations of host factor expression with viral load (zero-expression cells were excluded from regression analysis to reduce the effects of dropouts). (F) Expression levels of ISGs in cells with viral RNA detected. (G) Detection rates of SARS-CoV-2 genes in different cell types on both 10x Genomics 5′ and 3′ platforms. Given a viral gene g_v, the detection rate is defined as the ratio of the number of g_v⁺ cells to the total viral-RNA-positive cells of the specific cell type and then normalized by the gene length in the SARS-CoV-2 genome. (H) IHC staining of CD3 and SARS-CoV-2 spike protein in pulmonary tissue. Scale bar, 100 μM. See also Figure S4 and Tables S4 and S5.

Figure S4

Characteristics of SARS-CoV-2-RNA-positive epithelial and immune cells, related to Figure 4 (A) Associations of BSG with viral RNA load in neutrophils, plasma cells, T/NK cells, and ciliated epithelial cells (Person’s correlation). Grey points (no expression or dropouts) were excluded from the regression analysis to reduce the impacts of dropouts in scRNA-seq. (B) Violin plots showing the expression of ISGs in viral RNA-positive cells (from BALF) compared with viral RNA-negative cells from PBMC and BALF. Two-sided unpaired Wilcoxon test was used. (C) Pearson’s correlation between viral RNA load and the expression levels of ISGs. Grey points (no expression or dropouts) were excluded from the regression analysis to exclude the impacts of dropouts in scRNA-seq. (G) Detection rates of SARS-CoV-2 genes in different cell types on both 10x Genomics 5′ and 3′ platforms. (H) IHC staining of SARS-CoV-2 spike protein in lymphocytes in pulmonary tissue.

Figure 5

Impact of viral RNA presence on the expression and cell-cell interaction of epithelial subtypes (A) Volcano plot showing differentially expressed genes between squamous cells with or without viral RNA detected. Adjusted p value < 0.05, two-sided unpaired Wilcoxon test. ANXA1 is denoted in dark blue. (B) Enriched Gene Ontology (GO) terms in genes highly expressed in virus-positive squamous cells shown in (A). (C) Venn plot showing the intersection of genes upregulated in different epithelial cells with viral detection. (D) Cell-cell interaction networks of one severe COVID-19 patient (left) and one moderate COVID-19 patient (right) inferred by CSOmap-based on data from BALFs. Interactions with q values < 0.1 are shown. Significance: −log₁₀(q values). (E) Boxplots showing the self-distances among ciliated, secretory, and squamous cells with or without viral RNA detection in the pseudo-spaces predicted by CSOmap. Each dot represents an individual patient. Two-sided paired Wilcoxon test. (F) Violin plot showing the self-distances of three types of epithelial cells with viral detection (exemplified by one patient). Two-sided unpaired Wilcoxon test. (G) Boxplot showing the median self-distances of three type of epithelial cells with viral detection from all the patients with BALF samples. Each dot represents an individual patient. Two-sided unpaired Wilcoxon test. (H) Pie charts showing the ligand-receptor contribution to the interaction between virus-positive squamous cells and virus-positive neutrophils (top) and virus-positive squamous cells and virus-positive macrophages (bottom). (I) Boxplot showing the interactions between squamous cells (with and without viral detection) and macrophage (left) and neutrophils (left) with viral detection. Each dot represents an individual patient. Two-sided unpaired Wilcoxon test. Normalized connections: observed cell-cell interactions normalized by random expectation (n_A × n_B, where n is the cell number of type A or B). See also Figures S5 and S6.

Figure S5

Differences of various epithelial cells with viral RNA detected in the interaction potential with other cells, related to Figure 5 (A) Differential expression of ANXA1 in SARS-CoV-2 RNA-positive and negative squamous epithelial cells. (B) 2D visual view of the pseudo-space constructed by CSOmap with the location of ciliated cells highlighted. Each dot denotes a single cell and is colored by its cell type. (C and D) Self-distance of viral RNA-positive and negative ciliated and squamous cell groups in the pseudo-space shown in (B). Two-sided unpaired Wilcoxon test. (E) Comparison of interacting potentials of viral RNA-positive secretory epithelial cells with BALF Macro_c1-C1QC cells between moderate and severe patients. Spatial connections within the pseudo-space constructed by CSOmap were used for quantification, which were normalized by the cell numbers of both clusters. Error bar: s.e.m across different patients. (F) The ligand-receptor contribution between viral RNA-positive secretory epithelial cells and Macro_c6-VCAN cells. (G) Dot plot showing the expression level of MARCO in BALF samples. Pct, percentage of expressed cells.

Figure S6

The expression of selected genes in PBMC and BALF samples, related to Figure 5 (A) Dot plots showing the expression of S100A9 in cell clusters found in PBMCs. Each dot is colored by the mean expression and sized by the scaled mean (Z scores). The blue box highlights the expressions in patients belonging to the progression (severe) group. (B) Dot plots showing the expression of ANXA1 (top), FPR1 (middle) and TLR4 (bottom) in clusters found in PBMCs. (B) Dot plots showing the expression of ANXA1 (first panel), FPR1 (second panel), S100A9 (third panel), S100A8 (fourth panel) and TLR4 (bottom panel) in clusters found in BALFs. Each dot is colored by the means of the expression and sized by the scaled means (Z scores).

Figure 6

Mono_c1-CD14-CCL3 and megakaryocytes in peripheral blood appear as a dominant source for the inflammatory cytokine storm (A) t-SNE plots of PBMCs colored by major cell types (top left panel), inflammatory cell types (top right panel), cytokine score (middle panel), and inflammatory score (bottom panel). (B) Heatmap and unsupervised clustering of cell proportion of seven hyper-inflammatory cell subtypes (row normalized). (C) Boxplots of the cell proportion of Mono_c1-CD14-CCL3, Mega, and T_CD4_c08-GZMK-FOS^high clusters from healthy controls (n = 20), convalescent patients (moderate, n = 48), patients with progression (moderate, n = 18), convalescent patients (severe, n = 35), and patients with progression (severe, n = 38). Two-sided Wilcoxon test. (D) Ordinary least-squares model of age to cell proportion of Mono_c1-CD14-CCL3, Mono_c2-CD14-HLA-DPB1, and Mono_c3-CD14-VCAN clusters from healthy controls (blue, n = 20), convalescent patients (purple, n = 83), and patients with progression (red, n = 56). p value was assessed with the F-statistic for ordinary least-squares model. (E) Heatmap of cytokine expression among seven hyper-inflammatory cell subtypes (red font) and other clusters (gray font). (F) Boxplots of cytokine expression based on scRNA-seq and plasma profiling for healthy controls (n = 20 for scRNA-seq, and n = 5 for plasma), convalescent patients (severe, n = 5, for both scRNA-seq and plasma), and patients with progression (severe, n = 14, for both scRNA-seq and plasma). Two-sided Wilcoxon test. (G) Boxplots of the cytokine expression of Mono_c1-CD14-CCL3, Mega and T_CD8_c06-TNF clusters from healthy controls (n = 20), convalescent patients (moderate, n = 48), patients with progression (moderate, n = 18), convalescent patients (severe, n = 35), and patients with progression (severe, n = 38). Two-sided Wilcoxon test. (H) Ordinary least-squares model of age to cytokine expression of Mono_c1-CD14-CCL3, Mega, and T_CD8_c06-TNF clusters from healthy controls (blue, n = 20) and patients with progression (n = 18 + 38). p value was assessed with the F-statistic for ordinary least-squares model. In (C), (F), and (G), the box represents the second and third quartiles and median, whiskers each extend 1.5 times the interquartile range, and dots represent outliers. In (B) and (E), Mono_c1, Mono_c2, Mono_c3, T_CD4_c08, T_CD8_c09, T_CD8_c06, and Mega correspond to Mono_c1-CD14-CCL3, Mono_c2-CD14-HLA-DPB1, Mono_c3-CD14-VCAN, T_CD4_c08-GZMK-FOS^high, T_CD8_c09-SLC4A10, T_CD8_c06-TNF, and Mega, respectively. In (E), T_CD4_c11, T_CD8_c03, T_CD8_c04, T_CD8_c05, T_CD8_c07, T_gdT_c14, and T_CD8_c08, NK_c01 correspond to clusters of T_CD4_c11-GNLY, T_CD8_c03-GZMK, T_CD8_c04-COTL1, T_CD8_c05-ZNF683, T_CD8_c07-TYROBP, T_gdT_c14-TRDV2, T_CD8_c08-IL2RB, and NK_c01-FCGR3A, respectively. DC, dendritic cells; Mega, megakaryocytes; Mono, monocytes. See also Figure S7 and Tables S6 and S7.

Figure S7

Identification of hyper-inflammatory subtypes associated with cytokine storm in PBMCs, related to Figure 6 (A) t-SNE plots of PBMC cells colored by cytokine score (top panel) and inflammatory score (bottom panel). (B) The proportion of subtypes from healthy controls (n = 20), progression (severe, n = 38) and average of all samples (n = 159) (top panel); the inflammatory score (middle panel) and cytokine score (bottom panel) of subtypes from healthy controls (n = 20), convalescence (moderate, n = 48), progression (moderate, n = 18), convalescence (severe n = 35) and progression (severe, n = 38) patients. Significance was evaluated with Mann-Whitney rank test for each subtype versus all the other subtypes. ^∗∗∗∗p < 0.0001. (C) Boxplots of the proportion of inflammatory cell-types and other cell-types from healthy controls (n = 20), convalescence (moderate, n = 48), progression (moderate n = 18), convalescence (severe, n = 35) and progression (severe, n = 38) patients. Two-sided Wilcoxon rank-sum test. (D) Pie chart showing the proportion of 4 classified groups (named ‘both’, ‘Mono’, ‘Mega’ and ‘neither’) based on the proportion of Mono_c1-CD14-CCL3 and Mega cell-types in patients at the progression (severe) stage. (E) Boxplots of the inflammatory and cytokine score within 4 classified groups (named ‘both’, ‘Mono’, ‘Mega’ and ‘neither’). (F) Bar graphs showing cytokine concentration at the plasma levels of CCL3, IFNG, IL1RN and TNF from healthy controls (n = 5), convalescent (n = 7), non-severe (n = 4), severe (n = 4), death case (n = 7) patients. Shown are P values by Student’s t test. (G) Boxplot of CXCL8 expression of Mono_c1-CD14-CCL3 subtype and IFNG expression of T_CD8_c06-TNF subtype from the scRNA-seq datasets with influenza (n = 5). (H) Ordinary least-squares model of age to IFNG signal from array data (n = 310) with influenza. P value was assessed with F-statistic for ordinary least-squares model. In panel (B), (C), (E) and (G), the box represents the second, third quartiles and median, whiskers each extend 1.5 times the interquartile range; dots represent outliers. In panel (F), all points are shown and bars represent mean with the 95% confidence intervals. DC, dendritic cells. Mega, megakaryocytes. Mono, monocytes.

Figure 7

The interactions of hyper-inflammatory cell subtypes in lung and peripheral blood (A) t-SNE plots of BALF cells colored by major cell types (top panel), cytokine score (middle panel), and inflammatory score (bottom panel). (B) Boxplots of the inflammatory score (top panel) and cytokine score (bottom panel) of cell subtypes. Significance was evaluated with the Wilcoxon rank-sum test. ^∗∗∗∗p < 0.0001. (C) Heatmap and unsupervised clustering of cytokine expression of five hyper-inflammatory cell subtypes. (D) Circos plot showing the prioritized interactions mediated by ligand-receptor pairs between inflammation-related cell types from BALF and PBMCs, respectively. The outer ring displays color-coded cell types, and the inner ring represents the involved ligand-receptor interacting pairs. The line width and arrow width are proportional to the log fold change between severe and moderate progression groups in ligand and receptor, respectively. Colors and types of lines are used to indicate different types of interactions as shown in the legend. The bar plot at bottom indicates the interaction score for each interaction, which serves to measure the interaction strength. DC, dendritic cells; Epi, epithelial cells; Macro, macrophage cells; Mono, monocytes; Neu, neutrophils. See also Figure S8 and Table S6.

Figure S8

Intercellular interaction alterations among cell types between severe and moderate progression sample groups, related to Figure 7 (A). Circos plot showing the prioritized interactions mediated by ligand-receptor pairs between inflammation-related cell subtypes for each tissue, namely, PBMC (left panel) and BALF (right panel). The outer ring displays color coded cell types and the inner ring represents the involved ligand-receptor interacting pairs. The line width and arrow width are proportional to the log fold change between severe progression and moderate progression patient groups in ligand and receptor, respectively. Colors and types of lines are used to indicate different types of interactions as shown in the legend. The barplot at bottom indicates the interaction score for each ligand-receptor interaction which serves to measure the interaction strength. (B) Summary illustration depicting the potential cytokine/receptor interactions of hyper-inflammatory cell subtypes involved in the cytokine storm. DC, dendritic cells. Epi, epithelial cells. Macro, macrophage cells. Mono, monocytes. Neu, neutrophils. Mega, megakaryocytes.