Immune landscape of isocitrate dehydrogenase-stratified primary and recurrent human gliomas

Abstract

Background: Human gliomas are classified using isocitrate dehydrogenase (IDH) status as a prognosticator; however, the influence of genetic differences and treatment effects on ensuing immunity remains unclear.

Methods: In this study, we used sequential single-cell transcriptomics on 144 678 and spectral cytometry on over 2 million immune cells encompassing 48 human gliomas to decipher their immune landscape.

Results: We identified 22 distinct immune cell types that contribute to glioma immunity. Specifically, brain-resident microglia (MG) were reduced with a concomitant increase in CD8+ T lymphocytes during glioma recurrence independent of IDH status. In contrast, IDH-wild type-associated patterns, such as an abundance of antigen-presenting cell-like MG and cytotoxic CD8+ T cells, were observed. Beyond elucidating the differences in IDH, relapse, and treatment-associated immunity, we discovered novel inflammatory MG subpopulations expressing granulysin, a cytotoxic peptide that is otherwise expressed in lymphocytes only. Furthermore, we provide a robust genomic framework for defining macrophage polarization beyond M1/M2 paradigm and reference signatures of glioma-specific tumor immune microenvironment (termed GlioTIME-36) for deconvoluting transcriptomic datasets.

Conclusions: This study provides advanced optics of the human pan-glioma immune contexture as a valuable guide for translational and clinical applications.

Keywords: glioma; isocitrate dehydrogenase; microglia; tumor immune microenvironment.

Graphical Abstract

Figure 1.

Single-cell transcriptional landscape of the glioma TIME. (A) A scheme depicting the experimental workflow spanning sample preparation of resected brain tissues/tumors, scRNA-seq, spectral flow cytometry validation, and computational analysis. For scRNA-seq analyses, patients were stratified as the NGB (n = 3), IMP (n = 4), IMR (n = 6), IWP (n = 4), and IWR (n = 4) groups as a discovery patient cohort (hereafter collectively referred as glioma subtypes) (see details in Supplementary Table 1). The dissociated and CD45-APC-stained cells were sorted to obtain pure CD45⁺ tumor-associated leukocytes. Subsequently, matched scRNA-seq and bulk RNA-seq were performed followed by computational analysis. Cell types defined by scRNA-seq were further validated by a 40-plex protein marker spectral flow cytometry panel in an extended validation patient cohort comprising NGB (n = 3), IMP (n = 14), IMR (n = 9), IWP (n = 13), and IWR (n = 12). (B–D) UMAP visualization of unsupervised clustering analysis of leukocytes (n = 144 678) that passed quality filtering as shown in (B) all immune cells; (C) myeloid lineage clusters (n = 100 587) identified as microglia (MG), microglial-like cells (MG-like), interferon gene module-associated myeloid cells (Myeloid_IFNs), proliferative myeloid cells (Myeloid_Prolif), classical monocytes (c-Mo), nonclassical monocytes (nc-Mo), macrophages (MAC), monocyte-derived MAC (MDM), conventional dendritic cell 1 (cDC1), conventional dendritic cell 2 (cDC2), plasmacytoid DC (pDC), and neutrophil; (D) lymphoid lineage clusters (n = 44 091) identified as CD4 T, CD8 T, interferon-stimulated gene associated T (T_ISG), proliferative T (T_Prolif), gamma-delta T (γδ-T) lymphocytes, mucosal-associated invariant T (MAIT) cells, natural killer (NK), natural killer T (NKT), B lymphocytes (B cells), and plasma cells. Cells are color-coded for their inferred cell types (left) and the glioma subtypes of their corresponding tumors (middle). Stacked bar plots (right) show the percentage cell composition of each glioma subtype.

Figure 2.

