The CGGA RNA sequencing data [per kilobase of transcript per million mapped reads (RSEM)] and the relevant clinical information were downloaded from the CGGA data portal (http://www.cgga.org.cn/) as the training set. Similarly, the TCGA RNA sequencing data [fragments per kilobase of transcript per million mapped reads (FPKM)] and clinical information were downloaded from the TCGA data portal (https://www.cancergenome.nih.gov/) and used as a verification set. Moreover, the somatic mutation data from glioma patients were also downloaded from the TCGA data portal as MAF files. After the removal of samples with missing data on survival and WHO grade, we obtained 306 (CGGA dataset) and 616 (TCGA dataset) glioma patients. Thereafter, the CGGA and TCGA RNA sequencing data were normalized with Expectation-Maximization algorithm (log2 transformation) for the subsequent analysis. The clinical information for the CGGA and TCGA datasets was listed in Supplementary Table 1. Moreover, the seven RNA:m⁵C methyltransferases were screened according to the published literature (11, 17, 18, 20).

All human glioma tissue samples used in the study were obtained from patients who were operated on in The Second Affiliated Hospital of Nanchang University, Jiangxi, from 2015 to 2019. According to the World Health Organization classification, glioma samples were divided into four groups: Grade I, Grade II, Grade III, Grade IV. Twenty-five brain samples, including 9 high-grade glioma samples (WHO grade III, IV), 10 low-grade glioma samples (WHO grade II), and 6 non-neoplastic brain tissues (NBTs) were used as a control collected from the surgery of intractable epilepsy cases. The study was approved by the medical ethics committee of the Second Affiliated Hospital of Nanchang University and was performed in accordance with the approved guidelines. Informed consents were acquired from each glioma patient. The tumor tissues were immediately frozen and stored at −80°C.

To explore the function of gliomas' RNA:m⁵C methyltransferases-related genes, we identified three subtypes based on the RNA:m⁵C methyltransferase expression profiles of 306 patients with gliomas from the CGGA dataset using the R package “ConsensusClusterPlus” (50 iterations, 80% resampling rate Pearson correlation, http://www.bioconductor.org/) (27). The appropriate number of subtypes was calculated using cumulative distribution function (CDF) and consensus matrices. Thereafter, we used principal components analysis (PCA) to detect differential gene expression amongst the three subtypes using R package pca3d and rgl (28).

Function and interaction of the seven RNA:m⁵C methyltransferases were predicted using the String website (https://string-db.org/) and the R package “corrplot.” Gene Ontology (GO) analysis was performed with Metascape (http://metascape.org/) to annotate the function of differentially expressed genes in the different subtypes. Furthermore, the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment was performed using Database for Annotation, Visualization and Integrated Discovery (DAVID; https://david.ncifcrf.gov/; version 6.8) and visualized in the imageGP website (http://www.ehbio.com/ImageGP/). Gene set enrichment analysis (GSEA; version 4.0.3) was performed using the JAVA program and downloaded from the website (http://software.broadinstitute.org/gsea/). The hallmark gene set (h.all.v6.0.symbols.gmt) was also downloaded. Thereafter, the hallmark gene set was determined to be significantly enriched following normalization (p < 0.05) and a false discovery rate (FDR < 0.25).

Univariate Cox regression analysis was performed to identify the genes significantly related to survival (P < 0.05). Thereafter, six RNA:m⁵C methyltransferases (NOP2, NSUN2, NSUN4, NSUN5, NSUN6, NUSN7) were screened with the LASSO multivariate Cox regression algorithm using the R package “glmnet” (version 3.0). Finally, the signature genes and coefficients in the risk score signature were constructed based on the most suitable penalty parameter λ. The risk score formula we used was:

where Coef_i is the coefficient and Exp_i is the normalized expression of each signature gene. The risk score system of five RNA:m⁵C methyltransferases (NOP2, NSUN4, NSUN5, NSUN6, NUSN7) was constructed in the CGGA dataset, and evaluated with the TCGA dataset. Patients were ranked into high-risk and low-risk groups using the median risk score. Moreover, the genomic alterations of these glioma patients were analyzed using the R package “Maftools” (29).

The immunohistochemistry pathological specimen results from the prognostic RNA:m⁵C methyltransferases were downloaded from The Human Protein Atlas (https://www.proteinatlas.org/). Staining intensity, quantity and patients' information can be inquired about online.

To assess the mRNA expression levels of RNA:m⁵C methyltransferases in glioma patients we used trizol reagent (Invitrogen, Carlsbad, CA), extracted total RNA from brain tissues and the total RNA (2 μg) was reverse transcribed into cDNA, according to manufacturer's instructions. The qRT-PCR was performed on The LightCycler® 480 Real-Time PCR System. The expression level of genes was analyzed and normalized to GAPDH and the 2^−ΔΔCT method was performed to evaluate the fold change in gene expression. The following primer sequences (Shanghai GenePharma Co., Ltd.): NOP2 forward 5′-AAATGGGAGAAGGTGGCGTC-3′, and reverse 5′-CTCTCGGACATTAACCCGCA-3′; NSUN4 forward 5′-ATGTGCCCTGTACCACAGAC-3′, and reverse 5′-GGTAAGAGCTTCCAGTCCCC-3′; NSUN5 forward 5′-TGCTTCTCAGAGCCAACGAA-3′, and reverse 5′-CTTCCAGACTGTCTCACTGCAT-3′; NSUN6 forward 5′-CTCGGTGACAGGGAGGTCTA-3′, and reverse 5′-AGGCATCTCTCATGGTGAGG-3′; NSUN7 forward 5′-TCAGCTCCAGATTGCGATTC-3′, and reverse 5′-CCCTCGTGCCTTAGGATCAC-3′.

The One-way ANOVA and R packages “limma” were used to compare the expression level of RNA:m⁵C methyltransferases in gliomas with different WHO grades and IDH status. The t-tests were performed to contrast the expression levels in gliomas for clinicopathological features. A chi-square test was conducted to compare the distribution of clinical and molecular pathological features between the different groups stratified by risk scores. Univariate, multivariate, LASSO Cox regression and Kaplan-Meier analyses were performed to construct and evaluate the risk signature using the R packages “glmnet” and “survival” (30, 31). Roc curve analysis was performed to predict the OS of glioma patients using the R package “survivalROC.” All statistical analyses were two-sided and conducted using R software (version 3.6.1; https://www.r-project.org/) and SPSS (version 26). P < 0.05 was considered statistically significant.

m⁵C modifications of RNA play an important role in the occurrence and progression of tumors and influence many biological functions. We comprehensively studied the relationship between each RNA:m⁵C methyltransferase and the clinicopathological characteristics of gliomas, for example, WHO grades and isocitrate dehydrogenase (IDH)-mutant status. The relationship between aberrant expression of RNA:m⁵C methyltransferases and WHO grades are shown using heatmaps (Figures 1A,B). As shown in the heatmaps, the expression of almost all RNA:m⁵C methyltransferases was significantly correlated with WHO grade. The significant differentially expressed RNA:m⁵C methyltransferases include NOP2, NSUN2, NSUN4, NSUN5, NSUN6, and NSUN7. Thereafter we quantitatively analyzed the expression of these differentially expressed genes in the CGGA dataset (Figure 1C) and verified our findings using the TCGA dataset (Figure 1D). As shown in the figure, NOP2, NSUN2, NSUN4, NSUN5, and NSUN7 were up-regulated and NSUN6 was down-regulated with an increase of WHO grade.

Mutation of the IDH gene has been reported in gliomas, particularly in low-grade gliomas (LGGs), and the prognostic value of the mutation has been confirmed by many authors in the literature. Based on these findings, we investigated the relationship between aberrant expression of RNA:m⁵C methyltransferases and IDH-mutant status in LGGs. The results showed different expression of NSUN3, NSUN4, NSUN5, NSUN6, and NSUN7 between IDH-mutant and IDH-wildtype status in both CGGA (Figure 1E) and TCGA (Figure 1F) datasets. In patients harboring IDH mutants, the expression of NSUN4, NSUN5, and NSUN7 increased, while the expression of NSUN3 and NUSN6 decreased. We also studied the expression of RNA:m⁵C methyltransferases in glioblastomas (GBM) stratified by IDH mutant status, and findings indicated that NSUN5, NSUN6, and NSUN7 were still differentially expressed (Supplementary Figures 1A,B).

In addition, we predicted the gene mutation frequencies of seven RNA:m⁵ methyltransferases in the cBioPortal of the Cancer Genomics Database and verified these findings using the TCGA dataset, resulting in that mutations of these RNA:m⁵C methyltransferases were rare (Supplementary Figure 1C and Supplementary Table 4). Even for the top-ranked RNA:m⁵C methyltransferases like NOP2, the mutation frequencies were only 3%. It demonstrated that the aberrant expression of RNA:m⁵C methyltransferases may not be generated by a genetic mutation.

To investigate the close connection between RNA:m⁵C methyltransferases and the clinicopathological features of gliomas, we systematically investigated the function, interaction, and correlation of RNA:m⁵C methyltransferases. We found that all the RNA:m⁵C methyltransferase genes were involved in various types of methylation, among which NOP2, NSUN3, NSUN4, and NSUN5 were mainly involved in rRNA methylation and NSUN2, NSUN6, as well as NSUN3, participated in tRNA methylation. Their function and interaction were supported by text mining and co-expression (Figure 2A). Thereafter, a Pearson correlation analysis was performed to study the expression profile of seven RNA:m⁵C methyltransferases in the CGGA (Figure 2B) and TCGA (Supplementary Figure 2A) datasets. The expression of NOP2, NSUN2, NSUN4, NSUN5, and NSUN7 were positively correlated in gliomas. Whereas the expression of NSUN6 was significantly negatively correlated with NSUN4, NSUN5 as well as NSUN7. These results were consistent with the quantitative analysis of RNA:m⁵C methyltransferases expression profile in gliomas, indicating that the expression levels of NOP2, NSUN2, NSUN4, NSUN5, and NSUN7 were shown to be positively correlated with the malignant progression of glioma, while the expression of NSUN6 were negatively correlated with glioma.

Based on the RNA:m⁵C methyltransferase expression profiles of 306 patients with gliomas in the CGGA dataset, we used unsupervised consensus clustering analysis to identify three subtypes, namely MC1, MC2, and MC3 (Figure 2C). Clearly, k = 3 seemed to be a relatively stable distinction of the samples in the CGGA dataset with clustering stability increasing from k = 2 to k = 9 (Supplementary Figures 2B–E). Furthermore, the PCA analysis was used to compare the transcription profiles of these three subgroups. The results showed that they could be adequately divided into three distinct clusters (Figure 2E). Next, we investigated the relationship between the RNA:m⁵C methyltransferase expression profiles of these three subtypes and clinicopathological features of gliomas (Figure 2D). In these three subtypes, MC2 subtype compared to MC3 subtype and MC3 subtype compared to MC1 subtype were significantly correlated with a higher grade (P < 0.0001), IDH-wildtype status (P < 0.0001), 1p/19q-noncodel status (P < 0.0001), older average age at diagnosis (P < 0.0001), and receipt of additional chemotherapy (P < 0.01) (Supplementary Table 2). Moreover, we found significant differences in overall survival between the three groups (P < 0.0001). The survival of patients who fell into the MC2 subtype was obviously shorter than for the other two subtypes (Figure 2F). Moreover, the unsupervised consensus analysis of these RNA:m5C methyltransferases in the TCGA dataset were consistent with the results of the CGGA dataset (Supplementary Figures 2F–J). The results above indicating that consensus clustering of RNA:m⁵C methyltransferases could identify subtypes with different clinicopathological features and prognosis in gliomas.

To investigate the different clinicopathological features and overall survival rates of the three groups in gliomas, we annotated the biological processes of specific genes associated with the MC2 subtypes. A comparison was performed with the other two groups, in which 664 genes were up-regulated (log FC > 1.5, normalized P < 0.01 and FDR = 0.05) and 645 genes were down-regulated (log FC < −1.5 and normalized P < 0.01, FDR = 0.05) in MC2 relative to MC1 and MC3 subtypes. Go analysis of the up-regulated genes revealed that “extracellular matrix organization,” “vasculature development,” “epithelial cell proliferation,” “cell-substrate adhesion,” and “cellular response to tumor necrosis factor” were enriched in biological processes and pathways which might be highly correlated with malignant progression of gliomas. The top 20 significantly enriched biological processes were shown in the figure (Figure 3A). The KEGG pathway analysis further revealed that these genes were also notably associated with tumor-relevant signaling pathways, for example, ECM-receptor interaction, Jak-STAT signaling pathway, and P53 signaling pathway, amongst others (Figure 3B).

Moreover, gene set enrichment analysis (GSEA) was performed to investigate the hallmarks of tumors in the MC2 subtype. The results indicated that the tumor hallmarks, such as P53 pathway (NES = 1.720, P = 0.013), P13K/AKT/mTOR signaling (NES = 1.811, P = 0.008), DNA repair (NES = 2.050, P < 0.001), and MTORC1 signaling (NES = 1.894, P = 0.006) (Figures 3C,D) were enriched in MC2 subtype. Combined with the above analysis, the subtypes identified by RNA:m⁵C methyltransferases were significantly associated with the malignant progression of gliomas.

Based on the relationship between RNA:m⁵C methyltransferases and the malignant progression of gliomas, we further attempted to explore the prognostic role of RNA:m⁵C methyltransferases in gliomas by a univariate survival analysis using Cox proportional hazards models of expression levels in the CGGA dataset, which was set up as a training dataset. We obtained six genes associated with prognosis (P < 0.01), among which NOP2, NSUN2, NSUN4, NSUN5, and NSUN7 acted as risk factors (HR > 1), and NSUN6 played a protective role (HR < 1) in gliomas (Figure 4A). To improve the robustness of the six RNA:m⁵C methyltransferases, these genes were selected to conduct an additional analysis by the least absolute shrinkage and selection operator (LASSO) Cox regression algorithm in the CGGA dataset (Figure 4B). Five genes of RNA:m⁵C methyltransferase genes and coefficients (Figure 4C) were screened to construct the risk score signature, and the formula for the risk score is as follows: 0.884 × (expression value of NOP2) + 1.167 × (expression value of NSUN4) + 0.190 × (expression value of NSUN5) + (−0.161) × (expression value of NSUN6) + 0.014 × (expression value of NSUN7) both in the training dataset (CGGA) and the verification dataset (TCGA). In this risk score signature, four genes (NOP2, NSUN4, NSUN5, and NSUN7) were cancer-promoting and NSUN6 was a cancer-suppressor gene. To better understand the role of these five prognostic genes in glioma, K-M analysis was performed both in the CGGA dataset and TCGA dataset in which samples were classified by high or low expression according to the median gene expression level. All five of these RNA:m⁵C methyltransferase genes were significantly correlated with OS (P < 0.0001) (Figure 4D and Supplementary Figure 3A).

To obtain the predictive effect of the risk score signature for clinical outcomes in patients with gliomas, the median of all patients' scores was used as a standard, and the data were divided into high- and low-risk groups both in the CGGA and TCGA datasets. The analyses indicated that the number of patients who died increased significantly as the risk score increased (Figures 5A,D). In addition, there was a significant difference in overall survival (P < 0.0001) between the high-risk group and low-risk group (Figures 5B,E). Thereafter, the ROC curve analyses at 1-, 3-, and 5-years for prognostic risk scores were performed to test the predictive efficiency of the risk signature. The results showed that the risk score had high accuracy [the area under the curve (AUC) of all results in the ROC curve were >0.750] in distinguishing the OS of gliomas (Figures 5C,F). Moreover, we performed further studies on the prognostic value of this signature in glioma patients stratified by WHO grade and IDH-mutant status. The results revealed that the risk score signature could be used to divide glioma patients, in the CGGA dataset, into two distinct prognosis groups with different WHO grade subtypes (LGG and GBM) and IDH-mutant status subtypes (Figures 5G–J). The consensus results of the risk score signature were also obtained using the TCGA dataset (Supplementary Figures 3B–E). From the comprehensive analyses above, we concluded that the prognostic efficacy of the risk score was accurate and stable.

The expression of the five screened RNA:m⁵C methyltransferases in low- and high-risk patients in the CGGA dataset are represented by heatmaps (Figure 6A). We found statistically significant differences between the low-risk and high-risk groups both in the CGGA and TCGA datasets, based on WHO grade (P < 0.0001), histology (P < 0.0001), IDH status (P < 0.0001), 1p/19q status (P < 0.0001), age (P < 0.0001), and receipt of additional chemotherapy (P < 0.0001) (Supplementary Table 3). Thereafter, we quantitatively analyzed the risk scores in glioma to investigate the association between the risk scores and each clinicopathological characteristic. As shown in the figure, the risk scores were significantly different in these groups with the CGGA dataset (Figure 6B) compartmentalized by WHO grade, IDH status, 1p/19q status, histology, age and receipt of additional chemotherapy, and were verified in the TCGA dataset (Supplementary Figures 4A–G). The results demonstrated that the risk scores identified by five RNA:m⁵C methyltransferases were significantly correlated with malignancy of glioma. Moreover, considering the importance of reported glioma-associated driver-gene alterations in glioma initiation and progression, including ATRX, P53 pathway (TP53, MDM2, MDM4), RB pathway (CDK4, CDK6, CCND2, CDKNA/B, RB1), and P13K/RTK pathway (PIK3CA, PIK3R1, PTEN, EGFR, PDGFRA, NF1), we obtained the somatic mutation data from the TCGA database (32, 33). The mutational landscape of tumor driver-gene alterations between low- and high-risk patients with gliomas was rendered as a waterfall plot (Figure 6C). The high-risk patients were characterized by more frequent alterations in tumor driver-genes than low-risk patients. This means high-risk patients harbored more progressive cancer. Combining the above results, the risk scores can predict not only overall survival but also clinicopathological features.

Next, we investigated whether this risk score was an independent prognostic factor based on eight clinicopathological features. Univariate and multivariate Cox regression analyses were performed with the CGGA dataset. We observed that the risk score, age, WHO grade, IDH status, 1p/19q codel status, chemotherapy status, and radiotherapy status were significantly correlated with prognosis using the univariate analysis (Figure 7A). Multivariate analysis based on the above factors was performed and the risk score remained strongly associated with the OS (P < 0.001, Figure 7B). In the verification dataset (TCGA) we obtained similar results, including the same factors in the multivariate analysis; the risk score also remained strongly associated with the OS (P = 0.002, Figures 7C,D). The consensus results demonstrated that the risk score constructed by RNA:m⁵C methyltransferases was a powerful and independent prognostic factor for glioma outcome.

To assess the genes that are significant for construction of the risk score signature, we performed the qRT-PCR assay and acquired the immunohistochemistry pathological specimen data from the website of The Human Protein Atlas. As shown in the figure, all five RNA:m⁵C methyltransferase genes, NOP2, NSUN4, NSUN5, NSUN6, and NSUN7, were differentially expressed between normal brain tissues, low grade and high grade glioma tissues (Figures 8A–E). Furthermore, as shown in the figure, we found that the protein expression of NOP2 and NSUN4 in high-grade glioma tissues were much higher than those in low-grade glioma tissues (Figures 8F,G). However, the protein expression pattern of NSUN5, NSUN6, and NSUN7 was not presented here due to lack of immunohistochemistry images.