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 m6A 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 m6A 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 m6A 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.
Several studies have revealed the correlation between tumor microenvironment (TME) infiltrating immune cells and m 6 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 6 A was reduced in the high immunity subtype of lung adenocarcinoma, indicating that m 6 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.

Selection of m 6 A RNA Methylation Regulators
We used 12 m 6 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 6 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 6 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 6 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: 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 lowrisk 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.

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: Forward Reverse

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

Cell Viability and Cell Death Measurement
Cell viability was measured using the CellTiter-Glo R 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.

Consensus Clustering of m 6 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 6 A RNA methylation regulator in clinical prognosis, we performed consensus clustering based on gene expression of 12 key m 6 A RNA regulators in TCGA datasets. Due to the expression analogy of m 6 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 threecluster 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 6 A regulators, we assessed the interaction and correlation among these regulators (Supplementary Figure 3). 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.

Immune Landscape Was Significantly Associated With m 6 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).
In order to analyze the relationship between m 6 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 6 A RNA methylation regulators and immune cell infiltration, we assessed the relationships between the expressions of m 6 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 6 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 6 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 moduletrait 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 6 A RNA Methylation Regulators and Immune Infiltration
We explored the relationship between the expression levels of five hub genes and m 6 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 6 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 6 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 (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.

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).

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)(34)(35). However, no previous studies have investigated the efficacy of the combination of m 6 A and immune infiltration.
Here, we developed and validated a novel 4-gene prognostic model based on the combination of m 6 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 6 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)(45)(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 tumorinfiltrating 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)(71)(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)(75)(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 macrophagespecific PDPN knockout restrained lymphangiogenesis and reduced lymphatic cancer spread (79). PDPN-positive cancerassociated 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 userfriendly R shiny web app (http://www.houshixu.cn:3838/sampleapps/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.