Immunophenotyping of COVID-19 and influenza highlights the role of type I interferons in development of severe COVID-19

Abstract

Although most SARS-CoV-2-infected individuals experience mild coronavirus disease 2019 (COVID-19), some patients suffer from severe COVID-19, which is accompanied by acute respiratory distress syndrome and systemic inflammation. To identify factors driving severe progression of COVID-19, we performed single-cell RNA-seq using peripheral blood mononuclear cells (PBMCs) obtained from healthy donors, patients with mild or severe COVID-19, and patients with severe influenza. Patients with COVID-19 exhibited hyper-inflammatory signatures across all types of cells among PBMCs, particularly up-regulation of the TNF/IL-1β-driven inflammatory response as compared to severe influenza. In classical monocytes from patients with severe COVID-19, type I IFN response co-existed with the TNF/IL-1β-driven inflammation, and this was not seen in patients with milder COVID-19. Interestingly, we documented type I IFN-driven inflammatory features in patients with severe influenza as well. Based on this, we propose that the type I IFN response plays a pivotal role in exacerbating inflammation in severe COVID-19.

Fig. 1. Single cell transcriptomes of PBMCs from COVID-19 and influenza patients.

(A) tSNE projections of 59,572 PBMCs from healthy donors (HDs) (4 samples, 17,590 cells), severe influenza (FLU) patients (5 samples, 10,519 cells), COVID-19 patients (asymptomatic: 1 sample, 4,425 cells; mild COVID-19: 4 samples, 16,742 cells; severe COVID-19: 6 samples, 10,296 cells) colored by group information. (B) Normalized expression of known marker genes on a tSNE plot. (C) tSNE plot colored by annotated cell types. EM: effector memory, NK cell: natural killer cell, DC: dendritic cell, RBC: red blood cell. (D) Proportion of cell types in each group excluding ‘Uncategorized 1’, ‘Uncategorized 2’, ‘RBC’, and ‘Platelet’. The colors indicate cell type information. (E) Boxplots showing the fold enrichment in cell type proportions from mild COVID-19 (n=4), severe COVID-19 (n=6), and FLU (n=5) patients compared to the HD group (mild COVID-19 vs. HD: n=16, severe COVID-19 vs. HD: n=24, FLU vs. HD: n=20). For the boxplots, the box represents the interquartile range (IQR) and the whiskers correspond to the highest and lowest points within 1.5 × IQR. ‘Uncategorized 1’ (relatively high UMIs per cells and presence of multiple marker genes), ‘Uncategorized 2’ (B cell-like and high expression of ribosomal protein genes), ‘RBC’, and ‘Platelet’ were excluded. Two-sided Kolmogorov–Smirnov (KS) tests were conducted for each cell type between the disease and HD group. *p<0.05, **p<0.01, and ***p<0.001.

Fig. 2. Immune landscape of COVID-19.

(A) Hierarchical clustering using the Pearson correlation coefficient (PCC) of a normalized transcriptome between diseases in cell type resolution (n = 33). The color intensity of the heat map indicates the PCC values. The color bars above the heat map indicate the cell type and disease group. The black box indicates the cell types that highly correlate between the severe COVID-19 and FLU groups. (B) Illustration of the enrichment p-values for the select GO biological pathways (n = 49) of differentially expressed genes (DEGs) in COVID-19 and FLU patients (left 6 columns: DEGs for COVID-19 and FLU groups compared to HD, right 2 columns: DEGs between COVID-19 and FLU groups). (C) tSNE plot of representative gene expression patterns for GBP1 (FLU specific), CREM (COVID-19 specific), and CCL3 (COVID-19/FLU common). (D) Top, dendrogram from WGCNA analysis performed using relative normalized gene expression between the COVID-19 and FLU groups for the genes belonging to the select biological pathways in (B) (n=316). Bottom, heat map of relative normalized gene expression between the COVID-19 and FLU groups. The color bar (left) indicates cell type information clustered by hierarchical clustering based on the PCC for relative normalized gene expression. Modularized gene expression patterns by WGCNA are shown together (G1, n=10; G2, n=147; G3, n=27; G4, n=17; G5, n=12; G6, n=64; G7, n=34; G8, n=5).

Fig. 3. Subpopulation analysis of CD8⁺ T cells.

