The Non-N6-Methyladenosine Epitranscriptome Patterns and Characteristics of Tumor Microenvironment Infiltration and Mesenchymal Transition in Glioblastoma

Background An increasing number of RNA modification types other than N6-methyladenosine (m6A) modification have been detected. Nonetheless, the probable functions of RNA modifications beyond m6A in the tumor microenvironment (TME), mesenchymal (MES) transition, immunotherapy, and drug sensitivity remain unclear. Methods We analyzed the characteristics of 32 non-m6A RNA modification regulators in 539 glioblastoma (GBM) patients and the TME cell infiltration and MES transition patterns. Using principal component analysis, a non-m6A epitranscriptome regulator score (RM score) model was established. We estimated the association between RM score and clinical characteristics, TME status, GBM subtypes, and drug and immunotherapy response. Results Three definite non-m6A RNA modification patterns associated with diverse biological pathways and clinical characteristics were identified. The high RM score group was characterized by a poor prognosis, enhanced immune infiltration, and MES subtype. Further analysis indicated that the high RM score group had a lower tumor mutation burden as well as a weaker response to immunotherapy. The higher RM score group may benefit more from drugs targeting the EGFR and WNT signaling pathways. Conclusion Our study exposed the potential relationship of non-m6A RNA modification regulators with clinical features, TME status, and GBM subtype and clarified its therapeutic value.


