Expression of m6A Regulators Correlated With Immune Microenvironment Predicts Therapeutic Efficacy and Prognosis in Gliomas

Background N6-methyladenosine (m6A) RNA methylation and tumor immune microenvironment played crucial roles in cancer development. However, their association in gliomas remains to be fully elucidated. Methods A total of 2144 glioma patients from CGGA, TCGA, and Rembrandt databases were extracted in our study, in which 325 were set as the training cohort and 1819 were defined as the validation cohort. Survival differences evaluated by Kaplan–Meier analysis between groups. Patients were clustered into subgroups by consensus clustering. ESTIMATE algorithm was applied to calculate immune and stroma scores. The infiltration of immune cells was characterized by TIMER algorithm. The risk signature was constructed by multivariate Cox regression analysis. Results Nineteen m6A regulators were highly expressed in glioma tissues. The expression of m6A regulators was associated with prognoses, grade, isocitrate dehydrogenase (IDH) status, and 1p19q status of gliomas. Two subgroups were identified by consensus clustering, in which cluster 1 was associated with favorable prognosis, high stroma and immune scores, and high immune infiltration. When the patients were divided into high risk and low risk groups based on their risk scores, we found that patients in the high risk group had poor prognoses. Besides, patients in the high risk group had a higher stroma and immune scores, and higher abundance of immune infiltration. These results were further verified in the validation cohort, which contained three independent datasets. Moreover, patients in the low risk group enjoyed better prognoses without chemoradiotherapy or single chemotherapy. Conclusion Our study revealed that m6A regulators could predict the prognosis and therapeutic efficacy, and were also associated with the immune microenvironment in gliomas.


INTRODUCTION
Gliomas are the most common primary malignancies in the central nervous system. The morbidity of gliomas is approximately 7 cases per 100,000 people, which accounts for the majority of primary brain tumors (Ostrom et al., 2019). Based on the World Health Organization (WHO) classification, gliomas are classified into four grades, in which grade 1 and grade 2 gliomas are defined as low-grade glioma whereas grade 3 and grade 4 gliomas are termed as high-grade glioma (Louis et al., 2016). Typically, patients with a higher grade glioma suffer from a worse prognosis. The 10-year survival rate of patients with low-grade glioma is 47% with a median survival time of 11.6 years (Ohgaki and Kleihues, 2005;Smoll et al., 2012). Although the low-grade glioma is relatively more optimistic compared to high-grade glioma, almost 70% of low-grade glioma will progress to high-grade one in a few years (Maher et al., 2001). The median survival time of patients with grade 3 glioma is approximately 3 years and that of patients with grade 4 glioma is about 13 months (Bleeker et al., 2012). Since patients with gliomas have such a poor prognosis, there is a clear urgent to find novel biomarkers to predict the prognosis. Previous studies have found that gliomas with isocitrate dehydrogenase (IDH) mutation and 1p19q codeletion indicate a relatively favorable survival (Eckel-Passow et al., 2015). Besides, O-6-methylguanine-DNA methyltransferase (MGMT) promoter methylation has been found to increase the chemosensitivity of temozolomide treatment, and it is a strong prognostic biomarker in patients with glioblastoma (Wick et al., 2014). However, additional studies are needed to explore novel biomarkers to predict the prognosis of glioma patients.
In recent years, numerous studies have proved that tumor immune microenvironment plays a crucial role in cancer progression and therapeutic efficacy (Quail and Joyce, 2017). The brain was long considered to be an "immune privileged" organ. However, this viewpoint was challenged when the lymphatic vessels were discovered along with the dural sinuses in mice (Louveau et al., 2015). Immune cells could infiltrate into the brain and form the immune microenvironment with other components. Multiple studies have identified the immunosuppressive status in gliomas, in which tumorassociated macrophages and regulatory T (Treg) cells in glioma microenvironment suppress the activities of T cells, mediating the immune escape of gliomas (Colombo and Piconese, 2007). Besides, glioma cells increase the expression of immunosuppressive factors such as programmed cell death 1 ligand (PD-L1) to reduce the presentation of antigens (Bloch et al., 2013;Xu et al., 2020). Therefore, the content of the immune microenvironment was closely associated with the efficacy of immunotherapy.
A previous study has revealed that m6A regulators are associated with glioma progression and prognosis (Chai et al., 2019). However, the relationship between m6A regulators and the immune microenvironment in gliomas remains unclear. In this study, we extracted data from the Chinese Glioma Genome Atlas (CGGA), The Cancer Genome Atlas (TCGA), and the Rembrandt datasets to explore the expression pattern and prognostic value of m6A regulators in gliomas. Besides, clustering subgroups and risk models were established based on the expression of m6A regulators to validate the predictive value of m6A regulators in risk stratification and prognosis. Moreover, the association between m6A regulators and the immune microenvironment was explored using the constructed signature. Additionally, the predictive value of m6A signature in the efficacy of chemotherapy and radiotherapy was also investigated. Our study aims to comprehensively assess the correlation of m6A regulators with prognosis, immune microenvironment, and therapeutic efficacy in gliomas.

