Integrated analysis of genome-wide DNA methylation and cancer-associated fibroblasts identified prognostic biomarkers and immune checkpoint blockade in lower grade gliomas

Background: Cancer-associated fibroblasts (CAFs) are vital components of prominent cellular components in lower-grade gliomas (LGGs) that contribute to LGGs’ progression, treatment resistance, and immunosuppression. Epigenetic modification and immunity have significant implications for tumorigenesis and development.

Methods: We combined aberrant methylation and CAFs abundances to build a prognostic model and the impact on the biological properties of LGGs. Grouping based on the median CAFs abundances score of samples in the TCGA-LGGs dataset, differentially expressed genes and aberrantly methylated genes were combined for subsequent analysis.

Results: We identified five differentially methylated and expressed genes (LAT32, SWAP70, GSAP, EMP3, and SLC2A10) and established a prognostic gene signature validated in the CGGA-LGGs dataset. Immunohistochemistry (IHC) and in vitro tests were performed to verify these expressions. The high-risk group increased in tumor-promoting immune cells and tumor mutational burden. Notably, risk stratification had different ICB sensitivities in LGGs, and there were also significant sensitivity differences for temozolomide and the other three novel chemotherapeutic agents.

Conclusion: Our study reveals characteristics of CAFs in LGGs, refines the direct link between epigenetics and tumor stroma, and might provide clinical implications for guiding tailored anti-CAFs therapy in combination with immunotherapy for LGGs patients.

1. Introduction

Lower-grade gliomas (LGGs) including World Health Organization (WHO) grade II and III diffuse gliomas are slow-growing infiltrative brain tumors (1). Although the survival of LGGs patients after standardized treatment is better than that of glioblastoma (GBM), recurrent LGGs inevitably progress to GBM (2). With advances in the 2021 WHO Classification of Tumors of the Central nervous system, the understanding of molecular typing for glioma is gradually increasing. Exploration of epigenetics can help us better understand LGGs’ immunity and prognosis.

DNA methylation and gene expression are promising sources for identifying glioma’s molecular biomarkers. For instance, the promoter hypermethylation and epigenetic silencing of the O6-methylguanine-DNA methyltransferase (MGMT) gene have become a classical biomarker for temozolomide resistance glioma (3, 4). DNA demethylation and upregulation of IGF2BP3 can be involved in the malignant progression of glioma (5). Alternated DNA methylation in ZDHHC12 is associated with migration and invasion capabilities in glioma cell lines (6). GPX8 expression was correlated to the reduced DNA methylation at the promoter region and might be related to cancer-associated fibroblasts and immune infiltration levels in glioma (7). However, the clinical impact of these studies remains limited. Either due to the lack of drugs targeting these potential biomarkers or because of a breakthrough in immunotherapy.

Cancer-associated fibroblasts (CAFs) are the significant members of tumor stroma cells in the tumor microenvironment (TME) (8). Research on the significance of CAFs in cancer has recently gained momentum. Accumulating evidence has indicated that CAFs significantly affect tumor progression and migration, promote epithelial-mesenchymal transition (EMT), and induce chemoresistance and immunosuppression (9–12). On the other hand, CAFs and extracellular matrices constitute the tumor immune escape initiation mechanism (13). In response to this problem, harnessing CAFs-related immunosuppressive stromal environment has been proposed to ameliorate the response to immune checkpoint inhibitors (14, 15). However, whether CAFs are associated with the predictive value and immunotherapy of LGGs patients has not been elucidated.

We reasoned that LGGs samples with different CAFs scores broadly alter methylation levels and immune infiltration patterns. To verify the conjecture, we used a median of CAFs score as the grouping basis for the sample to gather genome-wide methylation and gene expression data to locate the altered methylations coupled with altered expression of the same genes. Then we constructed a risk score system containing five risk genes and validated them at tissue-level and cell-level. We found the risk score is an excellent predictive value for survival and a potential factor for immune checkpoint blockade (ICB) therapies. Applying this prognostic gene signature, the sensitivity of GDC0941, Bleomycin, and Axitinib showed a significant difference in sensitivity within the subgroups. These drugs may have different effects on patients with different levels of CAFs infiltration.

2. Material and methods

2.1. Data acquisition

