Immunologic Profiling of Mutational and Transcriptional Subgroups in Pediatric and Adult High-Grade Gliomas

High-grade gliomas are defined by aggressive and infiltrative growth, reducing the life expectancy in pediatric and adult patients (1). Despite multimodal treatment, including neurosurgical resection and combination of radiotherapy and chemotherapy, tumor recurrence is a certainty (1). Although treatment modalities such as immunotherapeutic approaches with checkpoint inhibitors have failed to show efficacy in a phase III study in recurrent glioblastomas (2, 3), ongoing studies are trying to harness the ability of the immune system to counteract tumor growth and induce long-term remission (4). To identify which patients might benefit from checkpoint inhibition or other immunotherapeutic approaches, such as tumor vaccines, chimeric antigen receptor (CAR) T-cell therapy, autologous T-cell transfer, or cytolytic virus injections, we must better understand the glioma-specific immune microenvironment (5). Large-scale immunologic profiling of gene-expression data in other cancer entities has yielded information to improve precision immunotherapy (6–8). Unfortunately, the intricate cellular, genetic, and transcriptional heterogeneity of malignant gliomas adds complexity to efforts to define targetable immunologic profiles (9, 10).

The 2016 World Health Organization (WHO) classification identified subcategories among high-grade gliomas (11). Molecular profiling, especially by isocitrate dehydrogenase (IDH) gene mutation analysis (12), improves understanding of subgroup-specific aggressiveness and prediction of overall survival (OS) in adult patients. Transcriptomic studies have identified additional subgroups. Phillips and colleagues described different molecular subgroups of astrocytoma and their prognostic relevance, and then Verhaak and colleagues refined this classification by introducing a gene-expression–based molecular classification of glioblastoma that stratifies tumors into mesenchymal, classic, neural, and proneural subgroups (13, 14). These groups share features of biological aggressiveness and response to therapy (13). The proneural subgroup was further subdivided into cGIMP and non-cGIMP methylated tumors (15). The neural subtype has been excluded in many subsequent studies due to suspected contamination of tumor samples with normal tissue (14, 16).

In addition to the groups defined by transcriptional subtypes or mutations in the IDH gene in adult tumors, genomic sequencing studies also identified lineage-defining aberrations, such as the somatic mutations of the histone 3 variants in the positions G34 and K27 (17–19). Furthermore, pediatric glioblastoma can be stratified by additional molecular subgroups according to their RTK or MYCN pathway aberrations (19). The K27 and G34 mutations in histone 3 were predominantly found in pediatric and adolescent glioblastoma, highlighting their biological distinctiveness and the need to stratify high-grade gliomas according to their mutational profile.

Accumulating data suggest that molecularly distinct glioma subgroups differ in their tumor microenvironment and their immunostromal profiles (20, 21). Investigations of The Cancer Genome Atlas (TCGA) data set, predominantly including IDH wild-type (WT) glioblastomas, revealed heterogeneity in the composition of the immune environment between mutational and transcriptional subgroups (22). Wang and colleagues demonstrated a subgroup-dependent macrophage/microglia infiltration according to the gene signature, with increased infiltration in NF1-mutated tumors (23). This was also confirmed by gene-enrichment analysis for transcriptional subtypes by Doucette and colleagues (24) and by our group by gene-expression profiling for different molecular subgroups in medulloblastoma (25).

Our analysis evaluates the subtype-specific immune microenvironment across a cohort of pediatric and adult high-grade gliomas, comparing mutational and transcriptional profiles in order to identify immune infiltration patterns and functional pathways mediating local immunosuppression.

