PDF
Citation
Print
Introduction
Glioblastoma (GBM) is the most lethal type of human glioma and median survival remains poor following treatment with the current gold standard therapy, which includes maximal surgical resection followed by concurrent chemoradiation and adjuvant temozolomide (TMZ) [1]. Among growing interest in identifying new treatment strategies to improve survival, researchers found that GBM creates an immunosuppressive tumor microenvironment (TME) that facilitates evasion of the host immune system and allows for disease progression. The TME is, in part, a result of angiogenesis and subsequent hypoxia that promote the recruitment and retention of regulatory T cells (Treg) [2,3]. Programmed death-ligand 1 (PD-L1) expression also contributes to the immunosuppressive TME and is the negative prognostic factor for GBM survival [4].
Immunotherapy that targets the programmed death receptor-1 (PD-1) or PD-L1 is an area of active investigation. As a co-inhibitory pathway to T-cell response, PD-1 can bind to PD-L1 that is often expressed on tumor cells, which subsequently downregulates the host immune response. Although pembrolizumab, one of the monoclonal antibodies that blocks PD-1 from binding to PD-L1, yielded unprecedented responses among patients with advanced stage cancers [5], the objective response rate was only 6% to 17%. To improve the response rate, combination with other treatment modalities, including radiation therapy (RT) and angiogenesis inhibitors, are being considered. Also, predictive biomarkers to select for patients that will benefit from immune-based therapy are being explored.
For newly diagnosed GBM, RT is the backbone of standard treatment regimens after maximal surgical resection [1]. RT exerts direct cytotoxic effects on tumors and downregulates immune cells that are sensitive to irradiation, which in turn suppresses the immune response [6]. The concept of radiation-induced tumor equilibrium was recently introduced such that RT might improve the efficacy of immune checkpoint blockade, thereby modulating the immunosuppressive TME towards an anti-tumor TME [7]. In addition, the combination of immunotherapy with an angiogenesis inhibitor could enhance the efficacy of immunotherapeutics because tumor-derived vascular endothelial growth factor A (VEGF-A) has an immunosuppressive function through inhibition of dendritic cell maturation.
In this study, we tested whether the radiosensitivity gene signature [8] and PD-L1 expression level were factors that could predict the clinical outcome of patients receiving RT. We hypothesized that GBM with high PD-L1 expression and radioresistant (RR) was associated with poor survival regarding to RT. Then, we identified a patient subgroup, defined as the PD-L1-high-RR, and characterized the differential expression/methylation pattern of genes and the differential expression of immune cells in the tumor and TME.
Materials and Methods
1. The Cancer Genome Atlas data acquisition
Clinical information from a total of 399 primary GBM patients was collected based on material from previously published research that used The Cancer Genome Atlas (TCGA) lower-grade glioma and GBM cohorts [9] (https://tcga-data.nci.nih.gov/docs/publications/lgggbm_2016/). We used normalized gene expression data from these patients. The detailed normalization process to correct for sample differences is described in the supplementary information of the previous study. In this process, the RALB (RAS Like ProtoOncogene B) gene among the 31 genes was eliminated, and a total of 30 gene signatures were used in the current study. Methylation data was retrieved from the Illumina Human Methylation 450 platform (Illumina, San Diego, CA) and matched into study cohorts using the 'TCGABiolinks' ver. 2.6.9 package [10] of 'R' statistical software (R Foundation for Statistical Computing, Vienna, Austria).
2. Gene signature clustering and PD-L1 grouping
For grouping by PD-L1 status, we determined the cutoff as the median value of CD274 gene expression in the overall study cohort. Patients with CD274 gene expression lower than the cutoff were clustered into the PD-L1-low group, whereas patients with CD274 gene expression higher than the cutoff were clustered into the PD-L1-high group. Previous studies [11,12] used this method to classify the PD-L1-high versus PD-L1-low groups in TCGA cohorts. Specifically, a study [12] that investigated the expression pattern of PD-L1 in 229 glioma samples revealed that patients showing PD-L1 positive expression accounted for 51% of all glioma patients. Thus, our approach that used the median value of CD274 expression among hundreds of patients in the TCGA cohort was acceptable.
For grouping by radiosensitivity, we classified all patients into two groups based on gene signature using consensus clustering (k=2). A total of 1,000 permutation tests, with a subsampling ratio of 0.9, were performed. The optimal number of groups was set to two in order to differentiate the RR and radiosensitive (RS) groups in the current study. Because the reliability of the clustering results and the median values depend on the total number of patients, we performed these grouping processes in the overall cohort, which included RT-treated and non-RT–treated patients. This process is summarized in S1 Fig.
Chi-square tests were used to compare clinical features between the “PD-L1-high-RR” and “others” groups. Clinical information, including baseline characteristics, overall survival (OS) data, and RT information, was obtained from the TCGAbiolinks package [10] in R software. The Kaplan-Meier method was used to compare OS rate between the PD-L1-high-RR and others groups, depending on whether RT had been performed. Cox proportional hazards models were established to identify factors that were significantly associated with the OS for all study cohorts in a univariate analysis. Significant factors were incorporated into multivariable models for RT-treated and non-RT–treated patients to prove the predictive values for receipt of RT. All analyses were performed using ‘R’ ver. 3.3.3 and STATA ver. 14 (StataCorp LP, College Station, TX) statistical software.
3. Transcriptional and epigenetic analysis
Differentially expressed genes (DEGs) for the PD-L1-highRR group were identified using edgeR embedded in the TCGAbiolinks package [10]. We performed this analysis by fitting a negative binomial generalized log-linear model to the read counts for each gene with a threshold to filter genes according to false discovery rate (FDR) cut as 0.05 and log fold change (FC) cut as 1.0. To understand the underlying biological process, gene set enrichment analysis (GSEA) was performed for the PD-L1-high-RR group to identify significant canonical pathways overrepresented by the DEGs. DEGs were mapped into the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, and their corresponding pathways were visualized.
The concept of transcriptional deconvolution of immune cells in tissue samples has been recently introduced. Regarding immune cell infiltration into the TME, a bioinformatics tool that uses gene signature, termed as CIBERSORT, has been developed and validated to predict the immune cell repertoire [13]. In the present study, the proportion of immune cells was enumerated by 'CIBERSORT' for each tumor sample within the PD-L1-high-RR and others groups. The mean proportion of various immune cells in each group was calculated and compared between the two groups.
We performed the entire process of 1,000 permutations and filtered samples that had p < 0.05. Furthermore, quanTIseq [14], which is a deconvolution tool that can quantify ten immune cell types from bulky RNA sequencing data, was used to support CIBERSORT results. Differentially methylated genes (DMRs) were searched for identifying functional regions that are involved in regulating gene expression. The beta-values between two groups were calculated, and the p-values were estimated with the Wilcoxon test adjusted by the Benjamin-Hochberg method. The minimum absolute value for delta beta-values was set to 0.2 and the FDR-adjusted Wilcoxon rank-sum p-value was set to 0.01 for detecting differences. Identified DMRs and DEGs were joined and plotted in the form of a starburst plot. Log10 values (FDR-corrected p-value) for each gene for DNA methylation and expression were used to identify which genes are epigenetically regulated in the PD-L1-high-RR group. All transcriptional and epigenetic analysis was performed using packages from TCGAbiolinks [10].
Results
1. Exploring the PD-L1-high-RR group
In total, 284 patients (28.8%) received RT and 115 patients (71.2%) did not. After we defined the PD-L1-high versus PD-L1-low and RR versus RS, we could obtain combinatorial four subgroups: PD-L1-low-RS, PD-L1-low-RR, PD-L1-high-RS, and PD-L1-high-RR. Significant difference in OS among those four subgroups was observed when RT was given (p=0.002) (Fig. 1A), however, this difference disappeared when RT was not given (p=0.895) (Fig. 1B). To test our hypothesis, we re-grouped initial four combinatorial groups into simply two groups: the PD-L1-high-RR (PD-L1-high-RR) versus the rest of other groups (others). Afterward, log-rank test revealed that OS was significantly different between the PD-L1-high-RR and others groups when those patients were primarily treated with RT (p=0.008) (Fig. 1C). However, this difference was not observed when RT was not given (p=0.572) (Fig. 1D). Thus, we performed following analyses in current study, comparing the PD-L1-high-RR and others groups.
Fig. 1.
Kaplan-Meier curves for overall survival among four combinatorial groups in RT-treated patients (A) and in no RTtreated ones (B). To test our hypothesis, four combinatorial groups were re-grouped into two groups (PD-L1-high-RR vs. others). Again, Kaplan-Meier curves for overall survival comparing two groups are depicted in RT-treated patients (C) and no RT-treated ones (D), respectively. p-value was estimated by the log-rank test. RT, radiation therapy; PD-L1, programmed death-ligand 1; RR, radioresistant; RS, radiosensitive.
2. Patient characteristics
Baseline characteristics from the PD-L1-high-RR and others group are summarized in Table 1. There were no significant differences in age, sex, race, and Karnofsky Performance Scale (KPS) between the two groups. In 2008, a study [15] characterized four transcriptional subgroups in GBM: proneural, neural, classical, and mesenchymal. Generally, the proneural type of GBM shows a survival benefit compared to the mesenchymal type. In the current study, the proneural type was found less in the PD-L1-high-RR group than in others group (4.2% vs. 25.4%, p < 0.001). Furthermore, molecular characteristics were significantly different between the two groups. Isocitrate dehydrogenase (IDH)–wild types were more observed in the PD-L1-high-RR compared to others group (81.5% vs. 69.3%, p=0.004). However, there is no significant difference in 1p19q non-codeletion status between the PD-L1-high-RR and the others groups (p=0.353). GBM with a methylated O6-methylguanine-DNA methyltransferase (MGMT) promoter was identified more in other groups than the PD-L1-high-RR group (33.2% vs. 21.8%, p=0.002).
Table 1.
Patient characteristics
3. Predictive value of radiosensitivity and programmed death-ligand 1 for RT
For all patients, the following factors were associated with OS (Table 2): age, sex, race, KPS, IDH mutational status, 1p19q codeletion status, MGMT promoter status, transcriptome subtype, Estimation of STromal and Immune cells in MAlignant Tumours using Expression data (ESTIMATE) immune score [16], and membership in the PD-L1-high-RR group.
Table 2.
Univariate analysis and multivariate models for overall survival in RT-treated and non-RT–treated patients
RT, radiation therapy; HR, hazard ratio; CI, confidence interval; KPS, Karnofsky Performance Scale; IDH, isocitrate dehydrogenase; MGMT, O6-methylguanine-DNA methyltransferase; ESTIMATE, Estimation of Stromal and Immune cells in MAlignant Tumours using Expression Data; PD-L1, programmed death-ligand 1; RR, radioresistant.
We next investigated whether membership in the “PD-L1- high-RR” group was an independent predictive factor when RT had been performed. A multivariate Cox model for RTtreated patients revealed that membership in the PD-L1- high-RR group (hazard ratio [HR], 1.70; 95% confidence interval [CI], 1.03 to 2.81; p=0.037), IDH-wild type of tumor (HR, 2.91; 95% CI, 1.24 to 6.81; p < 0.001), and male (HR, 1.71; 95% CI, 1.21 to 2.62; p < 0.001) were independently associated with inferior GBM survival. Conversely, membership in the PD-L1-high-RR group was not associated with OS for nonRT-treated patients. Collectively, these results suggested that radiosensitivity gene signature and PD-L1 status were potential predictive factors for survival of GBM patients who receive RT.
4. Differentially expressed and differentially methylated genes
Thirty-gene expression profiles according to the PD-L1- high-RR and the others groups are represented in a heatmap (Fig. 2A), showing the different pattern of transcriptomic profile between them. A total of 200 DEGs were identified between the PD-L1-high-RR and others group. All DEGs, log FC, and FDR are listed in S2 Table. Of these DEGs, the number of downregulated and upregulated genes in the PD-L1-high-RR group was 37 (18.5%) and 163 (81.5%), respectively. With respect to the biologic process, GSEA with these DEGs revealed that upregulated genes in the PD-L1-high-RR group were associated with the immune (common genes, n=19; FDR=1.01e-08) and inflammatory responses (common genes, n=11; FDR=6.55e-05), as shown in Fig. 2B. Furthermore, triggering receptor expressed on myeloid cells 1 (TREM1) signaling (common genes, n=8; FDR=1.97e-06) and pathways involved with chronic inflammation, multiple sclerosis granulocyte, rheumatoid arthritis, were enriched in the PD-L1-high-RR group (Fig. 2C). Other biologic process and pathways enriched in the PD-L1-high-RR are listed in S3 Table.
Fig. 2.
(A) Heatmap comparing the mean fraction of immune cells between the samples in the PD-L1-high-RR and the others groups, showing the differential 30 gene expression pattern. Rows list each of the 30 genes, and columns show each sample. The expression level was normalized between –10 and 10, represented by the colors (green to red). Gene set enrichment analysis and the relationship between differentially expressed genes and methylation-regulated genes. The top 20 enriched biologic processes (B) and pathways (C) are shown. The x-axis indicates a–log10 FDR, the red line represents the ratio of genes in that pathway, and the number in parenthesis is the number of common genes in that pathway. In the starburst plot (D), log10 (FDR-corrected p-value) is plotted for DNA methylation (x axis) and gene expression (y axis). The degree of change (log2 multiple) is indicated by color (red and green). RT, radiation therapy; PD-L1, programmed death-ligand 1; RR, radioresistant; MGMT, O6 -methylguanine-DNA methyltransferase; IDH, isocitrate dehydrogenase; FDR, false discovery rate.
A total of 1,855 DMRs were identified with 1,604 hypermethylated CpG sites (86.5%) and 251 hypomethylated sites (13.5%) in the “PD-L1-high-RR” group. To investigate the epigenetically down or upregulated genes, we constructed a starburst plot showing significant genes with corrected FDR values (Fig. 2D). Gene names are shown in Fig. 2D when they were matched to available methylation probes in TCGA dataset. Detail information including full list of those gene symbol, methylation probeID, difference of methylation beta-value, log FC, FDRs for differential expression/methylation, the status of expression/methylation in the PD-L1-high-RR group are presented in S4 Table. An epigenetically downregulated gene was BAI1 (Brain-specific Angiogenesis Inhibitor 1) (log FC=–0.332, FDR-adjusted p=0.001), while an epigenetically upregulated gene was MAP3K8 (also known as COT or TPL2) (log FC=0.406, p=0.001). In terms of MAP3K8, T-cell receptor signaling and mitogen-activated protein kinase signaling pathways were retrieved from the KEGG pathway database, and DEGs were mapped and visualized in S5A and S5B Fig., respectively. In the PD-L1-high-RR group, MAP3K8 genes were upregulated in both pathways.
5. Deconvolution of immune cells in the PD-L1-high-RR tumor samples
Transcriptional deconvolution of immune cells using CIB-ERSORT were performed in the PD-L1-high-RR and others groups, respectively. With statistically significant samples using Monte Carlo sampling for the deconvolution (p < 0.05), relative percentages of mean fraction of each immune cell are represented in Fig. 3A and listed in S6 Table. M2 type macrophages and Treg cells in tumors were more common for the PD-L1-high-RR group than others group (53.5% vs. 46.5% and 58.9% vs. 41.1%, respectively). On the other hand, the proportion of CD8 T cells was lower in the PD-L1-high-RR than the others group (40.1% vs. 59.8%). Furthermore, we performed deconvolution process using by quanTIseq, showing that M2 type macrophage was more observed in PD-L1-high-RR group compared with the others group (normalized fraction 24.0% vs. 12.1%) (Fig. 3B). However, Treg and CD8 T cells were not detected by quanTiseq. Furthermore, we performed deconvolution analyses in each four subgroup: PD-L1-high-RR, PD-L1-high-RS, PD-L1-low-RR, and PD-L1-low-RS. Detail fractions of immune cells according to four groups are shown in S7 Fig.
Fig. 3.
In the starburst plot (D), log10 (FDR-corrected p-value) is plotted for DNA methylation (x axis) and gene expression (y axis). The degree of change (log2 multiple) is indicated by color (red and green). RT, radiation therapy; PD-L1, programmed death-ligand 1; RR, radioresistant; MGMT, O6 -methylguanine-DNA methyltransferase; IDH, isocitrate dehydrogenase; FDR, false discovery rate. (A) From CIBERSORT, the normalized mean fraction of each of the immune cells found within samples in the PD-L1-high-RR and the others groups are demonstrated in the bar plot. (B) Results from the another deconvolutional tool, quanTIseq, are presented in bar graphs indicating the relative normalized mean fraction of immune cells between the PD-L1-high-RR and others groups. PD-L1, programmed death-ligand 1; RR, radioresistant.
Discussion
We showed that the PD-L1-high-RR is a predictor of inferior OS for GBM patients who receive RT. DEGs in this group were associated with an immune response that constitutes an immunosuppressive TME. Epigenetically regulated genes were related to angiogenesis and key oncogenic survival pathways. Also, immune cells representing the immunosuppressive TME were more enumerated with transcriptome analysis in the PD-L1-high-RR tumor samples compared to other tumor samples.
Higher PD-L1 expression was found to be an adverse prognostic factor for OS and the level of expression significantly correlated with the histologic grade of glioma [4]. A study with 976 glioma samples [17] revealed that a high PD-L1 expression level, as measured by immunohistochemistry, was correlated with immune responses, including T-cell activation and macrophage-related response. Compared to lower grade glioma, GBM showed an up-regulation of PD-L1 and a correlation with angiogenesis, suggesting a relationship between PD-L1 and malignancy. To date, there is no research to analyze clinical outcome in GBM patients with regards to the concept of PD-L1 expression, radiosensitivity, and their associated genomic profiles. Here, we showed that PD-L1 status and radiosensitivity could predict OS in GBM patients who received RT, and we correlated the immuno-suppressive microenvironment in the PD-L1-high-RR group with poor clinical outcome.
This conclusion was supported by the results of GSEA. TREM1 signaling was enriched in the PD-L1-high-RR group of current study. A previous TCGA-GBM cohort study [18] noted that high TREM1 expression was associated with shorter OS, and TREM1 expression was more found in the mesenchymal subtype. Furthermore, genes related to inflammatory and immune responses were differentially expressed in the PD-L1-high-RR group. This is consistent with other studies that reported enrichment of proinflammatory genes within the immunosuppressive mesenchymal subtype of GBM [19]. TCGA-GBM cohort analysis using mRNA sequencing data [20] revealed that immune infiltration, dominantly infiltrated by the tumor-associated macrophages (TAMs) and/or myeloid-derived suppressor cells (MDSCs), was associated with a poor clinical outcome. In the current study, the TME of the PD-L1-high-RR group had more M2 type macrophages, more Treg cells, and less CD8 T cells as compared with the other group. Thus, the immunosuppressive TME resulted from the immune repertoire of the PD-L1-high-RR group could potentially lead to tumor growth and poor OS.
Results
of epigenetic analysis also supported the existence of an immunosuppressive TME in the PD-L1-high-RR group. A variety of methylation-regulated genes related to tumor suppression, apoptosis, sensitivity to chemotherapy, and progression in GBM were found as novel targets [21]. In current study, hypomethylation-upregulated MAP3K8 and hypermethylation-downregulated BAI1 genes were identified in the PD-L1-high-RR group. In GBM, MAP3K8, as one of phosphoinositide 3-kinase pathway cluster genes, is linked with a hazard effect on survival in univariate analysis [22]. BAI1, an inhibitor of angiogenesis, was epigenetically down-regulated in the PD-L1-high-RR group. Angiogenesis is promoted as a consequence of downregulation of BAI1 which expression is epigenetically silenced in glioma [23]. In accord with these results, methylation status and expression level of BAI1 and MAP3K8 potentially contribute to the inferior OS in the PD-L1-high-RR.
GBM is regarded as immunologically cold tumors because its TME impedes infiltrating, proliferating, and activating of effector immune cells by various conditions including eminence of Treg and MDSCs, and hypoxia [24]. Overcoming these conditions is prerequisite for the good response to immunotherapy. One of attempts to reverse cold tumors into hot tumors is combination of immunotherapy and other anticancer treatments such as RT or VEGF inhibitor [25]. Given the correlation between BAI1 and VEGF [26], we speculate that a VEGF inhibitor could be beneficial in the treatment of the PD-L1-high-RR group. Angiogenesis results in the recruitment of immune-suppressive cells, such as M2 type TAMs, MDSCs, or Treg cells. Also, hypoxia-inducible factors such as hypoxia-inducible factor-1α, VEGF-A, or interleukin-10 can induce upregulation of PD-L1 [27]. Thus, control of tumor vasculature with an angiogenesis inhibitor might have a synergistic effect with immune checkpoint blockades [28]. Consequently, clinical trials to test combinational strategies that use an immune checkpoint blockade and an angiogenesis inhibitor are ongoing. For the treatment of recurrent GBM, a phase II trial of pembrolizumab in combination with bevacizumab (NCT-02337491) is underway. Preliminary results suggest that there is no dose limiting toxicity. Another ongoing phase II study is evaluating the efficacy of combination treatment of durvalumab with RT or bevacizumab relative to treatment with durvalumab alone in GBM (NCT02336165). Interim analysis suggested that durvalumab alone shows activity; however, this arm of the study was closed because prolonged OS was not shown.
Cumulating evidences indicate that irradiation can augment innate and adaptive immune response against tumors, thereby potentiating the radiation response. RT can induce an inflammatory microenvironment in the context of immunogenic cell death that is characterized by the release of tumor antigen and danger-associated molecular patterns, which in turn increases the production of proinflammatory cytokines and chemokines [7]. For example, Zeng et al. [29] reported the increased tumor-infiltrated cytotoxic T cells and Treg-induced decreased tumor-inflammation when stereotactic RT and anti–PD-1 blockade were applied to immunecompetent orthotopic murine glioma model. Based on these understandings, several clinical trials are ongoing in patients with GBM. A phase II study (NCT02667587) is currently comparing the combination of RT plus TMZ with nivolumab or placebo in primary GBM patients. A phase I study [30] testing the combination of pembrolizumab, bevacizumab, and hypofractionated RT (30 Gy/5 fractions) in recurrent GBM patients reported that this treatment strategy demonstrated acceptable toxicity. The rate of a durable response was 53%, and 1-year OS was 64%. Based on these results, we speculate that the PD-L1-high-RR group identified in this study may benefit from anti-PD-1/PD-L1 blockade and RT.
There are several limitations to the current study. First, since this was a retrospective study using TCGA cohorts, we sought to analyze the OS data. Given that GBM progresses locally in most patients, analyzing the local, and progression-free survival is also important to solidify the predictive value of the radiosensitivity and PD-L1 status. Nevertheless, we prioritized analyzing the OS over recurrence/progressionfree survival due to the potential bias that arises from variation and limitation with regards to follow-up information in the TCGA database. Regarding recurrences, further prospective studies are needed. Second, quantifying the degree of the radiosensitivity for each patient is difficult. Yet, classification of the PD-L1-high-RR and the others group is intuitive and can serve as a cornerstone for precision medicine. The research described here can form the basis for additional studies characterizing the factors predictive of RT and immunotherapy response.
In conclusion, we validated the predictive value of radio-sensitivity and PD-L1 status, and we identified the PD-L1- high-RR group, which showed a worse outcome following RT in TCGA-GBM dataset. This group demonstrated RT resistance and the immunosuppressive TME through a variety of mechanisms, suggesting that this group could potentially benefit from RT combined with PD-L1 blockade and angiogenesis inhibitor.
XML Download