Single-cell analyses and host genetics highlight the role of innate immune cells in COVID-19 severity

Abstract

Mechanisms underpinning the dysfunctional immune response in severe acute respiratory syndrome coronavirus 2 infection are elusive. We analyzed single-cell transcriptomes and T and B cell receptors (BCR) of >895,000 peripheral blood mononuclear cells from 73 coronavirus disease 2019 (COVID-19) patients and 75 healthy controls of Japanese ancestry with host genetic data. COVID-19 patients showed a low fraction of nonclassical monocytes (ncMono). We report downregulated cell transitions from classical monocytes to ncMono in COVID-19 with reduced CXCL10 expression in ncMono in severe disease. Cell-cell communication analysis inferred decreased cellular interactions involving ncMono in severe COVID-19. Clonal expansions of BCR were evident in the plasmablasts of patients. Putative disease genes identified by COVID-19 genome-wide association study showed cell type-specific expressions in monocytes and dendritic cells. A COVID-19-associated risk variant at the IFNAR2 locus (rs13050728) had context-specific and monocyte-specific expression quantitative trait loci effects. Our study highlights biological and host genetic involvement of innate immune cells in COVID-19 severity.

Fig. 1. Study design and single-cell transcriptional analysis of PBMCs from COVID-19 patients and healthy controls.

a, Overview of the study design. The image was created using BioRender.com. b, UMAP embedding of scRNA-seq data for all 895,460 cells. Thirteen cell types were defined by RNA expression of marker genes (Extended Data Fig. 1a). c, Graph representation of Nhoods identified by Milo. Nodes are Nhoods, colored by their log₂ FC between COVID-19 patients (n = 72) and healthy controls (n = 75) adjusted by age and sex. Nondifferential abundance Nhoods (FDR ≥ 0.1) are colored white, and sizes correspond to the number of cells in a Nhood. Graph edges depict the number of cells shared between adjacent Nhoods. d, Beeswarm plot showing the distribution of adjusted log₂ FC in abundance between COVID-19 patients and healthy controls in Nhoods according to 13 cell types. Colors are represented in the same way as in c. e,f, The module score of Type I IFN response and IFN-γ response in PBMCs. The score was calculated using a gene set termed ‘GOBP_RESPONSE_TO_TYPE_I_INTERFERON’ (GO:0034340) and ‘GOBP_RESPONSE_TO_INTERFERON_GAMMA’ (GO:0034341), respectively. Upper heatmaps depicting the difference between average scores of 13 cell types and that of all single cells. The module scores of cells in each cell type were compared with the average score of all PBMCs using a two-sided one-sample t-test. Lower heatmaps depicting the difference between average scores of moderate or severe disease group and those of the healthy group in each of 13 cell types (n = 75 healthy, n = 9 moderate, n = 64 severe). The module scores of cells of moderate or severe disease group were compared to those of healthy group in each cell type using a two-sided Welch’s t-test. *P_uncorrected < 1.0 × 10⁻⁵⁰, **P_uncorrected < 1.0 × 10⁻³⁰⁰. Mono, monocytes; Pro, proliferative; Nhood, neighborhood.

Fig. 2. Defective IFN-γ response and reduced transition potential to ncMono in monocytes of severe COVID-19.

a, UMAP embedding of 116,944 monocytes and DC. Seven cell types were defined by RNA expression of marker genes (Extended Data Fig. 2a). b, Graph representation of Nhoods identified by Milo. Nodes are Nhoods, colored by their log₂ FC between COVID-19 (n = 72) and healthy controls (n = 75) adjusted by age and sex. Nondifferential abundance Nhoods (FDR ≥ 0.1) are colored white. c, Beeswarm plot showing the distribution of adjusted log₂ FC in abundance between COVID-19 and healthy controls in Nhoods according to seven cell types. Colors are represented in the same way as in b. d, The top ten enriched biological processes by GO analysis of upregulated DEGs of moderate and severe disease compared to healthy group in five cell types. Dot color indicates the statistical significance of the enrichment (adjusted P values via the Benjamini–Hochberg method), and dot size represents gene ratio annotated to each term. e, The differential gene expression analysis between moderate (n = 8) and severe (n = 64) COVID-19 in ncMono. DEGs (FDR < 0.05 and FC > 2) are colored in light blue and labeled by gene symbols if log₂ FC > 1.5. f, Velocities derived from the dynamical model for monocytes and DC cluster from COVID-19 and healthy group are projected into a UMAP-based embedding. Colors indicate the same clusters as in a. g, Average unspliced ratio of each sample stratified by three monocytes clusters, colored by COVID-19 (n = 73) and healthy (n = 75) groups. Condition-specific regression lines are shown. P value for the interaction effect between three monocyte clusters and two clinical conditions is uncorrected and reflects two-sided test.

