Prognosis Analysis and Validation of m6A Signature and Tumor Immune Microenvironment in Glioma

Glioma is one of the most typical intracranial tumors, comprising about 80% of all brain malignancies. Several key molecular signatures have emerged as prognostic biomarkers, which indicate room for improvement in the current approach to glioma classification. In order to construct a more veracious prediction model and identify the potential prognosis-biomarker, we explore the differential expressed m⁶A RNA methylation regulators in 665 gliomas from TCGA-GBM and TCGA-LGG. Consensus clustering was applied to the m6A RNA methylation regulators, and two glioma subgroups were identified with a poorer prognosis and a higher grade of WHO classification in cluster 1. The further chi-squared test indicated that the immune infiltration was significantly enriched in cluster 1, indicating a close relation between m⁶A regulators and immune infiltration. In order to explore the potential biomarkers, the weighted gene co-expression network analysis (WGCNA), along with Least absolute shrinkage and selection operator (LASSO), between high/low immune infiltration and m⁶A cluster 1/2 groups were utilized for the hub genes, and four genes (TAGLN2, PDPN, TIMP1, EMP3) were identified as prognostic biomarkers. Besides, a prognostic model was constructed based on the four genes with a good prediction and applicability for the overall survival (OS) of glioma patients (the area under the curve of ROC achieved 0.80 (0.76–0.83) and 0.72 (0.68–0.76) in TCGA and Chinese Glioma Genome Atlas (CGGA), respectively). Moreover, we also found PDPN and TIMP1 were highly expressed in high-grade glioma from The Human Protein Atlas database and both of them were correlated with m6A and immune cell marker in glioma tissue samples. In conclusion, we construct a novel prognostic model which provides new insights into glioma prognosis. The PDPN and TIMP1 may serve as potential biomarkers for prognosis of glioma.

Introduction

Glioma is a common primary tumor in the central nervous system (CNS), accounting for about 80% of brain malignancies (1, 2). The lower-grade gliomas (LGGs) has a relatively favorable prognosis, consisting of the diffuse low-grade and intermediate-grade gliomas (World Health Organization [WHO] grades II and III), whereas glioblastoma (GBM) are generally high-grade gliomas (grade IV) (3, 4). Despite recent medical advances, patients with high-grade GBM are still associated with poor prognosis. Thus, identifying the difference in various gliomas may assist oncologists in finding the prognostic biomarkers and potential targets for glioma patients.

N6-Methyladenosine (m⁶A) is the most popular internal mRNA modification in diverse cell types and consists of the m6A methyltransferases, reverted by the demethylases and identified by m⁶A binding proteins (5–10). Generally, m⁶A modification has various regulatory functions in tumorigenesis, progression and immunity modulation (11–15). Meanwhile, tumor immune microenvironment also participates in tumor initiation and progression and influences the clinical outcomes of patients (16–18). Immune classification of cancers is crucial in therapeutic strategy establishing and prognosis assessment of patients with tumors (19, 20).

Several studies have revealed the correlation between tumor microenvironment (TME) infiltrating immune cells and m⁶A modification. In the gastric tumors, m6A modification patterns could predict the stages of tumor inflammation, TME stromal activity, genetic variation and patient prognosis. Lower m6A score indicated an inflamed TME phenotype and enhanced response to anti-PD-1/L1 immunotherapy (21). The high expression of WTAP, a m6A methyltransferase, was also associated with RNA methylation and its low expression was related to a high T cell-related immune response in gastric cancer (22). Additionally, m⁶A was reduced in the high immunity subtype of lung adenocarcinoma, indicating that m⁶A may mediate immune signatures and help to provide potential strategies (23). However, the potential roles of m6A modification in immune infiltration remain obscure, especially in glioma. Therefore, identification of immune infiltration characterizations mediated by multiple m6A regulators might be helpful for the survival prognosis of patients with glioma.

In this study, in order to investigate the novel prediction model and potential biomarkers for glioma, WGCNA and LASSO were applied to identify candidate genes that might take part in both m6A and immune infiltration in glioma based on TCGA database. Differentially expressed genes (DEGs) were identified, along with their prognostic values, and further validated by external datasets and tissue microarray. Besides, the constructed prediction model revealed a high efficacy for prognosis prediction. The potential predictive biomarkers were also identified to assist oncologists in clinic treatment.

Methods and Materials

Datasets Acquisition From TCGA Datasets

