doi: 10.1038/s41588-023-01375-1.
Epub 2023 Apr 24.
Single-cell analyses and host genetics highlight the role of innate immune cells in COVID-19 severity
Collaborators,
Affiliations
- PMID: 37095364
- PMCID: PMC10181941
- DOI: 10.1038/s41588-023-01375-1
Item in Clipboard
Single-cell analyses and host genetics highlight the role of innate immune cells in COVID-19 severity
Nat Genet.
2023 May.
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.
© 2023. The Author(s).
Conflict of interest statement
All authors declare no conflict of interest.
Figures
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 log2 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 log2 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. *Puncorrected < 1.0 × 10−50, **Puncorrected < 1.0 × 10−300. Mono, monocytes; Pro, proliferative; Nhood, neighborhood.
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 log2 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 log2 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 log2 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.
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 log2 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 log2 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 (*Puncorrected < 0.05, **Puncorrected < 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.
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 log2 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 log2 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 (*Puncorrected = 0.012, **Puncorrected = 1.3 × 10−7), 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.
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.
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.
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.
(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 log2 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 log2 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.
(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 log2 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 log2 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).
(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 log2 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 log2 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).
(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 log2 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 log2 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).
(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.
(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.
(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 (PInteraction < 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).
Similar articles
-
Integrating single-cell sequencing data with GWAS summary statistics reveals CD16+monocytes and memory CD8+T cells involved in severe COVID-19.Genome Med. 2022 Feb 17;14(1):16. doi: 10.1186/s13073-022-01021-1. Genome Med. 2022. PMID: 35172892 Free PMC article.
-
Single-cell transcriptomes of peripheral blood cells indicate and elucidate severity of COVID-19.Sci China Life Sci. 2021 Oct;64(10):1634-1644. doi: 10.1007/s11427-020-1880-y. Epub 2021 Feb 4. Sci China Life Sci. 2021. PMID: 33564978 Free PMC article.
-
Single-cell epigenomic landscape of peripheral immune cells reveals establishment of trained immunity in individuals convalescing from COVID-19.Nat Cell Biol. 2021 Jun;23(6):620-630. doi: 10.1038/s41556-021-00690-1. Epub 2021 Jun 9. Nat Cell Biol. 2021. PMID: 34108657 Free PMC article.
-
Variants in ACE2; potential influences on virus infection and COVID-19 severity.Infect Genet Evol. 2021 Jun;90:104773. doi: 10.1016/j.meegid.2021.104773. Epub 2021 Feb 17. Infect Genet Evol. 2021. PMID: 33607284 Free PMC article. Review.
-
Lessons from dermatology about inflammatory responses in Covid-19.Rev Med Virol. 2020 Sep;30(5):e2130. doi: 10.1002/rmv.2130. Epub 2020 Jul 12. Rev Med Virol. 2020. PMID: 32656939 Free PMC article. Review.
Cited by
-
Clonal landscape of autoantibody-secreting plasmablasts in COVID-19 patients.Life Sci Alliance. 2024 Sep 17;7(12):e202402774. doi: 10.26508/lsa.202402774. Print 2024 Dec. Life Sci Alliance. 2024. PMID: 39288992 Free PMC article.
-
Robust detection of infectious disease, autoimmunity, and cancer from the paratope networks of adaptive immune receptors.Brief Bioinform. 2024 Jul 25;25(5):bbae431. doi: 10.1093/bib/bbae431. Brief Bioinform. 2024. PMID: 39226888 Free PMC article.
-
Quantification of escape from X chromosome inactivation with single-cell omics data reveals heterogeneity across cell types and tissues.Cell Genom. 2024 Aug 14;4(8):100625. doi: 10.1016/j.xgen.2024.100625. Epub 2024 Jul 30. Cell Genom. 2024. PMID: 39084228 Free PMC article.
-
Immune Stimulation with Imiquimod to Best Face SARS-CoV-2 Infection and Prevent Long COVID.Int J Mol Sci. 2024 Jul 12;25(14):7661. doi: 10.3390/ijms25147661. Int J Mol Sci. 2024. PMID: 39062904 Free PMC article. Review.
-
A Novel Human Extravascular Monocyte Subset with Antiviral Functions Is Crucial for Resolving Lung Tissue Infection.bioRxiv [Preprint]. 2024 Aug 15:2024.03.08.583965. doi: 10.1101/2024.03.08.583965. bioRxiv. 2024. PMID: 38496468 Free PMC article. Preprint.
References
Publication types
MeSH terms
LinkOut - more resources
Full Text Sources
Medical