RNA-seq data and clinical data on LGGs were extracted from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) and Chinese Glioma Genome Atlas (CGGA, http://www.cgga.org.cn/). After transcripts per million (TPM) conversion, we held genes with expression levels larger than 0.1 TPM for analysis. Matched Normal data: 105 cortex tissues were obtained from the GTEx project (https://commonfund.nih.gov/GTEx/). We downloaded both methylation data (Illumina Infinium HumanMethylation450 BeadChip) and somatic mutation data generated by TCGA from the UCSC Xena browser (https://xenabrowser.net/hub/) (16). Fibroblast and glioma cell lines from The Cancer Cell Line Encyclopedia (CCLE) (https://portals.broadinstitute.org/ccle) (17). Infiltration Estimation for TCGA-LGGs was collected from TIMER2.0 (http://timer.comp-genomics.org/). The STRING (https://www.string-db.org) database produced the protein-protein interaction network and reconstructed them via Cytoscape software. A flowchart of the entire procedure can be shown in Figure 1.

FIGURE 1

Figure 1 Study Flow Chart.

2.3. Analysis of DNA methylation data

The Illumina HumanMethylation450 BeadChip array contains probes covering 99% of reference sequence genes and 96% of CpG islands. The raw methylation intensities for each probe were represented as β-values, which were converted into M-values with R package Lumi for statistics analysis (18). 5’-C-phosphate-G-3’ (CpG) methylation data between different groups were compared with R package limma to identify differentially methylated CpG sites. Benjamini-Hochberg (BH) method was used to adjust p-value as a false discovery rate (FDR). The CpG site and gene mapping files were downloaded from illumine (https://www.illumina.com/). CpGs span various gene regions, including 1500 bp and 200 bp upstream of the transcription start sites (TSS1500 and TSS200, respectively). The average β-values for each region were calculated according to all CpG sites at the corresponding region, and the average β-value was converted to an M-value. Average regional methylation data between different groups were compared with R package limma to identify differentially methylated regions (DMRs). The hypermethylated DMRs with a threshold of adjusted p-value< 0.05 combined a delta β-value > 0.2, and the hypomethylated DMRs with a threshold of adjusted p-value< 0.05 combined a delta β-value< −0.2. Genes harboring DMRs in any part of the gene features were differentially methylated genes (DMGs). The equations described above are listed below.

$\begin{array}{ll}{M}_i={\log }_2\left(\frac{{\beta }_i}{1-{\beta }_i}\right)& (1)\end{array}$

$\begin{array}{ll} {\beta }_{region}=2^{\left(\sum _{i=1}^k\log 2\left({\beta }_i\right)\right)∕k}& (2)\end{array}$

where M is the intensity of the methylated allele, i = 1, 2, 3,…, k, and k is the number of CpG sites in a region.

2.4. Cancer-associated fibroblasts (CAFs) infiltration estimation and immune score calculation

CAFs abundances were separately estimated via Estimate the Proportion of Immune and Cancer cells (EPIC) algorithm using the R package Immunedeconv (19). TCGA-LGGs were divided into a high-CAFs-score group and a low-CAFs-score group according to the median score. The estimated immune and stromal scores were computed using the R package ESTIMATE.

2.5. Analysis of DEGs and DMGs

Differential expression between the high-CAFs-score and low-CAF-score group samples was analyzed with the R package limma. False discovery rate (FDR) as adjusted p-value using the Benjamini-Hochberg (BH) method. The fold change was log2-transformed. Differentially expressed genes (DEGs) were calculated with a difference > 1.5-fold and p< 0.01. Methylation analysis results were carried out for joint analysis. Venn diagram analyses were performed to calculate the intersection of DMGs and DEGs and explored the differentially methylated and expressed genes (DMEGs). DMEGs were grouped according to four expression patterns: HypoUp, HypoDown, HyperUp, and HyperDown.

2.6. Functional enrichment analyses

To functionally annotate DMEGs of this study, Gene Ontology (GO) including biological process (BP), cellular component (CC), and molecular function (MF) analysis was performed in the R package ClusterProfiler (20). A p-value of< 0.05 and an FDR of< 0.05 were used for the cutoff value. The ClueGO Plugin version 2.5.8 in Cytoscape Version 3.8.2 was employed to identify hub genes and functional analysis (21).

2.7. Construction and validation of the risk score system

We selected HypoDown and HyperUp genes in TSS200 and TSS1500 for univariate cox regression analysis and filtrated the prognostic-related genes. Subsequently, we used the R package glmnet to conduct the least absolute shrinkage and selection operator (LASSO) Cox regression algorithm and develop a potential risk signature. The minimum value of lambda was derived from 1,000 cross-validations (‘1-se’ lambda), and which corresponding partial likelihood deviance value was the smallest for the risk model. Coefficients with regression were confirmed by the “cvfit” function with 1000 repeats. The risk score calculating equation, which contains five risk genes, is:

$\begin{array}{ll}R\text{i}skscor\text{e}={\displaystyle \sum _{i=1}^n}Co\text{e}{f}_{\text{i}}*x_i& (3)\end{array}$

where Coef_i means the coefficients, x_i is the expression value of each gene.

The predictive power of the prognostic signature was evaluated by the receiver-operating characteristic (ROC) curve. The independent clinical factor was validated by multivariate Cox regression analysis. Finally, a nomogram was constructed according to independent predictors. The calibration of the nomogram was evaluated by the calibration curve to assess the goodness of 1-, 3-, and 5-year overall survival.

2.8. Analysis of immunological characteristics

In this study, the mRNA expression matrix of LGGs was analyzed using the CIBERSORT R script downloaded from http://cibersort.stanford.edu. Based on deconvolution, we estimated the abundance of immune cell populations. The relationship between each immune cell and survival was measured by Kaplan-Meier (KM) survival analysis. We evaluated a total of 60 immune checkpoints (ICP) genes in two categories (Inhibitory ICP (22) and Stimulatory ICP (23)) from widely recognized literature (24). Then we assessed the expression and survival of these ICP in the TCGA-LGGs cohort for a comprehensive overview of the immunosuppressive landscape.

The Tumor Immunophenotype Profiling (TIP) was performed to quantify the extent of infiltrating immune cells and anticancer immunity (25). Assessment of antitumor immunity was conceptually divided into seven steps, including tumor cell antigen release (step 1), cancer antigen presentation (step 2), priming and activation (step 3), trafficking of immune cells to tumors (step 4), infiltration of immune cells (step 5), T cell recognition of cancer cells (step 6), and killing of cancer cells (step 7).

TIDE (http://tide.dfci.harvard.edu/), an excellent algorithm, was used to explore the prediction of clinical response to immune checkpoint blockade (ICB) therapy (22). The TIDE score was calculated to simulate two mechanisms of tumor immune evasion: the induction of T cell dysfunction with high infiltration of cytotoxic T lymphocytes and the retard of T cell infiltration in tumors with low cytotoxic T lymphocyte infiltration. The TIDE score is a good reflection of the responsiveness of the ICB. The SubMap (https://www.genepattern.org/) was carried out to validate the reliability of the prediction of TIDE.

The R package pRRophetic was used to predict chemotherapeutic response in LGGs patients (26). In addition to temozolomide, which patients with glioma widely used, we included three drugs with the therapeutic potential for glioma in this study: Axitinib, GDC-0941 (PI3K inhibitor), and Bleomycin.

2.9. Verification of gene expression at cellular level and tissue level

H4, SW1783, and HMF cell lines were purchased from ATCC. According to the manufacturer’s instructions, total RNA was isolated using Trizol reagent (Invitrogen, USA). 2μg of the total RNA was transcribed into cDNA. SYBR Green PCR kit (Takara, Japan) was used for qRT-PCR. We selected the 2−ΔΔCq method to calculate gene transcription level, with β-actin mRNA as control. Data represent the mean ± SD of triplicate real-time PCR. The primers were synthesized by Tsingke Biotechnology (Shanghai, China) and displayed in Supplementary Table S1. Immunohistochemistry (IHC) analyzed the protein expression levels. GSAP (ab106630), LATS2 (ab111054), SWAP70 (ab228846), and SLC2A10 (ab110528) antibody was purchased from Abcam. EMP3 (sc-81797) antibody was purchased from Santa Cruz Biotechnology. Clinical characteristics of LGGs patients cohort are displayed in Supplementary Table S2. All the patients and the hospital’s Ethics Committee approved this research. LGGs tissues were formalin-fixed, paraffin-embedded, and sectioned at 4 µm. Immune complexes were detected with the SP Kit (Solarbio, Beijing, China) and DAB Substrate Kit (Solarbio, Beijing, China). Signals were detected using an Olympus BX41 microscope. Quantification of Immunohistochemistry (IHC) staining was performed in a blinded fashion.

2.10. Statistical analysis

All the data were analyzed using the R software version 4.1.0. The overall survival (OS) between different groups was analyzed using Kaplan-Meier curves. Kruskal–Wallis tests were applied to compare gene expression in two groups. The Fisher test assessed different groups’ responses to ICB treatment. Somatic mutation data sorted in the form of Mutation Annotation Format (MAF) was analyzed using the R package maftools. ImageJ and ImageJ plugin IHC profiler was applied to quantify IHC staining analysis. IHC scoring data and qRT-PCR data were analyzed using GraphPad/Prism 9.0. In addition, tumor mutational burden (TMB) and mutation counts were computed from somatic mutation frequencies. p< 0.05 was marked as ‘*’, p< 0.01 was marked as ‘**’, p< 0.001 was marked as ‘***’. and p< 0.0001 was marked as ‘****’.

3. Results

3.1. Differentially methylated and expressed genes (DMEGs) in LGG

To identify DMEGs in LGGs, we first extracted the gene expression and DNA methylation data of TCGA-LGGs and performed a comparative analysis. Samples were divided into high- and low-CAFs-score groups according to the median CAFs score. From the summary estimate, 2393 statistically significant Differentially Expressed Genes (DEGs) were identified (difference > 1.5-fold, p-value< 0.01), including 263 upregulated and 2131 downregulated genes (Figure 2A, Table S3). Promoter regions (TSS200 and TSS1500) were enrolled in the primary study. As shown in volcano plots, 1276 DMGs were identified from two regions, including 896 DMGs in the TSS200 region (Figure 2B, Table S4) and 583 DMGs in the TSS1500 region (Figure 2C, Table S5).

FIGURE 2

Figure 2 Identification and functional enrichment of DMEGs. (A) Volcano plot presenting the DEGs (difference > 1.5-fold, p-value< 0.01) between two groups. (B) DMGs (difference > 1.5-fold, p-value< 0.01) volcano plot of TSS200 region. (C) DMGs (difference > 1.5-fold, p-value< 0.01) volcano plot of TSS1500 region. (D) Venn diagram exhibiting DMEGs expressed in the TCGA dataset. (E) The profiles of DMEGs and principal component analysis (PCA) between tumor and normal cortex tissues. (F) GO analysis shows significant GO terms in DMEGs.

We analyzed the relationship between methylation and gene expression by integrating DMGs and DEGs in two promoter regions (TSS1500 and TSS200). A total of 77 DMEGs were identified (Figure 2D), and then we performed a principal component analysis (PCA) of the DMEGs in normal tissue and LGGs (Figure 2E). The PCA profile revealed a clear separation of normal samples from LGGs. Afterwards, GO enrichment analysis was performed on these DMEGs, and the result indicated that DMEGs were significantly enriched in positive regulation of oligodendrocyte differentiation, negative regulation of canonical Wnt signaling pathway, and regulation of epithelial to mesenchymal transition (Figure 2F). These functions were closely related to the oncogenesis and progression of glioma.

3.2. DMEGs analysis in two promoter regions

To investigate the differences in DMEGs within each region, we classified these DMEGs into four groups (HypoUp, HyperUp, HyperDown, and HypoDown) for TSS200 and TSS1500, respectively (Figures 3A, B). The HyperDown group was the most prevalent in TSS200 and TSS1500 regions, and the HypoUp group had the second-highest proportion in the two regions (Figure 3C). After extracting the HyperDown and HypoUp DMEGs in TSS200 and TSS1500 regions separately, we used STRING to construct PPI networks. These genes were analyzed by GO enrichment using the “ClueGO” plugin for Cytoscape software (p< 0.01, Kappa score = 0.5). Functional enrichment analysis revealed that these DMEGs participated in critical biological processes. In the TSS200 region, DMEGs are mainly associated with response to mechanical stimulus, regulation of cell-substrate adhesion, and positive regulation of macrophage migration (Figure 3D). As shown in Figure 3D, MMP14 is critically involved in various immune-related functions. Meanwhile, DMEGs in the TSS1500 region are mainly related to regulating epithelial to mesenchymal transition, programmed necrotic cell death, and positive regulation of gliogenesis (Figure 3E). These simultaneous differential DEGs and DMGs pooled to DMEGs may be the main factor causing the altered biological function of LGGs.

FIGURE 3

Figure 3 Grouping and functional analysis of DMEGs. (A, B) Venn diagram showed four different groups (HypoDown, HyperUp, and HyperDown) of DMEGs in the TSS200 region (A) and TSS1500 region (B). (C) The bar chart shows the different groups of DMEGs in the two regions. (D, E) ClueGO Cytoscape network of statistically DMEGs in two regions.

3.3. Construct the DMEGs prognostic signature

Among the DMEGs in HyperDown and HypoUp groups, Univariate Cox regression analysis screened 203 DMEGs with prognostic values. Then a Lasso‐penalized Cox analysis was performed to shrink further the scope of DMEGs screening (Figure 4A) and lambda. Min was regarded as the optimal value in the cross-validation process (Figure 4B). Five DMEGs (GSAP, EMP3, LATS2, SWAP70, and SLC2A10) and corresponding coefficients (Table S6) were identified. We used TCGA-LGGs as the train set CGGA-325 and CGGA-693 as the test sets, and samples were split into high- and low-risk groups by the median value of the risk score. KM survival curves depicted that LGG patients with increased risk scores had worse clinical outcomes in both train set and test sets (Figures 4C, D, and Figure S1A). Statistical analysis was performed using a log-rank test (train set p< 0.001, test set p< 0.001). After that, we established ROC curves of the risk score model with 1-year, 3-year, and 5-year. The results revealed that the risk score could effectively distinguish LGGs patients with different survival statuses (Figure 4E 1-yer AUC = 0.86, 3-year AUC = 0.83, and 5-year AUC = 0.80). The results were similar and slightly lower in the test set (Figure 4F 1-yer AUC = 0.73, 3-year AUC = 0.79, and 5-year AUC = 0.77). Multivariate Cox regression analysis showed the independent prognostic value of this risk score (Figure 4G, p< 0.001, HR = 3.844). Then we examined the mRNA expression of five risk genes in the glioma cell lines and fibroblast cell lines with CCLE (Figure 4H). Although glioma lines are not representative of LGGs, we found the expression of GSAP was consistent with that of fibroblasts, and the expression of the other four risk genes was lower than that of fibroblasts in glioma cell lines. The mortality of patients increased with the increase in the risk score (Figure 4I). Finally, we established the Nomogram model, which contained risk score and age to assess the survival prediction in LGGs patients (Figure 4J) and calibration curves for nomogram predicted 1-, 3- and 5-year overall survival performed well with the risk model (Figure 4K).

FIGURE 4

Figure 4 Construction of prognostic signature. (A, B) The process of building the signature. The least absolute shrinkage and selection operator (LASSO) regression was performed, calculating the minimum criteria. (C, D) K-M curves showed that the high-risk subgroup had worse overall survival than the subgroup in the train set (p< 0.001) and test set (p< 0.001). (E, F) ROC curves showed the predictive efficiency of the risk signature on the 1-, 3- and 5-year survival rates of train set (E) and test set (F). (G) Independent prognostic factors were determined by the multivariate Cox regression analyses. (H) Expression of five risk genes in fibroblast cell lines and glioma cell lines in the CCLE dataset. (I) The heatmap was based on the expression of the five genes in the high- and low-risk group. (J, K) A nomogram (J) and decision curve analysis (K) of the risk score for predicting 1-, 3- and 5-year survival. (*p < 0.05, **p <0.01, ***p <0.001, ****p <0.0001).

3.4. Correlation between risk score and clinical characteristics

Clinical variables were introduced into the risk score system to analyze the relationship between risk scores and clinical characteristics. We initially compared the CAFs’ scores concerning the risk score. As shown in Figure 5A, the risk score was mildly positively associated with the CAFs score (r = 0.42, p< 0.01). In contrast, tumor purity was negatively correlated with the risk score (Figure 5B r = -0.56, p < 0.01). Other clinical features were then introduced. WHO grade III LGGs have a higher risk score than WHO grade II LGGs, which is consistent with their poor prognosis (Figures 5C, D, Figure S1B). The risk score in MGMT promoter methylated LGGs were significantly lower than that in MGMT promoter unmethylated LGGs samples (Figures 5E, F, Figure S1C). It is widely accepted that glioblastoma patients with MGMT promoter methylated are sensitive to temozolomide and suitable for TMZ chemotherapy. For another crucial molecular marker 1p19q codeletion status, the risk score was significantly higher in 1p19q-non-codeletion samples compared to codeletion samples in both datasets (Figures 5G, H, Figure S1D). Correspondingly, IDH1/2 wild-type cases showed a valid increased risk score compared to IDH1/2 mutant cases (Figures 5I, J, Figure S1E). The high-risk scores were seen in cases aged > 45 years. (Figures 5K, L, Figure S1F). Overall, some important molecular markers and clinical features of LGGs responded well to this risk score system.

FIGURE 5

Figure 5 Clinical characteristics and risk scores. (A) Correlation of risk scores with CAFs scores, r = 0.42 p< 0.001. (B) Correlation of risk scores with tumor purity, r = -0.56, p< 0.001. (C, D) Risk scores for different WHO-graded samples, (E, F) MGMT promoter status, (G, H) 1p19q codeletion status, (I, J) IDH mutation status, and (K, L) Age effect. (*p < 0.05, **p <0.01, ***p <0.001, ****p <0.0001).

3.5. Prognostic signature and immune landscape

To investigate the relationship between the prognostic signature and immune cell infiltration in LGGs. We evaluated immune scores, immune cell infiltration, and immune checkpoints separately. CAFs are the most prominent tumor stroma cell type in the TME. Comparing the stromal score and immune score of the LGGs datasets, we found a significantly positive correlation between the risk score and stromal score, and the same was true of the immune score (Figures 6A, B, stromal score p< 0.001 and Figures 6C, D immune score p< 0.001). CIBERSORT algorithm showed the proportions of distinct immune cell subpopulations in different risk groups (Figure 6E). The relative expression is shown in bar diagrams. The proportions of Macrophages M0 and Macrophages M2 in the high-risk group were significantly higher than in the low-risk group (Figure 6F). By contrast, the proportions of Monocytes, activated NK-cell and activated Mast-cells were higher in the low-risk group. Subsequently, we performed KM survival analysis to evaluate the OS with differing immunocytes infiltration samples (Figures 6G–I). High proportions of activated NK-cell and activated Mast-cell had a better OS (Figures 6G, I). Conversely, the level of macrophage M0 expression is inversely related to OS (Figure 6H). We next examined whether the risk score is associated with the expression of inhibitory and stimulatory immune checkpoint (ICP) molecules. Fortunately, the expression of many immune checkpoint molecules showed significant differences between high and low-risk groups (Figure 6J Inhibitory ICP and Figure 6K stimulatory ICP). The KM curves for survival analysis of these immune checkpoints were presented in Figure S2. In summary, the risk score was associated with the expression of these immune checkpoint molecules.

FIGURE 6

Figure 6 Analysis of immune infiltration in different risk groups. (A–D) Immune score and stromal score calculated by ESTIMATE in different groups (A, C TCGA data set, and B, D CGGA data set). (E, F) The relative infiltrating proportion of 22 immune cells in high- and low-risk groups. (G–I) KM curves of infiltrating immune cells associated with survival in LGGs patients (p< 0.05). (J, K) Immune checkpoints expression in the LGGs microenvironment, inhibitory ICP (J), and stimulatory ICP (K). (*p < 0.05, **p <0.01, ***p <0.001, ****p <0.0001).

3.6. Risk score-based stratification predicts the immune response and chemotherapy efficacy

To explore the risk score stratification and the associated characteristics of the antitumor immune response, we introduced the Tracking Tumor Immunophenotype (TIP) system, which analyzed the status of anticancer immunity and the proportion of tumor-infiltrating immune cells across seven-step Cancer-Immunity Cycle. As shown in Figure 7A, the scores of the high-risk group were significantly higher in tumor antigen release (step 1), immune cell recruitment (step 4), and immune cell infiltration into the tumor (step 5). Conversely, the low-risk group scored significantly higher in T cell priming and activation (step 3), T cell recognition of cancer cells (step 6), and killing of cancer cells (step 7). Then, we used the CGGA dataset to validate these results (Figure S3A). Next, we introduce the TIDE algorithm to assess the efficacy of risk score in predicting ICB responsiveness. We found there was a significant difference in response to ICB treatment between the two groups (p< 0.001), and the response to ICB treatment was more sensitive in the low-risk group (Figures 7B, C). SubMap was used to compare the prediction response to anti-PD1 and anti-CTLA4 therapy results with another dataset containing 47 patients with melanoma that responded to immunotherapies. Using this tool, we found that the high-risk group in the train and test sets showed comparable performance in predicting the LGGs’ response to anti-PD1 therapy (Figures 7D, E p< 0.05). Anti-CTLA4 therapy also showed a partial response in the TCGA train set (Figure 7D Bonferroni corrected p = 0.32).

FIGURE 7

Figure 7 Risk score-based analysis of the stratifiable immune response and chemotherapy efficacy. (A) The TIP system quantified seven steps of the antitumor immune response. (B, C) Predicted response of TIDE to immune checkpoint inhibitors. (D, E) Comparing the effectiveness of PD1 and CTLA4 in response to different risk groups. (F–I) The sensibility of chemotherapeutic drugs (F) Temozolomide, (G) Axitinib, (H) GDC0941, and (I) Bleomycin. (*p < 0.05, **p <0.01, ***p <0.001, ****p <0.0001).

We tried to analyze the response of two risk groups to chemo-drug efficacy. Here, in addition to TMZ, we selected three other chemo drugs from the literature that may have therapeutic potential for glioma. As expected, the sensitivity in the low-risk group was slightly better than that in the high-risk group (Figure 7F). That may be related to the level of MGMT promoter methylation in the low-risk group (Figure 5E). For the other three drugs, the estimated IC50 was significantly better in the low-risk group (Figures 7G–I). These results were validated using the CGGA325 dataset (Figure S3B–E). Although these drugs have not been used in clinical treatment on a large-scale, differences in sensitivity suggest that they have potential as novel therapeutic agents.

3.7. Genomic alterations of prognostic signature

Tumor genomic alterations have profound effects on immunity and drug therapy. We investigated the mutation frequencies of different risk groups and showed the top ten most frequently mutated genes (Figures 8A, B). As we expected, IDH1 had the highest mutation frequency, predominantly missense. Remarkably, some genes were associated with the immune microenvironment and immunotherapy. These genes were more prominent in the high-risk group, such as ATRX, EGFR, and PTEN. We targeted IDH1, the most frequently mutated of the two groups in our study and analyzed the relationship between the expression of five risk genes and IDH1 mutations. The results were shown in Figures 8C–G. The expression of all risk genes was higher in samples with wild-type IDH1 than in mutant IDH1. These risk genes and corresponding risk scores are consistent with the findings summarized in the preceding text (Figure 5I). In addition, we analyzed correlations between tumor mutational burden (TMB) and risk score. Like previous research, LGGs patients in the high TMB group have a poorer prognosis. Hence, we introduced the risk score and analyzed it jointly with TMB. Our results show that the risk score had a low positive correlation with TMB values (Figure 8H). Meanwhile, the high-risk group with a poorer prognosis corresponded to high TMB values (Figure 8I). Combining the two elements in the analysis, we found that LGGs samples with high-risk scores and high TMB values had the worst prognosis (Figure 8J). Finally, we compared the TIDE, dysfunction, and exclusion scores between the different risk groups. The TIDE score in patients with low-risk scores was significantly higher than those with high-risk scores (Figure 8K p< 0.01). In parallel, patients in the high-risk group had higher dysfunction scores constructed using dysregulated immune genes (Figure 8L p< 0.001). No difference was found between the two groups for exclusion scores constructed using immune rejection genes (Figure 8M). The mechanism of immune escape in these high-risk LGGs samples may be immune dysfunction rather than immune exclusion. These epigenetic alterations may affect the prognostic model for the therapeutic assessment of LGGs patients.

FIGURE 8

Figure 8 Epigenetic analysis of risk genes. (A, B) The mutation of high- and low-risk groups (Top 10 mutated genes). (C–G) Expression of 5 risk genes in the presence of different IDH1 mutations status. (H) Correlation analysis of risk score and TMB. (I) K-M curves of the high-TMB and low-TMB groups (p< 0.001). (J)The combined risk score and TMB analysis of K-M curves in LGGs patients. (K–M) TIDE algorithm to model tumor immune evasion, (K) TIDE score, (L) Dysfunction score, and (M) Exclusion score.

3.8. Experimental verification in cell lines and tissues

After obtaining the above five risk genes, we identified them at the cellular and protein levels. Given the rarity of fibroblasts in brain tissue, we chose a stable human mammary fibroblast (HMF) cell line as a control group. T98G and U251 cell lines are commonly used glioma cell lines for experiments. GSAP and SWAP70 expressed similar levels in different cell lines, and the expression of the remaining three genes was lower in glioma cell lines (Figures 9A–E). For protein expression, three patients were analyzed with IHC. LGGs tissues were obtained from the tumor center (TC) and tumor periphery (TP). We found the expression of five proteins was highly expressed in tumor periphery (Figure 9F). After a 4-step grading system was quantified, except for LATS2, other proteins showed high expression in the TP group (Figures 9G–K). This founding may be related to the research that tumor cells can influence the recruitment of CAFs precursors and induce the activation of normal fibroblasts into CAFs.

FIGURE 9

Figure 9 Experimental verification in cell lines and tissues. (A–E) The mRNA levels were determined by qRT-PCR in three cell lines. (F) The expression of five protein comparison of tumor center (TC) vs tumor periphery (TP) sites. (G–K) Bar graph showing the five protein levels obtained by quantification of immunohistochemical images. (*p<  0.05, **p < 0.01, and ***p < 0.001).

4. Discussion

Abnormal epigenetic alterations contribute to tumorigenesis and progression, as reflected in the latest guidelines for glioma classification (27). DNA methylation has been found to regulate microRNAs and predict overall survival in glioma (28, 29). Recurrence of LGGs occurs mainly within a few centimeters of the resected cavity, even in complete tumor resection and adjuvant chemotherapy (30, 31). Glioma recurrence and prognosis are closely related to alterations of TME (32, 33). Immune cells and CAFs are essential components of the TME. Here, we combined abnormal methylation and CAFs abundance for an in-depth analysis of LGGs, which is essential for a more comprehensive understanding of TME and developing stromal CAFs-targeted therapies.

Herein, we analyzed DMEGs in different CAFs abundance groups, and we found that the functional annotations of these DMEGs were enriched in pathways of tumorigenesis, progression, and malignant transformation. Compelling evidence shows that the extracellular matrix acts as the “soil” for malignant tumor progression and immune resistance (34). Among them, functions closely related to CAFs, such as the Wnt signaling pathway plays a vital role in the carcinogenic activity of LGGs (35). Oligodendrocyte differentiation reflects the stemness of glioma cells (36). EMT has been implicated in cancer stemness, invasiveness, and drug resistance (23). The genes with simultaneous alterations in gene expression and methylation levels may be factors that alter LGGs functions through CAFs. DNA methylation at these promoter regions is widely known to correlate negatively with gene expression levels (37). GO enrichment analysis revealed the main functions and hub genes involved in DMEGs. It is worth noting that in TSS200 and TSS1500 regions, MDK is a critical player in cancer progression and immune microenvironment (38). MMP14 regulates the activity of multiple extracellular and plasma membrane proteins, influencing cell-cell and cell-extracellular matrix communication (39).

After integrating clinical information, we constructed a prognostic signature based on five genes (EMP3, GSAP, LATS2, SLC2A10, and SWAP70). Compared to fibroblast cell lines and glioma cell lines using the CCLE database, we found that the expression of GSAP was similar. In the validation of in vitro cell experiments, in addition to GSAP, SWAP70 expression was also similar to fibroblasts. However, no evidence exists that any glioma cell line can represent LGGs or GBM. Its predictive value appears to be quite weak. In order to test protein expression, we performed sampling and IHC analysis in the center and periphery of LGGs. The expression of EMP3, GSAP, SLC2A10, and SWAP70 was higher in tumor periphery. This finding suggests that there may be more activated fibroblasts at the TP, and CAFs could function at tumor periphery. However, there is no significant difference in LATS2 expression between TC and TP. Despite LATS2 having been recognized as a target gene of CAFs-derived exosome microRNA-92 in breast cancer (40). SLC2A10 regulated fibroblasts in arterial tortuosity syndrome by encoding glucose transporter 10 (41). The PI3K-dependent recruitment of SWAP70 to the plasma membrane has been observed in growth factor-stimulated fibroblasts (42). EMP3 plays an important role in the regulation of membrane receptors associated with IDH-Wild type glioblastoma (43).

More and more research focused on the immuno-phenotype and immunotherapy of glioma cells. The high-risk scoring group showed increased antitumor immune cells macrophage M2 and M0. Despite glioma being defined as a cold tumor, proportions of macrophages can still constitute up to 30–50% of the TME (44–46). Surprisingly, mast cells activated were higher in the low-risk group. We quantified antitumor immunity in seven steps and further evaluated the antitumor immune process. Increased risk in the LGGs sample was accompanied by a decrease in the score of T cell priming and activation (step3) and destruction of tumor cells (step6, 7). It corresponds to a higher density of antitumor immune infiltration in the high-risk group. Immunotherapy, especially ICB, has brought paradigm shifts to cancer treatment. We found that the specific inhibitory immune checkpoints PD1 and CTLA4 were significantly overexpressed in the high-risk group (Figure 6J), while the Submap approach suggests that the high-risk group showed promising performance in predicting LGGs predicted response to anti-PD1 and anti-CTLA4 therapies (Figures 7D, E). Other immune checkpoints that are highly expressed in high-risk groups are also being studied in an expanding way, such as the inhibitory immune checkpoints CD276 (47) and IL-10 (48), and the stimulatory checkpoints ICOS (49) and CD40 (50) are involved in the regulation of T cell function. TGFB1 can increase endothelial and epithelial permeability. It is more inclined to promote GBM cell invasion (51). SLAMF7 is a cell surface receptor involved in natural killer cell activation that received approval for treating multiple myeloma (52). Treatment by blocking or stimulating these new checkpoints in LGGs holds the promise of going beyond traditional PD1 therapies. Although our risk score distinguished the expression of ICPs of LGGs and predicted anti-ICB therapy, it still lacks elucidation of the interaction mechanism between CAFs and ICB. We hope to provide new ideas on the relationship between the treatment of immune checkpoints and CAFs.

TMZ has become the conventional chemotherapeutic agent for glioma however, TMZ resistance is the main factor that leads to current studies aimed at expanding multiple chemotherapeutic agents for glioma (53). Besides TMZ, three promising drugs were introduced in this study (Figures 7G–I). Axitinib induces senescence-associated cell death and necrosis in glioma (54). The PI3K inhibitor GDC-0941 enhances radiosensitization and reduces chemoresistance to temozolomide in GBM (55). Bleomycin inhibits proliferation and promotes apoptosis of glioma via the TGF-β/Smad signaling pathway (56). All three drugs showed promising IC50 in the low-risk group, and we hope this study will provide potential directions for the relationship between new chemotherapeutic agents and CAFs. However, these speculations are still at the level of data analysis, and whether these drugs can be applied to LGGs still needs a lot of experimental verification, such as molecular docking.

Solid evidence suggests that TMB plays a vital role in tumor immune escape (57). TP53 was frequently mutated in the high-risk subtype, and its mutation was reported to be associated with a poorer prognosis. Mutations in IDH1 characterize the majority of lower-grade gliomas in adults and define a subtype associated with a favorable prognosis (58, 59). Glioma shows a markedly elevated mutation burden, referred to as TMB-H (60). A study suggests that some gliomas with high TMB may benefit from PD-1 blockade therapy (61). Interestingly, MMR deficiency gliomas with TMB-H also lack significant inflammatory CD8+ infiltrates (62). On the other hand, TIDE algorithm analysis revealed that the mechanism of immune escape in LGGs samples might be immune dysfunction (Figure 8L). Even in the presence of a high level of infiltration by cytotoxic T cells, immune escape is still inevitable (63). However, it is worth noting that our study also has some limitations. Although our study identified five risk genes, they could not directly serve as CAFs marker genes in LGGs. We should further confirm the role of these CAFs risk genes on glioma cells using in vitro co-culture or single-cell multi-omics. Bioinformatics has always been used to research CAFs, but how gliomas induce the production of CAFs and exercise their functions in TME still requires extensive in vivo or in vitro experimental studies. These efforts will be included in our future studies.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Materials. Further inquiries can be directed to the corresponding authors.

Ethics statement

The studies involving human participants were reviewed and approved by The Animal Ethical and Welfare Committee of Zhejiang Provincial People’s Hospital. The patients/participants provided their written informed consent to participate in this study.

Author contributions

SH, JWD, and HJ conceived and designed the study and revised the manuscript. JWD, HZ, JJ, XG, and XY provided analytical technical support. ZL, NW, and JZ participated in the production of charts and pictures. JWD and FW designed and completed qRT-PCR and IHC experiments. YB, JYD, and SM revised the manuscript. SH supervised the study. All authors contributed to the article and approved the submitted version.

Funding

This work was funded by the National Natural Science Foundation of China (No. 61575058).

Acknowledgments

The authors gratefully acknowledge TCGA, CGGA, and other public databases for their facilitation of data analysis.

Conflict of interest

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.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2022.977251/full#supplementary-material

Supplementary Figure 1 | Clinical characteristics and risk scores in CGGA-693 LGGs test set. (A) K-M curves showed that the high-risk subgroup had worse overall survival than the subgroup in the test set (p< 0.001). (B) K-M curves of different risk groups in WHO II and WHO III gliomas, (C) MGMT promoter status, (D) 1p19q codeletion status, (E) IDH mutation status, and (F) age effect.

Supplementary Figure 2 | Overall survival analysis on immune checkpoint for LGGs patients. Grouping by median expression of immune checkpoint genes.

Supplementary Figure 3 | Immunotherapy response and chemotherapy sensitivity in the CGGA. (A) TIP system analyzed the status of anticancer immunity and the proportion of tumor-infiltrating immune cells across the seven-step Cancer-Immunity Cycle. (B) Chemosensitivity of TMZ, p = 0.11. (C) Chemosensitivity of Axitinib, p< 0.001. (D) Chemosensitivity of GDC0941, p< 0.001. (E) Chemosensitivity of Bleomycin, p = 0.004.

Supplementary Table 1 | The primers used in this study.

Supplementary Table 2 | Clinical characteristics of patient cohort.

Supplementary Table 3 | DEGs in TCGA.

Supplementary Table 4 | DMGs in TSS200 region

Supplementary Table 5 | DMGs in TSS1500 region

Supplementary Table 6 | The formula of calculate risk score.

References

1. Ostrom QT, Cioffi G, Waite K, Kruchko C, Barnholtz-Sloan JS. CBTRUS statistical report: Primary brain and other central nervous system tumors diagnosed in the united states in 2014-2018. Neuro-oncology (2021) 23(12 Suppl 2):iii1–iii105. doi: 10.1093/neuonc/noab200

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Weller M, van den Bent M, Tonn JC, Stupp R, Preusser M, Cohen-Jonathan-Moyal E, et al. European Association for neuro-oncology (EANO) guideline on the diagnosis and treatment of adult astrocytic and oligodendroglial gliomas. Lancet Oncol (2017) 18(6):e315–e29. doi: 10.1016/s1470-2045(17)30194-8

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Yin D, Xie D, Hofmann WK, Zhang W, Asotra K, Wong R, et al. DNA Repair gene O6-methylguanine-DNA methyltransferase: promoter hypermethylation associated with decreased expression and G:C to A:T mutations of p53 in brain tumors. Mol carcinogenesis. (2003) 36(1):23–31. doi: 10.1002/mc.10094

CrossRef Full Text | Google Scholar

4. Hegi ME, Diserens AC, Godard S, Dietrich PY, Regli L, Ostermann S, et al. Clinical trial substantiates the predictive value of O-6-methylguanine-DNA methyltransferase promoter methylation in glioblastoma patients treated with temozolomide. Clin Cancer research: An Off J Am Assoc Cancer Res (2004) 10(6):1871–4. doi: 10.1158/1078-0432.ccr-03-0384

CrossRef Full Text | Google Scholar

5. Nomura M, Saito K, Aihara K, Nagae G, Yamamoto S, Tatsuno K, et al. DNA Demethylation is associated with malignant progression of lower-grade gliomas. Sci Rep (2019) 9(1):1903. doi: 10.1038/s41598-019-38510-0

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Lu F, Shen SH, Wu S, Zheng P, Lin K, Liao J, et al. Hypomethylation-induced prognostic marker zinc finger DHHC-type palmitoyltransferase 12 contributes to glioblastoma progression. Ann Trans Med (2022) 10(6):334. doi: 10.21037/atm-22-520

CrossRef Full Text | Google Scholar

7. Ren Z, He Y, Yang Q, Guo J, Huang H, Li B, et al. A comprehensive analysis of the glutathione peroxidase 8 (GPX8) in human cancer. Front Oncol (2022) 12:812811. doi: 10.3389/fonc.2022.812811

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Tommelein J, Verset L, Boterberg T, Demetter P, Bracke M, De Wever O. Cancer-associated fibroblasts connect metastasis-promoting communication in colorectal cancer. Front Oncol (2015) 5:63. doi: 10.3389/fonc.2015.00063

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Fiori ME, Di Franco S, Villanova L, Bianca P, Stassi G, De Maria R. Cancer-associated fibroblasts as abettors of tumor progression at the crossroads of EMT and therapy resistance. Mol cancer. (2019) 18(1):70. doi: 10.1186/s12943-019-0994-2

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Wu X, Tao P, Zhou Q, Li J, Yu Z, Wang X, et al. IL-6 secreted by cancer-associated fibroblasts promotes epithelial-mesenchymal transition and metastasis of gastric cancer via JAK2/STAT3 signaling pathway. Oncotarget (2017) 8(13):20741–50. doi: 10.18632/oncotarget.15119

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Lotti F, Jarrar AM, Pai RK, Hitomi M, Lathia J, Mace A, et al. Chemotherapy activates cancer-associated fibroblasts to maintain colorectal cancer-initiating cells by IL-17A. J Exp Med (2013) 210(13):2851–72. doi: 10.1084/jem.20131195

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Monteran L, Erez N. The dark side of fibroblasts: Cancer-associated fibroblasts as mediators of immunosuppression in the tumor microenvironment. Front Immunol (2019) 10:1835. doi: 10.3389/fimmu.2019.01835

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Gamradt P, De La Fouchardière C, Hennino A. Stromal protein-mediated immune regulation in digestive cancers. Cancers (2021) 13(1):146. doi: 10.3390/cancers13010146

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Miyai Y, Sugiyama D, Hase T, Asai N, Taki T, Nishida K, et al. Meflin-positive cancer-associated fibroblasts enhance tumor response to immune checkpoint blockade. Life Sci alliance. (2022) 5(6):e202101230. doi: 10.26508/lsa.202101230

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Bienkowska KJ, Hanley CJ, Thomas GJ. Cancer-associated fibroblasts in oral cancer: A current perspective on function and potential for therapeutic targeting. Front Oral Health (2021) 2:686337. doi: 10.3389/froh.2021.686337

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Goldman M, Craft B, Hastie M, Repečka K, McDade F, Kamath A, et al. The UCSC xena platform for public and private cancer genomics data visualization and interpretation. bioRxiv (2019), 326470. doi: 10.1101/326470

CrossRef Full Text | Google Scholar

17. Barretina J, Caponigro G, Stransky N, Venkatesan K, Margolin AA, Kim S, et al. The cancer cell line encyclopedia enables predictive modelling of anticancer drug sensitivity. Nature (2012) 483(7391):603–7. doi: 10.1038/nature11003

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Du P, Zhang X, Huang CC, Jafari N, Kibbe WA, Hou L, et al. Comparison of beta-value and m-value methods for quantifying methylation levels by microarray analysis. BMC Bioinf (2010) 11:587. doi: 10.1186/1471-2105-11-587

CrossRef Full Text | Google Scholar

19. Sturm G, Finotello F, List M. Immunedeconv: An r package for unified access to computational methods for estimating immune cell fractions from bulk RNA-sequencing data. Methods Mol Biol (Clifton NJ). (2020) 2120:223–32. doi: 10.1007/978-1-0716-0327-7_16

CrossRef Full Text | Google Scholar

20. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an r package for comparing biological themes among gene clusters. Omics: J Integr Biol (2012) 16(5):284–7. doi: 10.1089/omi.2011.0118

CrossRef Full Text | Google Scholar

21. Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, et al. ClueGO: A cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinf (Oxford England). (2009) 25(8):1091–3. doi: 10.1093/bioinformatics/btp101

CrossRef Full Text | Google Scholar

22. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med (2018) 24(10):1550–8. doi: 10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Zhang Y, Weinberg RA. Epithelial-to-mesenchymal transition in cancer: complexity and opportunities. Front Med (2018) 12(4):361–73. doi: 10.1007/s11684-018-0656-6

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Thorsson V, Gibbs DL, Brown SD, Wolf D, Bortone DS, Ou Yang TH, et al. The immune landscape of cancer. Immunity (2018) 48(4):812–30.e14. doi: 10.1016/j.immuni.2018.03.023

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Xu L, Deng C, Pang B, Zhang X, Liu W, Liao G, et al. TIP: A web server for resolving tumor immunophenotype profiling. Cancer Res (2018) 78(23):6575–80. doi: 10.1158/0008-5472.can-18-0689

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Geeleher P, Cox N, Huang RS. pRRophetic: an r package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PloS One (2014) 9(9):e107468. doi: 10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Louis DN, Perry A, Wesseling P, Brat DJ, Cree IA, Figarella-Branger D, et al. The 2021 WHO classification of tumors of the central nervous system: a summary. Neuro-oncology (2021) 23(8):1231–51. doi: 10.1093/neuonc/noab106

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Kang EM, Yin AA, He YL, Chen WJ, Etcheverry A, Aubry M, et al. A five-CpG signature of microRNA methylation in non-G-CIMP glioblastoma. CNS Neurosci Ther (2019) 25(9):937–50. doi: 10.1111/cns.13133

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Guo W, Ma S, Zhang Y, Liu H, Li Y, Xu JT, et al. Genome-wide methylomic analyses identify prognostic epigenetic signature in lower grade glioma. J Cell Mol Med (2022) 26(2):449–61. doi: 10.1111/jcmm.17101

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Lemée JM, Clavreul A, Menei P. Intratumoral heterogeneity in glioblastoma: don’t forget the peritumoral brain zone. Neuro-oncology (2015) 17(10):1322–32. doi: 10.1093/neuonc/nov119

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Niyazi M, Brada M, Chalmers AJ, Combs SE, Erridge SC, Fiorentino A, et al. ESTRO-ACROP guideline “target delineation of glioblastomas”. Radiotherapy oncology: J Eur Soc Ther Radiol Oncol (2016) 118(1):35–42. doi: 10.1016/j.radonc.2015.12.003

CrossRef Full Text | Google Scholar

32. Shen Y, Chi H, Xu K, Li Y, Yin X, Chen S, et al. A novel classification model for lower-grade glioma patients based on pyroptosis-related genes. Brain Sci (2022) 12(6):700. doi: 10.3390/brainsci12060700

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Schaafsma E, Jiang C, Nguyen T, Zhu K, Cheng C. Microglia-based gene expression signature highly associated with prognosis in low-grade glioma. Cancers (2022) 14(19):4802. doi: 10.3390/cancers14194802

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Kalluri R. The biology and function of fibroblasts in cancer. Nat Rev Cancer. (2016) 16(9):582–98. doi: 10.1038/nrc.2016.73

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Zhang M, Wang D, Su L, Ma J, Wang S, Cui M, et al. Activity of Wnt/PCP regulation pathway classifies patients of low-grade glioma into molecularly distinct subgroups with prognostic difference. Front Oncol (2021) 11:726034. doi: 10.3389/fonc.2021.726034

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Hashemi M, Abbasiazam A, Oraee-Yazdani S, Lenzer J. Response of human glioblastoma cells to hyperthermia: Cellular apoptosis and molecular events. Tissue Cell (2022) 75:101751. doi: 10.1016/j.tice.2022.101751

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Jones PA. Functions of DNA methylation: islands, start sites, gene bodies and beyond. Nat Rev Genet (2012) 13(7):484–92. doi: 10.1038/nrg3230

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Filippou PS, Karagiannis GS, Constantinidou A. Midkine (MDK) growth factor: a key player in cancer progression and a promising therapeutic target. Oncogene (2020) 39(10):2040–54. doi: 10.1038/s41388-019-1124-8

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Gonzalez-Molina J, Gramolelli S, Liao Z, Carlson JW, Ojala PM, Lehti K. MMP14 in sarcoma: A regulator of tumor microenvironment communication in connective tissues. Cells (2019) 8(9):991. doi: 10.3390/cells8090991

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Dou D, Ren X, Han M, Xu X, Ge X, Gu Y, et al. Cancer-associated fibroblasts-derived exosomes suppress immune cell function in breast cancer via the miR-92/PD-L1 pathway. Front Immunol (2020) 11:2026. doi: 10.3389/fimmu.2020.02026

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Zoppi N, Chiarelli N, Cinquina V, Ritelli M, Colombi M. GLUT10 deficiency leads to oxidative stress and non-canonical αvβ3 integrin-mediated TGFβ signalling associated with extracellular matrix disarray in arterial tortuosity syndrome skin fibroblasts. Hum Mol Genet (2015) 24(23):6769–87. doi: 10.1093/hmg/ddv382

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Kriplani N, Duncan RR, Leslie NR. SWAP70 undergoes dynamic conformational regulation at the leading edge of migrating cells. FEBS letters. (2019) 593(4):395–405. doi: 10.1002/1873-3468.13326

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Martija AA, Pusch S. The multifunctional role of EMP3 in the regulation of membrane receptors associated with IDH-Wild-Type glioblastoma. Int J Mol Sci (2021) 22(10):5261. doi: 10.3390/ijms22105261

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Huang L, Wang Z, Chang Y, Wang K, Kang X, Huang R, et al. EFEMP2 indicates assembly of M0 macrophage and more malignant phenotypes of glioma. Aging (2020) 12(9):8397–412. doi: 10.18632/aging.103147

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Grabowski MM, Sankey EW, Ryan KJ, Chongsathidkiet P, Lorrey SJ, Wilkinson DS, et al. Immune suppression in gliomas. J neuro-oncology. (2021) 151(1):3–12. doi: 10.1007/s11060-020-03483-y

CrossRef Full Text | Google Scholar

46. Guadagno E, Presta I, Maisano D, Donato A, Pirrone CK, Cardillo G, et al. Role of macrophages in brain tumor growth and progression. Int J Mol Sci (2018) 19(4):1005. doi: 10.3390/ijms19041005

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Picarda E, Ohaegbulam KC, Zang X. Molecular pathways: Targeting B7-H3 (CD276) for human cancer immunotherapy. Clin Cancer research: an Off J Am Assoc Cancer Res (2016) 22(14):3425–31. doi: 10.1158/1078-0432.ccr-15-2428

CrossRef Full Text | Google Scholar

48. Porichis F, Hart MG, Massa A, Everett HL, Morou A, Richard J, et al. Immune checkpoint blockade restores HIV-specific CD4 T cell help for NK cells. J Immunol (Baltimore Md: 1950) (2018) 201(3):971–81. doi: 10.4049/jimmunol.1701551

CrossRef Full Text | Google Scholar

49. Marin-Acevedo JA, Dholaria B, Soyano AE, Knutson KL, Chumsri S, Lou Y. Next generation of immune checkpoint therapy in cancer: new developments and challenges. J Hematol Oncol (2018) 11(1):39. doi: 10.1186/s13045-018-0582-8

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Tang T, Cheng X, Truong B, Sun L, Yang X, Wang H. Molecular basis and therapeutic implications of CD40/CD40L immune checkpoint. Pharmacol Ther (2021) 219:107709. doi: 10.1016/j.pharmthera.2020.107709

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Chen Y, Liu P, Sun P, Jiang J, Zhu Y, Dong T, et al. Oncogenic MSH6-CXCR4-TGFB1 feedback loop: A novel therapeutic target of photothermal therapy in glioblastoma multiforme. Theranostics (2019) 9(5):1453–73. doi: 10.7150/thno.29987

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Sherbenou DW, Mark TM, Forsberg P. Monoclonal antibodies in multiple myeloma: A new wave of the future. Clin lymphoma myeloma leukemia. (2017) 17(9):545–54. doi: 10.1016/j.clml.2017.06.030

CrossRef Full Text | Google Scholar

53. Liu B, Zhou J, Wang C, Chi Y, Wei Q, Fu Z, et al. LncRNA SOX2OT promotes temozolomide resistance by elevating SOX2 expression via ALKBH5-mediated epigenetic regulation in glioblastoma. Cell Death disease. (2020) 11(5):384. doi: 10.1038/s41419-020-2540-y

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Morelli MB, Amantini C, Nabissi M, Cardinali C, Santoni M, Bernardini G, et al. Axitinib induces senescence-associated cell death and necrosis in glioma cell lines: The proteasome inhibitor, bortezomib, potentiates axitinib-induced cytotoxicity in a p21(Waf/Cip1) dependent manner. Oncotarget (2017) 8(2):3380–95. doi: 10.18632/oncotarget.13769

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Shi F, Guo H, Zhang R, Liu H, Wu L, Wu Q, et al. The PI3K inhibitor GDC-0941 enhances radiosensitization and reduces chemoresistance to temozolomide in GBM cell lines. Neuroscience (2017) 346:298–308. doi: 10.1016/j.neuroscience.2017.01.032

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Jin H, Luo C. Bleomycin inhibits proliferation and promotes apoptosis of brain glioma cells via TGF-β/Smad signaling pathway. J BUON: Off J Balkan Union Oncol (2020) 25(2):1076–83.

Google Scholar

57. Chan TA, Yarchoan M, Jaffee E, Swanton C, Quezada SA, Stenzinger A, et al. Development of tumor mutation burden as an immunotherapy biomarker: Utility for the oncology clinic. Ann oncology: Off J Eur Soc Med Oncol (2019) 30(1):44–56. doi: 10.1093/annonc/mdy495

CrossRef Full Text | Google Scholar

58. Jiao Y, Killela PJ, Reitman ZJ, Rasheed AB, Heaphy CM, de Wilde RF, et al. CIC, FUBP1 and IDH1 mutations refine the classification of malignant gliomas. Oncotarget (2012) 3(7):709–22. doi: 10.18632/oncotarget.588

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Yan H, Parsons DW, Jin G, McLendon R, Rasheed BA, Yuan W, et al. IDH1 and IDH2 mutations in gliomas. New Engl J Med (2009) 360(8):765–73. doi: 10.1056/NEJMoa0808710

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Teter Z, Śliwczyński A, Brzozowska M, Świerkowski M, Jacyna A, Pinkas J, et al. The assessment of overall survival (OS) after adjuvant chemotherapy for patients with malignant endometrial cancer in Poland. Ginekologia polska. (2017) 88(6):296–301. doi: 10.5603/GP.a2017.0056

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Bouffet E, Larouche V, Campbell BB, Merico D, de Borja R, Aronson M, et al. Immune checkpoint inhibition for hypermutant glioblastoma multiforme resulting from germline biallelic mismatch repair deficiency. J Clin oncology: Off J Am Soc Clin Oncol (2016) 34(19):2206–11. doi: 10.1200/jco.2016.66.6552

CrossRef Full Text | Google Scholar

62. Touat M, Li YY, Boynton AN, Spurr LF, Iorgulescu JB, Bohrson CL, et al. Mechanisms and therapeutic implications of hypermutation in gliomas. Nature (2020) 580(7804):517–23. doi: 10.1038/s41586-020-2209-9

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Joyce JA, Fearon DT. T Cell exclusion, immune privilege, and the tumor microenvironment. Sci (New York NY). (2015) 348(6230):74–80. doi: 10.1126/science.aaa6204

CrossRef Full Text | Google Scholar