Fig. 3. Differential abundance analysis of T cells and NK cells and TCR analysis.

a, UMAP embedding of 628,715 T and NK cells. Thirteen cell types were defined by RNA expression of marker genes (Extended Data Fig. 3a). b, Graph representation of Nhoods identified by Milo. Nodes are Nhoods, colored by their log₂ FC between COVID-19 (n = 72) and healthy controls (n = 75) adjusted by age and sex. Nondifferential abundance Nhoods (FDR ≥ 0.1) are colored white. c, Beeswarm plot showing the distribution of adjusted log₂ FC in abundance between COVID-19 and healthy controls in Nhoods according to 13 cell types. Colors are represented in the same way as in b. d, The distribution of the clone state of T cells in each cluster according to disease status. Differences of average clonal expansion rate of each sample between clinical conditions were evaluated in each cluster using two-sided Welch’s t-test (*P_uncorrected < 0.05, **P_uncorrected < 0.005), where only cells mapped with TCRs were included in the analysis (n = 75 healthy, n = 9 moderate, n = 64 severe). e, UMAP embedding of T cells (nine cell types) colored by clonal expansion size. Left panel shows clonal expansion divided into three categories, and right panel shows clonal expansion sizes ranging from 0 to 500. f, Network plots showing similarity of TCRα and TCRβ CDR3 amino acid sequence for each sample, disease status and cell types. Clonotype clusters with clonal size ≥50 are selected. g, T cells that were suspected to be specific to SARS-CoV-2 based on CDR3 amino acid sequence were projected on UMAP. Ef, effector.

Fig. 4. Differential abundance, gene expression and clonal analysis of B cells.

a, UMAP embedding of four cell types of 123,728 B cells (Extended Data Fig. 4a). b, Graph representation of Nhoods identified by Milo. Nodes are Nhoods, colored by their log₂ FC between COVID-19 (n = 72) and healthy controls (n = 75) adjusted by age and sex. Nondifferential abundance Nhoods (FDR ≥ 0.1) are colored white. c, Beeswarm plot showing the distribution of adjusted log₂ FC in abundance between COVID-19 and healthy controls in Nhoods according to four cell types. d, The differential gene expression analysis between moderate (n = 8) and severe (n = 64) COVID-19 in B_naive, B_memory and B_activated. DEGs (FDR < 0.05 and FC > 2) are colored in light blue and labeled by gene symbols. e, UMAP embedding of B cells colored by clonal expansion size. Left panel shows clonal expansion divided into three categories, and right panel shows clonal expansion size from 0 to 50. f, The distribution of the clone state of B cells in each cluster according to disease status. Differences of average clonal expansion rate of each sample between disease status were evaluated in each cluster using two-sided Welch’s t-test (*P_uncorrected = 0.012, **P_uncorrected = 1.3 × 10⁻⁷), where only cells mapped with BCRs were included in the analysis (n = 75 healthy, n = 9 moderate, n = 64 severe). g, Network plots showing similarity of CDR3 amino acid sequence in BCR heavy and light chain for each sample, disease status and cell types. Clonotype clusters with clonal size ≥10 are selected.

Fig. 5. Differential cell–cell interactions between COVID-19 patients and healthy controls and within COVID-19 severity.