Data Extraction
All RNA-seq data and clinical characteristics of enrolled samples were extracted from CGGA 1 , TCGA 2 , and Rembrandt databases. A total of 2144 glioma samples were enrolled in our study, in which 325 samples extracted from the CGGA database (CGGA325) were defined as the training cohort; 693 samples extracted from the CGGA database (CGGA693), 651 samples extracted from the TCGA database, and 475 samples extracted from the Rembrandt database were defined as the validation cohort; 20 normal tissues extracted from the CGGA database were termed as the control group. Data of CGGA325 and CGGA693 datasets were obtained in fragments per kilobase of exon model per million mapped fragments (FPKM) format. Data of TCGA dataset were batch normalized counts format, whereas that of Rembrandt dataset was normalized microarray format. The characteristics of patients in the training and validation cohorts were summarized in Table 1.

Identification of m6A Regulators
A total of 22 m6A regulators were identified according to previous studies Yang et al., 2018;Yi et al., 2020). The expression level of these genes was assessed between glioma samples and normal pairs.

Consensus Clustering
The consensus clustering was performed using "ConsensusClusterPlus" R package to categorize patients with gliomas into subgroups. The clustering algorithm was partitioning around medoids and the distance was measured by the euclidean metric.

Stroma and Immune Scores Calculation
The stroma and immune scores were measured by Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE) analysis using "estimate" R package (Yoshihara et al., 2013). Tumor purity was calculated according to the algorithm.

Immune Cells Infiltration Analysis
The abundance of six immune cells including B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells was calculated by Tumor Immune Estimation Resource (TIMER) algorithm . The role of copy number alternations (CNAs) of m6A regulators on immune cell infiltration was evaluated using TIMER algorithm 3 .

Construction of Risk Signature
The selection of candidate risk m6A regulators was performed by the least absolute shrinkage and selection operator (LASSO) analysis. Multivariate Cox regression analysis was used to profile independent prognostic genes. Variance inflation factor (Vif) and hypothesis testing were used to filter out genes with high collinearity. Risk score for each patient in the training and validation cohort was calculated by the following algorithm: Risk score = 0.052 × YTHDF2 + 0.025 × ALKBH5 + 0.029 × KIAA1429 + 0.023 × IGF2BP3. The patients were divided into high risk and low risk groups based on the mean value of the risk score.

Statistical Analysis
Statistical analyses and visualization were mainly performed using R version 4.0.2, GraphPad Prism version 8.0.1, and TBtools (Chen et al., 2020). Kaplan-Meier and log-rank analysis were used to evaluate the survival differences between grouped patients. Subgroup analysis was used to assess the stability of the risk signature, in which patients were divided into two subgroups based on their age (≤41 and >41 years old) and gender (female and male). Time-dependent receiver operating characteristic (ROC) curve analysis was used to evaluate the predictive value of constructed risk model using "survivalROC" R package. Student's t-test and one-way ANOVA analysis were used to estimate the differences between two groups and more than two groups. The correlation of gene expression was calculated using "spearman" method. Two-sided p < 0.05 was regarded as statistically significant.

