Immune microenvironments differ in immune characteristics and outcome of glioblastoma multiforme

Abstract Understanding the interactions between tumors and the host immune system holds great promise to uncover biomarkers for targeted therapies, predict the prognosis of patients and guide clinical treatment. However, the immune signatures of glioblastoma multiforme (GBM) remain largely unstudied in terms of systematic analyses. We aimed to classify GBM samples according to immune‐related genes and complement the existing immunotherapy system knowledge. In this study, we designed a strategy to identify 3 immune subtypes representing 3 different immune microenvironments (M1‐M3) and associated with prognosis. The 3 subtypes were significantly different in terms of specific immune characteristics (immune cell subpopulations, immune responses, immune cells, and checkpoint gene interactions). In additional, copy number variations and methylation changes were identified that correlated with genes related to a worse prognosis subtype in the microenvironment. More importantly, in M3 (worst prognosis subtype) and M2 (best prognosis subtype), the interaction between immune cells and checkpoint genes was different, which had an important effect on the prognosis. Finally, we used risk scores of immune cells and checkpoint genes to evaluate the prognosis of GBM patients and validated the results with 3 independent datasets. Disordered interactions between immune cells and checkpoint genes result in a change in the immune microenvironment and affects the prognosis of patients. We propose that a better understanding of the immune microenvironment of advanced cancers may provide new insights into immunotherapy.

tumors, and targeted therapy will become an important means of tumor treatment in the future. 6 However, immunotherapy is only applicable to a substantial fraction of patients, while others either are not suitable candidates or fail to respond. This has become a difficult problem, so new standards are urgently needed to guide individualized treatment.
An imbalance between the host microenvironment and tumor cells can result in tumor invasion, progression, and metastasis. Immune responses play critical roles in carcinogenesis and the progression of solid tumors. 7,8 It has been reported that nascent transformed cells can be initially eliminated by the host immune system based on both innate and adaptive immunity. 9 The immune environment of a primary tumor is associated with the clinical response and benefit of immunotherapy. 10,11 For example, the immune checkpoint molecules programmed cell death protein 1 (PD1) and PD1 ligand 1 (PDL1) are targets of drugs in clinical practice. 12 Some studies have demonstrated that tumor-infiltrating CD8 + T lymphocyte cells and intratumoral TH1-type molecules are associated with positive responses to therapeutics by blockading PD1 and PDL1. However, therapies with an antibody targeting PD1 (anti-PD1) displayed response rates from 17% to 21%, with some responses being remarkably durable. 13,14 The reason for this phenomenon may be tumor mutational burden 15 or cytokine release syndrome. 16 However immune responses and diseases are not caused by a single immune cell or checkpoint, their interaction relationship changes. 17 Therefore, studies of the relationship between the tumor and host microenvironment have revealed that interactions between the immune microenvironment and molecules affect the therapeutic outcome of patients. Researchers have also explored promising candidate biomarkers for predictive and prognostic value to further guide the individualized treatment of patients.
This study aimed to investigate the overall immune landscape and its clinical relevance in GBMs. From The Cancer Genome Atlas (TCGA) GBM cohort, we identified 3 immune subtypes (M1-M3) of GBMs based on the expression of global immune genes and the subtypes were related to patients' survival. Each subtype was characterized by a different subpopulation of immune cells and immune responses. The abnormally expressed genes in M3 (the worst prognosis subtype) also exhibited different methylation and copy number variation (CNV) levels from the M1 and M2 subtypes. In M2 (the best prognosis subtype) and M3, the interactions between immune cells and checkpoint genes differed significantly, and we identified 3 checkpoints (CD27, PDL1, and CTLA4) with the most significant differences between the 2 subtypes. Then, we found that the total risk scores for these 3 checkpoints and 24 immune cells were related to patients' survival, furthermore, high risk scores predicted worse prognosis. In addition, we validated these findings in 3 independent datasets. We observed that the immune microenvironment of different patients was different, which would affect the immune treatment outcome of patients. Through this method, we hope to find target markers in different immune microenvironments that interact with immune cells to predict patients' prognosis in order to guide individualized treatment.