Immune contexture of human gliomas validated by spectral flow cytometry. Color-coded scatter bar plots represent the relative proportions and error bars represent an interval of maximum and minimum values ± standard deviation (SD) of the indicated immune cell types (of CD45⁺ leukocytes) across glioma subtypes: NGB (n = 3), IMP (n = 14), IMR (n = 9), IWP (n = 13), and IWR (n = 12). (A) P2RY12⁺CX3CR1⁺ MG (left), CD11c⁺ MG (right). (B) T lymphocytes (gated on CD3⁺ ΤCR-γδ⁻ ΤCRV-α7.2⁻) (left), CD8⁺ T cell (gated on T) (right). (C) NK cells (gated on CD56⁺) (left), CD56^lo NK (right). Statistical significance in A–C was determined by using the Kruskal–Wallis test at *P < .05, **P < .01 between NGB and glioma subtypes: IMP versus IMR, IWP versus IWR, and IMP versus IWP; n.s. = not significant. (D) Pie chart summarizing the proportion of indicated immune cell types in NGB and across glioma subtypes in TIME. (E) Circos plots showing the multiple correlation matrix between frequencies of leukocytes in human gliomas. P < .05. R > 0.6 is represented in red, and R < −0.6 is shown in green.

Figure 3.

Heterogeneity of MG, non-MG myeloid cells, DCs, and their functional states in the glioma TIME. (A) UMAP visualization of unsupervised clustering analysis of MG cell states (n = 63 332). Cells are color-coded for their inferred MG cell states (top) and their distribution across glioma subtypes (bottom). (B) Bubble plot showing the scaled expression (indicated by the color of the circle) and percentage of expression (indicated by the size of the circle) of selected inflammatory marker genes from the top 50 differentially expressed genes (DEGs) of inflammatory MG clusters. NGB (n = 3), IMP (n = 4), IMR (n = 6), IWP (n = 4), and IWR (n = 4). (C) Plot showing P-value and color-coded chi-squared test-derived standardized residuals as index of the statistical significance and strength of the association between GNLY gene expression and myeloid cell types from Antunes et al. (2021) in primary (top) and recurrent GBM (bottom). Red and blue indicate positive and negative enrichment for expression of GNLY in each myeloid cell type. P < 2.2e-16, χ² test. (D) Representative example of multiplex tissue images from brain tumor samples of GBM. Composite and unmixed images (20× and 40× magnification, scale bars represent 50 and 10 µm, respectively) showing a panoramic view of P2RY12 (green), GNLY (red), and TREM2 (white) and DAPI (4ʹ,6-diamidino-2-phenylindole) and high magnification showing cells co-expressing P2RY12, GNLY, and TREM2 and DAPI and individual marker expression. (E) Flow cytometry (FCM) dot plot showing the expression of TREM2 and GNLY on MG (in red), and GNLY fluorescence minus one (FMO) control (gray) (left). Corresponding color-coded scatter bar plots showing proportions of TREM2⁺ and inflammatory TREM2⁺GNLY⁺ cells (right). (F) FCM contour plot showing the expression of CD86 and HLA-DR on MG (left). Corresponding color-coded scatter bar plots showing proportions of HLA-DR⁺ and CD86⁺ HLA-DR⁺ cells (right). (G) UMAP visualization of unsupervised clustering analysis of non-MG myeloid cell states (n = 37 255). Cells are color-coded for their inferred cell states (top), and their distribution across glioma subtypes (bottom). (H) FCM contour plot showing the expression of LYZ and CD163 on MAC (top, left) and MDM (bottom, left). Color-coded scatter bar plots depicting the proportions of MAC (top, right) and MDM (bottom, right) expressing CD163⁺LYZ⁺. (I) UMAP visualization of unsupervised clustering analysis of cDC2 cells (n = 6066). Cells are color-coded for their inferred cell states (left), and their distribution across glioma subtypes (right). (J) Bubble plot showing the scaled expression (indicated by the color of the circle) and percentage of expression (indicated by the size of the circle) of selected cluster-specific genes from the top 50 DEGs of cDC2-associated clusters. (K) FCM pseudo-color plots showing the expression of CD1c and Clec9A (gated on CD11c⁺ HLA-DR⁺ cells, left). FCM pseudo-color plots showing the expression of CCR2 on Mo-cDC1 and Mo-cDC2 cells (middle). Corresponding color-coded scatter bar plots showing the proportions of bona fide Clec9A⁺CCR2⁻ b-cDC1 and CD1c⁺CCR2⁻ b-cDC2, and monocytic derivatives CCR2⁺Clec9A⁺ Mo-cDC1 and CCR2⁺CD1c⁺ Mo-cDC2 subsets (right). b-cDC and Mo-cDCs were gated on CD1c⁺ and Clec9A⁺ cells. (L) Heatmaps obtained from normalized APC score derived from median flourescence intensities (MFIs) of B2M, HLA-DR, CD86, and LYZ in indicated myeloid cell types. Statistical significance in D, E, G, and J was determined by using the Kruskal–Wallis test at *P < .05, **P < .01 between NGB and glioma subtypes, IMP versus IMR, IWP versus IWR, and IMP versus IWP; n.s. = not significant. Glioma subtypes used for spectral cytometry analysis: NGB (n = 3), IMP (n = 14), IMR (n = 9), IWP (n = 13), and IWR (n = 12).