Expression of m6A Regulators Was Involved in the Progression and Development of Gliomas
According to previous studies, we selected 22 m6A regulators for further investigation. The RNA expression of these genes was extracted from the training cohort. Results showed that 19 m6A regulators were significantly higher in gliomas compared with normal brain tissue (p < 0.05) ( Figure 1A); whereas no significant difference was detected between gliomas and normal tissue regarding the expression of METTL16, ALKBH5, and METTL5 (p > 0.05). Besides, Kaplan-Meier analysis revealed that 15 m6A regulators were independent prognostic genes and the other regulators including YTHDC1, FMR1, ZC3H13, METTL14, EIF3A, METTL3, METTL16, and IGF2BP1 were not associated the prognosis of glioma patients. Generally, patients with low expression of m6A regulators enjoyed favorable prognoses (p < 0.05); whereas patients with low expression of FTO had a less survival time compared to those with high expression of FTO (p < 0.05) ( Figure 1B). These results indicated that m6A regulators were highly involved in the progression and development of gliomas.

Expression of m6A Regulators Was Associated With Current Glioma Prognostic Markers
To further investigate the role of m6A regulators in gliomas, we explored the expression of m6A regulators in different subgroups of gliomas. Results showed that the expression of 15 m6A regulators, which had independent prognostic value, was significantly different in different grades of gliomas (Figure 2A). Generally, the high expression of m6A regulators indicated a higher grade of glioma; whereas the high expression of FTO indicated a lower grade of glioma. Additionally, the expression of 16 m6A regulators was significantly different in IDH mutant and IDH wild-type gliomas (p < 0.05), in which 6 genes were downregulated and 10 genes were up-regulated in IDH mutant gliomas compared with IDH wild-type gliomas ( Figure 2B). Moreover, the expression of 14 m6A regulators in 1p19q codeletion and non-codeletion gliomas was evident (p < 0.05), in which 11 genes were down-regulated and 3 genes were up-regulated in 1p19q codeletion gliomas compared with non-codeletion gliomas ( Figure 2C). A total of 11 m6A regulators had a consistent expression pattern in IDH and 1p19q subgroups, in which four genes were up-regulated (FTO, FMR1, EIF3A, and ZC3H13) and seven genes were down-regulated (ALKBH5, IGF2BP2, IGF2BP3, RBM15, WTAP, YTHDF1, and YTHDF3) in IDH mutant and 1p19q codeletion gliomas. Given that patients with the lower grade, IDH mutant, and 1p19q codeletion gliomas had better prognoses, our findings suggested that m6A regulators had promising values in predicting the prognosis of glioma patients.

Consensus Clustering for m6A Regulators Correlated With Glioma Prognosis and Immune Microenvironment
The unsupervised clustering method, consensus clustering, was performed to classify patients in the training cohort into subgroups based on the expression of m6A regulators. K = 2 was identified with optimal clustering stability ( Figure 3A-C). A total of 325 patients in the training cohort were clustered into two subgroups, 181 patients in cluster 1 and 144 patients in cluster 2. The expression pattern of m6A regulators in cluster 1 and cluster 2 was shown by the heatmap (Figure 3D). The expression level of m6A regulators (except FTO) were lower in cluster 1 compared with cluster 2. Besides, the stroma (p < 0.05), immune (p < 0.05), and ESTIMATE scores (p < 0.05) were significantly higher whereas the tumor purity was markedly lower (p < 0.05) in cluster 1 compared with cluster 2 ( Figure 3E). Moreover, the overall survival of patients in cluster 1 was significantly longer than those in cluster 2 (p < 0.05) ( Figure 3F). Additionally, the abundance of neutrophil, macrophage, and dendritic cell was significantly higher in cluster 1 compared with cluster 2 (p < 0.05), whereas no difference was detected regarding B cell, CD4+ T cell, and CD8+ T cell ( Figure 3G). These results indicated that the clustering subgroups based on m6A regulators were closely related to prognosis and the immune microenvironment in gliomas.