a, Heatmaps depicting number of ligand–receptor pairs connecting cell–cell pairs for COVID-19 (n = 73) and healthy controls (n = 75), respectively. Rows indicate cells expressing the ligands and columns indicate cells expressing the receptors. An asterisk indicates cell–cell interactions with a number of ligand–receptor of more than 100, and cell types that share such interactions with at least one cell type are highlighted in red. b, Heatmap depicting the cell-connectivity-summary networks based on mean expression weight between COVID-19 (n = 73) and healthy controls (n = 75). An asterisk indicates the cell–cell interactions with FC of mean expression ≥1.5, and cell types that share cell–cell interactions with FC of mean expression ≥1.5 with five or more cell types are highlighted in red. c, Heatmap depicting the cell-connectivity-summary networks based on mean expression weight between moderate (n = 9) and severe (n = 64) COVID-19. An asterisk and red highlight mean the same as in b. d, Cell–cell interaction of IFNG/IFNGR and CXCL10/CXCR3 around nonclassical monocytes. Heatmaps depicting the cell-connectivity-summary networks based on mean expression weight of IFNG/IFNGR (left) and CXCL10/CXCR3 (right) according to three conditions, respectively (n = 75 healthy, n = 9 moderate, n = 64 severe). The image was created using BioRender.com. Ef, effector.

Fig. 6. Associations of PBMC cell types with host genetic risk of COVID-19.

a, UMAP embedding of COVID-19 PBMCs (n = 72) colored by scDRS disease score calculated from GWAS summary statistics of three phenotypes from COVID-19 HGI (round 6). b, Differences of disease score in individual cell-level among three phenotypes. c, Heatmaps depicting each cell type-disease association for three phenotypes, respectively. Heatmap colors denote uncorrected P value of cell-type-disease association evaluated using scDRS. Squares denote significant cell type-disease associations (FDR < 0.05), and cross symbols denote significant heterogeneity in association with disease across individual cells within a given cell type (FDR < 0.05). FDR was calculated via the Benjamini–Hochberg method across all pairs of cell types and three phenotypes.

Fig. 7. Context-specific and cell type-specific eQTL analysis of COVID-19-associated risk variants.

a, Heatmaps depicting eQTL effect of COVID-19-associated risk variants on V2G, highest gene prioritized by the V2G score of Open Target Genetics, by six major cell types, separately for COVID-19 (n = 67) and healthy controls (n = 75). Heatmap colors denote uncorrected P value of eQTL effect. b, Heatmaps depicting eQTL effect of COVID-19-associated risk variants by three cell types of monocytes, separately for COVID-19 (n = 67) and healthy controls (n = 75). c, rs13050728 eQTL effect on IFNAR2 expression in monocytes. The box plot shows the eQTL effect in COVID-19 (red) and healthy controls (blue). P values are uncorrected and reflect two-sided tests. Boxes denote the IQR, and the median is shown as horizontal bars. Whiskers extend to 1.5 times the IQR. All samples are shown as individual points (n = 67 COVID-19, n = 75 healthy). In a and b, an asterisk indicates FDR < 0.05, and FDR was calculated via the Benjamini–Hochberg method across all pairs of the cell types and five variants, separately for COVID-19 and healthy controls. IQR, interquartile range.

Extended Data Fig. 1. Cell type annotation of PBMC, differential abundance analysis and IFN responses.

(a) Violin plots showing the expression distribution of selected canonical cell markers in the 13 clusters. The rows represent selected marker genes and the columns represent clusters. (b) Tile plot showing percentage concordance between the manually annotated 25 clusters and Azimuth annotation. (c) A bar plot of the proportion of cell types shown in Fig. 1b, separated by three conditions (n = 75 healthy, n = 9 moderate, n = 64 severe). (d) Graph representation of neighborhoods identified by Milo in COVID-19 patients. Nodes are neighborhoods, colored by their log₂ fold change between moderate (n = 8) and severe (n = 64) COVID-19 patients adjusted by age, sex, time since symptom onset and duration of systemic steroids treatment. Nhood, neighborhood. (e) Box plot showing the distribution of adjusted log₂ fold change in abundance between moderate and severe COVID-19 in neighborhoods according to 13 cell types. Boxes denote the interquartile range (IQR), and the median is shown as horizontal bars. Whiskers extend to 1.5 times the IQR, and outliers are shown as individual points. (f,g) UMAP embedding of PBMCs colored by Type I IFN response score and IFN-γ response score. The score was calculated using a gene set termed ‘GOBP_RESPONSE_TO_TYPE_I_INTERFERON’ (GO:0034340) and ‘GOBP_RESPONSE_TO_INTERFERON_GAMMA’ (GO:0034341), respectively.