| Patients and samples
The GBM samples used in this study were obtained from Beijing Tiantan Hospital between January 2005 and December 2009 (detailed information about the tissue samples is presented in Table S1). Overall survival (OS) was calculated from the date of diagnosis until death or the end of follow-up, while progression-free survival was defined as the time between diagnosis and the first unequivocal clinical or radiological sign of disease progression. All samples were rinsed with normal saline after surgical resection and divided into 2 portions: one was placed in liquid nitrogen immediately and then stored at −80°C until use, and the other was placed in formaldehyde solution and hematoxylin-eosin stained to assess the percentage of tumor cells. All patients signed informed consent forms. The use of human tissue samples and the experimental procedures for this study were reviewed and approved by the Ethics Committee of the Cancer Institute and Hospital, Chinese Academy of Medical Sciences.

| RNA isolation and gene expression microarray
Total RNA from frozen samples was extracted using the mirVana miRNA Isolation kit (Ambion 1561) according to the manufacturer's protocol. RNA concentrations were determined by an ND-1000 UV-VIS Spectrophotometer (NanoDrop Technologies, Wilmington, DE), and RNA integrity was evaluated using the RNA 6000 LabChip kit in combination with the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA). The RNA samples used in this study all exhibited RNA concentrations >40 ng/μL and RNA integrity numbers >7.0. After RNA concentration and integrity analysis, the samples were analyzed using Agilent 4 × 44K Whole Human Genome Oligo Microarrays at the Cancer Hospital, Chinese Academy of Medical Sciences, according to the manufacturer's specifications. In brief, 500 ng of purified total RNA was reversed transcribed in vitro using the Low RNA Input Linear Amplification Kit PLUS (Agilent) and then transcribed into cRNA labeled with Cy3. In total, 1.65 μg of cRNA was hybridized to each microarray. After hybridization, the slides were washed and scanned with an Agilent G2505B Microarray Scanner System. The fluorescence intensities of the scanned images were extracted and preprocessed using Agilent Feature Extraction Software (v9.1). The raw data were normalized using the GeneSpring GX software program, version 11.5 (Silicon Genetics, Redwood City, CA). The raw and processed data are publicly available at the Gene Expression Omnibus (GEO) website under the accession number GSE122586.

