Dictionary of immune responses to cytokines at single-cell resolution

Cytokines mediate cell–cell communication in the immune system and represent important therapeutic targets1–3. A myriad of studies have highlighted their central role in immune function4–13, yet we lack a global view of the cellular responses of each immune cell type to each cytokine. To address this gap, we created the Immune Dictionary, a compendium of single-cell transcriptomic profiles of more than 17 immune cell types in response to each of 86 cytokines (>1,400 cytokine–cell type combinations) in mouse lymph nodes in vivo. A cytokine-centric view of the dictionary revealed that most cytokines induce highly cell-type-specific responses. For example, the inflammatory cytokine interleukin-1β induces distinct gene programmes in almost every cell type. A cell-type-centric view of the dictionary identified more than 66 cytokine-driven cellular polarization states across immune cell types, including previously uncharacterized states such as an interleukin-18-induced polyfunctional natural killer cell state. Based on this dictionary, we developed companion software, Immune Response Enrichment Analysis, for assessing cytokine activities and immune cell polarization from gene expression data, and applied it to reveal cytokine networks in tumours following immune checkpoint blockade therapy. Our dictionary generates new hypotheses for cytokine functions, illuminates pleiotropic effects of cytokines, expands our knowledge of activation states of each immune cell type, and provides a framework to deduce the roles of specific cytokines and cell–cell communication networks in any immune response.

a-k, Upregulated GPs following cytokine treatment with respect to PBS controls for a, IFN-γ, b, IL-18, c, IL-33, d, IL-36α, e, IL-2, f, IL-4, g, IL-7, h, IL-15, i, IL-3, j, GM-CSF, and k, TNF-α.GPs that are significantly upregulated between cytokine and PBS treatment (FDR < 0.01 and effect size > 1) in any cell type are shown in the figure.Significant GPs (FDR < 0.05 and effect size > 0) for each cell type are represented as circles, with circle size indicating statistical significance and color representing effect size (capped at 10).The top-weighted genes in each upregulated GP, which are co-expressed genes that may be related to cell type identity or cellular response to cytokines, are shown on the right in blue.0.0 0.5 1.0   Tgd-e IL-2,-7,-15,-18,-36α Eif5a, Ppa1, Hspa5, Igtp, Ifi47       pDC-e IL-36α Bcl2a1d, Bcl2a1b, Il4i1, Tmem176a, Wfdc21 Most highly expressed genes in each subcluster  Ifitm3  Top genes in GP   MigDC-e GM-CSF; IL-18 Cish, Jaml, Nr4a3, Pkib, Coro2a LC-e GM-CSF; IL-3,-18; TSLP Plek2, Gpbp1, Pkib, Nr4a3, Tcaf2 NK cells in each cytokine treatment h Expression of selected IL-18 regulated genes in NK cells with comparisons to other cytokine treatments 12 | Further characterization of IL-18 induced NK-f polarization state and IL-1α/β induced NK-c polarization state in NK cells.See next page for caption.Supplementary Fig. 13 | Cell type properties in the cell-cell interactome.Properties of Extended Data Fig. 11.For each cell type (a circle), showing the number of cytokines expressed and the total number of cell types targeted through any cytokine.Size of the circle is proportional to the number of paths in the cell-cell interactome.cell-cell interactome Supplementary Fig. 14 | Additional IREA software output on human COVID-19 blood samples.a-b.IREA analysis on peripheral blood cells collected from severe COVID-19 patients relative to healthy volunteers.IREA compass plots showing enrichment scores for each of 86 cytokines in CD8+ T cells and CD4+ T cells in ventilated COVID-19 patients relative to healthy controls.Bar length represents enrichment score, shade represents FDR adjusted Pvalue (two-sided Wilcoxon rank-sum test), with darker colors representing more significant enrichment (red: enriched in ventilated COVID-19 patients, blue: enriched in healthy control).Cytokines with receptors expressed are indicated by black filled boxes.