Construction of Risk Signature Based on m6A Regulators Expression in the Training Cohort
Then, the risk signature was established to evaluate the predictive value of m6A regulators. LASSO analysis filtered out seven m6A regulators with the minimum lambda value ( Figure 4A). After hypothesis testing, four m6A regulators including ALKBH5, IGF2BP3, KIAA1429, and YTHDF2 were selected with p-value of less than 0.05 and Vif of less than 2. The risk signature was constructed by multivariate Cox analysis (Supplementary Figure 1). Patients in the training cohort were divided into high risk and low risk groups based on risk scores ( Figure 4B). Multivariate Cox analysis revealed that risk score was the independent risk factor for patients in the training cohort ( Table 2). Time-dependent ROC analysis revealed that the predictive accuracy of risk score was highest in predicting 5-year survival ( Figure 4C). The area under curve (AUC) of 1-year, 3year, and 5-year survival was 0.801, 0.871, and 0.887, respectively. Moreover, patients in the low risk group had a longer overall survival compared with those in the high risk group (p < 0.05) ( Figure 4D); when patients were divided into subgroups based on their age and gender, those in the low risk group still had a longer survival time (p < 0.05) (Supplementary Figure 2A). The expression of four candidate m6A regulators was higher in high risk group compared with the low risk group (Figure 4E). The stroma, immune, and ESTIMATE scores were significantly higher (p < 0.05) whereas tumor purity was lower in the high risk group compared with the low risk group (p < 0.05) ( Figure 4F). In cluster 1, which had a better prognosis, the risk score was notably lower than cluster 2 (p < 0.05) ( Figure 4G). Besides, the risk score was elevated in the higher grade, IDH wild-type, and 1p19q noncodeletion subtype of glioma (p < 0.05) (Figures 4H-J). As for the histological subtype of glioma, the risk score was increased in the relatively malignant glioma (e.g., anaplastic astrocytoma vs astrocytoma; anaplastic oligodendroglioma/oligoastrocytoma vs oligodendroglioma/oligoastrocytoma); the glioblastoma, which was the most malignant glioma, had the highest risk score ( Figure 4K). Moreover, the abundance of CD8+ T cell, neutrophil, macrophage, and dendritic cell was significantly higher in the high risk group compared with the low risk group (p < 0.05) ( Figure 4L). These results indicated that the constructed risk signature exhibited a potent value in predicting the prognosis of glioma patients. The risk score was closely related with the immune microenvironment in gliomas.

Risk Signature Was Associated With Glioma Prognosis and Immune Microenvironment in the Validation Cohort
To further verify the predictive value of risk signature, we applied the risk score algorithm in the validation cohort. The risk score of patients in the validation cohort was calculated. All patients were divided into high risk and low risk groups based on the mean value of risk score (Figures 5A-C). Multivariate Cox analysis revealed that risk score was the independent risk factor for patients in the validation cohort (Supplementary Table 1). Besides, Kaplan-Meier analysis revealed that patients in the low risk group had a better prognosis compared with those in the high risk group (p < 0.05) (Figures 5D-F); when patients were divided into subgroups based on their age and gender, those in the low risk group still had a better prognosis (p < 0.05) (Supplementary Figure 2B-D). Time-dependent ROC analysis showed that predictive accuracy of risk score was highest in predicting long-term survival (Figures 5G-I). The AUC of 1-year, 3-year, and 5-year survival in CGGA693 dataset was 0.596, 0.659, and 0.682, respectively; that in TCGA dataset was 0.719, 0.762, and 0.787, respectively; that in Rembrandt dataset was 0.620, 0.704, and 0.727, respectively. Besides, the risk score was elevated in the high-grade, IDH wild-type, and 1p19q non-codeletion subtype of glioma in the validation cohort (p < 0.05), which was consistent with the training cohort (Figures 5J-L). Moreover, the risk score was relatively high in the malignant subtype of glioma and glioblastoma had the highest risk score, which was consistent with the training cohort (p < 0.05) ( Figure 5M). As for the immune microenvironment of gliomas, the stroma, immune, and ESTIMATE scores were significantly higher whereas the tumor purity was lower in the high risk group compared with the low risk group in the validation cohort (p < 0.05) (Figures 6A-C), which was consistent with the training cohort. Moreover, the abundance of CD8+ T cell, neutrophil, macrophage, and dendritic cell was significantly higher in the high risk group compared with the low risk group in the validation cohort (p < 0.05), which was consistent with the training cohort (Figures 6D-F). These results confirmed the establishment of the risk signature and its association with the prognosis and immune microenvironment of gliomas.
Additionally, the effect of CNAs of four candidate m6A regulators in risk signature on immune cell infiltration was further investigated to provide a novel insight into the association between risk score and immune cell infiltration. The arm-level gain of ALKBH5 and IGF2BP3 was closely associated with the infiltration of six types of immune cells (Supplementary  Figures 3A,B); whereas CNAs of KIAA1429 did not exhibit an association with immune cell infiltration (Supplementary Figure 3C). Besides, the high amplication and arm-level deletion of YTHDF2 were significantly related to immune cell infiltration (Supplementary Figure 3D). These results suggested that the four candidate genes of the risk signature were closely associated with the infiltration of immune cells.