All data analyses were performed using the statistical programming language R (26), including the packages survival, caret, randomForest, igraph, gplots, and Rtsne (https://CRAN.R-project.org/package=survival, https://CRAN.R-project.org/package=caret, http://igraph.org, https://CRAN.R-project.org/package=gplots, http://www.jmlr.org/papers/v9/vandermaaten08a.html, and https://cran.r-project.org/package=randomForest).

Raw gene-expression data (.CEL files) from high-grade gliomas including eight series on the Affymetrix Human Genome U133 Plus 2.0 Array and one series on the Affymetrix HT_HG-U133A array were downloaded from Gene Expression Omnibus (GEO; Fig. 1A; Supplementary Table S1; ref. 27).

To minimize batch effects, the fRMA normalization, which includes a batch-effect correction, was used as implemented in the R-package fRMA (28) using the input vectors from the R-packages hthgu133ahsentrezgfrmavecs and hgu133plus2-hsentrezgfrmavecs (http://bioconductor.org/packages/release/data/annotation/html/hthgu-133afrmavecs.html, and http://bioconductor.org/packages/release/data/annotation/html/hgu133-plus2frmavecs.html) and the custom CDF files hthgu133ahsentrezgcdf and hgu133plus2hsentrezgcdf (version 19; ref. 29). All data originating from the same center/study were transformed to median 0, as there were still remaining batch effects (Supplementary Fig. S1A). Because the data were not homogeneous with respect to mutational profile, localization, and grade, no adjustment of the variability of each gene was performed between series measured on the same microarray. However, the mean absolute deviation, which is less influenced by outliers than the standard deviation, was aligned between the two microarray platforms for each gene separately. After batch-effect removal, clustering of the data revealed concordance of similar tumors (with respect to diagnosis, mutational profile, or transcriptional subgroup) from different data sets.

Duplicates defined by a Pearson correlation > 0.999 were removed, yielding a final data set of 79 anaplastic astrocytomas (AA), 28 diffuse intrinsic pontine gliomas (DIPG), and 1,028 glioblastomas. Analyses involving both arrays were restricted to the 12,061 common genes.

Regarding the information on the IDH subgroup, the study by Gravendeel and colleagues performed sequencing of the IDH1 locus only (30), whereas TCGA data sets included sequencing information of the complete IDH gene (13, 31).

To obtain specific signatures of tumor-infiltrating immune and stromal cells, we first defined candidate genes from the signatures published by Becht and colleagues (32) and Danaher and colleagues (33) for 10 microenvironmental cell populations (MCP): T cells, CD8⁺ T cells, cytotoxic lymphocytes, B lineage cells, natural killer (NK) cells, monocytic lineage cells, myeloid dendritic cells, neutrophils, fibroblast-like cells, and endothelial cells (Supplementary Tables S2 and S3). Second, all candidate genes were manually reviewed, and those likely to overlap with glioma signatures were removed (Supplementary Table S2). In a final step, a correlation network was built by linking all genes with Pearson correlation over 0.4 for each signature separately. Only those genes belonging to the largest correlation module were retained as final candidate genes, because uncorrelated genes are unlikely to be specifically expressed by a unique cell population (33).

Scores of tumor-infiltrating immune and stromal cell signatures were defined as the average values of the corresponding marker genes. Signatures for different immunologic processes were obtained from Doucette and colleagues (24) and Thorsson and colleagues (34) and computed analogously (Supplementary Table S3). This analysis is based on gene-expression data, and no histologic or flow-cytometric analysis was performed. Conclusions of the presence of different immune cells were drawn based on expression of genes or gene signatures.

Inference of the mutational status (IDH, K27, G34, and WT) and the transcriptional group (classic, mesenchymal, proneural, and neural) were performed with a random forest classifier coupled with a confidence estimation using the R-package randomForest. The classifier was built on the 1,000 most variable probes in the training sets with 1,000 trees. Default settings were used for the remaining parameters. Selection of the most variable genes was based on the standard deviation computed for each training fold separately.

As the frequency of IDH, K27, and G34 mutations was age dependent in our data set, separate classifiers were built for young patients [reference set: all WT (n = 42), K27 (n = 36), and G34 (n = 13) patients under 20 years, as well as 20 IDH-mutated patients under 30 years] and older patients [reference set: all IDH-mutated (n = 79) and all WT (n = 485) patients over 20 years]. Patients with a known transcriptional subgroup from the TCGA cohort were chosen as reference for the transcriptional subgroup classifier. To consider imbalanced class sizes, each group size was reduced to at most twice the number of samples of the smallest class for the WT/K27/G34/IDH classifier by downsampling. Group size was reduced to the size of the smallest class for the two other classifiers.

Estimates of the classification accuracy were obtained from 5-fold cross-validation. As no parameter optimization was performed, nested cross-validation was not required. Because raw classification results yielded accuracies under 90% for each task, we used tail probabilities to increase accuracy, and samples with a low proportion of votes for the predicted class were excluded and considered unclassifiable. To this end, the lowest threshold value was determined for each classification task such that samples from the validation fold with prediction score over the threshold were correctly assigned with an accuracy of over 91%. This allowed increasing accuracies over 91% for each of the classification tasks. The random forest classifier was compared with a radial basis function kernel support vector, which did not yield superior classification accuracy in this setting.

Data were trimmed at z-scores ± 3 prior to heat map analysis. Heat maps were generated using the average-linkage method and Pearson correlations as similarity measure. Significance of differential expression between subgroups was assessed using the Kruskal–Wallis test. Multiple testing correction was performed with the Benjamini–Hochberg (BH) method.

Immunostromal clusters were defined using k-means clustering of the previously defined marker genes. The number of clusters (n = 6) was selected using the elbow method. For visualization, t-distributed stochastic neighbor embedding (t-SNE) was used with Euclidian distance, perplexity 50, 1,000 iterations, and a learning rate η = 100. Default values were chosen for the remaining parameters. Survival analysis was performed using the proportional hazards model as implemented by the coxph function from the R-package survival. Expression values for the different signatures were transformed to quartiles and included as numerical values. Multiple testing correction for the cell populations/signatures was performed using the BH method where appropriate if not otherwise indicated.

Our aim in the present study was to combine gene-expression profiling data from the Affymetrix U133A and U133p2 gene array chips and analyze immune cell and immune pathway expression profiles in pediatric and adult high-grade gliomas according to age as well as mutational and transcriptional subgroups. We identified nine studies with appropriate gene-expression data sets available on GEO (Fig. 1A; refs. 21, 30, 31, 35–40). After exclusion of doublets and other tumor entities not representing gliomas, overall 1,135 cases were included in this study (Supplementary Table S1). To compare gene-expression profiles from different centers as well as platforms, we corrected for batch effects using fRMA preprocessing and additional normalization (Supplementary Fig. S1A; Materials and Methods).

Next, the combined data set was characterized in more detail to allow correlative analyses between different glioma subgroups. However, transcriptional subtype information was available only in 495 samples [not available (NA): n = 640], and the mutational status was available only in 663 samples (NA: n = 472). We therefore used existing gene-expression data sets with known mutational and transcriptional subtype to train a machine-learning algorithm, which was then able to predict the mutational and transcriptional subtype of gene-expression data sets with unknown mutational and transcriptional subtype information. To do so, we applied a random forest approach to separately predict the mutational and transcriptional subgroups of the missing samples (Fig. 1B). Machine learning was able to predict the mutational subgroup, defined by either histone H3 mutations at positions G34 or K27, as well as IDH mutations, or tumors harboring a wild-type for IDH (WT) and no histone H3 mutations (17–19). Next, machine learning predicted the transcriptional subgroup, defined by classic, mesenchymal, neural, and proneural signatures, by extrapolating from the known expression profiles. In an additional in silico validation step, the model accuracy was determined by internal cross-validation (Fig. 1C). Using tail probabilities and applying a threshold of >91% prediction accuracy (Fig. 1D), we were able to predict the mutational subgroup in an additional 344 cases and the transcriptional subgroup in 560 additional cases based on the known gene-expression profiles (Fig. 1E). The age distribution of the individual groups is shown in Supplementary Fig. S1B. We kept the neural phenotype in our prediction model despite the controversial literature in this regard (16). We did this, although data suggest that a neural profile might be contaminated with large amounts of nonneoplastic brain cells. Including the neural group therefore enabled us to investigate the immune microenvironment of high-grade gliomas that are particularly diffuse or infiltrated by normal brain cells.

To analyze the abundance of tumor-infiltrating immune cell populations by inferring from gene-expression data, several gene signatures have been published using different tools (32, 33, 41). These approaches are based on gene-expression signatures of the cell types under consideration and are mostly derived from pancancer analyses. However, some of the known signatures also include genes that are expressed during neuronal development or tumorigenesis of malignant brain tumors. We therefore manually reviewed these signatures and tested the correlation of the published immune cell population signatures in our data set for internal correlation to extract a glioma-specific immune subset signature, which is independent of tissue-specific interferences (Fig. 2; Supplementary Table S2).

Inspection of diverse marker signatures revealed several genes (n = 22/144) that were not immune cell specific. Some of these genes are involved in neural development or glioma biology [e.g., KLRC3 (42), EOMES (43), GREM1 (44); Supplementary Table S2]. Furthermore, many signatures included genes with low or negative correlation in our data set, which is in contradiction with the hypothesis that these genes are specifically expressed by a unique type of cell in high-grade glioma (Fig. 2, gray dots; ref. 33). As genes with different coexpression patterns are unlikely to be specifically expressed by a single-cell population, only those genes belonging to the largest connected component were retained as final markers for the following analyses (Fig. 2, orange dots).

After optimization of the gene-expression signatures of infiltrating immune cell populations, we analyzed the composition of immune and stromal cell gene-expression clusters in correlation to age, transcriptional subgroup, and mutational subgroup (Fig. 3A). Unsupervised hierarchical clustering identified four major clusters, which segregated according to their cellular signature into four clusters: vascular (immune 1), monocytic and stromal (immune 2), monocytic and T-cell dominant (immune 3), and an immune cluster with a gene-expression profile indicative of infiltration by antigen-presenting cells (APC), NK cells, and T cells (immune 4; Fig. 3A). Hierarchical clustering revealed that certain immune clusters correlated with the transcriptional subtype. Immune 1 primarily correlated with the classic transcriptional subtype, whereas immune 2 and 3 were dominated by a mesenchymal signature. Immune 4 demonstrated a mixed picture with no association to a transcriptional subgroup. Furthermore, we observed an increase of immune cluster 4 in G34-mutated tumors, as well as the expected increase of immune cluster 2 in mesenchymal tumors (Supplementary Fig. S2A and S2B). There was no obvious correlation of the immune clusters with age (Fig. 3A).

Next, we analyzed the prevalence of each immune/stromal cell population, as defined by their abovementioned gene signature (Supplementary Table S3), individually according to mutational subgroup, age, transcriptional subgroup, and histologic diagnosis (Fig. 3B). Age did not have a significant influence on the gene expression of classic immune cell subsets (Fig. 3B, second row). Age did, however, affect the stromal component, as tumors of pediatric and young patients demonstrated significantly lower gene-expression profiles of fibroblast-like/stroma (P = 1.9e⁻⁵) and in part endothelial signatures (P = 0.0016). Low expression of the fibroblast-like signature was most frequent in IDH-mutated tumors (Fig. 3B). Association of stromal expression with age was therefore most likely due to the age distribution of IDH-mutated tumors, which were most prevalent in young adults (Supplementary Fig. S1) and not to intrinsically age-related mechanisms, as the stromal expression signature was not observed in WT gliomas. The stromal component was associated with monocytic cells (Supplementary Fig. S2C). This association was primarily seen in the immune 2 cluster, which was increased in tumors with a mesenchymal transcriptional subtype (Supplementary Fig. S2). T-cell signatures were most prevalent in WT glioma and in tumors with mesenchymal signatures and were significantly less expressed in histone H3–mutated (K27 and G34) tumors (P = 0.0022; Fig. 3B, first and third rows). A similar pattern was observed for cytotoxic T-lymphocyte (CTL) genes. B lineage cell gene expression was the strongest in G34-mutated tumors but was not influenced by the transcriptional subtype. Monocytic lineage clusters, which include microglia and tumor-infiltrating macrophages, were underrepresented in G34-mutated tumors (P = 2.5e⁻⁷) but overrepresented in gliomas with a mesenchymal signature (P = 7.9e⁻⁹²). Signatures for myeloid dendritic cells (mDC), neutrophils, and NK cells, on the other hand, were expressed in G34-mutated and IDH-mutated tumors, although no correlations to a transcriptional subgroup were found. The endothelial component was increased in IDH-mutated and WT tumors (P = 0.0018) with a classic signature (P = 2.7e⁻⁵), indicating a higher vascularization in these tumors compared with histone H3–mutated gliomas. Some of the strongest differences were observed with the stroma/fibroblast-like gene signature, which was reduced in IDH-mutated tumors but expressed in WT tumors (P = 5.2e⁻²⁰) and tumors with a mesenchymal signature (P = 1.1e⁻⁶⁶; Fig. 3B, right column). Taken together, these data show an impact of the mutational and transcriptional subgroups on the immune cell composition, demonstrating that the tumor biology, driven by the transcriptional and mutational phenotypes rather than the age of the immune system or the age at tumor occurrence, shapes the immune landscape.

Next, we analyzed the immune microenvironment in more detail by investigating both gene signatures for immune cell populations and functional pathways involved in immune activation or inhibition, which have been defined by glioma or pancancer studies by Thorsson and colleagues and Doucette and colleagues (24, 34). First, we performed a t-SNE analysis using individual gene-expression signatures for immune cell populations that have been optimized for glioma samples (Fig. 2). We included published signatures for pathways involved in immune activation/inhibition (Supplementary Table S3; Fig. 4A). K-means clustering for all samples (n = 1,135) was undertaken to identify microenvironmental subtypes. A Kaplan–Meier plot of the OS within these clustered subtypes is shown in Fig. 4B. Here, the different clusters show significant differences (P = 2.1e⁻¹⁰), which may, in part, be explained by the fact that clusters 1 (green) and 6 (gray) contained mostly IDH-mutated tumors, which have a better prognosis than K27, G34, or WT tumors (Fig. 4A and B). The progression-free survival (PFS) for clusters 1 to 6 (P = 0.0043) and tumor samples that have been collected after previous treatment, i.e., chemotherapy or radiotherapy, are shown in Supplementary Fig. S3A to S3D. No significant impact of previous treatment on the immune gene-expression landscape, except for a decrease of endothelial cell signatures on pretreated samples (P = 2.3e⁻⁵), was observed in this analysis.

In analogy to the data presented in Fig. 4A, t-SNE plots were drawn for the transcriptional and mutational subgroups, as well as for age and histologic diagnosis (Fig. 4C–F). Gene-expression profiles for the t-SNE plotted samples of the individual signatures of immune cell populations and functional pathways from the pancancer (34) and glioma analysis (24) are shown in Fig. 4G–I. The gene expression of the individual signatures is illustrated in Supplementary Fig. S4A and S4B. Immune gene t-SNE clusters correlated best with the transcriptional subgroup. Here, the mesenchymal and classic subtypes demonstrated the strongest expression of immune genes, contributing to polarization of the samples. Cluster 4 was mainly defined by genes indicating an infiltration of regulatory T cells (Treg), fibroblast-like cells, and cells of monocytic and myeloid lineage (Fig. 4G, second row). A strong signature of CTL or the corresponding CD8⁺ T-cell markers was detected only in a minority of these samples, indicating that a stroma-rich environment correlates with strong Treg and monocytic infiltration and precludes infiltration of potentially tumor-specific cytotoxic CD8⁺ T cells. No clear polarization to any cluster was found for other gene signatures, such as mDCs, neutrophils, B cells, or endothelial cells, within the combined analysis of the microenvironment. Although almost all IDH-mutated gliomas are known to express a proneural transcriptional subtype and therefore have a rather good prognosis, the proneural type is associated with poor prognosis in WT glioblastomas (45). When stratifying by their immunostromal gene signatures, these two biologically different types group together, indicating that the transcriptional phenotype has a stronger impact on the immune microenvironment than the mutational IDH status.

Functional immune inhibition or activation, described by pathways of stimulatory and inhibitory receptor/ligand interaction, was primarily associated with clusters 4 and 2. Although the stimulatory receptor/ligand pathways were both expressed in the same tumor samples, the inhibitory receptor/ligand axis was discordantly expressed (Fig. 4H, middle and bottom). Inhibitory ligands, including TGFB1, VEGFA/B, and IL10, among others, were predominantly found in tumors with a mesenchymal transcriptional subgroup and extended beyond the monocytic and regulatory T-cell gene environment. Samples with strong inhibitory ligand gene expression were especially low in CTL signatures. Although the expression of inhibitory receptors, such as PDCD1, CTLA4, LAG3, and KIR2DL1/3, can usually be attributed to CD3⁺ T cells, their expression did not correlate with tumors rich in T-cell gene-expression clusters and was mainly associated with tumors with low immune cell signatures and low stimulatory receptor/ligand pathways (Fig. 4H). The inhibitory receptor expression represented one of the most distinct clusters (Fig. 4H; red circle, cluster 6), distant from immunosuppressor and immunosuppressive signaling signatures defined by Doucette and colleagues (Fig. 4I) or other glioma antigen and ARG1 + IDO signatures (Fig. 4J; ref. 24). Solely the coinhibitory pathways, which include PDCD1LG2 and SLAMF7, and glioma antigen–rich samples colocalized to the inhibitory receptor area. These samples, which express high amounts of glioma antigen and harbor cells with strong expression of PDCD1 and CTLA4, might identify tumors amenable to checkpoint inhibition treatment. Moreover, within this region, most IDH-mutated tumors can be found. Although they do not express large amounts of antigen, the tumors seem to depend on the inhibitory receptor pathway.

Finally, we wanted to determine the impact of individual immune cell population signatures and functional immune pathways on the OS of high-grade gliomas. First, we stratified all tumor samples with available OS information (n = 1,057) by diagnosis, age, transcriptional subgroup, and mutational subgroup and depicted the hazard ratio of individual gene signatures as heat map representations (Fig. 5A; Supplementary Table S4A). Across all samples, functional immune pathways, including immunosuppressors, immunoactivators, immunosuppressive cytokines, and checkpoints or immune effectors, among others, are all associated with worse OS. These effects are pronounced in tumors with glioblastoma and AA histology and in patients >40 years of age. The negative impact of functional immune pathways on the OS is strongest in tumors with proneural transcriptional phenotype, which also include IDH-mutated glioma, whereas the OS of mesenchymal tumors is independent of their immune microenvironment. This also holds true for WT tumors, in which no association with worse survival and negative immune regulators was found. Improved survival is observed primarily in association with gene-expression profiles of cytotoxic (T and NK cells) and APCs (mDCs and B cells) or inhibitory receptors (PDCD1, CTLA4, and LAG3; Fig. 5A). These data show that the immune microenvironment in tumors with mesenchymal transcriptional subgroup or patients at young age (<19 years) does not appear to affect OS, potentially indicating that these tumors are characterized by an immunosuppressed environment.

To identify the influence of immune cell signatures on OS (Fig. 5B) as well as PFS (Supplementary Fig. S3D) and to determine a marker profile that could be used in clinical diagnostics to classify the different immune microenvironments, we defined four immune clusters (Fig. 3A). When we plotted the OS of each immune cluster, we were able to confirm that immune cluster 2 (monocytic/stromal) was associated with a significantly worse prognosis (Fig. 5B). Other immune clusters with an immune microenvironment rich in APC and cytolytic cells (NK and T cells) present the best OS.

Because of continuing study of immune-checkpoint inhibitors in clinical trials, we analyzed the association of major immune checkpoints with the different high-grade glioma groups (Fig. 5C; Supplementary Table S4B). Here, CD274 (PD-L1), CTLA4, PDCD1, and LAG3 demonstrated significant correlations. As expected, PD-L1, a marker that is associated with immune escape, was associated with a worsened survival (hazard ratio >1), whereas the other markers (LAG3 and CTLA4) correlated with better survival. PD-L1 was informative in tumors with AA histology and with K27 mutation. Some checkpoints appeared to have a contrary impact on OS, as TGFB1 expression resulted in better prognosis in tumors with classic transcriptional phenotype and a worse survival in patients under the age of 20 years. Taken together, these data demonstrate that immunologic phenotype refines response of high-grade gliomas to immunologic treatment approaches.

Molecular stratification of high-grade gliomas is paving the way to a better understanding of tumor heterogeneity. Integrated molecular analyses allow stratification into different subtypes and prediction of clinical response to therapy (13). Analyses in adult glioblastoma patients suggest that genomically diverse subtypes display distinct microenvironmental profiles (20, 22, 24). It is therefore hypothesized that genetic driver mutations promote unique immunologic landscapes (20, 23). Given the diverse genomic alterations between pediatric and adult high-grade glioma, i.e., K27 and G34 mutations, as well as the potential clinical implications for the translational use of therapeutic regimens, including chemotherapy and immune-checkpoint inhibition, we compared the immunologic microenvironment of high-grade gliomas across age groups and stratified by clinically relevant genetic subgroups.

In our study, we demonstrate that different mutational and transcriptional subgroups shape the composition of the immune subset signatures, independent of age or major histologic diagnoses. We identified four distinct immunologic patterns of immune cell profiles, which are defined by either vascular (immune 1), monocytic/stromal (immune 2), monocytic/T-cell (immune 3), and APC/NK cell/T-cell (immune 4)–dominant immune clusters. These clusters correlated with the different transcriptional subgroups of high-grade gliomas. Based on a TCGA gene-expression data set of 544 patients, previous work by Doucette and colleagues highlighted that the mesenchymal subgroup displays a gene-expression profile reflective of a proinflammatory antitumor response (24). Our cohort of 1,135 gene-expression profiles further enabled stratification of the mesenchymal subgroup into two immunologically distinct infiltration clusters, as the prevalence of monocytic markers together with a stromal and endothelial composition (immune 2) could be distinguished from a monocytic lineage profile with reduced stromal component but increased cytotoxic T-cell signature (immune 3), respectively. This observation suggests that the mesenchymal subgroup may harbor a distinct immunologic subtype with increased gene expression indicative for higher CTL infiltration, which might be amenable to immune-checkpoint inhibition.

Work published by Luoto and colleagues analyzed the immune cell composition of 156 glioblastoma RNA sequencing (RNA-seq) samples by TCGA to investigate the impact of genomic alterations on the immune microenvironment in IDH WT and IDH-mutated glioblastoma (22). That study identifies three immune response groups, a “negative” macrophage, granulocyte, and macrophage-dominated group in IDH1-mutated tumors, a “humoral,” lymphocyte-enriched signature in EGFR-gained samples, and a “cellular-like,” macrophage-infiltrated subgroup in EGFR-amplified tumors (22). We confirm the observation of a lower monocyte/macrophage signature component in IDH-mutated tumors. However, we did not see an IDH-associated immune phenotype, as in our cohort, K27- and G34-mutated gliomas correlated more with an immune cluster 4 phenotype (immune 4: APC/NK cell/T cell). We showed that IDH WT and mutated tumors displayed the most distinct differences in the stromal/fibroblast-like cell signature. Yet, only IDH-mutated tumors seemed to leverage certain functional pathways, as they expressed “inhibitory receptors” such as PD-1, CTLA4, LAG3, or KIRs, which concurs with the finding of a “negative regulation of lymphocyte function” in IDH tumors by Luoto and colleagues (22).

In a global pancancer analysis, glioblastoma was classified as a lymphocyte-depleted, macrophage-enriched tumor entity with a suppressed T_H1 axis (34). Although the monocytic macrophage lineage signature was also evident in our analysis, we describe a subset of glioblastoma with low expression of monocytic lineage gene cluster, as the vascular (immune 1: vascular) and APC/NK cell/T-cell (immune 4) clusters were depleted of monocyte/macrophage signatures. The K27 and G34 tumors do not seem to follow the glioblastoma pattern described by Thorsson and colleagues (34). We show that K27 and G34 tumors displayed lower gene-expression signatures of monocytic lineage cells and harbored lower cytotoxic and overall T-cell gene expression compared with IDH-mutated and WT gliomas. The depletion of lymphocytic clusters in the K27 and G34 tumors is comparable with the “immunologically quiet” immune subtype described for low-grade gliomas in the pancancer analysis (34). However, the increased B cell gene-expression prevalence in G34-mutated tumors in our study points to a primary signature indicative for CD4⁺ T_H2 cell infiltration, which therefore deviates from the observed dominance of a T_H1 response in low-grade glioma by Thorsson and colleagues (34).

The abovementioned studies demonstrate that the molecular and transcriptional heterogeneity in malignant gliomas correlates with high variability of the immunologic landscape (20, 22, 24, 34). These studies emphasize the need for further stratification to individualize upcoming clinical studies. Our data demonstrate that age and in part histologic diagnosis are not relevant for the immune subset composition. However, immune cell signatures and functional immune pathways are associated with the transcriptional subgroups (13, 16). Especially the mesenchymal subgroup demonstrated immune activity and association with monocytic lineage profiles and immune clusters 2 and 3. The strongest correlation of immune-relevant cell signatures and functional pathways was, however, observed in the proneural subgroup. Here, immunosuppressive and tumor-associated macrophage–related pathways, as well as stromal/fibroblast-like components correlated with a decreased OS. Although the mesenchymal subgroup displayed an immunologic component, almost no correlation with OS was present. In addition, immune gene expression did not correlate with either survival in IDH WT gliomas or with patient age below 20 years, implying an immunologic “inertness” in these groups.

We initially hypothesized that high-grade gliomas in children and young adults might demonstrate distinct immunologic profiles due to a still developing or young immune system. Except for the above-described lack of correlation of immunologic gene expressions with OS, we did not find any association with a specific immune cell pattern. However, when dissecting high-grade gliomas of the young and adolescent patients into their molecular subgroups, i.e., K27- and G34-mutated tumors, opposing effects of immunologic genes on OS became apparent. In the K27-mutated tumors, CD274 (PD-L1) and CTLA4 correlated with a worse prognosis, whereas the opposite was true for G34-mutated tumors. Furthermore, G34 tumors appear to rely primarily on TGFB1 and HAVCR2 as pathways for immune escape. This finding implicates that checkpoint inhibition of the PD-1/PD-L1 pathway might be beneficial for patients with diffuse midline gliomas harboring histone H3 K27M mutations.

Our study shows that immunologic heterogeneity in high-grade gliomas is associated with distinct transcriptional and mutational subgroups, independent of patient age. The immune cell signatures correlate with the transcriptional subgroup and can be clustered into four immunologic groups. The mesenchymal subgroup shows an immunologic activation pattern, which does not correlate with OS, reflecting the inert immunobiology of this transcriptional subtype as well as IDH WT glioblastoma. The proneural subgroup, on the other hand, which includes IDH-mutated tumors, displays expression on “inhibitory immune receptors” that correlates with OS. K27 and G34 gliomas exhibit diverse immune pathway expression, with K27 tumors especially deprived of lymphocyte, NK cell, and antigen-presenting mDC signatures and presumably other mechanisms of immune escape, such as the PD-L1 pathway. Taken together, our study demonstrates that future immunotherapeutic strategies should take into account the distinct immunologic landscapes that may modulate therapeutic efficacy, especially for rare mutational subgroups of high-grade gliomas.

M. Westphal is a consultant/advisory board member for AbbVie. No potential conflicts of interest were disclosed by the other authors.

Conception and design: M. Bockmayr, C.L. Maire, U. Schüller, M. Mohme

Development of methodology: M. Bockmayr, F. Klauschen, M. Mohme

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): M. Bockmayr, C.L. Maire, S. Rutkowski, K. Lamszus, M. Mohme

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): M. Bockmayr, M. Westphal, K. Lamszus, M. Mohme

Writing, review, and/or revision of the manuscript: M. Bockmayr, F. Klauschen, C.L. Maire, S. Rutkowski, K. Lamszus, U. Schüller, M. Mohme

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): M. Bockmayr, M. Mohme

Study supervision: F. Klauschen, K. Lamszus, M. Mohme

This study was supported by the Werner-Otto-Stiftung (to U. Schüller) and the Fördergemeinschaft Kinderkrebszentrum Hamburg (to M. Bockmayr and U. Schüller), the Anni-Hofmann Stiftung (to K. Lamszus), and the Else Kröner-Fresenius Stiftung (to M. Mohme).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.