2 | CD4+ T cell responses to cytokines: polarization states, subclusters, marker gene expression, and gene programs. a,
Top, UMAP visualization of CD4+ T cells for all cytokines, colored by polarization states; bottom, table with cell type polarization states (left column), single cytokine drivers (middle column), and top marker genes (right column); reproduced from Fig.3for ease of reference.b, Pairwise Pearson correlation coefficients between polarization states.c, UMAP visualization of CD4+ T cells shown independently for each cytokine treatment, colored by cytokine treatment (blue) or PBS treatment control (gray).d, UMAP visualization of CD4+ T cells for all cytokine or PBS treatment; cells colored for CD4+ T cell Louvain subclusters.e, Top overexpressed genes in each Louvain subcluster in d; color, columnscaled average expression; size of circle, percentage of cells in the subcluster expressing each gene.f-i, CD4+ T cell responses to each cytokine stimulation.f, Fraction of cells per subcluster in each cytokine treatment.Colors represent subclusters defined in d. g, Enrichment of each subcluster in each cytokine treatment; size of circle, Bonferroni-adjusted P-value of hypergeometric test relative to PBS; black fills, P < 0.01.h, Row-normalized relative expression of representative marker genes of each polarization state in cytokine-treated vs. PBS-treated cells.i, Enrichment of CD4+ T cell gene programs obtained from NMF analysis of all CD4+ T cells in cytokine-treated cells relative to PBS-treated cells; Supplementary Fig. 2 | CD4+ T cell responses to cytokines: polarization states, subclusters, marker gene expression, and gene programs.See next page for caption.Supplementary Fig. size of circle, FDR-adjusted P-value from two-sided Wilcoxon rank-sum test; shade, effect size representing the mean difference in gene program weight.j, Top weighted genes in each gene program in i. k, Average gene program weight in each subcluster.Rows and columns were hierarchically clustered using the complete-linkage method on Euclidean distances.0.0 0.5 1.0 Correlation Subcluster S01

Supplementary Fig. 3 | CD8+ T cell responses to cytokines: polarization states, subclusters, marker gene expression, and gene programs. See next page for caption. Supplementary Fig. 3 | CD8+ T cell responses to cytokines: polarization states, subclusters, marker gene expression, and gene programs. a, Top
, UMAP visualization of CD8+ T cells for all cytokines, colored by polarization states; bottom, table with cell type polarization states (left column), single cytokine drivers (middle column), and top marker genes (right column); reproduced from Fig.3for ease of reference.b, Pairwise Pearson correlation coefficients between polarization states.c, UMAP visualization of CD8+ T cells shown independently for each cytokine treatment, colored by cytokine treatment (blue) or PBS treatment control (gray).d, UMAP visualization of CD8+ T cells for all cytokine or PBS treatment; cells colored for CD8+ T cell Louvain subclusters.e, Top overexpressed genes in each Louvain subcluster in d; color, columnscaled average expression; size of circle, percentage of cells in the subcluster expressing each gene.f-i, CD8+ T cell responses to each cytokine stimulation.f, Fraction of cells per subcluster in each cytokine treatment.Colors represent subclusters defined in d. g, Enrichment of each subcluster in each cytokine treatment; size of circle, Bonferroni-adjusted P-value of hypergeometric test relative to PBS; black fills, P < 0.01.h, Row-normalized relative expression of representative marker genes of each polarization state in cytokine-treated vs. PBS-treated cells.i, Enrichment of CD8+ T cell gene programs obtained from NMF analysis of all CD8+ T cells in cytokine-treated cells relative to PBS-treated cells; size of circle, FDR-adjusted P-value from two-sided Wilcoxon rank-sum test; shade, effect size representing the mean difference in gene program weight.j, Top weighted genes in each gene program in i. k, Average gene program weight in each subcluster.Rows and columns were hierarchically clustered using the complete-linkage method on Euclidean distances.

4 | γδ T cell responses to cytokines: polarization states, subclusters, marker gene expression, and gene programs. a,
Top, UMAP visualization of γδ T cells for all cytokines, colored by polarization states; bottom, table with cell type polarization states (left column), single cytokine drivers (middle column), and top marker genes (right column); reproduced from Fig.3for ease of reference.b, Pairwise Pearson correlation coefficients between polarization states.c, UMAP visualization of γδ T cells shown independently for each cytokine treatment, colored by cytokine treatment (blue) or PBS treatment control (gray).d, UMAP visualization of γδ T cells for all cytokine or PBS treatment; cells colored for γδ T cell Louvain subclusters.e, Top overexpressed genes in each Louvain subcluster in d; color, column-scaled average expression; size of circle, percentage of cells in the subcluster expressing each gene.f-i, γδ T cell responses to each cytokine stimulation.f, Fraction of cells per subcluster in each cytokine treatment.Colors represent subclusters defined in d. g, Enrichment of each subcluster in each cytokine treatment; size of circle, Bonferroni-adjusted P-value of hypergeometric test relative to PBS; black fills, P < 0.01.h, Row-normalized relative expression of representative marker genes of each polarization state in cytokine-treated vs. PBS-treated cells.i, Enrichment of γδ T cell gene programs obtained from NMF analysis of all γδ T cell in cytokine-treated cells relative to PBS-treated cells; size of circle, FDR- 1 2 3 4 Supplementary Fig. 4 | γδ T cell responses to cytokines: polarization states, subclusters, marker gene expression, and gene programs.See next page for caption.Supplementary Fig. adjusted P-value from two-sided Wilcoxon rank-sum test; shade, effect size representing the mean difference in gene program weight.j,Top weighted genes in each gene program in i. k, Average gene program weight in each subcluster.Rows and columns were hierarchically clustered using the complete-linkage method on Euclidean distances.