Figure 3.

Heterogeneity of MG, non-MG myeloid cells, DCs, and their functional states in the glioma TIME. (A) UMAP visualization of unsupervised clustering analysis of MG cell states (n = 63 332). Cells are color-coded for their inferred MG cell states (top) and their distribution across glioma subtypes (bottom). (B) Bubble plot showing the scaled expression (indicated by the color of the circle) and percentage of expression (indicated by the size of the circle) of selected inflammatory marker genes from the top 50 differentially expressed genes (DEGs) of inflammatory MG clusters. NGB (n = 3), IMP (n = 4), IMR (n = 6), IWP (n = 4), and IWR (n = 4). (C) Plot showing P-value and color-coded chi-squared test-derived standardized residuals as index of the statistical significance and strength of the association between GNLY gene expression and myeloid cell types from Antunes et al. (2021) in primary (top) and recurrent GBM (bottom). Red and blue indicate positive and negative enrichment for expression of GNLY in each myeloid cell type. P < 2.2e-16, χ² test. (D) Representative example of multiplex tissue images from brain tumor samples of GBM. Composite and unmixed images (20× and 40× magnification, scale bars represent 50 and 10 µm, respectively) showing a panoramic view of P2RY12 (green), GNLY (red), and TREM2 (white) and DAPI (4ʹ,6-diamidino-2-phenylindole) and high magnification showing cells co-expressing P2RY12, GNLY, and TREM2 and DAPI and individual marker expression. (E) Flow cytometry (FCM) dot plot showing the expression of TREM2 and GNLY on MG (in red), and GNLY fluorescence minus one (FMO) control (gray) (left). Corresponding color-coded scatter bar plots showing proportions of TREM2⁺ and inflammatory TREM2⁺GNLY⁺ cells (right). (F) FCM contour plot showing the expression of CD86 and HLA-DR on MG (left). Corresponding color-coded scatter bar plots showing proportions of HLA-DR⁺ and CD86⁺ HLA-DR⁺ cells (right). (G) UMAP visualization of unsupervised clustering analysis of non-MG myeloid cell states (n = 37 255). Cells are color-coded for their inferred cell states (top), and their distribution across glioma subtypes (bottom). (H) FCM contour plot showing the expression of LYZ and CD163 on MAC (top, left) and MDM (bottom, left). Color-coded scatter bar plots depicting the proportions of MAC (top, right) and MDM (bottom, right) expressing CD163⁺LYZ⁺. (I) UMAP visualization of unsupervised clustering analysis of cDC2 cells (n = 6066). Cells are color-coded for their inferred cell states (left), and their distribution across glioma subtypes (right). (J) Bubble plot showing the scaled expression (indicated by the color of the circle) and percentage of expression (indicated by the size of the circle) of selected cluster-specific genes from the top 50 DEGs of cDC2-associated clusters. (K) FCM pseudo-color plots showing the expression of CD1c and Clec9A (gated on CD11c⁺ HLA-DR⁺ cells, left). FCM pseudo-color plots showing the expression of CCR2 on Mo-cDC1 and Mo-cDC2 cells (middle). Corresponding color-coded scatter bar plots showing the proportions of bona fide Clec9A⁺CCR2⁻ b-cDC1 and CD1c⁺CCR2⁻ b-cDC2, and monocytic derivatives CCR2⁺Clec9A⁺ Mo-cDC1 and CCR2⁺CD1c⁺ Mo-cDC2 subsets (right). b-cDC and Mo-cDCs were gated on CD1c⁺ and Clec9A⁺ cells. (L) Heatmaps obtained from normalized APC score derived from median flourescence intensities (MFIs) of B2M, HLA-DR, CD86, and LYZ in indicated myeloid cell types. Statistical significance in D, E, G, and J was determined by using the Kruskal–Wallis test at *P < .05, **P < .01 between NGB and glioma subtypes, IMP versus IMR, IWP versus IWR, and IMP versus IWP; n.s. = not significant. Glioma subtypes used for spectral cytometry analysis: NGB (n = 3), IMP (n = 14), IMR (n = 9), IWP (n = 13), and IWR (n = 12).