| Public clinical and molecular data collection
In TCGA datasets, RNAseq data (level 3, RSEM-normalized data), methylation array data (Illumina Human Methylation 450) and CNV data (Affymetrix SNP 6.0) of GBMs were downloaded from the NIH National Cancer Institute GDC Data Portal (https ://portal.gdc.cancer.gov/). Two external independent transcriptome datasets of GBMs were used for validation, namely, the Chinese Glioma Genome Atlas RNA sequencing dataset (CGGA data) and GEO microarray datasets GSE16011, as well as their corresponding clinical information.
A total of 404 immuno-related human genes, including 4 immune response types, 24 immune cells, and 22 immune response categories, that were curated from the nCounter® PanCancer Immune Profiling Panel (NanoString) were implemented as candidate genes in this study. Detailed annotations for these 404 genes are listed in Table S2.

GBMS by immune genes
The expression of immune-related genes was used to identify immune subtypes. Through T-Distributed Stochastic Neighbor Embedding (t-SNE), the TCGA GBM cohort was divided into 3 immune subtypes, representing 3 different immune microenvironments. The optimal cluster number was determined by NbClust, an R package that determines the best number of clusters in a dataset.

| Immune cell subpopulations and responses in GBM subtypes
We used gene set enrichment analysis to identify immune cell types that were enriched in each immune subtype. The expression level of each gene was log2-transformed for subsequent analysis. For each patient, genes were ranked in descending order according to their expression values, and the association was represented by a normalized enrichment score. An immune cell type was considered enriched in a patient when the P-value was <0.1. The percentage of each immune cell with significant enrichment in each immune subtype was calculated to compare the proportion of immune cell types in different classifications. 4 The enrichment scores of immune responses in each immune subtype were determined by single sample enrichment analysis in the R package GSVA. Gene set variation analysis (GSVA) is used to estimate the variation in gene set enrichment through samples of expression datasets, so the enrichment scores of immune responses in each subtype to be computed can be compared. 18

| Regulation of abnormally expressed genes in the M3 subtype
For each subtype, the featured genes were identified by comparing the samples in this subtype with the remaining samples using Student's t test. Through comparing the differentially expressed genes of M3 subtype with M1 subtype and M2 subtype, respectively, we selected significantly differentially expressed genes in the M3 subtype, and then to explore the mechanisms of changes differentially expressed genes of M3 subtype via epigenetic variations and CNVs.
As for mRNA expression, median expression levels (used to summarize expression in each subtype) were computed using only samples with nonmissing values.
To study the relationship between the expression and DNA methylation of those genes, we mapped DNA methylation probes to the genes. The methylation level of a particular gene was defined as the mean value of all probes mapping to that gene. For a given gene, the beta value was evaluated within each immune subtype. In addition, we used Student's t test to examine whether these genes were differentially methylated in M3 subtype compared with M1 subtype and M2 subtype.
We used the GISTIC2 method to estimate the thresholded gene-level CNV of GBMs to examine CNV changes in M3 abnormal genes across subtypes. 19

| Reconstruction of the immune response interaction network
Based on the gene expression data, we reconstructed immune cell-gene (immune checkpoint inhibitors and co-inhibitory and co-stimulatory markers of the innate and adaptive immune systems) networks; these immune checkpoint genes were selected from a previous report. 20 Through GSVA to estimate the enrichment scores of each immune cells, this score can represent the expression or infiltration of immune cells. Then we used Pearson correlation to computer the correlation of the enrichment score of immune cells and the expression level of the target checkpoint genes. The network was visualized using Cytoscape. 21 The edge weights of the network were based on the Pearson correlation coefficient between the immune cells and the target checkpoint genes. The cyan elliptical nodes represent immune cells, and the salmon octagonal nodes represent immune checkpoint genes.

| Prognostic analysis
To identify the infiltration of immune cells and checkpoint genes that could predict GBM patient prognosis, a risk factor score was calculated to assess the survival of patients. In brief, we used univariate Cox regression analysis to evaluate the association between survival time and the infiltration of each immune cell type and the expression levels of checkpoint genes. Regression coefficients with a plus sign indicated that increased expression was associated with decreased survival, that is risky factors; conversely, a minus sign indicated that increased expression was associated with increased survival, that is protective factors. A mathematical formula was constructed to predict survival, and we assigned a risk score to each patient by the regression coefficients from the univariate Cox regression analysis. 22 The risk score of each patient was calculated as follows: i and β represent the regression coefficients of gene expression values and the infiltration degree of immune cells, respectively. All patients in the dataset were thus assigned to high-risk and low-risk groups using the median risk score as the cut-off point. Patients with higher risk scores were expected to have poor survival outcomes. The Kaplan-Meier method was used to estimate the OS time for the two subgroups, and differences in survival time were analyzed using the log rank test (R package "survival").

| Statistical analysis
All statistical analyses in this study were performed using R software (http://www.r-proje ct.org). T-SNE was implemented in the R package Rtsne. Heatmaps and Circos plots were generated by the R packages pheatmap and OmicCircos, respectively. All statistical tests were 2-sided, and a P value less than 0.05 was considered statistically significant.

| Immune subtypes in cancer
To characterize the immune characteristics of all TCGA GBM samples, immune signatures according to the nCoun-ter® PanCancer Immune Profiling Panel classification into 4 immune response types (adaptive response, inflammation response, humoral response and innate response) containing a panel of 404 immune-related genes were used to describe the immune landscape in GBM samples ( Figure 1A). We observed that different immune responses in different samples were associated with differences in the immune microenvironment, additionally, the interactions between tumors and hosts were diverse. To identify common immune subtypes and evaluate whether tumor microenvironment features can predict outcomes, we analyzed the microenvironments across the landscape of all GBM samples. Four immune response types, that is 404 genes were used to perform a cluster analysis, and 3 immune subtypes were obtained, representing 3 different immune microenvironments (referred to as M1, M2, and M3, Figure 1B). Furthermore, the immune microenvironment subtypes were associated with OS. M2 had the best prognosis, while M3 had the least favorable outcome ( Figure  1C). These findings suggested that the immune microenvironment affects the prognosis of patients. Therefore, we adopted the perspective of immune microenvironments to explore immunotherapy effects, which provided us with fresh clues.

| Composition of the 3 immune subtypes
The immune cell proportions of the tumor stromal fraction varied across immune subtypes. Using a transcriptome expression dataset, we performed GSVA of 24 immune cell subpopulations, including 11 innate and 13 adaptive immune cell subpopulations. The results showed that the proportion of subpopulations of immune cells was different in each immune subtype. In M1, the immune cell subpopulations with the top 3 highest enrichments included Tcm (central memory T) cells (79%), TFH (T follicular helper) cells (42%), and B cells (39%). In M2, the top 3 highest enrichments included Tcm cells (84%), macrophages (80%), and TFH cells (31%); however, in M3, the top 3 highest enrichments were Tcm cells (79%), macrophages (63%), and B cells (51%). These results, which are shown in Figure 2A-C, indicated that the proportion of immune cells was associated with the immune microenvironment, and different immune microenvironments affected the prognosis of patients.
In addition to immune cell subpopulation differences, through GSVA, we observed that the 3 subtypes demonstrated different immune responses. The M3 subtype (with the worst prognosis) showed high enrichment scores in immunosuppression, phagocytosis, leukocyte migration, TNF superfamily members and receptors, while the M2 subtype (with the best prognosis) exhibited high enrichment scores in CD8-positive response, B-cell activation, B-cell receptor signaling pathway, transcription factors, and so on ( Figure 2D).

| Regulation of abnormally expressed genes in the M3 subtype
The abnormally expressed genes in M3 were not conducive to the prognosis of patients; therefore, we speculated that these genes might be critical for cancer occurrence and development. To further explore this theory, understanding of their expression in different states of the immune microenvironment was needed. We examined the expression and expression control of those genes via epigenetic and CNV mechanisms.
The expression profiles of 179 abnormally expressed genes in M3 immune subtype varied across different immune subtypes, perhaps indicating their role in shaping the immune microenvironment. The methylation levels of these genes in M3, there were 70 (39.1%) genes had lower methylation levels than in the M1 and M2 subtypes, and, 57 (31.8%) genes had higher methylation levels than in the M1 and M2 subtypes, indicating that the methylation levels affected the expression of multiple genes and varied across immune subtype ( Figure 3A). However, the CNVs of most genes remained unchanged; yet, for the CNVs in M3 compared with the M1 and M2 immune subtypes, the amplification frequency (82 genes, 45.8%) was the highest, while the deletion frequency (51 genes, 28.5%) was higher than that of M2, but lower than that in the M1 subtype ( Figure 3B). These findings suggest that changes in copy number may affect gene expression levels to some extent. The changes in the transcriptome expression, methylation and CNV levels of these genes are shown F I G U R E 1 Differentially expressed immune signatures in GBM. A, The expression of an immune signature containing 404 immunerelated genes in GBM samples is shown as heatmap. Generally, these genes were classified into 4 groups: adaptive, inflammation, humoral, and innate immune response-related genes. B, GBM samples were divided into 3 groups. The blue shading represents immune microenvironment 1 (M1, 38 samples), the magenta shading represents immune microenvironment 2 (M2, 45 samples), and the red shading represents immune microenvironment 3 (M3, 68 samples). C, Kaplan-Meier survival plots of each subtype. The blue line represents the M1 subtype, the magenta line represents the M2 subtype, and the red line represents the M3 subtype on chromosome. In the circos plot, an ideogram of a normal karyotype is shown in the outermost ring. The next outermost ring is the heatmap of gene expression at corresponding genomic coordinates, red represents high gene expression and blue represents low gene expression; the next ring represents the CNVs, multiple lines in stair steps show different individuals, red lines represent copy number amplification and blue lines represent copy number deletion; the innermost ring is heatmap illustrating DNA methylation β values, red represents high methylation β values and blue represents low methylation β values; and the middlemost lines showing the interaction between genes and genes ( Figure 3C). Overall, these marked differences in M3 abnormally expressed genes may be reflective of modulation of the immune microenvironment by immune cells and cancer cells and may affect the therapeutic outcome of patients. The observed differences in the regulation of those genes might have implications for therapeutic development and combination immune therapies, and the multiple mechanisms at play in evoking them further highlight their biological importance.

| Networks modulating immune response interactions
A number of immune checkpoint inhibitors and co-inhibitory and co-stimulatory markers of the innate and adaptive immune systems are currently under investigation for immunotherapy in various cancers. In an attempt to identify promising candidates for GBM immunotherapy, we reconstructed a network of immune cells and immunomodulatory molecules to explore which molecules are available. The interaction network between immune cells and these checkpoint genes was then filtered for molecules that were significantly associated with prognosis. We found that in the M2 (good prognosis) and M3 (bad prognosis) subtypes, the interactions between immune cells and checkpoint molecules were significantly different, and the interaction coefficient was higher in M2 subtype than M3 subtype, that is the stronger interaction between immune cells and checkpoint genes, the better the prognosis of patients. CD27, PDL1, and CTLA4 were at the center of the network, at the same time, the 3 checkpoint genes were the most significant changes in their interactions between M2 subtype and M3 subtype, we hypothesize that these three genes in different immune microenvironments significantly affect the prognosis of patients ( Figure 3D-E). It has been suggested that the prognosis of patients or the occurrence of disease is determined not by a single molecule but by the interaction of genes or molecules. Because different patients have significant differences in medicine sensitivity, their prognosis also differs significantly. This requires that individual treatment be adopted for different individuals in order to achieve the best therapeutic effect. Studying the interactions between target molecules and other molecules or cells and then connecting them with the prognosis of patients is F I G U R E 3 Regulation of abnormally expressed genes in the 3 subtypes and the interaction network between 24 immune cells and 18 checkpoint genes in the M2 and M3 subtypes. A, From left to right: mRNA expression (median normalized expression levels); methylation expression (median DNA methylation beta value). B, The ratio of amplification and deletion frequencies in 3 immune subtypes, respectively. C, Circos plot displaying the distribution of gene expression, DNA methylation, CNV, and interactions between genes on chromosomes. D-E, The interaction networks between 24 immune cells and 18 checkpoint genes in the M2 and M3 subtypes, respectively. The cyan elliptical nodes represent immune cells, and the salmon octagonal nodes represent checkpoint genes. The edges represent correlation coefficients, and the colors of the edges represent the Pearson correlation coefficients, where red represents a high correlation coefficient and blue represents a low correlation coefficient of great importance for determining patients' prognosis and guidelines for clinical therapy.

| Correlation between interaction relationship and overall survival
To explore the influence of network relationships on survival, we tested the prognostic value of CD27, PDL1, and CTLA4 and 24 immune cells. We performed a univariate Cox regression analysis to evaluate the association between survival time and the 3 checkpoint genes expression levels and 24 immune cells enrichment scores. In our datasets (GSE122586) and 2 external independent datasets, survival analysis indicated that the relationship could predict the prognosis of GBM patients, and high risk scores were associated with poor survival (Figure 4A-D). These results suggested that the interaction relationship was able to predict the prognosis of GBM patients, assisting with postsurgical treatment management and providing a priori guidance for individualized treatment and targeted treatment.

| DISCUSSION
We have developed an integrated strategy to characterize the interaction relationship between the immune microenvironment and the clinical outcome of human GBM. Our approach to deeply mine large datasets enabled us not only to disentangle tumor-immune interactions but also to devise strategies for guiding cancer immunotherapy in GBM.
There is a growing body of evidence suggesting that the interaction between cancer cells and the host immune microenvironment plays a critical role in the occurrence and development of tumors. 23 Research into immune checkpoints marks the beginning of a new era in cancer immunotherapy. However, immunotherapy works in only a small number of cancers, and only a fraction of patients have good treatment results. 24 Further research into how to guide the clinical treatment and predict the prognosis of patients offers considerable potential.
GBM is the most malignant tumor of the nervous system, but the function of the immune response in GBM progression and prognostication remains completely unknown. High-throughput sequencing technology provides objective data for this purpose. In this study, we divided GBM samples into 3 subtypes utilizing immune signatures. The 3 subtypes differed significantly in terms of specific immune characteristics (immune cell subpopulations, immune response, interactions between immune cells, and checkpoint genes). In additional, CNV and methylation changes were identified that correlated with genes related to a worse prognosis in the microenvironment. These findings may provide new insight into strategies for immunotherapy of GBM.
Previous studies have reported some molecular subtypes of glioma based on genome-wide profiles. 25,26 In our analysis, we focused only on global immune-related gene profiles, which could provide more detailed information about the immune landscape of GBM. The M2 subtype had the most favorable prognosis, suggesting that the CD8-positive response and an increase in B-cell infiltrates are needed for cancer control, consistent with previous reports. [27][28][29] In contrast, the M3 subtype conferred the worst prognosis and displayed composite signatures reflecting immunosuppression response, leukocyte migration, and so on.
Possible impacts of methylation changes and CNVs in abnormally expressed genes in M3 were seen. Most of these genes were methylation changes, indicating that the methylation levels of these genes affected the expression profiles of multiple genes. Regarding the CNVs in M3 compared with the M1 and M2 immune subtypes, the amplification frequency was the highest, and the deletion frequency was higher than that in M2, but lower than that in the M1 subtype. These findings suggest that changes in CNV may affect gene expression levels to some extent. Further work is needed to determine the functional aspects of these associations.
The immune response is determined by the collective states of intracellular molecular networks in tumor, immune, and other stromal cells via soluble proteins such as cytokines, which mediate interactions among those cells. 30 Therefore, studying the interactions between immune cells and molecules is important for the immunotherapy of various cancers. We found that the interactions between immune cells and immune checkpoint genes differed between the M2 and M3 subtypes. In M2 compared with the M3 subtype, the interaction coefficients between immune checkpoint genes and immune cells were greater, the disordered interactions between immune cells and checkpoint genes result in changes in the immune microenvironment and affect the prognosis of patients. Based on immune cell-gene interactions, we identified the top 3 genes significantly different interactions with immune cells in the M2 and M3 subtypes, and then we used risk scores of 24 immune cells and 3 checkpoint genes (CD27, PDL1, and CTLA4) in 4 independent datasets to predict the prognosis of patients and to guide their clinical treatment. Of note, using these methods, it is not always possible to fully ascertain whether a particular interaction works in the tumor, immune, or stromal cell compartment, but this could be improved by incorporating additional cell type-specific knowledge.
There are some limitations in our study. First, the histological samples were too small and not sufficiently representative to evaluate immune cell infiltration. Second, for most tumor types in TCGA, samples with fewer than 60% tumor cell nuclei according to a pathologist review were excluded from study, 31 thus potentially excluding the most immune-infiltrated tumors from the analysis. The degree to which this biased the results, relative to the general population of cancer patients, is difficult to ascertain. In addition, our analyses were limited by restricting them to data from molecular assays in the absence of targeted classical cellular immunology assays for confirming cell phenotype distribution, as those types of data have not been collected from TCGA patients.
In conclusion, using the gene expression profile of global immune genes, we identified 3 immune subtypes in GBM samples. These subtypes were associated with prognostic, genetic, and immune-modulatory alterations that may shape the specific types of immune environments we observed. With our increasing understanding that the tumor immune environment plays an important role in prognosis as well as response to therapy, defining the immune subtype of a tumor may play a critical role in predicting the disease outcome as opposed to relying solely on features specific to individual cancer types. These findings regarding the intratumoral immune microenvironment may shed new light on immunotherapy strategies for advanced gliomas.
The cancer microenvironment has high prognostic value. The analysis of the immune context of tumors revealed a set of cellular and molecular immune markers that could be used to effectively and reproducibly classify patients according to their survival. Further studies are needed to validate our conclusions and to elucidate the underlying mechanisms of these phenomena.