Risk Stratification Was Associated With the Efficacy of Chemoradiotherapy and Immunotherapy
The association between therapeutic efficacy and risk stratification was also explored. In the training cohort, patients in the high risk group who received chemoradiotherapy had longer survival compared with those receiving other therapies (p = 0.0006), and radiotherapy exhibited a potent efficacy (p = 0.0006) ( Figure 7A); patients in the low risk group who received single radiotherapy had better prognoses compared with those without radiotherapy (p = 0.0028), and chemotherapy did harm to the survival rates (p = 0.0007) (Figure 7B). In the validation cohort, monotherapy of chemotherapy and radiotherapy had no significant effect on improving the survival rate of patients in the high risk group ( Figure 7C); however, patients in the low risk group who received single chemotherapy or radiotherapy had better prognoses compared with those receiving other therapies (p = 0.0074), and chemotherapy decreased the survival rates (p = 0.0244) ( Figure 7D). Besides, the risk score was correlated with the expression of immune checkpoints (i.e., PD-1, PD-L1, CTLA-4, and B7H3) in the training and validation cohorts, in which the risk score had a notable relationship with the expression of B7H3 (Figures 7E-H). These results showed that for patients in the low risk group, radiotherapy exhibited a potent efficacy whereas chemotherapy even decreased the survival rates. Therefore, the risk stratification was associated with the efficacy of chemoradiotherapy and might predict the efficacy of immunotherapy in gliomas.