(A) tSNE plot of the non-EM-like CD8⁺ T cell subpopulations in all groups (left, n=6,253), COVID-19 (top right, n=2,653), FLU (middle right, n=1,452), and HD (bottom right, n=2,148) colored by cluster information. (B, C) Boxplots showing the proportion of individual sub-clusters from the non-EM-like CD8⁺ T cell cluster within each group (COVID-19, n=10; FLU, n=5; HD, n=4). The proportions follow normal distribution as tested by the Shapiro-Wilk normality test except the proportion of cluster 3 in the COVID-19 group (p=0.04). Cluster 1 and cluster 3 were highly enriched in the FLU and COVID-19 group, respectively. Two-sided Welch’s t test p-values were 4.4E-03 between COVID-19 and FLU in cluster 1, 3.5E-02 between FLU and HD donor in cluster 1, 8.6E-03 between COVID-19 and FLU in cluster 3, and 5.8E-3 between COVID-19 and HD in cluster 3. *p<0.05, **p<0.01. (D) STRING analysis using the top 30 up-regulated genes in cluster 1 (left) and cluster 3 (right). (E) Bar plots showing enrichment p-values of eight representative GO biological pathways for pro-inflammation and interferon in cluster 1 or cluster 3-specific up-regulated genes (cluster 1, n=66; cluster 3, n=183).

Fig. 4. Transcriptome of classical monocytes in COVID-19 patients.

(A) Venn diagram of differentially expressed genes (DEGs) in COVID-19 and FLU compared to HD. The representative genes are shown together. (B) K-means clustering of DEGs between all pairs of FLU, mild COVID-19, and severe COVID-19 (n=499). The color indicates the relative gene expression between the diseases and HD. The representative genes are shown together. (C) Bar plots showing the average –log10(p-value) values in enrichment analysis using the perturbed genes of four different cell lines listed in L1000 LINCS for up-regulated genes in cluster 2 (C2, left) and cluster 3 (C3, right). Error bars indicate standard deviation. (D) Combined enrichment scores were compared between C2 and C3 for the gene sets of the type I IFN response (left; GSE26104) and TNF response (right; GSE2638, GSE2639). **p<0.01. Each dot indicates an individual subject. (E) Bar plots showing the average –log10(p-value) values in the enrichment analysis using the perturbed genes listed of four different cell lines in L1000 LINCS for up-regulated genes in cluster 4 (C4, left) and cluster 5 (C5, right). Error bars indicate standard deviation (C and E).

Fig. 5. Trajectory analysis of classical monocytes.

(A) Volcano plot showing DEGs between mild and severe COVID-19 groups. Each dot indicates individual gene, colored by red when a gene is significant DEG. (B) Bar plot showing the average –log10(p-value) values in enrichment analysis using the perturbed genes of four different cell lines listed in L1000 LINCS for up-regulated genes in the severe COVID-19 group. Error bars indicate standard deviation. (C) Trajectory analysis of classical monocytes from specimens obtained at two different time points in a single COVID-19 patient (mild: C7-2, 1,197 cells; severe: C7-1, 631 cells). The color indicates cluster information (left) or the severity of COVID-19 (right). (D) Relative expression patterns of representative genes in the trajectory analysis are plotted along the Pseudotime. The color indicates the relative gene expression calculated by Monocle 2. (E) Bar plots showing the average –log10(p-value) values in the enrichment analysis using the perturbed genes of four different cell lines in L1000 LINCS for up-regulated genes in cluster 3 (left) and cluster 1 (right). Error bars indicate standard deviation. (F) Comparison of combined enrichment scores between cluster 3 and cluster 1 for the gene sets from systemic lupus erythematosus (SLE) (n=16) and rheumatoid arthritis (RA) (n=5). ***p<0.001; ns, not significant. (G) GSEA of up-regulated genes in cluster 3 (left) and cluster 1 (right) to the ‘class 1’ gene module of monocyte-derived macrophages by Park et al. (2017). NES: normalized enrichment score, FDR: false discovery rate.

Fig. 6. Validation of the combined IFN-I and inflammatory responses in the transcriptome of post-mortem lung tissues from lethal COVID-19.

(A) UCSC Genome Browser snapshots of representative genes. (B) Bar plot showing the average –log10(p-value) values from the enrichment analysis using the perturbed genes of four different cell lines in L1000 LINCS for up-regulated genes (n= 386) in post-mortem lung tissues compared to biopsied healthy lung tissue. Error bars indicate standard deviation. (C) GSEA of significantly up- and down-regulated genes in post-mortem lung tissues for gene sets originated from up-regulated genes in C2 (n=96), C3 (n=143), C4 (n=218), and C5 (n=30) of Fig. 4B. (D and E) GSEA of significantly up- and down-regulated genes in post-mortem lung tissues for gene sets originated from the top 200 up-regulated genes in cluster 3 (left) and cluster 1 (right) from the trajectory analysis in Fig. 5C (D), and from gene sets originated from the top 200 up-regulated genes in classical monocytes of mild (left) and severe (right) COVID-19 (E).