Figure 4.

Diversity of glioma-associated cytotoxic lymphoid cells. (A) UMAP visualization of unsupervised clustering analysis of CD8 T cell clusters defined by scRNA-seq. Cells are color-coded for their inferred CD8 T lymphocyte subpopulations (left) and the glioma subtypes of their corresponding tumors (right). (B) Color-coded scatter bar plots represent the proportions of activated and effector CD8 T cell subsets as defined by Ki-67⁺, GZMB⁺, GNLY⁺, and GZMB⁺GNLY⁺ cells (top panels), and PRF1⁺GZMB⁺, PRF1⁺GZMB⁺GNLY⁺ and PRF1⁺GZMB⁻GNLY⁻ (bottom panels) cells. (C) UMAP visualization of unsupervised clustering analysis of NK cell clusters defined by scRNA-seq. Cells are color-coded for their inferred NK cell subpopulations (left), and the glioma subtypes of their corresponding tumors (right). (D) Color-coded scatter bar plots represent the proportions of activated and effector NK cell subsets as defined by Ki-67⁺, GZMB⁺, GNLY⁺, and GZMB⁺GNLY⁺ cells (top panels), and PRF1⁺GZMB⁺, PRF1⁺GZMB⁺GNLY⁺, and PRF1⁺GZMB^-GNLY⁻ cells (bottom panels) expression. (E) Pseudo-color FCM plots showing the expression of Ki-67, GZMB, and GNLY on NKT cells, and corresponding color-coded scatter bar plots represent the proportions of activated and effector NKT cell subsets as defined by Ki-67⁺, GZMB⁺, GNLY⁺, and GZMB⁺GNLY⁺ cells. (F) Pseudo-color FCM plots showing the expression of Ki-67, GZMB, and GNLY on MAIT cells, and corresponding color-coded scatter bar plots represent the proportions of activated and effector MAIT cell subsets as defined by Ki-67⁺, GZMB⁺, GNLY⁺, and GZMB⁺GNLY⁺ cells. Statistical significance in B, D, E, and F was determined by using the Kruskal–Wallis test at *P < .05, **P < .01, and ***P < .001 between NGB and glioma subtypes, IMP versus IMR, IWP versus IWR, and IMP versus IWP; n.s. = not significant. Glioma subtypes used for spectral cytometry analysis: NGB (n = 3), IMP (n = 14), IMR (n = 9), IWP (n = 13), and IWR (n = 12).

Figure 5.