Extended Data Fig. 2. Immunological features of monocytes and dendritic cells.

(a) Violin plots showing the expression distribution of selected canonical cell markers in the seven clusters. The rows represent clusters and the columns represent selected marker genes. (b) A bar plot of the proportion of cell types, separated by three conditions (n = 75 healthy, n = 9 moderate, n = 64 severe). (c) Graph representation of neighborhoods identified by Milo in COVID-19 patients. Nodes are neighborhoods, colored by their log₂ fold change between moderate (n = 8) and severe (n = 64) COVID-19 patients adjusted by age, sex, time since symptom onset and duration of systemic steroids treatment. Nhood, neighborhood. (d) Box plot showing the distribution of adjusted log₂ fold change in abundance between moderate and severe COVID-19 in neighborhoods according to seven cell types. Boxes denote the interquartile range (IQR), and the median is shown as horizontal bars. Whiskers extend to 1.5 times the IQR, and outliers are shown as individual points. (e) The top 20 enriched biological processes by GO analysis of upregulated DEGs of COVID-19 (n = 72) compared to healthy controls (n = 75) in five cell types. Dot color indicates the statistical significance of the enrichment (adjusted P-values via the Benjamini-Hochberg method), and dot size represents gene ratio annotated to each term. (f) Heatmaps depicting the average unspliced ratio of seven cell types stratified by COVID-19 (n = 73) and healthy controls (n = 75).

Extended Data Fig. 3. Immunological features of T cells and NK cells and TCR analysis.

(a) Violin plots showing the expression distribution of selected canonical cell markers in 13 clusters. The rows represent clusters and the columns represent selected marker genes. (b) A bar plot of the proportion of cell types, separated by three conditions (n = 75 healthy, n = 9 moderate, n = 64 severe). (c) Graph representation of neighborhoods identified by Milo in COVID-19 patients. Nodes are neighborhoods, colored by their log₂ fold change between moderate (n = 8) and severe (n = 64) COVID-19 patients adjusted by age, sex, time since symptom onset and duration of systemic steroids treatment. Nhood, neighborhood. (d) Box plot showing the distribution of adjusted log₂ fold change in abundance between moderate and severe COVID-19 in neighborhoods according to 13 cell types. Boxes denote the interquartile range (IQR), and the median is shown as horizontal bars. Whiskers extend to 1.5 times the IQR, and outliers are shown as individual points. (e) UMAP embedding of TCR detection. (f) Repertoire overlap according to clonotype size between COVID-19 (n = 73) and healthy controls (n = 75) in all clonotypes (left) and in clonotypes with size < 50 (right). (g) Repertoire overlap according to clonotype size between moderate (n = 9) and severe (n = 64) group in all clonotypes (left) and in clonotypes with size < 50 (right). (h) The distribution of the clonal state of T cells which were suspected to be specific to SARS-CoV-2 based on CDR3 amino acid sequence in each cluster according to disease status (n = 75 healthy, n = 9 moderate, n = 64 severe).

Extended Data Fig. 4. Immunological features of B cells and BCR analysis.

