Immune Cytolytic Activity Is Associated With Genetic and Clinical Properties of Glioma

This study retrospectively enrolled RNA sequencing data of 325 glioma patients, ranging from WHO grade II–IV, and corresponding clinical information from the Chinese Glioma Genome Atlas (CGGA) database (www.cgga.org.cn↗) as training cohort. The RNA sequencing data from the Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/↗) of 699 glioma patients, ranging from grade II–IV, and corresponding clinical information were collected as validation cohort. Moreover, we collected another 20 patients for independent clinical validation. The OS was calculated from the data of initial diagnosis until death or the last follow-up. Thus, in total 1,024 glioma patients were evaluated in this study. Information on these patients are available from corresponding data portal.

IDH mutations typically leading to the replacement of an arginine with a histidine at amino acid residue 132 of the protein (R132H). The mutation results in a loss of normal enzymatic function and the abnormal production of 2-hydroxyglutarate (2-HG) (13). In CGGA cohort, IDH1/2 mutations were detected by DNA pyro-sequencing as previous reported (14). And the IDH1/2 mutations information were downloaded from TCGA website in TCGA cohort.

We obtained Fragments Per Kilobase Million (FPKM) values of glioma samples from two datasets. The FPKM values was transformed into Transcript Per Million (TPM) values with following formula:

The CYT was calculated as the geometric mean of the GZMA and PRF1 expression in TPM (10). Significance analysis of microarrays (SAM) was applied to identify the differentially expressed genes based on the threshold of fold change more than 2 and false discovery rate (FDR) value < 0.01 adjusted by 1,000 times permutation test. Gene ontology (GO) analysis was performed in DAVID (http://david.abcc.ncifcrf.gov/home.jsp↗) for functional annotation of the genes.

MR images of patients from the CGGA and TCGA databases were obtained from the CGGA imaging database (http://www.cgga.org.cn↗) and the Cancer Imaging Archive (http://www.cancerimagingarchive.net↗), respectively. Tumor contrast enhancement was assessed by two experienced neuroradiologists blinded to the patients' clinical information. A third senior neuroradiologist re-examined the images and determined which should be used if the types of enhancement identified by the first two neuroradiologists were inconsistent. Briefly, a small (or no) region of edema (-) was defined as edema extending ≤1 cm from the margin of the tumor based on T2-weighted images; otherwise edema was graded as moderate to severe (+) (15).

Difference in variables between groups were tested by the Student's t-test or chi-squared test, as appropriated. The Kaplan-Meier curves were employed to estimate overall survival distribution. Log-rank test was applied to compare the statistic difference of survival curves between two groups. All figures and statistical analysis were conducted using R software (version 3.4.2; http://www.R-project.org↗) and GraphPad Prism 7.0 software. A p < 0.05 was considered as statistically significant. All statistical tests were two-sided.

To assess the stability and quality of CYT, we profiled the distribution of GZMA and PRF1 using their RNA expression data in CGGA and TCGA datasets, respectively. As shown in Figure S1, GZMA and PRF1 exhibited significantly positively correlation relationship in CGGA (r = 0.57, p < 0.0001) and TCGA (r = 0.7, p < 0.0001) datasets. Due to the intratumor heterogeneity and plasticity of glioma, we investigated the relationship of CYT and previous widely accepted prognostic and predictive factors, including grade, IDH status, and 1p/19q status. In CGGA dataset, cytolytic activity was significantly increased with increasing grade and the CYT exhibited the highest value in glioblastoma (Figure 1A). The similar result was well validated in TCGA dataset (Figure 1B). Moreover, IDH wildtype glioblastoma showed the highest CYT values (p < 0.001) when comparing to other four subtypes of glioma in both CGGA (Figure 1C) and TCGA (Figure 1D) datasets.

In addition, we also explored that CYT was highly enriched in IDH wild-type glioma group in both CGGA (Figure S2A, p < 0.0001) and TCGA datasets (Figure S2B, p < 0.0001). Meanwhile, we found CYT value was significantly upregulated in 1p/19q intact glioma group than 1p/19q co-deleted group in both CGGA (Figure S2C, p < 0.0001) and TCGA datasets (Figure S2D, p < 0.0001). These findings indicated that cytolytic T-cell-related immune activities were more prevalent in these relatively malignant subtypes of glioma.

Glioma can be classified into four molecular subtypes: classical, proneural, neural, and mesenchymal (16). The glioma patients in mesenchymal subtype group always confer aggressiveness, treatment resistance and worse prognosis (17). The core genes to identify the mesenchymal subtype glioma was shown in Table S1. To investigate the relationship of glioma molecular subtype and cytolytic T-cell activity, we compared the difference of CYT value between four molecular subgroups. In CGGA cohort, CYT values were significantly increased (p < 0.0001) in mesenchymal subgroup (Figure 2A). Meanwhile, this trend of CYT distribution was validated in TCGA cohort significantly (p < 0.0001) (Figure 2B). Then Receiver Operating Characteristic (ROC) curves were used to estimate the diagnostic ability of CYT to discriminate the true status of mesenchymal subtype. The area under the curve (AUC) was 92.1 and 91.5% in CGGA (Figure 2C) and TCGA (Figure 2D), respectively. Those results emphasized that CYT was closely connected with malignant progression of glioma and could serve as a biomarker for mesenchymal subtype.

To reveal the molecular mechanism underlying the cytolytic T-cell activity of glioma, we divided gliomas of TCGA dataset into CYT-high cohort (tumors in the top 25th percentile) and CYT-low cohort (tumors in the bottom 25th percentile). We identified the diversity mutational profiles characteristic of gliomas (Figure 3A). Overall, considering all 326 glioma samples, we confirmed 21 genes (Chi-square test, p < 0.05) significantly mutated in CYT-high or CYT-low cohorts. Generally, 91% glioma samples harbored IDH1 mutation, which mostly concurrent with either TP53 and ATRX mutations in astrocytoma samples or with CIC and FUBP1 mutations in oligodendroglioma samples in CYT-low cohort. Meanwhile, 16 genes were founded with significantly higher frequencies of nonsense, splice-site and frameshift mutations (p < 0.05) in CYT-high cohort, which could affect RTK-RAS-PI3K pathway, blood cell hereditary, ATPase activity, ERK pathway, immune system, SWI/SNF complex, and chromatin remodeling.

Besides coding mutations, we also estimated the copy number alterations information of gliomas between CYT-high and CYT-low cohorts (Figure 3B). The most commonly recurrent deletion events on chromosome 1p (ARID1A, CDKN2C, FUBP1) and chromosome 19q (CIC) were observed in gliomas of CYT-low cohort. Meanwhile, the recurrent focal amplifications of chromosome 1q32.1 (MDM4), 4q12 (FIP1L1, PDGFRA), 7p11.2 (EGFR), and 12q14.1(CDK4), and deletion regions of 9p21.3 (CDKN2A, CDKN2B) and 10q23.31(PTEN) were the most frequently identified in CYT-high cohort. In general, glioma samples in CYT-low cohort possessed fewer somatic copy number alterations (SCNA events) and copy number alterations than CYT-high cohort.

To depict the biological features correlated with cytolytic activity, we compared the gene expression profiles of glioma samples between CYT-low group and CYT-high group. Differentially expressed genes were selected with conditions of | log2 (fold-change) | >1 and adjusted p < 0.001. We found 354 down-regulated genes and 1,166 up-regulated genes strongly correlated with CYT in CGGA dataset (Figure S3A), while 534 down-regulated genes and 2,034 up-regulated genes strongly correlated with in TCGA dataset (Figure S3B). To yield an accuracy analysis, the differentially expressed genes which were shared by two datasets (847 up-regulated genes and 188 down regulated genes) were selected for Gene Ontology analysis with online methods (DAVID, https://david.ncifcrf.gov/↗). We found that genes positively correlated with cytolytic T-cell activity were more involved in immune response, regulation of angiogenesis and cell adhesion, dendritic cell chemotaxis, leukocyte migration, monocyte chemotaxis, neutrophil chemotaxis, T cell activation, and negative regulation of B cell activation (Figures 4A,B). Meanwhile, genes negatively correlated with cytolytic T-cell activity were related to normal biology activities, such as chemical synaptic transmission and regulation of ion transmembrane transport, etc. These results indicated that samples in CYT-high group always possessed a more complex and strongly immune response system.

To our knowledge, immune checkpoints were crucial for suppress unwanted and harmful self-directed activities and modulating the duration of immune system. However, in cancer, immune checkpoints could shut down the T-cell powerful fighting activities, allowing tumor to grow undisturbed. Then we explore the relationship of CYT and checkpoints expression (CTLA4, IDO1, LAG3, PD1, PD-L1, PD-L2, TIGIT, and TIM-3). The findings were intriguing, both in CGGA and TCGA datasets, all checkpoints were significantly positively associated with T-cell cytolytic activity (Figures 4C,D). These results indicated that checkpoints were induced more as immune suppressors in CYT-high group in tumor microenvironment when immune cell and immune response were more active. These data provided significantly evidence that CYT may be the representation of immune checkpoints synergistically effect in glioma and might play an important role in the immunosuppressive response through a series of immune checkpoints.

Edema, a reflection of inflammation, was a common pathophysiological entity surrounding gliomas. The secretory products of macrophages were reported might contribute to the formation of glioma edema. Then, we compared the edema extending differences between CYT-low and CYT-high group to investigate the potential association of CYT and peritumoral edema. As shown in Figure 5, we demonstrated that gliomas in CYT-high group significantly possessed more-severe edema than other patients with lower cytolytic T-cell activity in CGGA (p = 0.0003) and TCGA dataset (p = 0.0002), respectively. And we validated this intriguing phenomenon in an independent validation group (p = 0.04). Those results suggested peritumoral brain edema probably was a reflection of CYT and tumor immune response.

To further analysis the prognostic value of CYT, we employed dichotomization to separate cases for depicting the Kaplan Meier survival curves according to the cutoff value (median CYT value). In CGGA dataset, we found that patients with higher CYT value had a significant shorter overall survival than their counterparts in glioma (p < 0.0001) and GBM (p = 0.021) (Figures 6A,B). Simulation results were shown in TCGA dataset that CYT value was a dismal biomarker in glioma (p < 0.0001) and GBM (p = 0.036) (Figures 6C,D). Next, we conducted cox regression model to evaluate the prognostic value of CYT and other prognostic factors. In CGGA dataset (Figure S4A), the result confirmed that CYT value was significantly associated with overall survival (p < 0.0001) after adjusting other factors, such as age at diagnosis (p = 0.54), grade (p < 0.001), IDH status (p = 0.09), and 1p19q status (p < 0.001). In TCGA dataset (Figure S4B), the result suggested that CYT value (p = 0.04) was an independent prognostic indicator. These results demonstrated that CYT value was an independent unfavorable predictor for patients with glioma.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.