Transcriptomic utility signatures defined by scRNA-seq for characterizing spectral polarization modules of glioma-associated MG/MAC and defining GlioTIME-36 gene signatures to infer leukocyte abundance. (A) Circos plot showing overrepresented stimulus-specific polarization gene expression modules for MG (upper) and MAC/MDM (lower) as pseudo-bulk versus respective multispectral macrophage polarization references. Bar heights and color indicate statistical significance (P < .05–.0001). Gene expression modules are represented as glucocorticoid (GC), interferon-gamma (IFN-γ), lauric acid (LA), linoleic acid (LiA), oleic acid (OA), palmitic acid (PA), prostaglandin E2 (PGE2), standard lipopolysaccharide (sLPS), tumor necrosis alpha (TNF-α), TPP (TNF-α + PGE2 + Pam3CysSerLys4), and interleukin-4 (IL-4). (B) Circos plot (upper panel) showing overrepresented stimulus-specific polarization gene expression modules for MG cluster (MG_inflam 1) defined by scRNA-seq versus respective multispectral macrophage polarization references. Bubble plot (lower panel) showing gene set enrichment analysis of MG and MAC/MDM subsets against stimulus-specific polarization gene expression modules. Each bubble represents a gene expression module. Bubble size corresponds to gene ratio, and color indicates statistical significance. Only modules with ≥5 overlapping genes and an adjusted P-value < .05 are shown here. (C) A schematic of calculating the normalized enrichment score (NES) of the GlioTIME-36 and LM22 signatures derived from bulk mRNA-seq datasets of CD45⁺ leukocytes (immune pure) from glioma subtypes; NGB (n = 1), IMP (n = 3), IMR (n = 4), IWP (n = 3), and IWR (n = 4), and GLASS dataset (n = 299). (D) Heatmap of GlioTIME-36 with differential NES activity in different glioma subtypes. Statistical significance was determined using the MWW test at P < .05 and effect size > 0.3. The top annotation tracks define glioma subtypes, histology, and 1p 19q codeletion status. (E) Box plots represent the mean NES ± SD of indicated tumor-associated leukocyte subpopulation across glioma subtypes inferred from the GlioTIME-36 signatures in the GLASS dataset. P-values derived by the Kruskal–Wallis H test are shown on the x axis. P-values on the top were derived from post hoc correction by Nemenyi’s test for multiple subgroups comparison of continuous variables.

Figure 5.

Transcriptomic utility signatures defined by scRNA-seq for characterizing spectral polarization modules of glioma-associated MG/MAC and defining GlioTIME-36 gene signatures to infer leukocyte abundance. (A) Circos plot showing overrepresented stimulus-specific polarization gene expression modules for MG (upper) and MAC/MDM (lower) as pseudo-bulk versus respective multispectral macrophage polarization references. Bar heights and color indicate statistical significance (P < .05–.0001). Gene expression modules are represented as glucocorticoid (GC), interferon-gamma (IFN-γ), lauric acid (LA), linoleic acid (LiA), oleic acid (OA), palmitic acid (PA), prostaglandin E2 (PGE2), standard lipopolysaccharide (sLPS), tumor necrosis alpha (TNF-α), TPP (TNF-α + PGE2 + Pam3CysSerLys4), and interleukin-4 (IL-4). (B) Circos plot (upper panel) showing overrepresented stimulus-specific polarization gene expression modules for MG cluster (MG_inflam 1) defined by scRNA-seq versus respective multispectral macrophage polarization references. Bubble plot (lower panel) showing gene set enrichment analysis of MG and MAC/MDM subsets against stimulus-specific polarization gene expression modules. Each bubble represents a gene expression module. Bubble size corresponds to gene ratio, and color indicates statistical significance. Only modules with ≥5 overlapping genes and an adjusted P-value < .05 are shown here. (C) A schematic of calculating the normalized enrichment score (NES) of the GlioTIME-36 and LM22 signatures derived from bulk mRNA-seq datasets of CD45⁺ leukocytes (immune pure) from glioma subtypes; NGB (n = 1), IMP (n = 3), IMR (n = 4), IWP (n = 3), and IWR (n = 4), and GLASS dataset (n = 299). (D) Heatmap of GlioTIME-36 with differential NES activity in different glioma subtypes. Statistical significance was determined using the MWW test at P < .05 and effect size > 0.3. The top annotation tracks define glioma subtypes, histology, and 1p 19q codeletion status. (E) Box plots represent the mean NES ± SD of indicated tumor-associated leukocyte subpopulation across glioma subtypes inferred from the GlioTIME-36 signatures in the GLASS dataset. P-values derived by the Kruskal–Wallis H test are shown on the x axis. P-values on the top were derived from post hoc correction by Nemenyi’s test for multiple subgroups comparison of continuous variables.