Integrating single-cell sequencing data with GWAS summary statistics reveals CD16+monocytes and memory CD8+T cells involved in severe COVID-19

Abstract

Background: Understanding the host genetic architecture and viral immunity contributes to the development of effective vaccines and therapeutics for controlling the COVID-19 pandemic. Alterations of immune responses in peripheral blood mononuclear cells play a crucial role in the detrimental progression of COVID-19. However, the effects of host genetic factors on immune responses for severe COVID-19 remain largely unknown.

Methods: We constructed a computational framework to characterize the host genetics that influence immune cell subpopulations for severe COVID-19 by integrating GWAS summary statistics (N = 969,689 samples) with four independent scRNA-seq datasets containing healthy controls and patients with mild, moderate, and severe symptom (N = 606,534 cells). We collected 10 predefined gene sets including inflammatory and cytokine genes to calculate cell state score for evaluating the immunological features of individual immune cells.

Results: We found that 34 risk genes were significantly associated with severe COVID-19, and the number of highly expressed genes increased with the severity of COVID-19. Three cell subtypes that are CD16+monocytes, megakaryocytes, and memory CD8+T cells were significantly enriched by COVID-19-related genetic association signals. Notably, three causal risk genes of CCR1, CXCR6, and ABO were highly expressed in these three cell types, respectively. CCR1⁺CD16+monocytes and ABO⁺ megakaryocytes with significantly up-regulated genes, including S100A12, S100A8, S100A9, and IFITM1, confer higher risk to the dysregulated immune response among severe patients. CXCR6⁺ memory CD8+ T cells exhibit a notable polyfunctionality including elevation of proliferation, migration, and chemotaxis. Moreover, we observed an increase in cell-cell interactions of both CCR1⁺ CD16+monocytes and CXCR6⁺ memory CD8+T cells in severe patients compared to normal controls among both PBMCs and lung tissues. The enhanced interactions of CXCR6⁺ memory CD8+T cells with epithelial cells facilitate the recruitment of this specific population of T cells to airways, promoting CD8+T cell-mediated immunity against COVID-19 infection.

Conclusions: We uncover a major genetics-modulated immunological shift between mild and severe infection, including an elevated expression of genetics-risk genes, increase in inflammatory cytokines, and of functional immune cell subsets aggravating disease severity, which provides novel insights into parsing the host genetic determinants that influence peripheral immune cells in severe COVID-19.

Keywords: COVID-19; GWAS; Immune cells; Inflammatory storm; Single-cell sequencing.

Fig. 1

The workflow for this integrative genomic analysis. A Combination of single-cell RNA sequencing data and GWAS summary statistics on severe COVID-19 based on two independent methods. One method is the regression-based polygenic model based on whole scRNA-seq profiles, and another is the generalized linear regression model based on the top 10% most specific genes for each cell type. B An increase in genetics-risk genes and cytokines for severe COVID-19. C Cellular interaction analysis of genetics-influenced immune cell subsets with epithelial cells

Fig. 2