INTRODUCTION
Glioblastoma (GBM) is the most ordinary malignant intracranial cancer (1). Patients with GBM have an approximately 14-16-month median survival time, and GBM is highly resistant to standard therapies because of its extraordinary immunosuppressive microenvironment and mesenchymal (MES) transition (2)(3)(4).
Detected at an approximately 10-fold lower m 6 A level, m 1 A is involved in the methylation of the adenine base on the first nitrogen atom as well as carries positive electricity (22). The m 1 A methylation process is mediated by methyltransferases (writers) consisting of TRMT61B, TRMT6, and TRMT61A, while the removal process is catalyzed via demethylases (erasers) containing ALKBH3, ALKBH1, and FTO (23)(24)(25). Various specific RNA-binding proteins (readers), including YTHDC1, YTHDF3, YTHDF2, and YTHDF1, can recognize the m 1 A motif, thus affecting m 1 A functions (26). m 1 A is necessary for tRNA structure stabilization, proper rRNA biogenesis, and methylated mRNA translation (22,27). If 2′-O-methyladenosine is the first ribotide behind the m 7 G cap of mRNA, it could be methylated at the N 6 position to render m 6 A m , which is related to mRNA metabolism (28,29). PCIF1 is the writer that creates the m 6 A m modification, and FTO is the eraser of m 6 A m (25,30). APA cleaves mRNA at more than one site and adds poly (A) tails to generate different transcripts of the 3′-untranslated region (3′ UTR) or coding regions, which regulates the function, stability, and translation efficiency of RNAs (31). PABPN1, NUDT21, CLP1, PCF11, and CPSF6 can regulate the APA process (32).
At levels similar to those of m 1 A, m 7 G, a modification carrying a positive charge at the 5′ cap, is necessary for mRNA export, splicing, and translation (33). In addition, the METTL1-catalyzed m 7 G tRNA methylome is important for the translation of mRNA (34). N m is relatively widespread in rRNA, tRNA, snRNA, and microRNA and can be installed by FBL (35,36). In addition, N m is essential for adjusting the structure and function of ribosomes (37).
Appearing in both coding and non-coding regions, m 5 C sites are chiefly CG-rich regions that accumulate in the UTRs of mRNA (38). The writer of m 5 C is NSUN2, while the eraser is TET2 (14,39,40). In addition, YBX1 has been reported to be an m 5 C reader (41). m 5 C modification plays an essential role in mRNA export, tRNA structure stabilization, and rRNA translational fidelity (39,42). Catalyzed by NAT10 and THUMPD1, ac 4 C was subsequently discerned in serine and leucine 18S rRNA and tRNAs. ac 4 C, the sole acetylation modification in eukaryotic RNA, regulates mRNA translation and ribosome biogenesis (43,44).
RNA editing is a kind of programmed posttranscriptional mechanism altering nucleotides in selected transcripts (45). Catalyzed by ADAR, ADARB1, and ADARB2, RNA editing can transform the sequence as well as alter transcriptional procedures (46). Accumulating in 3′UTRs, U RNA editing alters the protein level (47). UPP1 can catalyze the phosphorolysis of uridine to uracil (48,49).
To thoroughly comprehend the significance of RNA modifications beyond m 6 A, investigation of the crosstalk among diverse patterns of RNA alterations is urgently needed. Ten kinds of non-m 6 A RNA modification "regulators" may build a complex and significant regulatory network in GBM, which might help elucidate GBM tumorigenesis mechanisms.
Immune checkpoint blockade therapy (ICB therapy), also known as immunotherapy, has delivered promising clinical outcomes for various cancers; nevertheless, it commonly exhibits a poor response due to the tumor microenvironment (TME) (50,51). Unfortunately, a large number of patients with GBM have not experienced outstanding survival benefits; that is, immunotherapy for GBM is far from reaching clinical expectations (52). Accordingly, to enhance the effect of immunotherapy for GBM, it is crucial to thoroughly investigate the immunosuppressive TME. Recent studies have revealed that m 6 A plays a significant role in complex TME formation (53,54). Estimating m 6 A patterns of a type of tumor to characterize TME infiltration could guide more effective immunotherapy strategies (55,56). Non-m 6 A RNA modification as well as correlative regulators are highly related to the microenvironment of immune cells as well as tumor cells. TET2 deficiency in Treg cells results in T-cell activation (57). NAT10 regulates the function of CD4 + T cells (58). However, few studies have systematically analyzed the non-m 6 A RNA modification patterns of individual tumors, while these RNA modifications play a non-negligible role in tumorigenesis (16,59).
Based on genetic transcription signatures, GBM can be categorized into three subtypes (MES; classical, CL; as well as proneural, PN) (3). MES transition promotes radiochemotherapy resistance in GBM (60,61). The MES subtype is particularly aggressive among these three subtypes, while the PN subtype yields the best prognosis (62,63). In addition, the MES-subtype GBM promotes the formation of an immunosuppressive TME, while the TME advances the MES transition in GBM (62,64). An increasing number of studies have demonstrated that non-m 6 A RNA modification regulators are significantly correlated with MES transition. It has been reported that ADAR and YBX1 promote MES transition by regulating oncogenic microRNA maturation (65,66). Nevertheless, research on RNA modifications beyond m 6 A is still not as mature as that for m 6 A for many reasons, including technology limitations. Accordingly, it is significant to investigate the non-m 6 A epitranscriptome patterns and characteristics of TME infiltration and MES transition in GBM.
In this study, we explored somatic mutation data for 390 GBM cases from The Cancer Genome Atlas (TCGA, http://portal.gdc.cancer.gov/; http://xena.ucsc.edu/) and integrated gene expression data from 539 GBM samples from TCGA and Chinese Glioma Genome Atlas (CGGA, https://www. cgga.org.cn/) cohorts to evaluate the RNA modification patterns. We discovered that non-m 6 A RNA modification patterns were associated with TME cell-infiltrating characteristics as well as MES transition. Next, according to differentially expressed genes (DEGs), we established a non-m 6 A epitranscriptome "regulator" score (RM score) to evaluate the effect on the non-m 6 A epitranscriptome in individual patients. Eventually, we verified the pertinence of the RM score in distinguishing the posttranscriptional and transcriptional events as well as evaluated its importance in predicting the response to targeted and ICB therapy.

Data Collection and Processing
The study workflow chart is shown in Figure S1A. Genetic transcription data as well as clinical information on GBM patients were obtained from TCGA and the CGGA. Somatic mutation counts and copy number variation (CNV) were obtained from the TCGA database. In total, 539 GBM samples (CGGA_batch_1:139, CGGA_batch_2:249, TCGA:151) were gathered in this study for further analyses. We used sva package's "ComBat" algorithm to correct non-biotechnology deviations causing batch effects. R Bioconductor packages and R (version 4.10) were used to analyze the data.

Unsupervised Clustering for 32 Non-m 6 A RNA Modification Regulators
We extracted 32 non-m 6 A RNA modification regulators and their expression from the TCGA and CGGA databases. These non-m 6 A RNA modification regulators included 10 RNA modification types, which are listed in Table S1. Using unsupervised cluster analysis, different non-m 6 A RNA modification patterns concerning 32 non-m 6 A RNA modification regulatory factors were identified, and the patients were classified into different clusters. We applied the Consensus Cluster Plus package to execute the steps above with1,000 repetitions to guarantee classification stability.

TME Cell Infiltration Estimation
First, single-sample gene set enrichment analysis (ssGSEA) was used to identify tumor immune-infiltrating cell abundance in GBM samples. The markers of 28 immune-related cells and types were obtained from the dataset of Bindea et al. (67). Using an unsupervised hierarchical clustering algorithm, 539 GBM samples were assigned to two groups based on immune infiltration. Next, ESTIMATE was used to calculate the immune score as well as stromal score (68). Then, we used CIBERSORT (http://cibersort.stanford.edu/) to measure the infiltration levels of 22 distinct immune cells among the GBM samples.

Estimation of the MES/PN Score
Verhaak et al. identified three prognostic subtypes by integrated genomic analysis (69). There are significant prognostic differences between different subtypes, and accurate identification of subtypes can also guide treatment strategies. We downloaded the gene sets "VERHAAK_GLIOBLASTOMA_MESENCHYMAL" and "VERHAAK_GLIOBLASTOMA_PRONEURAL" from the MSigDB database v7.4 and quantified the MES and PN scores of the GBM samples using the ssGSEA algorithm.

Identification of DEGs Among Non-m 6 A RNA Modification Patterns
To identify non-m 6 A RNA modification-related genes, DEGs among three different non-m 6 A RNA modification characteristics were arranged via the limma R package. The significance criterion for DEGs was p-value <0.01. To evaluate the function of DEGs, Gene Ontology (GO) analysis was performed to calculate biological process (BP), molecular function (MF), and cellular component (CC) terms using DAVID (https://david.ncifcrf.gov/) with a p-value cutoff of <0.05.

Small Interfering RNA and Plasmid Transfection
Small interfering RNAs (siRNAs) targeting UPP1 and plasmid overexpression of FTO were synthesized (GenePharma, China). SiRNAs and plasmids were transfected with Lipofectamine ™ 3000 reagent (Thermo Fisher Scientific) according to the protocol of the manufacturer.

Western Blotting
Protein was extracted from GBM cell lines. The blots were incubated with primary antibodies against UPP1 (Abcam, UK), FTO (Cell Signaling Technology, USA), and ACTIN (Cell Signaling Technology).

Transwell Assay
To evaluate migratory ability, GBM cells were added to the top chamber in DMEM without FBS, and the bottom chamber was filled with 10% FBS DMEM. To measure the ability to recruit macrophages, macrophages were added to the top chamber in RPMI-1640 without FBS, and the bottom chamber was filled with 10% FBS RPMI-1640 and growth media of different groups of GBM cells. After 24 h of incubation, the membrane was fixed in 4% paraformaldehyde and stained with crystal violet.

Constructing the RM Scoring System to Evaluate Individual GBM Samples
First, overlapping DEGs distinguished from non-m 6 A RNA modification clusters were chosen for analysis. Then, we analyzed the prognostic value for each overlapping DEG via a univariate Cox regression model, and significant prognostic genes were identified. Next, the PCA algorithms were used to generate the RM score. Both principal component 2 and principal component 1 were selected to describe the RM scoring system.

Immunofluorescence
Tumor tissues were obtained from 12 patients treated for GBM at Qilu Hospital. RNA quantity and quality were estimated via an Agilent 2100 bioanalyzer (Agilent Technologies) and a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific). The mRNA library was prepared by using the NEBNext ® Ultra ™ RNA Library Prep Kit (Beijing Novel Bioinformatics Co., Ltd.). The HiSeqX platform (Illumina) on the Illumina standard protocol was used for RNA sequencing. Next, we performed a PCA algorithm based on the expression of DEGs to calculate the RM scores of the 12 GBM patients. We selected three cases with the highest as well as three cases with the lowest scores and performed immunofluorescence (IF) staining. After blocking with 5% goat serum (Gibco, USA) for 30 min, samples were incubated with primary antibodies at 4°C overnight. The following primary antibodies were used: rabbit SOX2 antibody (3579, Cell Signaling Technology, 1:400, USA), mouse CD44 antibody (3570, Cell Signaling Technology, 1:400), mouse CD163 antibody (ab156769, Abcam, 1:100, UK), and rabbit FoxP3 antibody (12653, Cell Signaling Technology, 1:400). Next, the samples were incubated with fluorophore-conjugated secondary antibodies at 37°C for 1 h. Alexa Fluor 594 anti-rabbit antibody (A11037; Invitrogen; 1:400) and Alexa Fluor 488 antimouse antibody (A11029; Invitrogen; 1:400) were used. DAPI was used to stain the cell nuclei. All immunostained samples were analyzed using Leica Application Suite Software and a Leica TCS SP8 confocal system. ImageJ was used to analyze the fluorescence intensity.

Summary of Genomic and Clinical Information of the Immunotherapy Cohort
We analyzed gene expression-related clinical information of the immunotherapy cohort and identified two independent immunotherapy cohorts: advanced urothelial tumor with atezolizumab treatment (IMvigor210 cohort, http://researchpub.gene.com/IMvigor210CoreBiologies/packageVersions/) and metastatic melanoma with pembrolizumab treatment (GSE78220 cohort). We obtained the expression data and clinical annotations of the GSE78200 cohort from the Gene Expression Omnibus (GEO) database.

Association Analysis Between Drug Sensitivity and RM Score
Cancer Cell Line Encyclopedia (CCLE) RNA-seq data were obtained from https://portals.broadinstitute.org/ccle/. Cancer cell line drug responses, measured as area under the curves (AUCs) and drug pathways, were obtained from Genomics of Drug Sensitivity in Cancer (GDSC; https://www.cancerrxgene. org/downloads). Spearman correlation analysis was applied to estimate the association between drug sensitivity and RM score.

Comparison of the m 6 A Score and the RM Score
The Bayesian information criterion (BIC) as well as Akaike information criterion (AIC) was applied to compare the two models. A model with lower BIC and AIC values was identified as a better model.

RM Score Among Clinical Traits in Pan-Cancer
Univariate Cox regression analysis was used to investigate the time-dependent prognostic value of RM score in Pan-cancer. Overall survival (OS) and progression-free survival (PFS) were selected to study the relationship between RM score and prognosis. In addition, the correlation between RM score and 22 immune cell infiltration was calculated in 33 cancers, respectively.

Statistical Analysis
The statistical analyses were generated by R-4.0.1. Distance and Spearman correlation analyses were applied to calculate 32 RNA modification expression correlation coefficients. Kruskal-Wallis and one-way ANOVA tests were applied to appraise difference comparisons of more than three collections. To assess the relationship between patient survival and the RM score, the optimal cutoff point for survival information for each dataset was determined using the survminer package. The log-rank test and Kaplan-Meier analysis were applied to assess survival among different clusters as well as RM score groups. The hazard ratio (HR) of DEGs was calculated via a univariate Cox regression model. To evaluate whether the RM score serves as an independent predictor for survival, sex, age, and isocitrate dehydrogenase (IDH) status were considered variables in the multivariate Cox regression model analysis. Statistical analysis was two-sided, as well as p <0.05 was regarded as statistically significant.

Landscape of Genetic Variations in Nonm 6 A RNA Modification Regulators in GBM
A total of 32 non-m 6 A RNA modification regulators (Table S1) were included in the current study. Figure 1A shows that these nonm 6 A RNA modification regulators can dynamically and reversibly add, remove, and recognize non-m 6 A RNA-modified sites and potentially change biological functions, such as splicing, stability, export, translation, and degradation of RNA. We first summarized the somatic mutation prevalence of 32 non-m 6 A RNA modification regulators among GBM. Among 390 samples, 24 showed genetic alterations of non-m 6 A RNA modification regulators, with a frequency of 6.67% ( Figure 1B). These mutations were multifarious and included splicing-related mutations, missense mutations, and deletions ( Figure 1B). The mutation frequency of ADARB2 was the highest, while PABPN1, TRUB1, TRUB2, NAT10, PUS1, and TRMT61A showed no mutations among the 390 GBM patients ( Figure 1B). Further analyses revealed mutation co-occurrence relationships among YTHDF2, YTHDF1, TRMT61B, THUMPD1, METTL1, FTO, FBL, YTHDF3, and UPP1 ( Figure S1B). Moreover, the investigation of 32 non-m 6 A RNA modification regulators exhibited that CNV-related mutations were widespread. METTL1, CPSF6, PUS1, ADARB2, ADAR, and UPP1 showed widespread CNV amplification, while PABPN1, FBL, ALKBH1, and TRUB1 had CNV deletions ( Figure 1C). The CNV alteration locations of 32 non-m 6 A RNA modification regulators on chromosomes are presented in Figure 1D. To determine whether genetic variations affected nonm 6 A RNA modification regulator expression among GBM cases, we analyzed the expression of these regulators between GBM and normal samples to determine whether CNV alterations could result in perturbations in non-m 6 A RNA modification regulator expression. Non-m 6 A RNA modification regulators with CNV amplification demonstrated substantially higher expression in GBM tissues than in normal brain tissues, while those with deletions exhibited the opposite trend ( Figures 1C, E). Moreover, we determined that the regulator expression patterns significantly varied among the three subgroups ( Figure 1F). O 6 -methylguanine DNA methyltransferase (MGMT) methylation, isocitrate dehydrogenase 1 (IDH1) mutation, and chromosomal 1p/19q codeletion have been found to confer a favorable prognosis in GBM (70). The same phenotypes were investigated in groups classified based on molecular subtypes ( Figures S1D-F). Based on the expression of 32 non-m 6 A RNA modification regulators, we distinguished patients with GBM from normal controls ( Figure  S1G). The analyses revealed high genetic and transcriptomic landscape heterogeneity among non-m 6 A regulators between GBM and normal cases, indicating that genetic variations and expression alterations among non-m 6 A RNA modification genes play an essential role in GBM progression and occurrence.

Non-m 6 A RNA Modification Patterns Mediated by 32 Regulators
To obtain an overall understanding of the expression pattern, 539 GBM cases from the TCGA and CGGA datasets that contained clinical information and OS data were enrolled. A univariate Cox regression model revealed that 10 of 32 non-m 6 A regulators were associated with the survival of GBM patients ( Figure S1C).
To investigate the relationships among regulators, we calculated pairwise correlations among the expression of the 32 regulators in the GBM cases ( Figure S2A). We identified that the expression of DKC1, YBX1, TRMT6, and PUS7 was positively correlated with that of many other regulators, while the expression of UPP1 and ALKBH3 was negatively correlated with that of other regulators ( Figure S2A). In addition, the protein-protein interactions between non-m 6 A RNA modification regulators are shown in Figure S1H. The comprehensive landscape of the intricate associations between non-m 6 A RNA modification regulators and the prognostic value for GBMs was visualized with the non-m 6 A RNA modification regulator network (Figure 2A; Figure S2A). We identified significant correlations among 32 non-m 6 A RNA modification regulators, which revealed that the crosstalk between the regulators of erasers, writers, and readers might participate in the formation of various non-m 6 A RNA modification patterns.
Next, we used ConsensusClusterPlus to sort patients with three non-m 6 A RNA modification patterns based on the expression profiles of 32 selected non-m 6 A RNA modification regulators. After unsupervised clustering, 144 cases were termed Cluster_A, 227 cases were identified in Cluster_B, and the other 168 cases were termed Cluster_C ( Figures 2B, S2B). The expression of the 32 regulators in different clusters is shown in Figure 2C. In the survival analysis of non-m 6 A RNA modification patterns, a dominant survival disadvantage was found for the Cluster_B non-m 6 A RNA modification pattern ( Figure 2D). To identify the biological functions of the three non-m 6 A RNA modification patterns, GSVA was performed. As exhibited in Figure 2E, Cluster_B was markedly enriched in pathways associated with the formation of the immunosuppressive TME and MES transition, such as extracellular matrix (ECM) receptor interaction, TGF b, JAK STAT, WNT, and NOTCH signal paths. In addition, Cluster_A was prominently enriched in pathways regulating the cell cycle, cell aging, and DNA damage, while Cluster_C was markedly related to apoptosis processes ( Figure S2C).

TME Cell Infiltration and MES Transition Characteristics in Distinct Non-m6A RNA Modification Patterns
SsGSEA was used to assess TME-infiltrating cells in GBM tissues. According to the abundance of infiltrating immune cells, we classified the GBM samples into high and low infiltration groups ( Figure 3A). Patients with 1p19q codeletion status or IDH mutation status were mainly enriched in the low infiltration group ( Figure 3A). Analyses of TME cell infiltration via the CIBERSORT algorithm confirmed that Cluster_B was obviously enriched in innate immune cell infiltration, including the infiltration of mast cells, natural killer (NK) cells, macrophages, and dendritic cells (DCs) ( Figure 3B; Table S2). In addition, abundant stromal elements in the TME are associated with immunosuppression (71). We uncovered that the stromal score of Cluster_B was the highest among the three clusters via the ESTIMATE algorithm ( Figure 3B). When we investigated the MES and PN scores, we found that the MES score was the highest and the PN score was the lowest in Cluster_B ( Figure 3C). The MESsubtype GBM tended to exhibit an immunosuppressive TME. Previous studies demonstrated that there were abundant immune cells in tumors with an immune-excluded group, and immune cells participated in the formation of an immunosuppressive TME (72,73). We found that three non-m 6 A patterns possessed definite TME infiltration characteristics. Cluster_B was termed an immuneexcluded group, characterized by abundant innate immune cell infiltration, stromal activation, and a high MES score; Cluster_A was termed an immune-inflamed group, marked by abundant immune cell infiltration; and Cluster_C was termed an immunedesert phenotype, characterized by rare immune cell infiltration ( Figures 3A-C and S3A). In addition, GBMs with an immuneexcluded subtype had the worst prognosis ( Figure 2D). Furthermore, to analyze the immune tolerance as well as activity condition of each pattern, we selected LGALS9 , LAG3, ICOS,  IL23A, LDHA, CTLA4, CD274, PTPRC, PDCD1, IL12A,  ICOSLG, CD28, TNFRSF4, VTCN1, PDCD1LG2, TNFRSF18,  TNFSF18, TNFSF4, CD80, CD86, CD8A, B2M, TNFSF9,  YTHDF1, PVR, FGL1, CD40, SIGLEC15, LAMA3, CD40LG,  HAVCR2, TNFRSF9, JAK2, JAK1, and LDHB as immune checkpoint-related molecules (74,75). We uncovered that diverse immune checkpoint-related genes were significantly overexpressed in Cluster_B and Cluster_A ( Figure S3C). In addition, we chose CD44, CHI3L1, FN1, SERPINE1, and TIMP1 as MES markers and OLIG2, DLL3, NCAM1, ASCL1, and SOX2 as PN markers (69,76,77). The expression of MES markers was high in Cluster_B, while the expression of PN markers was high in Cluster_C ( Figure S3D).
In addition, ssGSEA revealed that the immunosuppressive TME and MES transition pathways, such as the NF-kB and STAT3 signal paths, were obviously enriched in Cluster_B, whereas the monocyte and macrophage chemotaxis pathways were enriched in Cluster_A ( Figure 3D). Figure S3B illustrates the proportion of GBM patients with the indicated clinical status and molecular traits in the three groups with different non-m 6 A RNA modification patterns. Patients over 60 years old and those who were alive at the last follow-up were primarily in Cluster_C, while patients with IDH wild-type status and no MGMT methylation were mainly found in Cluster_B ( Figure S3B). We analyzed the specific relationship between the TMEinfiltrating cell type as well as non-m 6 A RNA modification regulators ( Figure S3E). Patients were classified into low or high immune infiltration, MES, and PN score groups. Thirty-two non-m 6 A RNA modification regulator expression levels in different groups are shown in Figures S3F-H. We focused on UPP1, a U regulator, which had a significant positive correlation with numerous TME-infiltrating immune cells ( Figure S3E). Patients were classified into low and high UPP1 expression subgroups, and a beneficial prognosis was presented in patients with low UPP1 ( Figure S4A). We used ESTIMATE to assess the overall levels of immune cell infiltration as well as the stromal score. The results showed that patients with a high expression of UPP1 exhibited high immune and stromal scores ( Figure S4B). The specific difference in TME-infiltrating levels between patients with low and high UPP1 expression was explored. We observed that tumors with high UPP1 showed obviously increased infiltration of Treg T cells and M2 macrophages ( Figure S4B). We uncovered that increased UPP1 resulted in the overexpression of immune checkpoint-related genes ( Figure  S4C). In addition, the expression of MES markers was higher, while that of PN markers was lower in patients with high UPP1 expression ( Figure S4D). Subsequent pathway enrichment analyses showed that high UPP1 samples exhibited an obvious enhancement of the immunosuppressive TME and MES transition-related pathways, including the NF-kB, STAT3, and epithelial-mesenchymal transition (EMT) signal paths ( Figure  S4E). From the above results, we speculated that UPP1-mediated U may promote the formation of an immunosuppressive TME and facilitate MES transition. In addition, we investigated the characteristics of one m 1 A and m 6 A m eraser, FTO, which was negatively associated with various TME-infiltrating immune cells ( Figure S3E). Patients with a high FTO exhibited a better prognosis ( Figure S5A). We then found that a high expression of FTO was related to low immune and stromal scores ( Figure S5B). We also discovered that a low FTO was related to the overexpression of a few immune checkpointrelated genes ( Figure S5C). MES marker expression was higher, while PN marker expression was lower in patients with low FTO expression ( Figure S5D). In addition, patients with low FTO expression showed prominent enrichment of the hypoxia pathway, STAT3 signaling pathway, and EMT signaling pathway ( Figure S5E). Thus, FTO may inhibit the formation of an immunosuppressive TME and the process of MES transition.
To assess the role of UPP1 and FTO in diverse cellular processes of GBM cells, several experiments were performed. Downregulation of UPP1 or overexpression of FTO resulted in a significant decrease in the EdU-positive cell percentage ( Figures 4A, B). In addition, the results of the Transwell assays revealed that UPP1 silencing and FTO overexpression reduced the number of GBM cells migrating to the membrane ( Figure 4C). Interestingly, we found that UPP1 knockdown or FTO overexpression significantly inhibited the capability of conditioned medium from GBM cells to recruit THP1differentiated macrophage ( Figure 4D).

Non-m 6 A RNA Modification Phenotype-Related DEGs in GBM
To characterize the potential biological behavior of the three non-m 6 A patterns, we confirmed 243 RNA phenotype-related DEGs ( Figure S6A). Then, a univariate Cox regression model was applied to perform prognostic analysis for each DEG in the signature ( Table S3). The DEGs with significant prognostic value were extracted for further analysis. To further corroborate this regulatory mechanism, unsupervised clustering methods were executed according to 150 DEGs ( Figure S6B). The marked enriched biological functions of these DEGs are presented in Figure S6C. The unsupervised clustering algorithm sorted the cases into two subtypes: gene.cluster_1 and gene.cluster_2 ( Figure 5A). The two patterns were identified in gene.cluster_1 and gene.cluster_2, marked by different signature genes ( Figure 5A). Patients with 1p19q codeletion status or IDH mutation status were mainly concentrated in gene.cluster_2 ( Figure 5A). In addition, in the two non-m 6 A RNA modification DEG clusters, noticeable differences in non-m 6 A RNA modification regulator expression were analyzed ( Figure 5B). A total of 242 patients with GBM were termed gene.cluster_2, which was validated to be correlated with a favorable prognosis ( Figure 5C). To our surprise, gene.cluster_1 was remarkably abundant in several immune cell infiltrates (Figures 5D, E). In addition, the immune cell infiltration of gene.cluster_1 was stronger than that of gene.cluster_2 ( Figures 5D, E). However, we found that the stromal score of gene.cluster_1 was higher than that of gene.cluster_2 ( Figure 5E). In addition, gene.cluster_1 possessed a higher MES score and lower PN score than gene.cluster_2 ( Figure 5F). These results indicate that patients in gene.cluster_1 had a highly immunosuppressive TME, which is related to a poor prognosis. Moreover, we determined that diverse immune checkpoint-related genes were significantly overexpressed in gene.cluster_1 ( Figure S6D). In addition, the expression of MES markers was high in gene.cluster_1, while the expression of PN markers was high in gene.cluster_2 ( Figure  S6E). Moreover, ssGSEA revealed that the immunosuppressive TME and MES pathways were significantly enriched in gene.cluster_1, whereas T-cell receptor signaling pathways were enriched in gene.cluster_2 ( Figure 5G).

Generation of the RM Score and Analysis of Clinical Traits
The above results showed that non-m 6 A plays an essential role in shaping different TME landscapes as well as facilitating MES transition. However, these analyses only assessed the characteristics of the population and could not correctly evaluate the detailed information of non-m 6 A RNA modification of each GBM case. Taking into account the complexity and individual heterogeneity of non-m 6 A RNA modifications, we generated a scoring model to assess the non-m6A RNA modification patterns of individual GBM patients according to phenotype-related DEGs; this system was termed the RM score. The relationships among the non-m 6 A RNA modification patterns, MES score, gene.cluster, and RM score were visualized in an alluvial diagram and in Table S4 ( Figure 6A). The Wilcoxon test revealed an obvious difference in the RM score between the gene.clusters. Gene.cluster_2 showed a lower median score than gene.cluster_1, which indicated that a high RM score could be closely linked to the immunosuppressive TME and MES transition signatures ( Figure 6B). Moreover, Cluster_B showed a significantly increased RM score compared with Cluster_A and Cluster_C ( Figure 6B). In addition, patients with stronger immune cell infiltration possessed higher RM scores ( Figure S7A). We investigated the value of the immune, MES, and PN scores in predicting patient prognosis. A high immune score, high MES score, or low PN score indicated a poor prognosis ( Figure S6F). In addition, we observed a significantly negative relationship between the MES score and PN score ( Figure S6G). Next, we evaluated the significance of the RM score in predicting GBM outcomes. With a cutoff value of −10.21, GBMs were classified into high or low RM score groups. GBMs with lower RM scores presented a significantly prolonged survival time ( Figure 6C). To evaluate the stability of the RM scoring system, we used the RM score in other independent GBM databases to verify its prognostic significance, and the results indicated that a low RM score was correlated with better clinical benefit ( Figures S8A, B). We analyzed whether the RM score could be regarded as an independent factor of survival for GBM. Multivariate Cox regression analysis was performed, confirming that the RM score was an independent and robust factor for predicting GBM patient outcomes [ Figure S8C; HR (low RM score) = 0.53 (0.39-0.73)]. We also examined the survival prediction efficiency for the combination of RM score with immune score, RM score and MES score, and RM score and PN score. The results presented that the low RM and low immune score groups showed the best OS among the groups ( Figure S7B). Furthermore, we observed that patients with low RM and MES scores or low RM scores and high PN scores had the most significant survival benefits (Figures S8D, E). We assessed the value of the RM scoring feature to evaluate the efficacy of radiotherapy or chemotherapy in GBMs. Among the patients who received adjuvant radiotherapy or chemotherapy at the same time, those with low RM scores exhibited the most significant therapeutic benefit ( Figure S7C). Furthermore, patients over 60 years old, patients alive at the last follow-up, and patients with IDH mutations, 1p/19q codeletion, and MGMT methylation were mainly found in the low RM score group ( Figures S8F, G).
Further analysis showed that the immune cell infiltration of patients with high RM scores was significantly higher than that of patients with low RM scores, while patients with high RM scores had higher stromal scores ( Figure S7D). In addition, patients with high RM scores had higher MES scores and lower PN scores ( Figure S7E). GBM patients with the MES subtype were characterized by higher TME levels, whereas GBM patients with the PN subtype were associated with lower TME levels. To exclude the influence of phenotypic and tumor immune infiltration correlations on the analysis results, we analyzed the correlation between RM scores and TME levels in patients with PN and MES subtypes ( Figure S9). Moreover, to reveal the significance of the RM score in estimating the level of TME immune infiltration and MES transition, we studied the expression of immune checkpointrelated genes and MES/PN marker genes in different RM score groups. We found that most immune checkpoint-related genes were upregulated in the high RM score group ( Figure S7F). In addition, the expression of MES markers was higher in high RM score patients, while the expression of PN markers was higher in low RM score patients ( Figure S7G). We examined the relationship between the RM score and known biological signatures using Spearman analysis. A heatmap of the association matrix exhibited that the RM score was positively related to immunosuppressive TME and MES transition signatures ( Figure S7H).
To verify the role of the RM score in evaluating the characteristics of TME infiltration and MES transition, IF staining was performed. The proportions of M2 macrophages (CD163 + ) and FoxP3 + Tregs infiltrated in tumors with high RM scores were higher than those in tumors with low RM scores ( Figure 6D). Moreover, CD44 expression was also markedly increased in the high RM score group, while SOX2 expression was decreased ( Figure 6D). These results further validated that tumors with high RM scores showed an immunosuppressive phenotype.

The Correlation Between the RM Score and TMB
Increasing evidence has demonstrated a significant correlation between responsiveness to immunotherapy and tumor mutation burden (TMB) (78). Taking the clinical significance of TMB into account, we further analyzed the relationship between the RM score and the TMB. By comparing the TMB values for patients with high and low RM scores, we concluded that the RM score was negatively associated with the TMB value ( Figure 7A). Patients with higher RM scores exhibited a lower TMB than patients with lower RM scores ( Figure 7B). We divided the patients into different groups according to the immune setting of the TMB, as described previously (79). As shown in Figure 7C, we observed that patients with a high TMB had prolonged survival times compared with those with a low TMB. Considering the contrary survival significance of the TMB and RM score, we analyzed the cross-influence of these scores on the prognosis of patients with GBM. Stratified survival analyses revealed that patients with low RM scores together with high TMB scores had the most significant survival benefits ( Figure 7D). Furthermore, we evaluated the somatic variant distribution of GBM driver genes among the high and low RM score groups using the maftools package. As shown in Figures 7E, F, the low RM score group presented more extensive TMB than the high RM score group. An increasing number of studies have shown that patients with a high TMB are more sensitive to anti-PD-1/PD-L1 immunotherapy (80). In addition, the alteration frequency of various genes was different between the low and high RM score groups ( Figures 7E, F). In conclusion, the above findings corroborated that the difference in non-m 6 A RNA modification patterns could be an important factor that predicted response efficiency to anti-PD-1/PD-L1 immunotherapy. These results might lead to novel strategies for exploring the mechanism of non-m 6 A RNA modification composition and gene mutation in ICB therapy.

The Predictive Value of the RM Score in Predicting ICB Therapy Response
Significant efforts have been made to screen biomarkers that predict ICB treatment response; some previously studied biomarkers include the TMB as well as the PD-L1 expression level (81). Considering that the RM score is related to the TME, we examined the value of the RM score in predicting the response of GBMs to immunotherapy. The exploration was based on two independent immunotherapy studies ( Figures 8A-H). We observed that patients who had lower RM scores exhibited obviously prolonged OS in the anti-PD-L1 research (IMvigor210; Figure 8A) and anti-PD-1 research (GSE78220; Figure 8G) (74,82). The patients in the IMvigor210 cohort exhibited diverse degrees of response to anti-PD-L1 blockers, including partial response (PR), progressive disease (PD), complete response (CR), and stable disease (SD). We explored the differences in RM scores for patients with diverse responses to ICB therapy and concluded that the proportion of patients in the response groups (CR and PR) was significantly lower in the high RM score group than in the low RM score group, while the proportion of patients in the no/limited response groups (SD and PD) showed the opposite trend, indicating that the RM score could reveal the response of patients to ICB therapy ( Figure 8B). The expression of PD-L1 on tumor cells (TCs) was substantially related to the RM score. The TC1 group exhibited the lowest RM score and was obviously different from the TC0 as well as TC2 + groups ( Figure 8C). Moreover, we observed that the PD-L1 level of immune cells (ICs) was positively related to the RM score, with IC0 exhibiting the lowest RM score and IC2 exhibiting the highest RM score ( Figure 8D). By analyzing the relationship of the RM score with tumor neoantigen burden, we observed that patients with low RM scores together with a high neoantigen burden exhibited the most prolonged survival time ( Figure 8E). In addition, patients with high RM scores showed significantly high levels of PD-L1 ( Figure 8F). Finally, we confirmed the significant clinical response and therapeutic advantages of anti-PD-1 immunotherapy in patients with low RM scores versus those with high RM scores ( Figures 8H, S8H). The RM score was also negatively correlated with the TMB in the GSE78220 cohort ( Figure S8I).

Exploration of the Therapeutic Significance of the RM Score
To investigate the value of the RM score in predicting drug sensitivity, we calculated the correlation between the RM score and the AUC for drugs in multiple cancer cell lines. According to Spearman correlation analysis, we selected 68 drugs for which the RM score and drug sensitivity were significantly correlated in the GDSC database ( Figure 8I) (83). The RM score was positively correlated with sensitivity to 12 drugs, including picolinic acid, chromatin histone methylation inhibitor EPZ5676, EGFR inhibitor AZD3759, gefitinib, erlotinib, and WNT inhibitors MN.64 and XAV939; the RM score was negatively correlated with drug sensitivity (and this positively correlated with drug resistance) for 56 drugs (Figures 8I, J). Furthermore, we explored the signal path targeted by the selected drugs. We observed that drugs whose sensitivity was positively related to the RM score targeted the EGFR and WNT signaling pathways. In contrast, drugs whose sensitivity was negatively related to the RM score targeted the cell cycle and apoptosis signal path ( Figure 8J). In conclusion, these results revealed that the RM score model might serve as a potential factor for exploring proper therapy strategies.

Construction of the m 6 A Score and Comparison of the m 6 A and RM Scores
A total of 20 m 6 A RNA modification regulators were enrolled in this study. Next, we classified patients into different m 6 A RNA modification patterns according to the expression profiles of 20 selected m 6 A RNA modification genes. After unsupervised clustering, 199 cases were termed m 6 A Cluster_A, 197 cases were termed m 6 A Cluster_B, and the other 143 cases were identified in m 6 A Cluster_C ( Figure 9A). In the survival analysis of m 6 A RNA patterns, m 6 A Cluster_B exhibited a particularly outstanding survival benefit ( Figure S10A; p = 0.0018). In addition, we identified 295 DEGs using the limma package. Then, univariate Cox regression analysis was applied to determine the prognostic value for each DEG (Table  S5), and the DEGs with significant prognostic value were extracted for further analysis. To further corroborate this regulatory information, an unsupervised clustering method was used according to the 172 DEGs. The patients were divided into two groups, m 6 A gene.cluster_1 and m 6 A gene.cluster_2, using the unsupervised clustering algorithm ( Figure 9B). A total of 272 patients with GBM were termed m 6 A gene.cluster_2, which was validated to be associated with a better prognosis ( Figure S10B). Based on these DEGs, we generated a scoring system to assess the m 6 A RNA modification pattern of each GBM patient. The relationships among m 6 A RNA modification patterns, MES score, m 6 A gene.cluster, and m 6 A score were visualized in an alluvial diagram and in Table S6 ( Figure 9C). Patients with low m 6 A scores demonstrated a prolonged survival time ( Figure 9D). The Wilcoxon test revealed a significant difference in the m 6 A score between the m 6 A gene.clusters. M 6 A gene.cluster_2 showed a lower median score than gene.cluster_1, which indicated that a high m 6 A score could be closely linked to poor prognosis ( Figure S10C). Moreover, compared with the other clusters, m 6 A Cluster_B showed a significantly decreased m 6 A score ( Figure S10D). In addition, patients with high RM scores exhibited higher immune, stromal, and MES scores and lower PN scores ( Figures 9E, F).
The value of the RM and m 6 A scores in predicting patient outcomes and the TME was compared using the BIC and AIC values. With lower BIC and AIC values, the RM scoring model exhibited a better description of TME infiltration and MES transition than the m 6 A scoring model ( Figure 9G). There was no obvious difference between the RM and m 6 A scores in predicting prognosis ( Figure S10E).

RM Score Correlates With Immune Infiltration and Survival Prognosis in Pan-Cancer
According to the forest plots, a positive association was obvious between RM score and OS in KIRC, LGG, and LAML ( Figure  S11A). In addition, the PFS forest plot confirmed the role of RM score as a risk factor in KIRC and LGG ( Figure S11A). As presented in Figure S11B, RM score is positively related to the TMB in BRCA, KIRC, LIHC, and THCA, whereas a negative association was observed in LAML and THYM. UCEC, and UCS ( Figure S11C). In BLCA, CESC, KICH, LAML, and TGCT, RM score was positively associated with M2 macrophage infiltration ( Figure S11C).

DISCUSSION
Increasing evidence has demonstrated that various RNA epigenetic modifications have an important effect on infection, the TME, and antitumor proliferation by interacting with several regulators. However, most studies have focused on a single regulator or a single kind of RNA modification, such as m 6 A, and the mutual relationships and effects of regulators of various non-m 6 A RNA modifications in cancer are not fully understood. Few studies have demonstrated that a given non-m 6 A RNA modification influences tumorigenesis in GBM, probably because non-m 6 A RNA modifications are not as abundant as m 6 A modifications.
Here, we revealed global non-m 6 A RNA alterations at the genetic as well as transcriptional levels and showed mutual correlations in GBMs. Surprisingly, there are complicated associations among 32 non-m 6 A RNA alteration regulators. We constructed three non-m 6 A RNA modification patterns based on 32 regulators. Cluster_A was characterized by abundant immune cell enrichment and classified as an immune-inflamed subtype; Cluster_B was characterized by the suppression of immune function as well as MES transition, corresponding to the immune-excluded subtype; and Cluster_C was marked by rare immune cell infiltration, corresponding to the immune-desert subtype. Although the immune-excluded subtype was related to abundant immune cell infiltration, the immune cells were inactive because of the distribution of effective immune cells in the tumor stroma. The stroma could permeate the tumor itself or might be confined to the tumor envelope, making immune cells appear to be inside the tumor (84). In addition, the MES score was the highest and the PN score was the lowest in Cluster_B. The immunosuppressive TME and MES transition pathways, such as the NF-kB and STAT3 signaling pathway, were markedly enriched in Cluster_B. The MES-subtype GBM tends to exhibit an immunosuppressive TME, while an immunosuppressive TME promotes GBM malignant progression via MES transition (62,64). Our data demonstrated that there may be a positive feedback loop between the immunosuppressive TME and MES transition, resulting in a poor prognosis.
Furthermore, the mRNA expression differences among nonm 6 A RNA modification patterns have been confirmed to be significantly related to RNA modification-, immune-, and MESrelated biological signaling pathways. The DEGs were considered a non-m 6 A RNA modification-related gene set. Similarly, two nonm 6 A-related gene subtypes were identified based on the DEGs, which were obviously related to the immunosuppressive TME as well as MES transition signaling pathways. The MES subtype is particularly aggressive among these three subtypes, while the PN subtype has the best prognosis (85). This demonstrated again that the various non-m 6 A RNA modifications were of great significance in the formation of TME patterns. In conclusion, a systematic evaluation of the non-m 6 A RNA modification patterns will improve our understanding of TME cell infiltration and MES transition characterization.
Considering the heterogeneity of non-m 6 A RNA modifications across individual tumors, a method for quantifying the non-m 6 A RNA modification patterns of individual tumors is urgently needed. Accordingly, we developed a scoring model, the RM score, to evaluate the nonm 6 A RNA modification pattern of individual GBMs. Patients with low RM scores demonstrated a prominent survival benefit. Patients in Cluster_B, characterized by an immune-excluded phenotype, exhibited a higher RM score. In addition, the RM scor e was signifi cantl y positively as sociated with immunosuppressive TME and MES transition signatures.
Immunotherapy is a developing field but is far from reaching clinical expectations in GBM patients because of the immunosuppressive TME. The MES transition not only promotes radiochemotherapy resistance but also shapes the immunosuppressive TME in GBM. Numerous studies have emphasized the non-negligible interaction between non-m 6 A RNA regulators, the TME, and the MES transition. The RM scoring model could predict patient prognosis and sensitivity to radiochemotherapy as well as immunotherapy. Patients who were sensitive to radiochemotherapy and ICB therapy were mainly enriched in the low RM score group. In addition, the RM score was markedly positively correlated with immunosuppressive TME and MES transition signatures.
Our findings also revealed a significantly negative association between the RM score and TMB. Increasing evidence has demonstrated that patients with a high TMB present acceptable responses to immunotherapy. A high TMB status was correlated with a favorable prognosis. In addition, RM score could predict the response of a patient to anti-PD-1/L1 immunotherapy. Patients receiving ICB therapy in the IMvigor210 and GSE78220 cohorts were assessed, and we identified that the RM score was markedly decreased in patients responding to ICB, which validated the predictive significance of the RM score model. Overall, this study showed that patients with low RM scores had more obvious benefits from immunotherapy.
In addition, we investigated the possible treatment outcome of non-m 6 A RNA modification regulators in GBMs. The RM score was positively correlated with the resistance of drugs that targeted the cell cycle and apoptosis signaling pathways and negatively correlated with the sensitivity of drugs that targeted the EGFR and WNT signaling pathways. These findings indicated that patients with higher RM scores might be more suitable for drugs targeting the EGFR and WNT signal paths instead of drugs targeting the cell cycle or apoptosis signaling pathways. Thus, non-m 6 A RNA modification patterns could be used as a qualified predictor that can be used to predict the clinical outcome of targeted therapies or chemotherapy. Our findings reveal novel possibilities for improving the efficacy of ICB treatment, revealing different TME phenotypes as well as MES subtypes and suggesting the potential for personalized and precise immunotherapy for GBM.
Increasing evidence has demonstrated that the m 6 A modification pattern is closely associated with the TME, prognosis, and other clinical characteristics. We constructed an RM scoring system to evaluate the potential roles of RNA modifications beyond m 6 A in the TME, MES transition, immunotherapy, and drug sensitivity. To further evaluate the superiority of the RM scoring model compared with the m 6 A modification pattern, we constructed an m 6 A scoring model and compared two models using the BIC and AIC algorithms. Finally, we confirmed that the RM scoring system is better than the m 6 A scoring system in assessing TME and GBM subtypes.
In this study, we systematically assessed the non-m 6 A RNA modification patterns of 539 GBMs on the basis of 32 regulator genes and associated these patterns with TME infiltration as well as MES transition characteristics. We constructed the RM score model to predict patient prognosis and the response to immunotherapy and targeted therapy. The systematic assessment of non-m 6 A RNA modification patterns could not only improve our knowledge of the crosstalk of RNA modifications but also contribute to the development of more personalized and precise immunotherapy regimens.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethics Committee of the Qilu Hospital (Jinan, China) and performed in accordance with the relevant guidelines and regulations; the patient data were acquired from the publicly available datasets, and informed consent was given.

AUTHOR CONTRIBUTIONS
JX, ZG, and GL designed this work. JX, ZG, HX, and XG integrated and analyzed the data. JX, ZG, KL, ZZ, PZ, and LD wrote this manuscript. JX, ZG, YF, SW, HW, QW, and RZ edited and revised the manuscript. All authors contributed to the article and approved the submitted version.