5 | Treg responses to cytokines: polarization states, subclusters, marker gene expression, and gene programs. a,
Top, UMAP visualization of Tregs for all cytokines, colored by polarization states; bottom, table with cell type polarization states (left column), single cytokine drivers (middle column), and top marker genes (right column); reproduced from Fig.3for ease of reference.b, Pairwise Pearson correlation coefficients between polarization states.c, UMAP visualization of Tregs shown independently for each cytokine treatment, colored by cytokine treatment (blue) or PBS treatment control (gray).d, UMAP visualization of Tregs for all cytokine or PBS treatment; cells colored for Treg Louvain subclusters.e, Top overexpressed genes in each Louvain subcluster in d; color, column-scaled average expression; size of circle, percentage of cells in the subcluster expressing each gene.f-i, Treg responses to each cytokine stimulation.f, Fraction of cells per subcluster in each cytokine treatment.Colors represent subclusters defined in d. g, Enrichment of each subcluster in each cytokine treatment; size of circle, Bonferroni-adjusted P-value of hypergeometric test relative to PBS; black fills, P < 0.01.h, Row-normalized relative expression of representative marker genes of each polarization state in cytokine-treated vs. PBS-treated cells.i, Enrichment of Treg gene programs obtained from NMF analysis of all Tregs in cytokine-treated cells relative to PBS-treated cells; size of circle, FDR-adjusted P-value from twosided Wilcoxon rank-sum test; shade, effect size representing the mean difference in gene program weight.j, Top weighted genes in each gene program in i. k, Average gene program weight in each subcluster.Rows and columns were hierarchically clustered using the complete-linkage method on Euclidean distances.

6 | pDC responses to cytokines: polarization states, subclusters, marker gene expression, and gene programs. a,
Top, UMAP visualization of pDCs for all cytokines, colored by polarization states; bottom, table with cell type polarization states (left column), single cytokine drivers (middle column), and top marker genes (right column); reproduced from Fig.3for ease of reference.b, Pairwise Pearson correlation coefficients between polarization states.c, UMAP visualization of pDCs shown independently for each cytokine treatment, colored by cytokine treatment (blue) or PBS treatment control (gray).d, UMAP visualization of pDCs for all cytokine or PBS treatment; cells colored for pDC Louvain subclusters.e, Top overexpressed genes in each Louvain subcluster in d; color, column-scaled average expression; size of circle, percentage of cells in the subcluster expressing each gene.f-i, pDC responses to each cytokine stimulation.f, Fraction of cells per subcluster in each cytokine treatment.Colors represent subclusters defined in d. g, Enrichment of each subcluster in each cytokine treatment; size of circle, Bonferroni-adjusted P-value of hypergeometric test relative to PBS; black fills, P < 0.01.h, Row-normalized relative expression of representative marker genes of each polarization state in cytokine-treated vs. PBS-treated cells.i, Enrichment of pDC gene programs obtained from NMF analysis of all pDCs in cytokine-treated cells relative to PBS-treated cells; size of circle, FDR-adjusted P-value from twosided Wilcoxon rank-sum test; shade, effect size representing the mean difference in gene program weight.j, Top weighted genes in each gene program in i. k, Average gene program weight in each subcluster.Rows and columns were hierarchically clustered using the complete-linkage method on Euclidean distances.
1 2 3 Supplementary Fig. 6 | pDC responses to cytokines: polarization states, subclusters, marker gene expression, and gene programs.See next page for caption.Supplementary Fig.