Risk genes and pathways associated with hospitalized COVID-19 from meta-GWAS summary data. A Manhattan plot and quantile-quantile (QQ) plot of meta-GWAS analysis highlighting eight risk genetic loci for hospitalized COVID-19. The red horizontal line represents the genome-wide significance threshold of P < 5×10⁻⁸. The genomic inflation factor λ = 1.02. B Nine index SNPs within eight genomic loci associated with hospitalized COVID-19. The left panel shows the P value of each index SNP, and the right panel shows the odds ratio with 95% confidence interval. C Circus plot showing the results of the S-MultiXcan-based analysis. The inner ring demonstrates the 22 autosomal chromosomes (Chr1-22). In the outer ring, a circular symbol represents a specific gene and color marks the statistical significance of the gene for hospitalized COVID-19 (red marks FDR < 0.05, orange indicates 6.96×10⁻⁵ ≤ P < 0.001, light blue marks 0.001 ≤ P ≤ 0.05, and dark blue indicates P > 0.0). D PPI network of these 34 identified risk genes based on the STRING database (v11.0, https://string-db.org/). The orange ring represents druggable genes targeted by at least one known drug. E Network module constructed by using the Jaccard distance showing the connectivity of 10 significant pathways enriched by 34 risk genes. F Heatmap showing the results of hierarchical clustering analysis of 27 risk genes on COVID-19 severity. Seven risk genes did not express in dataset #1, and the expression level of each gene was scaled. G The proportion of highly expressed genes among 27 risk genes in normal controls and in the three phases of COVID-19 (mild, moderate, and severe patients). Using 10,000 times of permutation analysis to calculate the significance of the observation (permuted P = 0.023). H Plot showing an increase of the significantly enriched pathways in the network module with elevated COVID-19 severities. Orange color represents a significantly enriched pathway (FDR ≤ 0.05) and gray color represents a non-significantly enriched pathway (FDR > 0.05)

Fig. 3

Integrative analysis identifies genetic associations between peripheral immune cells and severe COVID-19. A Bar graph showing the results of the combination of scRNA-seq data and GWAS summary statistics on severe COVID-19 based on the RolyPoly among normal controls and patients with different severities (i.e., mild, moderate, and severe). The y-axis shows the 13 cell types, and the x-axis shows the mean negative log-transformation P value (-Log2(P)). Orange color indicates a cell type showing a significant association, and light blue represents there is no significant association. B UMAP projections of peripheral immune cells colored by annotated cell types. The plot showing the region of CD16+monocytes, megakaryocytes, and memory CD8+T cells. The red dot represents positive gene expressions of CCR1⁺, ABO⁺, and CXCR6⁺, and gray stands for negative cells

Fig. 4

CCR1⁺ CD16+momocytes contribute higher risk to cytokine storms among severe COVID-19 patients. A Boxplot showing the difference in inflammatory cytokine score between CCR1⁺ and CCR1⁻ CD16+ monocytes. A two-side Wilcoxon sum-rank test was used. B Volcano plot showing differentially expressed genes between CCR1⁺ and CCR1⁻ CD16+ monocytes. C Significantly enriched pathways by 351 highly expressed genes among CCR1⁺ CD16+ monocytes. Color legend represents the log-transformed FDR value (-Log10(FDR)). D Bar graph showing the proportion of CCR1⁺ CD16+ monocytes among normal, mild, moderate, and severe groups. E Boxplot showing the inflammatory cytokine score of CCR1⁺ CD16+ monocytes among normal, mild, moderate, and severe groups. The Mann-Kendall trend analysis was used. F Bar graph showing the differentially up-DEGs among different COVID-19 patients compared with normal controls. Namely, mild COVID-19 vs. normal, moderate COVID-19 vs. normal, and severe COVID-19 vs. normal. Venn plot on top of bar showing the overlapped up-DEGs between moderate and severe patients. G The correlation of up-DEGs between moderate and severe patients. Pearson correlation analysis was used to calculate the correlation coefficient and P value. H–J Representative up-DEGs among CCR1⁺ CD16+ monocytes showing significantly elevated expressions with increased COVID-19 severities. HS100A8, IS100A9, and JIFITM1. K Disease-term enrichment analysis on 190 up-DEGs based on the GLAD4U database. The y-axis shows -Log10(FDR), and the x-axis shows the enrichment ratio

Fig. 5

Multi-functionality of CXCR6⁺ memory CD8+T cells for severe COVID-19. A–D Boxplots showing the difference in A cytokine score, B chemokine score, C IFN-ɑ/β response score, and D T cell activation score between CXCR6⁺ and CXCR6⁻ memory CD8+T cells. A two-side Wilcoxon sum-rank test was used. E Volcano plot showing differentially expressed genes between CXCR6⁺ and CXCR6⁻ memory CD8+T cells. F Bar graph showing the proportion of CXCR6⁺ memory CD8+T cells among normal, mild, moderate, and severe groups. G–I Boxplots showing the G chemokine score, H T cell activation score, and I migration score of CXCR6⁺ memory CD8+T cells among normal, mild, moderate, and severe groups. The Mann-Kendall trend analysis was used. J Venn plot showing the overlapped up-DEGs between pairwise comparisons: mild vs. normal, moderate vs. normal, and severe vs. normal. K Representative gene of GZMH among CXCR6⁺ memory CD8+T cells showing significantly elevated expressions with increased COVID-19 severities. L Heatmap showing up-DEGs in CXCR6⁺ memory CD8+T cells from pairwise comparisons: mild vs. normal, moderate vs. normal, severe vs. normal. The up-DEGs listed in the green panel were from mild vs. normal, in the yellow panel were from moderate vs. normal, and in the orange panel were from severe vs. normal. M Scatter plot showing the enriched GO biological processes by 108 up-DEGs among CXCR6⁺ memory CD8+T cells. The x-axis shows -Log10(FDR), and the y-axis shows the enrichment ratio

Fig. 6

Cell-to-cell interactions of CCR1⁺ CD16+momocytes and CXCR6⁺ memory CD8+T cells with other cells in PBMC and BALF. A, B Boxplot showing the number of cellular interactions of ACCR1⁺ CD16+ monocytes and BCCR1⁻ CD16+ monocytes with other immune cells in PBMC between normal controls and patients with increased COVID-19 severities. C Predicted cellular interactions of CCR1⁺ CD16+ monocytes with other immune cells in PBMC, comparing severe COVID-19 vs. normal control. D, E Boxplot showing the number of cellular interactions of DCXCR6⁺ memory CD8+T cells and ECXCR6⁻ memory CD8+T cells with other immune cells in PBMC between normal controls and patients with increased COVID-19 severities. F Predicted cellular interactions of CXCR6⁺ memory CD8+T cells with other immune cells in PBMC, comparing severe COVID-19 vs. normal control. G Boxplot showing an increase in cellular interactions with other cells in BALF for CCR1⁺ CD16+ monocytes than CCR1⁻ CD16+ monocytes. H Predicted cellular interactions with other cells in BALF, comparing CCR1⁺ CD16+ monocytes with CCR1⁻ CD16+ monocytes. I Boxplot showing an increase in cellular interactions with other cells in BALF for CXCR6⁺ memory CD8+T cells than CXCR6⁻ memory CD8+T cells. J Predicted cellular interactions with other cells in BALF, comparing CXCR6⁺ memory CD8+T cells with CXCR6⁻ memory CD8+T cells. The circular size represents the significance of each ligand-receptor axis, and color represents the communication probability