DISCUSSION
According to WHO grade classification of gliomas, glioma was categorized into four grades, in which a higher grade indicated a worse survival. The past 30 years of research into glioma biology led to the discovery of hundreds of molecular alterations in grade II, III, and IV gliomas. Molecular exploration is in need of improved outcomes and the value of prognosis. Turkalp et al. (2014) found that IDH1/2, as an early event in the development of glioma, provided some reference in predicting glioma. Akagi et al. (2018) also showed that oligodendroglioma patients with IDH-mutant and 1p/19q co-deleted, rather than the WHO grade, demonstrated a better overall survival.
The previous studies inspired us to explore the changes in the expression of other genes and proteins in relation to the occurrence and development of glioma and to guide our intensive study on the role of glioma genes. With the booming development of high-throughput sequencing, various databases about genomic profiles have been established by researchers, which makes us to clearly acknowledge the genomic changes. Previous studies have revealed that immune-related gene signature was associated with prognosis and immune infiltration in gliomas (Qin et al., 2020;Zhang et al., 2020).
Their risk scores exhibited potent predictive values in the diagnosis and prognosis of gliomas. Huang et al. (2020c) found that ATP metabolism-related signature was associated with prognosis and immune microenvironment in gliomas. Besides, the risk score had a potent predictive accuracy identifying mesenchymal subtype glioma. Tao et al. (2020) revealed that epithelial-mesenchymal transition (EMT) based signature was also correlated with prognosis and clinicopathological features of glioma. However, few research explored the signature based on m6A regulators in gliomas.
Our study would fill the blank of m6A-related signature in the prediction of gliomas. The m6A regulator-based signature was closely associated with the immune microenvironment of glioma. Consensus clustering based on m6A regulator expressions was performed to divide the training cohort into two clusters. Survival analysis revealed that patients in cluster 1 had a favorable survival. Besides, the immune and stroma scores and immune cell infiltration was higher in cluster1. However, previous studies revealed that high immune and stroma scores as well as high infiltration of macrophages were associated with poor prognosis, which was reversed with our results in consensus clustering (Deng et al., 2020;Ni et al., 2020;Tian et al., 2020). Therefore, we further constructed a risk model to explore and validate the association between m6A regulators and the immune microenvironment. This study constructed the risk signature based on the expression of m6A regulators and also revealed that risk score was associated with prognosis and immune infiltration, which was similar to other studies. The risk stratification based on risk signature might be used to predict the efficacy of chemoradiotherapy and immunotherapy.
The risk signature correlates with glioma types, instead of other factors such as age and gender. In different histological subtypes of glioma, the risk score had a elevated tendency in the relatively malignant subtype. For example, the risk score was significantly higher in anaplastic astrocytoma compared to astrocytoma. Although the risk score of anaplastic oligodendroglioma or oligoastrocytoma was not significantly higher than astrocytoma, it is notably higher than oligodendroglioma or oligoastrocytoma. Moreover, the most malignant glioma, glioblastoma, had the highest risk score. Therefore, the risk score was closely associated with the histological subtype of glioma. In order to exclude the interference of other factors such as age and gender, we conducted subgroup analysis to further validate the predictive value of risk score. Results showed that patients with high risk score in the age and gender subgroups still had poor prognoses in four datasets. Besides, multivariate Cox analysis confirmed that the risk score was independent risk factor for glioma patients in four datasets. These results suggested that the risk signature was not affected by other factors such as age and gender.
This risk signature in compliance with current molecular pathology based survival prediction. IDH mutation, MGMT promoter methylation, and 1p19q codeletion in gliomas were associated with more favorable prognoses. A total of 11 m6A regulators had similar expression patterns in IDH and 1p19q subgroups, in which four were up-regulated and seven were down-regulated in IDH mutant and 1p19q codeletion gliomas. Unexpectedly, FMR1, EIF3A, and ZC3H13, whose high expressions indicated poor prognosis, were up-regulated in the IDH mutant and 1p19q codeletion gliomas since the two genomic mutations indicated favorable prognoses. However, IDH mutant and 1p19q codeletion only stand for a large proportion of patients with relatively good prognoses but not absolutely all patients have longer survival time. Therefore, those genes up-regulated in the IDH mutant and 1p19q codeletion subgroups may indicate worse prognoses within the subgroup. By the combination of the risk score and current indicators, we can predict the prognosis of glioma patients more accurately and efficiently.
The risk signature was established in multiple types of cancers with potent predictive values (Huang et al., 2020d;Wang et al., 2020;Zhang et al., 2020). On one hand, the risk score was calculated based on the coefficient and expression of candidate genes. Typically, a higher risk score indicated a poor survival. In our study, four candidate genes (ALKBH5, IGF2BP3, KIAA1429, and YTHDF2) were selected to construct the risk signature. The risk score exhibited promising value in predicting 5-year survival in the training and validation cohorts (AUC = 0.887,0.682,0.787,and 0.727). Although the predictive value in CGGA693 dataset was relatively low, its value was further validated by the other two datasets. Therefore, we believed that the risk score could predict long-term survival with a relatively high accuracy. On the other, the immune and stroma scores, as well as immune cell infiltration, were significantly higher in the high risk group. Moreover, the risk score was notably lower in cluster 1 compared with cluster 2, indicating the risk signature could replace the consensus clustering with a better evaluation efficiency. Further, the constructed risk signature was verified in the validation cohort, which exhibited consistent results with the training cohort. Four candidate genes were highly expressed in the high risk group and exhibited as risk factors for glioma patients. ALKBH5 was found to play an oncogenic role in epithelial ovarian cancer (Zhu et al., 2019). Besides, the inhibition of ALKBH5 could reduce the proliferation of glioblastoma stem-like cells via FOXM1 axis . Previous studies revealed that IGF2BP3 promoted tumor proliferation and progression in colorectal cancer, gastric cancer, and glioblastoma (Suvasini et al., 2011;Zhou et al., 2017;Xu et al., 2019). In the end, KIAA1429 was proved to facilitate cancer progression in breast and live cancer (Lan et al., 2019;Qian et al., 2019). Similarly, YTHDF2 was shown to be involved in the progression and angiogenesis in hepatocellular carcinoma, and its inhibition exhibited promising efficacy in acute myeloid leukemia Li et al., 2018;Hou et al., 2019;Paris et al., 2019). These findings revealed that the dysregulation of candidate m6A regulators was involved in the tumor progression in various kinds of cancer.
The low expression of FTO were associated with poor prognosis in gliomas, although there have been a contentious and divisive. Among 22 m6A regulators, 19 were highly expressed in glioma tissues compared with normal tissues, indicating that these m6A regulators might play crucial roles in gliomas. Survival analysis revealed that 15 m6A regulators were identified as prognostic indicators, in which the high expression of 14 m6A regulators and the low expression of FTO were associated with poor prognosis, whereas the high expression of FTO indicated a favorable prognosis. These results suggested that although FTO was highly expressed in glioma tissues, it might play a protective role for patients with gliomas.  found that the expression of m6A regulators was dysregulated in osteosarcoma, and low expression of FTO and high expression of YTHDF3 was associated with poor prognosis, which was consistent with our findings. Besides, FTO could suppress the self-renewal of ovarian cancer stem cells via inhibiting cyclic adenosine monophosphate (cAMP) signaling pathway (Huang et al., 2020b). However, Su et al. (2020) revealed that the inhibition of FTO also exhibited potent anti-tumor effects in multiple types of cancers. FTO inhibition enhanced the cytotoxicity of T cell and reduced immune evasion by suppressing the expression of immune checkpoint genes. These findings suggested that FTO played a complicated role in different types of cancer, and its role in glioma needed further investigation.
Those m6A regulators were involved in the immune microenvironment of gliomas. When patients in the training and validation cohorts were divided into high risk and low risk groups based on the risk score, the high risk group had a higher stroma and immune scores and a lower tumor purity. Previous studies found that stroma and immune scores calculated by ESTIMATE algorithm were meaningful for the classification and prognosis of glioblastoma and lower-grade gliomas (Jia et al., 2018;Deng et al., 2020). Besides, the abundance of immune cells such as CD8 + T cell, neutrophil, macrophage, and dendritic cell was highly infiltrated in the high risk group. Further, we identified that CNAs (i.e., arm-level gain, arm-level deletion, and high amplication) of four candidate genes were markedly associated with immune cell infiltrations.
The risk stratification could facilitate the determination of therapeutic options to improve prognoses. Current management for gliomas includes surgery, chemotherapy, and radiotherapy, which cannot reach optimal remission although great advancement has been achieved. In our study, patients in the high risk group who received chemoradiotherapy exhibited the most favorable survival, whereas radiotherapy showed great efficacy. However, this result was not consistent in the validation cohort. Notably, patients in the low risk group who received chemoradiotherapy showed worse prognosis compared with those receiving chemotherapy or radiotherapy. In contrast, single radiotherapy exhibited promising efficacy in the training and validation cohorts although radiotherapy had no effect in improving the survival rates of patients. However, it should be noted that a part of patients with or without radiotherapy would receive chemotherapy, which was demonstrated to impair the survival rate of patients. Therefore, single radiotherapy should be considered for patients in the low risk group. Although patients received single chemotherapy exhibited favorable prognosis in the validation cohort, the finding was not consistent with the training cohort. Since chemotherapy was found to impair the survival rates of patients with low risk in training and validation cohorts, chemotherapy should be deliberated for patients in the low risk group. Additionally, the risk score was positively correlated with the expression of B7H3. The novel chimeric antigen receptor T (CAR-T) product targeting B7H3 exhibited promising efficacy in the treatment of glioblastoma both in vitro and in vivo. Currently, the clinical studies were undergoing to evaluate the efficacy of B7H3 CAR-T cell in the treatment of glioblastoma (NCT04077866) and pediatric gliomas (NCT04185038) (Xu et al., 2020). The risk stratification based on the risk score might help predict the efficacy of B7H3 CAR-T therapy.
There are still limitations in our study. Firstly, our findings are based on open accessed databases without our cohort. Secondly, the interactions between m6A regulators and immune cells were not validated by experiments. Last but not least, the regulatory mechanism of m6A regulators in glioma immune microenvironment is not elucidated, which warrants further investigation to provide a better understanding.

CONCLUSION
To sum up, our study comprehensively assessed the expression pattern and prognostic value of m6A regulators in gliomas. The expression of m6A regulators was associated with the classification of glioma subtypes. Besides, the consensus clustering and risk signature based on m6A regulator expressions could be used to predict prognosis and were associated with the immune microenvironment in gliomas. Additionally, the risk stratification could facilitate the prediction of the efficacy of chemoradiotherapy and might be associated with the efficacy of immunotherapy. These findings indicated that m6A regulators might be potential biomarkers indicating the prognosis and therapeutic efficacy for patients with gliomas and were associated with glioma immune microenvironment. Further studies are needed to explore regulatory mechanisms of m6A regulators in glioma progression and therapeutic efficacy.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
ZL and CL conceived, designed, and supervised the study. SX and LT drafted the manuscript. SX, LT, and GD collected the data. SX performed all data analysis. All authors reviewed and approved the final manuscript.