The Cancer Genome Altas (TCGA) GBMLGG datasets (n = 665) were downloaded from the University of California Santa Cruz (UCSC) Xena browser (https://xenabrowser.net/datapages/). The gene expression data were presented as FPKM values derived from TCGA level 3 data. Batch effects were removed before analyzing (24). Clinical data of TCGA datasets were downloaded from the UCSC Xena browser, including clinical information (age, gender), tumor information (subtypes) and survival information (overall survival) for patients with gliomas (Table 1). The RNA-seq transcriptome data and corresponding clinicopathological information of 420 LGG patients and 237 GBM patients were obtained from CGGA (www.cgga.org.cn) as a validation set. The RNA-seq transcriptome data were transformed as FPKM values. GSE16011 (25) expression data was downloaded from GEO database. Robust multi-array average (RMA) normalized files were used in this study. The probe was converted into gene symbol by median gene expression. The microarray data were estimated as log2(x+1) normalized expression value.

TABLE 1

Table 1. Summary table of the TCGA clinical information.

Selection of m⁶A RNA Methylation Regulators

We used 12 m⁶A RNA methylation regulators from published literature. Then, the expression of these m6A RNA methylation regulators in gliomas were systematically compared with different clinical outcomes using Gliovis (http://gliovis.bioinfo.cnio.es/) (26).

Unsupervised Analysis With ConsensusClusterPlus

In order to investigate the function of m⁶A RNA methylation regulators in glioma, we divided patients with glioma into different groups with “ConsensusClusterPlus” (50 iterations, resample rate of 80%). The principal component analysis was then performed with the R package “PCA” for R v3.5.1 to study the gene expression patterns in different glioma clusters. In order to determine the optimal K, Average Silhouette method and Gap Statistic method were applied, the results showed that the two groups were the best grouping number (Supplementary Figure 1). Wilcoxon signed rank test was used to compare the tumor mutation burden of cluster 1 and cluster 2.

Function Analysis of m⁶A Cluster Subgroups and Immune Infiltration Analysis Based on Single-Sample Gene Set Enrichment Analysis (ssGSEA)

Gene Set Variation Analysis (GSVA) was performed with the R package “gsva” to evaluate pathway enrichment for different clusters. To investigate the immune infiltration landscape of glioma, ssGSEA was performed to assess the level of immune infiltration (recorded as ssGSEA score) in a sample according to the expression levels of immune cell-specific marker genes with R package “gsva.” Most immune cell types related marker genes were obtained from the article published by Bindea et al. (27).

Cox Regression Analysis

We assessed the impact of immune cell types on clinical survival data and survival time by Cox proportional hazards regression analysis based on the R package “survival” and “forestplot.” Cell types with a high hazard ratio were considered to be risk factors to OS.

Hub Genes Correlated With m⁶A RNA Methylation Clusters and Immune Infiltration Based on Weighted Correlation Network Analysis (WGCNA)

We extracted all the DEGs (according to adj. p-value < 0.01, |logFC| ≥ 2, total = 729) from limma analysis with expression data retrieved from TCGA GBM/LGG datasets to perform Weighted correlation network analysis (WGCNA) using R package “limma.” We applied R package “WGCNA” to find clinical traits-related modules and hub genes among them (28). The adjacency matrix was then transformed into topological overlap matrix (TOM). Genes were divided into different gene modules according to the TOM-based dissimilarity measure. We set soft-thresholding power as 9 (scale free R2 = 0.85), cut height as 0.2, and minimal module size as 30 to identify key modules. Those with gene significance (GS) > 0.5 and module membership (MM) > 0.9 were defined as hub genes.

Validation of Prognostic Values of Hub Genes

To predict the clinical outcomes of glioma patients with the hub genes, we applied LASSO Cox regression algorithm to the 5 hub genes in the TCGA datasets. We selected four genes to build the risk signature based on the minimum criteria, and the coefficients obtained from the LASSO algorithm were used to calculate the risk score for each patient as follows:

$\begin{array}{l}Riskscore={\displaystyle \sum _{i=1}^n}*{\beta }_i\end{array}$

where n was the number of prognostic genes, exp_i the expression value of gene i, and β_i the regression coefficient of gene i in the LASSO algorithm. Using the median risk score as a cutoff value, glioma patients were divided into high- and low-risk score groups. Moreover, the relation between the prognosis signature and OS was investigated based on the external cohort CGGA datasets.

The Kaplan-Meier method was used to assess the differences of overall survival (OS) between low- and high-risk score glioma patients with R package “survival”.

The time-dependent receiver operating characteristic (ROC) curve was used to measure the prognostic performance by comparing the areas under the ROC curves (AUC) using R package “pROC.” 10-fold cross method was applied for ROC validation and AUC value calculation.

All the scripts were uploaded at Github website (https://github.com/mvpsc30/FIO-m6A-immune).

Assessment of Immunohistochemistry Data

The PDPN and TIMP1 immunohistochemistry results were acquired from the Human Protein Atlas (HPA, https://www.proteinatlas.org/) database (29). The EMP3 and TAGLN2 protein levels of selected genes were evaluated through commercially glioma tissue-microarrays and H-scores between Low-grade gliomas and High-grade gliomas.

Real-Time RT-PCR

Total RNA was extracted from tissue samples and cells using TRIzol reagent (Invitrogen) after washing with PBS. cDNA was synthesized from purified RNA using a SuperScript III First-Strand cDNA synthesis system (18080051, Life Technologies) according to the manufacturer's instruction. SYBR Green PCR Master Mix (Applied Biosystems, CA, USA) was used for PCR amplification and a real-time PCR machine (iQ5, Bio-Rad Laboratories) was used to quantify the expression of mRNAs. β-actin was used as endogenous control and the expression levels were quantified using the methods of 2–ΔΔCt.

Primers:

Cell Culture and Transfection

Human glioma cell line U87 and A172 were acquired from the American Type Culture Collection (ATCC) and cultured in DMEM medium (Gibco, Life Technologies, Grand Island, NY) supplemented with 10% fetal bovine serum (Gibco) and 100 U/ml penicillin/streptomycin (Gibco). According to the manufacturer's instructions, the Lipo 2000 transfection reagent was applied for the transfection. The siRNAs against TIMP1 (siRNA ID: s14143, ThermoFiher), PDPN (EHU119431, Sigma) and negative control (SIC001, Sigma) were purchased.

Western Blotting

Western blot (WB) assays was performed as previously described (30). Briefly, we prepared cell extracts for Western blotting in RIPA buffer. Then, lysates were separated by SDS-PAGE and were transferred to PVDF membranes (Millipore, Billerica, MA). Primary antibodies PDPN (Abcam, ab236529,1:1000), TIMP1 (Abcam, ab109125,1:1000), EMP3 (Santa cruz, sc-81797, 1:100), TAGLN2 (Proteintech, 10234-2-AP, 1:200), and GAPDH (Abcam, ab181602, 1:10000) were used along with HRP-labeled secondary antibody (1:10000, Sigma) in Western blot. The immune complex was detected by chemiluminescence (GE Healthcare, Wauwatosa, WI).

Cell Viability and Cell Death Measurement

Cell viability was measured using the CellTiter-Glo^®luminescent cell viability assay (Promega) based on the manufacturer's instructions. For phosphatidylserine exposure, cells were stained with annexin V-PE as instructed by the manufacturer (BD Biosciences, San Jose, CA), and assayed by flow cytometry (CyAn ADP, Beckman Coulter, Brea, CA, USA).

Statistical Analysis

Experimental results were analyzed with a Student's t-test and graphed using Graphpad Prism application (GraphPad Software, Inc., La Jolla, CA). Data are expressed as mean ± SD. A p < 0.05 was considered with statistical significance. The correlation between the expression profiles of TIMP1 and PDPN with immune and macrophage marker was analyzed using Spearman's rank test.

Results

Consensus Clustering of m⁶A RNA Methylation Regulators Identified Two Clusters of Gliomas With Distinct Immune Infiltration

A flowchart of this study is shown in Supplementary Figure 1. Based on biological functions of each m⁶A RNA methylation regulator in clinical prognosis, we performed consensus clustering based on gene expression of 12 key m⁶A RNA regulators in TCGA datasets. Due to the expression analogy of m⁶A regulators, the clustering analysis would classify the samples into different clusters. After evaluating the relative change in the area under the cumulative distribution function (CDF) curve and consensus heatmap, we selected a three-cluster solution (K = 2), which has no obvious increase in the area under the CDF curve (Supplementary Figures 2A–D). To further determine the optimal K, two methods (Average Silhouette method and Gap Statistic method) were applied. Based on these methods, two subgroups clustered by k = 2, namely, cluster 1 and cluster 2 subgroups were found (Supplementary Figures 2E,F). Most parts of m6A RNA methylation regulators' expressions showed clear distinction and significant difference in two cluster subgroups (Figures 1A,B). In order to better understand the interaction among the 12 m⁶A regulators, we assessed the interaction and correlation among these regulators (Supplementary Figure 3).

FIGURE 1

Figure 1. Identification of consensus clusters by m⁶A RNA methylation regulators overall survival of gliomas in the cluster 1/2 subgroups. (A) Violin plot of the two clusters (cluster1/2) defined by the m6A RNA methylation regulators consensus expression. (B) Principal component analysis of the total RNA expression profiles in the TCGA-GBM/LGG datasets. Gliomas in the cluster1 subgroup are marked with red. (C) Kaplan–Meier overall survival (OS) curves for 665 TCGA glioma patients of different cluster. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0. 0001.

The Kaplan-Meier survival analysis revealed a significant shorter OS in cluster 1 subgroup than the cluster 2 subgroup (Figure 1C). Moreover, we analyzed the DEGs between cluster1 and cluster2, and annotated their function Gene Set Variation Analysis (GSVA) for biological processes. The results indicated that DEGs are enriched in immune-related biological processes, including IL2/STAT5, IL6/JAK/STAT3, and Interferon-γ response signaling (Table 2) and the two categories identified by consensus clustering are correlated with immune infiltration of glioma.

TABLE 2

Table 2. Differences in pathway activities scored per sample by GSVA between cluster 1 and cluster 2, cluster 2 vs. cluster 1.

Immune Landscape Was Significantly Associated With m⁶A RNA Methylation Regulators

To explore the roles of immune cells in the malignant progression of gliomas, the RNA-seq data of 665 patients with gliomas from TCGA-GBM/LGG datasets were analyzed to evaluate the immune landscape. The high and low immune infiltration were defined by Euclidean distance and the ssGSEA scores of immune cells. The results indicated that B cells, Tcm cells, and T helper cells were enriched in high immune infiltration glioma. Relatively, gliomas with low infiltration were characterized for macrophages, eosinophils, neutrophils, and aDC cells (Figure 2A).

FIGURE 2

Figure 2. Immune landscape of glioma. (A) Heatmap of ssGSEA scores of TCGA-GBM/LGG and Table of cluster and immune infiltration subgroups (Chi-square test: X-squared = 116.63, p < 2.2e−16). (B) Kaplan–Meier overall survival (OS) curves for 665 TCGA glioma patients of different immune infiltration subgroups. (C) Forest plot for immune cells. The hazard ratios (HR), 95% confidence intervals (CI) calculated by univariate Cox regression are shown. (D) Mine plot of relationships between 12 m⁶A methylation regulators and 24 immune cells (*P < 0.05, **P < 0.01).

In order to analyze the relationship between m⁶A cluster group and immune infiltration, Chi-squared test was carried out (p < 2.2 × 10−16, Figure 2A). Moreover, we compared the immune infiltration score between cluster 1 and cluster 2, indicating that the proportion of most immune cells types was significantly different between clusters 1 and 2 (Supplementary Figure 4). Then Kaplan-Meier survival curve analysis was performed to explore the roles of immune cell infiltration on the prognosis of patients with glioma. The results revealed that patients with low immune infiltration had worse OS compared with patients with high immune infiltration (Figure 2B). We also applied a univariate Cox regression analysis on the immune cells of TCGA datasets, and found that 23/24 cell types were significantly correlated with OS (P < 0.05). Among these 23 immune cells, aDC, DC, iDC, cytotoxic cells, Eosinophils, Macrophages, Neutrophils, NK.CD56dim cells, NK cells, T cells, Th17 cells, and Th2 cells are risky immune cells with HR > 1, while CD8 T cells, B cells, Mast cells, NK.CD56bright cells, pDC, Tem, Tcm, T helper cells, TFH, Tgd, and Th1 cells were protective immune cells with HR < 1 (Figure 2C).

To further determine the relationship between m⁶A RNA methylation regulators and immune cell infiltration, we assessed the relationships between the expressions of m⁶A RNA methylation regulators and immune cells infiltration subgroups. The results indicated that high immune infiltration was strongly related to higher expressions of FTO, MELLT14, METTL3, YTHDC1, YTHDC2, and ZC3H13. Correspondingly, low immune infiltration with higher expressions of ALKBH5, HNRNPC, WTAP, YTHDF1, and YTHDF2 (Supplementary Figure 5). Then we calculated the relationships between each m⁶A RNA methylation regulators and immune cells, revealing that FTO, ZC3H13, and YTHDC1 had a significant positive correlation with Tcm cells. Meanwhile, macrophages had a negative relationship with FTO and ZC3H13 (Figure 2D). These data indicated that m⁶A clusters were highly associated with immune infiltration.

WGCNA and Identification of the Key Module

In order to explore the key genes that were mostly associated with m6A and immune cell infiltration subtypes in glioma, we performed WGCNA on the TCGA-GBM/LGG datasets. Glioma sample information such as age, m6A cluster subgroups, immune infiltration subgroups, OS and OS status were retrieved from TCGA-GBM/LGG (Supplementary Figure 6A). Eventually identified 6 modules by setting soft-thresholding power as 9 (scale-free R2 = 0.85) and cut height as 0.2 (Supplementary Figures 6B,C). From the heatmap of module-trait correlations, we evaluated that the black module was the most highly related to clinical traits (Supplementary Figure 6D), especially the immune infiltration and outcomes (correlation coefficient = −0.86 and 0.5, P = 4E-206 and 1E-39; respectively, Supplementary Figures 6E–G). Lastly, we selected 5 hub genes (TAGLN2, PDPN, TIMP1, EMP3, CHI3L1) from the black module by setting module membership (MM) >0.9 and gene significance (GS) >0.5. These genes were closely related to each other Supplementary Figure 6H).

Association of Hub Genes With m⁶A RNA Methylation Regulators and Immune Infiltration

We explored the relationship between the expression levels of five hub genes and m⁶A RNA methylation regulators to elucidate the underlying mechanisms of abnormal up-regulation in glioma. The correlation analysis showed that the expression of many hub genes was significantly correlated with m⁶A RNA methylation regulators (Supplementary Figure 7). Additionally, we found that TAGLN2, PDPN, EMP3, and CHI3L1 were positively associated with WTAP (Supplementary Figure 7), while TIMP1 was negatively correlated with YTHDC1.

Then we utilized the Spearman method to study the potential relationship between the expression of glioma hub genes and infiltration of immune cells. Interestingly, hub genes were all positively associated with Macrophages (Supplementary Figure 8). Conversely, negative relationship was observed between these five genes and the infiltration of B cells, Tcm cells and Tem cells (Supplementary Figure 8). These data indicated that the selected five hub genes were highly correlated with m⁶A RNA regulators and immune infiltration.

Validation of Hub Genes in Datasets

To predict the clinical outcomes of glioma with the hub genes, we applied the LASSO Cox regression algorithm to the five hub genes in the TCGA datasets (Supplementary Figures 9A,B). Four genes were highly associated with clinical features, such as grades, transcriptome subtype and IDH status (Supplementary Figures 10A–C). Moreover, these four genes were used to set up the risk signature based on the minimum criteria. Next, to assess the differences of survival time between low- and high-risk glioma patients, the Kaplan-Meier method was performed. Meanwhile, the log-rank test was also used to determine the statistical significance between groups. The time-dependent ROC curve was employed to measure the prognostic performance by comparing the AUC. Compared with those in the low-risk group, we illustrated that the glioma patients in the high-risk group had shorter OS, (Figures 3A,B, TCGA: HR = 1.07, 95% CI = 1.06–1.08, P < 0.01; CGGA: HR = 1.19, 95% CI = 1.16–2.23, P < 0.01). The time-dependent ROC curves revealed that the AUC for the 4-gene signature achieved 0.80 (0.76–0.83) and 0.72 (0.68–0.76) for the OS in TCGA and CGGA datasets, respectively (Figures 3C,D). Furthermore, the risk score exhibited a higher prognostic accuracy for OS than clinical histology, grade, IDH status and age (Figures 3E,F). These findings suggested an effective performance for predicting OS for glioma patients.

FIGURE 3

Figure 3. Validation of hub genes highly correlated with clinical traits. In each dataset, the risk score distribution, gene expression profiles, and patients' survival status are displayed (A, TCGA; B, CGGA). The black-dotted line represents the median cut-off, dividing patients into high- and low-risk groups. Kaplan-Meier and ROC curves with 95% confidence interval for the 4-gene signature in the four datasets. Patients with high risk scores had poor outcome in terms of overall survival (C, TCGA; D, CGGA). ROC curves comparing prognostic accuracy of risk score with clinical histology, grade, IDH status, and age in internal validation, and external validation cohorts (E, TCGA; F, CGGA).

Validation the Expression and Function of TIMP1 and PDPN

To further validate the expression of four genes in gliomas, we next detected their expressions in The Human Protein Atlas database, and the results revealed the PDPN and TIMP1 were higher expression in high-grade gliomas (Figure 4A). In addition, TAGLN2 and EMP3 were performed in commercially glioma tissue-microarrays. The H-score of both proteins was not statistically significant between low and high-grade gliomas (Supplementary Figure 11). Moreover, in the correlation analysis, we uncovered that TIMP1 and PDPN were positively correlated with marker genes of macrophage (Figure 4B, Table 3 and Supplementary Figure 8). TIMP1 was negatively related with YHDC1, while PDPN was positively related with WTAP (Figure 4C and Supplementary Figure 7). By knockdown the expression of PDPN or TIMP1, the cell proliferation was decreased, and the apoptosis and necrosis were increased in U87 and A172 (Supplementary Figure 12).

FIGURE 4

Figure 4. TIMP1 and PDPN reveal its higher expression in GBM and malignant biological phenotype in vitro. (A) The translation expression level of the TIMP1 and PDPN were positively correlated with disease status as they were upregulated in gliomas samples. (B) The relationship of TIMP1 and PDPN with macrophage marker CD68. (C) The relationship of TIMP1 and PDPN with indicated m6A gene.

TABLE 3

Table 3. Clinical data of patients.

Discussion

As the most aggressive primary brain tumor, glioma is considered as an enigma in neurosurgery (31, 32). Advanced knowledge of its genomic changes has promoted the discovery of prognostic signatures to facilitate the personalized treatment decisions (33–35). However, no previous studies have investigated the efficacy of the combination of m⁶A and immune infiltration. Here, we developed and validated a novel 4-gene prognostic model based on the combination of m⁶A RNA methylation and landscape of immune microenvironment. The developed 4-gene signature was able to identify the glioma patients with different risk levels for prognosis, which may compensate the already known prognostic indicators, such as age, tumor grade or histology. Additionally, we confirmed that PDPN and TIMP1 were higher expressed in high-grade glioma, and the Pearson correlation validated that PDPN and TIMP1 were correlated with marker gene of macrophage and indicated m6A gene.

m⁶A, the most prevalent intra-mRNA modification, is required for post-transcriptional regulation of mRNA in various cell types (11, 12, 36). Previous studies have shown that m6A could be a signature for predicting the prognosis in different type of cancers, such as renal cell carcinoma, hepatocellular carcinoma, bladder cancer and head and neck squamous cell carcinoma (37–40). We found that WTAP and HNRNPC were significantly increased in cluster 1 than cluster 2 (Figure 2). In the GBM, WTAP was found to be overexpressed and regulate migration and invasion in vitro (41). Its high expression was associated with poor postoperative survival (42). In addition, HNRNPC could also control the aggressiveness of GBM cells and be regarded as the potential prognostic biomarker and therapeutic targets of GBM (43).

With the high-speed development of omics, high-throughput tumor databases have been established, including TCGA and CGGA, which provided a solid foundation for analyzing the RNA modification and microenvironments of glioma (3, 44–46). One of the emerging strategies of management is based on the roles of immune cells in the growth and maintenance of tumors (47). According to the recent studies, myeloid-derived suppressor cells (MDSC) and tumor-associated macrophages (TAMs) have been identified as promising targets for anti-cancer treatment (48, 49). Neoantigen-targeting vaccines have also increased tumor-infiltrating T cells and altered the immune milieu of glioblastoma (50). According to the TCGA database, Jia et al. has drawn a list of 44 tumor microenvironment related genes and proved them in an independent GBM cohort as potential biomarkers for GBM (51). However, the outcomes may lead to the discordance generally based on only one factor (51, 52). In our current study, we integrated m6A and immune infiltration in TCGA to build a model to improve the overall prediction of outcome for patients with glioma. Four survival-related genes (TAGLN2, PDPN, TIMP1, and EMP3) were identified and verified by four external datasets. These combination of these four genes provided a more reliable signature, relative to that extracted from a single dataset. Furthermore, PDPN and TIMP1 were confirmed that they were higher expression in high-grade glioma and knockdown their expression decreased the glioma cell proliferation in vitro.

TAGLN2 is considered as a smooth muscle cytoskeletal protein (53). It has been proposed to be associated with growth and migration in bladder cancer (54, 55), esophageal squamous cell carcinoma (56), and gliomas (57). Moreover, it's up-regulation is associated with tumorigenesis and tumor progression (54, 58). Silence of TAGLN2 in gliomas cell lines significantly inhibited invasion and tumor growth (57). Increased expression of TAGLN2 was correlated with deteriorative tumor grade, and the function and regulation made it as a candidate prognostic biomarker (57). Jin et al. has also shown TAGLN2 as a potential biomarker of tumor-derived lung-cancer endothelial cells (59). Another study demonstrated that TAGLN2 could be a prospective tumor tissue marker for diagnosis and evaluating lymph node metastasis in bladder cancer patients (60).

EMP3 belongs to the PMP-22/EMP/MP20 family, which is thought to be involved in cell proliferation, cell-cell interactions and function as a tumor suppressor. Alaminos et al. have suggested that EMP3 was associated with poor survival (61). EMP3 overexpression in breast cancer was related to stronger HER-2 expression that may indicate a novel therapeutic target (62). Ma et al. have demonstrated that EMP3-mediated miR-663a inhibits the gallbladder cancer progression via the MAPK/ERK pathway (63). Recently, the bioinformatics analysis also found that EMP3 was one of the validated gene panel independently and was correlated with the GBM survival (64, 65). Another bioinformatics analysis though significant analysis of microarray (SAM) identified that EMP3 could be used to estimate glioma patient prognosis (66).

TIMP metallopeptidase inhibitor 1 (TIMP1) is a glycoprotein which antagonized mostly known MMPs. The encoded protein can promote cell proliferation in many cell types and may also have an anti-apoptotic function. A high serum level was found as a poor prognostic indicator in GBMs (67). TIMP1 has been suggested to interact with P75NTR in metastatic carcinoma and glioma cells (68), and silence of TIMP1 or inhibition of NF-kappa B activity led to slower tumor growth in vivo (69). Several studies have shown that TIMP1 was an important part of prognosis model and could be a biomarker for diagnosis (70–72). Furthermore, Jackson et al. have reviewed that TIMP1 overexpression is consistently correlated with cancer progression or poor prognosis (73).

Podoplanin (PDPN) is a transmembrane receptor that participates in various physiological and pathological processes, such as cell motility, tumor metastasis and angiogenesis (74–76). It regulated mammary stem cell function that reduced mammary tumor formation in breast cancer and could be a new regulator of Wnt/β-catenin signaling (77). PDPN receptor are upregulated in cancer cells, immune cells, synoviocytes, and fibroblasts that increase tissue inflammation and invasion to promote both arthritis and cancer (78). PDPN-expressing macrophages (PoEMs) stimulated local matrix remodeling, and macrophage-specific PDPN knockout restrained lymphangiogenesis and reduced lymphatic cancer spread (79). PDPN-positive cancer-associated fibroblasts (CAFs) contributed to an essential role in primary resistance to epidermal growth factor receptor tyrosine kinase inhibitors (EGFR-TKI) (80). Moreover, PDPN has been considered as a novel biomarker, chemotherapeutic target and a target for CAR T-cell therapy that may be a potential adoptive immunotherapy to treat GBM (81, 82).

Our finding provides a novel insight into the relationship between m6A and immune infiltration, and we laid a solid foundation for four genes that could be a new prognosis indicator for gliomas patients. In addition, we also developed a user-friendly R shiny web app (http://www.houshixu.cn:3838/sample-apps/fio/) for easier usage. Remarkably, several limitations should be noted. In this study, prognostic factors were found by combining m6A and immune microenvironment. However, we do not have large quantities of samples to verify them and the clustering of glioma by m6A regulators is probably skewed by the grade of glioma. Whether TAGLN2 and EMP3 modulate cell proliferation were unclear. Moreover, the signature requires further validation in prospective studies and multicenter clinical trials.

Conclusions

We construct a novel prognostic model that provides new insights into glioma prognosis. The PDPN and TIMP1 may serve as potential biomarkers for prognosis of glioma.

Data Availability Statement

Publicly available datasets were analyzed in this study, these can be found in The Cancer Genome Atlas via the University of California Santa Cruz (UCSC) Xena browser (https://xenabrowser.net/datapages/).

Ethics Statement

The studies involving human participants were reviewed and approved by Shanghai General Hospital, Shanghai Jiao Tong University School of Medicine, Shanghai 200080, China. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

This work was carried out in collaboration with all authors. ML, MW, and TM contributed to the conception and design of the study. SL and HX carried out the experiments. SL, HX, and AZ contributed to all figures and tables. SL, TM, MW, and YN revised the manuscript. SL, HX, AZ, and YX contributed to data collection and analysis. All authors have read and approved the final manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by National Natural Science Foundation of China (No. 81701359), Natural science foundation of Shanghai (18ZR1430400), and the Cross Research Fund of Medicine and Engineering of Shanghai Jiaotong University (Nos. YG2016QN32 and YG2019QNA67).

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.

Supplementary Material

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

Supplementary Figure 1. An overall flowchart of this work.

Supplementary Figure 2. Identification of consensus clusters by m6A RNA methylation regulators. Consensus clustering matrix for k = 2 to K = 5 (A–D). (E) Consensus clustering cumulative distribution function (CDF) for k = 2–9. (F) Relative change in area under CDF curve for k = 2–5. Silhouette analysis (G) and Gap analysis (H) showed that 2 clusters were appropriate classification for the data.

Supplementary Figure 3. Spearman correlation analysis of the 12 m⁶A modification regulators.

Supplementary Figure 4. The immune infiltration score between cluster 1 and cluster 2.

Supplementary Figure 5. Relationship between m⁶A RNA methylation regulators and immune infiltration (A–J) violin plot for 12 m⁶A regulators and immune infiltration subgroups.

Supplementary Figure 6. Identification of key modules correlated with clinical traits in the TCGA-GBM/LGG datasets through WGCNA. (A) Clustering dendrograms of genes. Color intensity varies positively with age, m⁶A cluster subgroups, immune infiltration subgroups, overall survival and overall survival status. Analysis of the scale-free fit index (B) and the mean connectivity (C) for various soft-thresholding powers. (D) Heatmap of the correlation between module eigengenes and clinical traits of diffuse gliomas. Each cell contains the correlation coefficient and P-value. (E) Dendrogram of all DEGs clustered based on a dissimilarity measure (1- TOM). (F) Clustering of module eigengenes. The red line indicates cut height (0.2). (G) Scatter plot of module eigengenes in the black module. (H) Hub genes show strong associations with each other. Red and blue colors indicate positive and negative coefficients and labels from−1 to 1 indicate correlation strength.

Supplementary Figure 7. Association of hub genes' expression with 12 m6A RNA methylation regulators in gliomas. (A) TAGLN2 (B) PDPN (C) TIMP1 (D) EMP3 (E) CHI3L1.

Supplementary Figure 8. Association of hub genes' expression with immune infiltration cells in gliomas. (A) TAGLN2 (B) PDPN (C) TIMP1 (D) EMP3 (E) CHI3L1.

Supplementary Figure 9. (A,B) The process of building the risk scores containing 4 hub genes and the coefficients calculated by least absolute shrinkage and selection operator (LASSO) Cox regression algorithm are shown.

Supplementary Figure 10. Expression of 4 hub genes in gliomas with different clinicopathological features, from right to left, TIMP1, TAGLN2, PDPN, and EMP3, respectively. (A) The expression levels of 4 hub genes in gliomas with different WHO grades. (B) The expression levels of 4 hub genes in gliomas with different transcriptome subtypes. (C) The expression levels of 4 hub genes in gliomas with different IDH status.

Supplementary Figure 11. H-score of EMP3 (A) and TAGLN2 (B) of glioma tissue-microarrays.

Supplementary Figure 12. (A,D) Western blot analysis validated the knockdown of TIPM1 or PDPN in U87 and A172 cells. (B,E) Cell proliferation was determined by ATP assay. (C,F) Flow cytometric analysis of Annexin V/PI staining in U87 and A172 cells after transfection with 50 nmol/L siTIMP1 or siPDPN and siCtrl for 72 h.

References

1. Ostrom QT, Gittleman H, Truitt G, Boscia A, Kruchko C, Barnholtz-Sloan JS. CBTRUS statistical report: primary brain and other central nervous system tumors diagnosed in the United States in 2011–2015. Neuro Oncol. (2018) 20:iv1–86. doi: 10.1093/neuonc/noy131

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Gittleman H, Sloan AE, Barnholtz-Sloan JS. An independently validated survival nomogram for lower grade glioma. Neuro Oncol. (2019) 22:665–74. doi: 10.1093/noajnl/vdz007

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Brat DJ, Verhaak RG, Aldape KD, Yung WK, Salama SR, Cooper LA, et al. Comprehensive, integrative genomic analysis of diffuse lower-grade gliomas. N Engl J Med. (2015) 372:2481–98. doi: 10.1056/NEJMoa1402121

PubMed Abstract | CrossRef Full Text

4. Perry JR, Laperriere N, O'callaghan CJ, Brandes AA, Menten J, Phillips C, et al. Short-course radiation plus temozolomide in elderly patients with glioblastoma. N Engl J Med. (2017) 376:1027–37. doi: 10.1056/NEJMoa1611977

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Wojtas MN, Pandey RR, Mendel M, Homolka D, Sachidanandam R, Pillai RS. Regulation of m(6)A transcripts by the 3'–>5' RNA helicase YTHDC2 is essential for a successful meiotic program in the mammalian germline. Mol Cell. (2017) 68:374–87.e312. doi: 10.1016/j.molcel.2017.09.021

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Deng X, Su R, Weng H, Huang H, Li Z, Chen J. RNA N(6)-methyladenosine modification in cancers: current status and perspectives. Cell Res. (2018) 28:507–17. doi: 10.1038/s41422-018-0034-6

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Ding C, Zou Q, Ding J, Ling M, Wang W, Li H, et al. Increased N6-methyladenosine causes infertility is associated with FTO expression. J Cell Physiol. (2018) 233:7055–66. doi: 10.1002/jcp.26507

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Tang C, Klukovich R, Peng H, Wang Z, Yu T, Zhang Y, et al. ALKBH5-dependent m6A demethylation controls splicing and stability of long 3'-UTR mRNAs in male germ cells. Proc Natl Acad Sci USA. (2018) 115:E325–e333. doi: 10.1073/pnas.1717794115

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Yang Y, Hsu PJ, Chen YS, Yang YG. Dynamic transcriptomic m(6)A decoration: writers, erasers, readers and functions in RNA metabolism. Cell Res. (2018) 28:616–24. doi: 10.1038/s41422-018-0040-8

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Yao QJ, Sang L, Lin M, Yin X, Dong W, Gong Y, et al. Mettl3-Mettl14 methyltransferase complex regulates the quiescence of adult hematopoietic stem cells. Cell Res. (2018) 28:952–4. doi: 10.1038/s41422-018-0062-2

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Desrosiers R, Friderici K, Rottman F. Identification of methylated nucleosides in messenger RNA from Novikoff hepatoma cells. Proc Natl Acad Sci USA. (1974) 71:3971–5. doi: 10.1073/pnas.71.10.3971

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Wang X, Lu Z, Gomez A, Hon GC, Yue Y, Han D, et al. N6-methyladenosine-dependent regulation of messenger RNA stability. Nature. (2014) 505:117–20. doi: 10.1038/nature12730

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Cui Q, Shi H, Ye P, Li L, Qu Q, Sun G, et al. mA RNA methylation regulates the self-renewal and tumorigenesis of glioblastoma stem cells. Cell Rep. (2017) 18:2622–34. doi: 10.1016/j.celrep.2017.02.059

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Zheng Q, Hou J, Zhou Y, Li Z, Cao X. The RNA helicase DDX46 inhibits innate immunity by entrapping mA-demethylated antiviral transcripts in the nucleus. Nat Immunol. (2017) 18:1094–103. doi: 10.1038/ni.3830

CrossRef Full Text | Google Scholar

15. Visvanathan A, Patil V, Arora A, Hegde AS, Arivazhagan A, Santosh V, et al. Essential role of METTL3-mediated mA modification in glioma stem-like cells maintenance and radioresistance. Oncogene. (2018) 37:522–33. doi: 10.1038/onc.2017.351

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Cooper LA, Gutman DA, Chisolm C, Appin C, Kong J, Rong Y, et al. The tumor microenvironment strongly impacts master transcriptional regulators and gene expression class of glioblastoma. Am J Pathol. (2012) 180:2108–19. doi: 10.1016/j.ajpath.2012.01.040

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Galon J, Pages F, Marincola FM, Thurin M, Trinchieri G, Fox BA, et al. The immune score as a new possible approach for the classification of cancer. J Transl Med. (2012) 10:1. doi: 10.1186/1479-5876-10-1

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Senbabaoglu Y, Gejman RS, Winer AG, Liu M, Van Allen EM, De Velasco G, et al. Tumor immune microenvironment characterization in clear cell renal cell carcinoma identifies prognostic and immunotherapeutically relevant messenger RNA signatures. Genome Biol. (2016) 17:231. doi: 10.1186/s13059-016-1092-z

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Hanahan D, Weinberg RA. The hallmarks of cancer. Cell. (2000) 100:57–70. doi: 10.1016/S0092-8674(00)81683-9

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Hanahan D, Coussens LM. Accessories to the crime: functions of cells recruited to the tumor microenvironment. Cancer Cell. (2012) 21:309–22. doi: 10.1016/j.ccr.2012.02.022

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Zhang B, Wu Q, Li B, Wang D, Wang L, Zhou YL. m(6)A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. Mol Cancer. (2020) 19:53–53. doi: 10.1186/s12943-020-01170-0

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Li H, Su Q, Li B, Lan L, Wang C, Li W, et al. High expression of WTAP leads to poor prognosis of gastric cancer by influencing tumour-associated T lymphocyte infiltration. J Cell Mol Med. (2020) 24:4452–65. doi: 10.1111/jcmm.15104

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Xu F, Chen J-X, Yang X-B, Hong X-B, Li Z-X, Lin L, et al. Analysis of lung adenocarcinoma subtypes based on immune signatures identifies clinical implications for cancer therapy. Mol Therapy Oncolyt. (2020) 17:241–9. doi: 10.1016/j.omto.2020.03.021

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. (2012) 28:882–3. doi: 10.1093/bioinformatics/bts034

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Gravendeel LaM, Kouwenhoven MCM, Gevaert O, De Rooi JJ, Stubbs AP, Duijm JE, et al. Intrinsic gene expression profiles of gliomas are a better predictor of survival than histology. Cancer Res. (2009) 69:9065–72. doi: 10.1158/0008-5472.CAN-09-2307

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Bowman RL, Wang Q, Carro A, Verhaak RGW, Squatrito M. GlioVis data portal for visualization and analysis of brain tumor expression datasets. Neuro Oncol. (2017) 19:139–41. doi: 10.1093/neuonc/now247

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Bindea G, Mlecnik B, Tosolini M, Kirilovsky A, Waldner M, Obenauf AC, et al. Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity. (2013) 39:782–95. doi: 10.1016/j.immuni.2013.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinform. (2008) 9:559. doi: 10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Uhlen M, Fagerberg L, Hallstrom BM, Lindskog C, Oksvold P, Mardinoglu A, et al. Proteomics. Tissue Based Map Human Proteome Sci. (2015) 347:1260419. doi: 10.1126/science.1260419

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Lin SJ, Wu ZR, Cao L, Zhang Y, Leng ZG, Guo YH, et al. Pituitary tumor suppression by combination of cabergoline and chloroquine. J Clin Endocrinol Metab. (2017) 102:3692–703. doi: 10.1210/jc.2017-00627

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Shergalis A, Bankhead A 3rd, Luesakul U, Muangsin N, Neamati N. Current challenges and opportunities in treating glioblastoma. Pharmacol Rev. (2018) 70:412–45. doi: 10.1124/pr.117.014944

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Van Tellingen O, Yetkin-Arik B, De Gooijer MC, Wesseling P, Wurdinger T, De Vries HE. Overcoming the blood-brain tumor barrier for effective glioblastoma treatment. Drug Resist Updat. (2015) 19:1–12. doi: 10.1016/j.drup.2015.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Davis B, Shen Y, Poon CC, Luchman HA, Stechishin OD, Pontifex CS, et al. Comparative genomic and genetic analysis of glioblastoma-derived brain tumor-initiating cells and their parent tumors. Neuro Oncol. (2016) 18:350–60. doi: 10.1093/neuonc/nov143

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Wang J, Cazzato E, Ladewig E, Frattini V, Rosenbloom DI, Zairis S, et al. Clonal evolution of glioblastoma under therapy. Nat Genet. (2016) 48:768–76. doi: 10.1038/ng.3590

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Zhou D, Alver BM, Li S, Hlady RA, Thompson JJ, Schroeder MA, et al. Distinctive epigenomes characterize glioma stem cells and their response to differentiation cues. Genome Biol. (2018) 19:43. doi: 10.1186/s13059-018-1420-6

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Dominissini D, Moshitch-Moshkovitz S, Schwartz S, Salmon-Divon M, Ungar L, Osenberg S, et al. Topology of the human and mouse m6A RNA methylomes revealed by m6A-seq. Nature. (2012) 485:201–6. doi: 10.1038/nature11112

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Chen M, Nie ZY, Wen XH, Gao YH, Cao H, Zhang SF. m6A RNA methylation regulators can contribute to malignant progression and impact the prognosis of bladder cancer. Biosci Rep. (2019) 39:2892. doi: 10.1042/BSR20192892

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Zhao X, Cui L. Development and validation of a m(6)A RNA methylation regulators-based signature for predicting the prognosis of head and neck squamous cell carcinoma. Am J Cancer Res. (2019) 9:2156–69.

PubMed Abstract | Google Scholar

39. Qu N, Qin S, Zhang X, Bo X, Liu Z, Tan C, et al. Multiple m(6)A RNA methylation modulators promote the malignant progression of hepatocellular carcinoma and affect its clinical prognosis. BMC Cancer. (2020) 20:165. doi: 10.1186/s12885-020-6638-5

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Wang J, Zhang C, He W, Gou X. Effect of m(6)A RNA methylation regulators on malignant progression and prognosis in renal clear cell carcinoma. Front Oncol. (2020) 10:3. doi: 10.3389/fonc.2020.00003

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Jin D-I, Lee SW, Han M-E, Kim H-J, Seo S-A, Hur G-Y, et al. Expression and roles of Wilms' tumor 1-associating protein in glioblastoma. Cancer Sci. (2012) 103:2102–9. doi: 10.1111/cas.12022

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Xi Z, Xue Y, Zheng J, Liu X, Ma J, Liu Y. WTAP expression predicts poor prognosis in malignant glioma patients. J Mol Neurosci. (2016) 60:131–6. doi: 10.1007/s12031-016-0788-6

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Park YM, Hwang SJ, Masuda K, Choi K-M, Jeong M-R, Nam D-H, et al. Heterogeneous nuclear ribonucleoprotein C1/C2 controls the metastatic potential of glioblastoma by regulating PDCD4. Mol Cell Biol. (2012) 32:4237–44. doi: 10.1128/MCB.00443-12

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Li J, Lu Y, Akbani R, Ju Z, Roebuck PL, Liu W, et al. TCPA: a resource for cancer functional proteomics data. Nat Methods. (2013) 10:1046–7. doi: 10.1038/nmeth.2650

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Ceccarelli M, Barthel FP, Malta TM, Sabedot TS, Salama SR, Murray BA, et al. Molecular profiling reveals biologically discrete subsets and pathways of progression in diffuse glioma. Cell. (2016) 164:550–63. doi: 10.1016/j.cell.2015.12.028

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Wang Q, Hu B, Hu X, Kim H, Squatrito M, Scarpace L, et al. Tumor evolution of glioma-intrinsic gene expression subtypes associates with immunological changes in the microenvironment. Cancer Cell. (2017) 32:42–56.e46. doi: 10.1016/j.ccell.2017.06.003

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Maximov V, Chen Z, Wei Y, Hope Robinson M, Herting CJ, Shanmugam NS, et al. Tumour-associated macrophages exhibit anti-tumoural properties in Sonic Hedgehog medulloblastoma. Nat Commun. (2019) 10:2410. doi: 10.1038/s41467-019-10458-9

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Kohanbash G, McKaveney K, Sakaki M, Ueda R, Mintz AH, Amankulor N, et al. GM-CSF promotes the immunosuppressive activity of glioma-infiltrating myeloid cells through interleukin-4 receptor-α. Cancer Res. (2013) 73:6413–23. doi: 10.1158/0008-5472.CAN-12-4124

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Pyonteck SM, Akkari L, Schuhmacher AJ, Bowman RL, Sevenich L, Quail DF, et al. CSF-1R inhibition alters macrophage polarization and blocks glioma progression. Nat Med. (2013) 19:1264–72. doi: 10.1038/nm.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Keskin DB, Anandappa AJ, Sun J, Tirosh I, Mathewson ND, Li S, et al. Neoantigen vaccine generates intratumoral T cell responses in phase Ib glioblastoma trial. Nature. (2019) 565:234–9. doi: 10.1038/s41586-018-0792-9

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Jia D, Li S, Li D, Xue H, Yang D, Liu Y. Mining TCGA database for genes of prognostic value in glioblastoma microenvironment. Aging. (2018) 10:592–605. doi: 10.18632/aging.101415

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Chai RC, Wu F, Wang QX, Zhang S, Zhang KN, Liu YQ, et al. m(6)A RNA methylation regulators contribute to malignant progression and have clinical prognostic impact in gliomas. Aging. (2019) 11:1204–25. doi: 10.18632/aging.101829

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Jang SH, Jun CD, Park ZY. Label-free quantitative phosphorylation analysis of human transgelin2 in Jurkat T cells reveals distinct phosphorylation patterns under PKA and PKC activation conditions. Proteome Sci. (2015) 13:14. doi: 10.1186/s12953-015-0070-9

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Zhang HJ, Jiang MJ, Liu QJ, Han ZX, Zhao YQ, Ji SQ. miR-145-5p inhibits the proliferation and migration of bladder cancer cells by targeting TAGLN2. Oncol Lett. (2018) 16:6355–60. doi: 10.3892/ol.2018.9436

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Zhao F, Zhou LH, Ge YZ, Ping WW, Wu X, Xu ZL, et al. MicroRNA-133b suppresses bladder cancer malignancy by targeting TAGLN2-mediated cell cycle. J Cell Physiol. (2019) 234:4910–23. doi: 10.1002/jcp.27288

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Du YY, Zhao LM, Chen L, Sang MX, Li J, Ma M, et al. The tumor-suppressive function of miR-1 by targeting LASP1 and TAGLN2 in esophageal squamous cell carcinoma. J Gastroenterol Hepatol. (2016) 31:384–93. doi: 10.1111/jgh.13180

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Han MZ, Xu R, Xu YY, Zhang X, Ni SL, Huang B, et al. TAGLN2 is a candidate prognostic biomarker promoting tumorigenesis in human gliomas. J Exp Clin Cancer Res. (2017) 36:9. doi: 10.1186/s13046-017-0619-9

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Pei J, Li P, Zhang ZY, Zhang HL, Gao YH, Wang DY, et al. Effect of TAGLN2 in the regulation of meningioma tumorigenesis and development. Eur Rev Med Pharmacol Sci. (2018) 22:307–13. doi: 10.26355/eurrev_201801_14173

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Jin H, Cheng X, Pei Y, Fu Y, Lyu Z, Peng H, et al. Identification and verification of transgelin-2 as a potential biomarker of tumor-derived lung-cancer endothelial cells by comparative proteomics. J Proteom. (2016) 136:77–88. doi: 10.1016/j.jprot.2015.12.012

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Chen C-L, Chung T, Wu C-C, Ng K-F, Yu J-S, Tsai C-H, et al. Comparative tissue proteomics of microdissected specimens reveals novel candidate biomarkers of bladder cancer. Mol Cell Proteom. (2015) 14:2466–78. doi: 10.1074/mcp.M115.051524

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Alaminos M, Davalos V, Ropero S, Setien F, Paz MF, Herranz M, et al. EMP3, a myelin-related gene located in the critical 19q13.3 region, is epigenetically silenced and exhibits features of a candidate tumor suppressor in glioma and neuroblastoma. Cancer Res. (2005) 65:2565–71. doi: 10.1158/0008-5472.CAN-04-4283

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Wang Y-W, Cheng H-L, Ding Y-R, Chou L-H, Chow N-H. EMP1, EMP 2, and EMP3 as novel therapeutic targets in human cancer. Biochim Biophys Acta Rev Cancer. (2017) 1868:199–211. doi: 10.1016/j.bbcan.2017.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Ma Q, Zhang YJ, Liang HB, Zhang F, Liu FT, Chen SL, et al. EMP3, which is regulated by miR-663a, suppresses gallbladder cancer progression via interference with the MAPK/ERK pathway. Cancer Lett. (2018) 430:97–108. doi: 10.1016/j.canlet.2018.05.022

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Shu C, Wang Q, Yan XL, Wang JH. Whole-genome expression microarray combined with machine learning to identify prognostic biomarkers for high-grade glioma. J Mol Neurosci. (2018) 64:491–500. doi: 10.1007/s12031-018-1049-7

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Guo XX, Su J, He XF. A 4-gene panel predicting the survival of patients with glioblastoma. J Cell Biochem. (2019) 120:16037–43. doi: 10.1002/jcb.28883

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Gao Y-F, Zhu T, Mao C-X, Liu Z-X, Wang Z-B, Mao X-Y, et al. (2016). PPIC, EMP3, and CHI3L1 are novel prognostic markers for high grade glioma. Int J Mol Sci. 17:111808. doi: 10.3390/ijms17111808

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Sreekanthreddy P, Srinivasan H, Kumar DM, Nijaguna MB, Sridevi S, Vrinda M, et al. Identification of potential serum biomarkers of glioblastoma: serum osteopontin levels correlate with poor prognosis. Cancer Epidemiol Biomarkers Prev. (2010) 19:1409–22. doi: 10.1158/1055-9965.EPI-09-1077

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Rojiani M, Brem S, Murphy S, Obadia M, Rojiani A. TIMP1 interactions with P75(NTR) in metastatic carcinoma and glioma cells. Cancer Res. (2009) 69. Available online at: https://cancerres.aacrjournals.org/content/69/9_Supplement/1129

Google Scholar

69. Friedmann-Morvinski D, Narasimamurthy R, Xia YF, Myskiw C, Soda Y, Verma IM. Targeting NF-kappa B in glioblastoma: a therapeutic approach. Sci Adv. (2016) 2:1501292. doi: 10.1126/sciadv.1501292

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Resovi A, Bani MR, Porcu L, Anastasia A, Minoli L, Allavena P, et al. Soluble stroma-related biomarkers of pancreatic cancer. EMBO Mol Med. (2018) 10:e8741. doi: 10.15252/emmm.201708741

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Zhong C, Wang G, Xu T, Zhu Z, Guo D, Zheng X, et al. Tissue inhibitor metalloproteinase-1 and clinical outcomes after acute ischemic stroke. Neurology. (2019) 93:e1675–85. doi: 10.1212/WNL.0000000000008564

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Loon EV, Gazut S, Yazdani S, Lerut E, de Loor H, Coemans M, et al. Development and validation of a peripheral blood mRNA assay for the assessment of antibody-mediated kidney allograft rejection: a multicentre, prospective study. EBioMed. (2019) 46:463–72. doi: 10.1016/j.ebiom.2019.07.028

PubMed Abstract | CrossRef Full Text | Google Scholar

73. Jackson HW, Defamie V, Waterhouse P, Khokha R. TIMPs: versatile extracellular regulators in cancer. Nat Rev Cancer. (2017) 17:38–53. doi: 10.1038/nrc.2016.115

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Hwang YS, Zhang X, Park KK, Chung WY. Functional invadopodia formation through stabilization of the PDPN transcript by IMP-3 and cancer-stromal crosstalk for PDPN expression. Carcinogenesis. (2012) 33:2135–46. doi: 10.1093/carcin/bgs258

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Sikorska J, Gawel D, Domek H, Rudzinska M, Czarnocka B. Podoplanin (PDPN) affects the invasiveness of thyroid carcinoma cells by inducing ezrin, radixin and moesin (E/R/M) phosphorylation in association with matrix metalloproteinases. BMC Cancer. (2019) 19:5239. doi: 10.1186/s12885-018-5239-z

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Suzuki-Inoue K, Inoue O, Ozaki Y. Novel platelet activation receptor CLEC-2: from discovery to prospects. J Thrombosis Haemostasis. (2011) 9:44–55. doi: 10.1111/j.1538-7836.2011.04335.x

PubMed Abstract | CrossRef Full Text

77. Bresson L, Faraldo MM, Di-Cicco A, Quintanilla M, Glukhova MA, Deugnier M-A. Podoplanin regulates mammary stem cell function and tumorigenesis by potentiating Wnt/β-catenin signaling. Development. (2018) 145:160382. doi: 10.1242/dev.160382

PubMed Abstract | CrossRef Full Text | Google Scholar

78. Krishnan H, Miller WT, Blanco FJ, Goldberg GS. Src and podoplanin forge a path to destruction. Drug Disc Today. (2019) 24:241–9. doi: 10.1016/j.drudis.2018.07.009

PubMed Abstract | CrossRef Full Text | Google Scholar

79. Bieniasz-Krzywiec P, Martín-Pérez R, Ehling M, García-Caballero M, Pinioti S, Pretto S, et al. Podoplanin-expressing macrophages promote lymphangiogenesis and lymphoinvasion in breast cancer. Cell Metabol. (2019) 30:917–36.e910. doi: 10.1016/j.cmet.2019.07.015

PubMed Abstract | CrossRef Full Text | Google Scholar

80. Yoshida T, Ishii G, Goto K, Neri S, Hashimoto H, Yoh K, et al. Podoplanin-positive cancer-associated fibroblasts in the tumor microenvironment induce primary resistance to EGFR-TKIs in lung adenocarcinoma with EGFR mutation. Clin Cancer Res. (2015) 21:642–51. doi: 10.1158/1078-0432.CCR-14-0846

PubMed Abstract | CrossRef Full Text | Google Scholar

81. Krishnan H, Ochoa-Alvarez J, Shen Y, Nevei E, Kephart D, Nguyen A, et al. Podoplanin (PDPN): novel biomarker and chemotherapeutic target. Cancer Res. (2015) 75:4375. doi: 10.1158/1538-7445.AM2015-4375

CrossRef Full Text | Google Scholar

82. Shiina S, Ohno M, Ohka F, Kuramitsu S, Yamamichi A, Kato A, et al. CAR T cells targeting podoplanin reduce orthotopic glioblastomas in mouse brains. Cancer Immunol Res. (2016) 4:259–68. doi: 10.1158/2326-6066.CIR-15-0060

PubMed Abstract | CrossRef Full Text | Google Scholar