(a) Violin plots showing the expression distribution of selected canonical cell markers in four clusters. (b) A bar plot of the proportion of cell types, separated by three conditions (n = 75 healthy, n = 9 moderate, n = 64 severe). (c) Graph representation of neighborhoods identified by Milo in COVID-19 patients. Nodes are neighborhoods, colored by their log₂ fold change between moderate (n = 8) and severe (n = 64) COVID-19 patients adjusted by age, sex, time since symptom onset and duration of systemic steroids treatment. Nhood, neighborhood. (d) Box plot showing the distribution of adjusted log₂ fold change in abundance between moderate and severe COVID-19 in neighborhoods according to four cell types. Boxes denote the interquartile range (IQR), and the median is shown as horizontal bars. Whiskers extend to 1.5 times the IQR, and outliers are shown as individual points. (e) The top 20 enriched biological processes by GO analysis of upregulated DEGs of COVID-19 (n = 72) compared to healthy controls (n = 75) in four cell types. Dot color indicates the statistical significance of the enrichment (adjusted P-values via the Benjamini-Hochberg method), and dot size represents gene ratio annotated to each term. (f) The top 10 enriched biological processes by GO analysis of upregulated DEGs of moderate (n = 8) and severe (n = 64) compared to healthy group (n = 75). Dot color and dot size mean the same as in (e). (g) Repertoire overlap according to clonotype size between COVID-19 (n = 73) and healthy controls (n = 75) in all clonotypes (left) and in clonotypes with size < 10 (right). (h) Repertoire overlap according to clonotype size between moderate (n = 9) and severe (n = 64) group in all clonotypes (left) and in clonotypes with size < 10 (right).

Extended Data Fig. 5. Cell-cell interaction analysis using CellPhoneDB and NATMI.

(a) Heatmaps depicting cell-cell communications quantified by CellPhoneDB in COVID-19 (n = 73) and healthy controls (n = 75), respectively. (b) Heatmaps depicting cell-cell communications quantified by NATMI in moderate (n = 9) and severe (n = 64) group, respectively. Rows indicate cells expressing the ligands and columns indicate cells expressing the receptors. Asterisk indicates cell-cell interactions with a number of ligand-receptor more than 100, and cell types which share more than 100 cell-cell interactions with at least one cell type are highlighted in red.

Extended Data Fig. 6. Assessing the association of PBMC cell types with host genetic risk.

(a) Flowchart of integrated analysis of scRNA-seq data and COVID-19 HGI (round6) GWAS summary statistics using scDRS. The figure design is based on the original paper of scDRS. (b) Differences of scDRS disease score between COVID-19 phenotypes in six major cell types. Differences of average disease scores of each sample were evaluated using two-sided paired t-test (adjusted for multiple comparisons using Bonferroni’s correction across all pairs of six cell-types and three GWAS comparisons, *P < 1×10^-10, **P < 1×10^-20, ***P < 1×10^-30). Boxes denote the interquartile range (IQR), and the median is shown as horizontal bars. Whiskers extend to 1.5 times the IQR, and all COVID-19 patients (n　= 72) are shown as individual points. (c) Heatmaps depicting cell type-disease association for three phenotypes in the manually annotated 25 clusters of COVID-19 scRNA-seq dataset. Heatmap colors denote uncorrected P-value of cell type-disease association evaluated using scDRS. Squares denote significant cell type-disease associations (FDR < 0.05), and cross symbols denote significant heterogeneity in association with disease across individual cells within a given cell type (FDR < 0.05). FDR was calculated via the Benjamini-Hochberg method across all pairs of 25 cell types and three GWAS phenotypes.

Extended Data Fig. 7. Single-cell eQTL analysis of rs13050728 and differential expression analysis of IFNAR2.

(a) Flowchart of how to select COVID-19-associated variants in single-cell eQTL analysis. (b) rs13050728 eQTL for IFNAR2 in classical monocytes. The box plot shows the eQTL in COVID-19 (red) and healthy controls (blue). Boxes denote the interquartile range (IQR), and the median is shown as horizontal bars. Whiskers extend to 1.5 times the IQR. All samples are shown as individual points (n = 67 COVID-19, n = 75 healthy). (c) Co-plots of the eQTL effect sizes of rs13050728 and 95% confidence intervals between COVID-19 and healthy controls. Dots represent the effect sizes and bars represent the 95 % confidence intervals. Cell types with significant interaction (P_Interaction < 0.05) in eQTL effect between COVID-19 and healthy controls are colored in red, and the rest are colored in blue. (d) Expression changes of IFNAR2 in differential expression (DE) analysis between COVID-19 (n = 72) vs healthy controls (n = 75) are indicated for each cell cluster. Cell types with P < 0.05 in DE analysis are colored in green, and the rest is colored in gray. P-values are uncorrected and reflect two-sided tests in (b-d).