The mRNA expression profile of 105 normal brain tissues of GTEx was obtained from UCSC Xena (https://xena.ucsc.edu/). Transcriptomic expression profiles (FPKM normalized) and corresponding clinicopathological data (including IDH status, 1p/19q status, O(6)-methylguanine DNA methyltransferase (MGMT) promoter status, WHO classification, age, gender, and survival information) were retrieved from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). Samples without survival information were excluded, and 504 LGG samples were retained for further analysis. In addition, the expression profile of 513 LGG samples as well as corresponding demographics (Table 1) were downloaded from the Chinese Glioma Genome Atlas (CGGA) database (http://www.cgga.org.cn/). Datasets were subjected to the R package ‘limma’ for removing batch effect. The mRNA-seq datasets (FPKM normalized) of 33 types of cancers were also retrieved from the TCGA GDC project in UCSC Xena data portal. LGG somatic mutation data in UCSC Xena were analyzed and visualized using the R package ‘maftools’ (19).

13 confirmed m5C regulatory factors were selected from the available literature, including DNMT1, DNMT3A, DNMT3B, ALYREF, NOP2, NSUN2, NSUN3, NSUN4, NSUN5, NSUN6, NSUN7, TET2, TRDMT1 (20, 21). The R package ‘limma’ was used to integrate and remove batch effects from TCGA and GTEx data, as well as to analyse differences in 13 regulators. The p value < 0.05 was considered significant.

The Search Tool for Interacting Genes (STRING) database (https://string-db.org/) was used to construct PPI interaction networks based on methods including scientific text mining, extraction of experimental data, and calculation of genomic features (22). The potential interactions between 13 m5C regulators were explored by setting the minimum required interaction scores on the STRING online website to a high confidence level (> 0.7). The cytoscape software (version 3.7.2) was used to visualize the PPI network. Furthermore, hub genes were identified using maximal clique centrality (MCC) computing method with the cytohubba plugin (23).

The Genome Reference Consortium Human Build 38 (GRCh38) annotation file was obtained from the GENCODE website (https://www.gencodegenes.org/human/) to annotate the TCGA and CGGA datasets, and 14,142 and 1,004 lncRNAs were annotated, respectively. Pearson correlation analysis was used to screen out 222 m5C-related lncRNAs from TCGA and 81 m5C-related lncRNAs from CGGA, based on 13 m5C regulatory factors. Absolute value of correlation coefficient > 0.5 and p value < 0.01 were set as cutoff.

Firstly, all lncRNAs associated with m5C were screened by univariate Cox regression analysis, and 107 and 46 lncRNAs with prognostic value were obtained in the TCGA and CGGA datasets, respectively. The least absolute shrinkage and selection operator (LASSO) regression analysis was performed on six lncRNAs co-expressed in both datasets using the package ‘glmnet’ to screen out biomarkers for further development m5CrLS. The following equation was used to calculate the m5CrLS score for LGG:

The samples were split into high and low m5CrLS score groups by the median value.

Differentially expressed genes (DEGs) between low and high m5CrLS score groups were calculated using the R package ‘limma’ based on the TCGA dataset with the criteria of | log2(FC) | > 1, p < 0.05, and FDR < 0.05. DEGs up-regulated [log2(FC) > 0] in the high m5CrLS score group were selected as candidates for Gene Ontology (GO) analysis, including cellular component (CC), molecular function (MF), biological process (BP), and Kyoto Encyclopedia of Genes and Genomes (KEGG) based on the R package ‘clusterProfiler’ and ‘enrichplot’ for visualization.

The relative activation of classical tumor-associated pathways in high and low m5CrLS score groups was characterized using single-sample gene set enrichment analysis (ssGSEA). The signature genes used for the calculations were derived from a summary of recent researches (24, 25). In the analysis, ssGSEA scores were normalized to a unit distribution, where 0 was the minimum value of activation for each pathway and 1 represented the maximum.

‘ESTIMATE’ algorithm was used to assess the content of immune and stromal cells in the TME of each sample, which was expressed in the form of four scores: immune score, matrix score, estimated score and tumor purity (26).

The fraction of B-cell, CD4⁺ T-cell, CD8⁺ T-cell, dendritic cell, macrophage, and neutrophil infiltration was estimated using the TIMER2.0 webtool (http://timer.cistrome.org/) based on mRNA expression profiles (27). The R package ‘limma’ was used for differential analysis and Spearman correlation analysis was performed to calculate the correlation between the m5CrLS score and the fraction of 6 immune infiltrations.

The ssGSEA algorithm was performed to further quantify the infiltration of 28 immune cells in the TME. The characteristic genes for each type of immune cell were summarized by JIA et al. (28). The ssGSEA scores distributed from 0 and 1 represented the minimum and maximum of each immune cell infiltration abundance, respectively.

Tumor Immunophenotype Profiling (TIP) (http://biocc.hrbmu.edu.cn/TIP/) is a web-based tool that allows convenient and rapid analysis and visualization of the extent of tumor-infiltrating immune cells and the 7 phase events of anti-cancer immunity, including tumor cell antigen release (step 1), cancer antigen presentation (step 2), stimulation and activation (step 3), immune cell transfer to the tumor (step 4), immune cell infiltration (step 5), T cell recognition of cancer cells and killing of cancer cells (steps 6, 7) (29). In view of this, the antitumor immune function activity of LGG was quantified by TIP and further analyzed for differences between high and low m5CrLS score groups.

The Tumor Immune Dysfunction and Exclusion (TIDE) (http://tide.dfci.harvard.edu) web platform integrates published immune checkpoint blockade (ICB) trials to predict clinical response to ICB (30). Both datasets were assessed by the TIDE algorithm and samples with TIDE scores below a threshold (default value of 0) were set to be responsive to ICB treatment.

All statistical analysis was performed on R (version 4.1.1). Kaplan-Meier (K-M) analyses and log-rank statistical tests were used to compare overall survival (OS) between groups. Survival analysis grouping based on 13 m5C regulators and four lncRNAs was divided according to the median of each gene expression. In addition, the median value of m5CrLS scores from TCGA dataset was used as node for survival analysis between m5CrLS groups as well as determination of sensitivity to TMZ and radiotherapy treatment. Receiver operating characteristic (ROC) curves and the corresponding area under the curve (AUC) values were used to evaluate the predictive power of m5CrLS scores established based on the expression of four biomarkers for prognosis. Univariate and multifactorial Cox regression analyses were used to assess the independent prognostic value of the m5CrLS and 13 m5C regulators. Wilcoxon test was performed to compare intra-subgroup differences in clinicopathological subtype, immune cell infiltration, ssGSEA score, the expression of immune checkpoint, and the TIP score. Kruskal test was used to compare differences between histopathological subtypes (astrocytoma, oligodendroglioma, oligoastrocytoma). A chi-square test was used to compare gender, age, WHO grade, IDH status, 1p/19q status, O(6)-methylguanine DNA methyltransferase (MGMT) promoter methylation status, and distribution characteristics of histopathological subtypes between low and high m5CrLS score groups. The difference in response to ICB treatment between low and high m5CrLS score groups was assessed using the fisher test.

To investigate the expression of 13 m5C regulatory factors in LGG, normal and tumor samples from the GTEx and TCGA databases were integrated. As a result, DNMT3A, DNMT3B, ALYREF, NSUN3, NSUN6, NSUN7, and TET2 were significantly up-regulated in tumor samples (p < 0.05), and NOP2, NSUN2, NSUN4, NSUN5, and DNMT1 regulators were significantly down-regulated (p < 0.05), while TRDMT1 was the only regulator without significant change in expression (p = 0.276) (Figure 1A). Then, K-M analysis was performed to explore the expression of these m5C regulators on LGG survival. The results showed that patients with decreased expression of DNMT3A, DNMT3B, DNMT1, NOP2, NSUN4, and NSUN7 had improved OS, which was also validated by the CGGA dataset (Figures S1A, B). The univariate and multivariate Cox regression analysis showed that NSUN4 and NSUN7 were independent risk factors for LGG based on both datasets (Figures 1B, S1C).

To further characterize the interactions between these m5C regulators, a protein-protein interaction (PPI) network was constructed (Figure 1C). As a result, TRDMT1 scored highest and was the hub gene in the network, and genes with their corresponding centrality weights were collated as Table S1. Subsequently, correlation analysis found that TET2 had the highest negative correlation with NSUN5 (cor = -0.44, p < 0.05) and the highest positive correlation with NSUN3 (cor = 0.68, p < 0.05) (Figure S1D). Overall, these results demonstrated the LGG-related m5C regulators and their interactions.

To further explore the impact of lncRNAs on the m5C regulatory patterns of LGG, we firstly conducted the Pearson correlation analysis between the expression of lncRNAs and m5C regulators. With | cor | > 0.5 and p < 0.01 as the threshold, a total of 222 and 81 m5C-related lncRNAs were identified in TCGA and CGGA samples. With univariate Cox analysis, 107 and 46 m5C-related lncRNAs retained prognostic significance in TCGA and CGGA samples (Tables S2, S3). We noted that LINC00265, CIRBP-AS1, GDNF-AS1, ZBTB20-AS4, NNT-AS1 and TRAF3IP2-AS1 were all identified as m5C-related lncRNAs in both datasets (Figure 1D). Also, the correlations between the expression of these six lncRNAs and m5C regulators were strong (cor > 0.5) (Figures 1E, S1E and Tables 2, S4).

To identify biomarkers, we identified four of the six lncRNAs with significant prognostic significance by LASSO regression analysis, namely LINC00265, CIRBP-AS1, GDNF-AS1 and ZBTB20-AS4 (Figure S2A). Particularly, GDNF-AS1 and ZBTB20-AS4 were protective factors, while LINC00265 and CIRBP-AS1 were risk factors (Figures 2A, B). Taking into account the clinicopathological features of LGG, the expression of GDNF-AS1 and ZBTB20-AS4 was significantly increased in samples of 1p/19q co-deletion and IDH mutation in both datasets (p < 0.001), and the expression of CIRBP-AS1was decreased in samples with the IDH mutation (Figures 2C, S2B). K-M curves showed that LGG with increased expression of GDNF-AS1 and ZBTB20-AS4 had improved survival, whereas patients with elevated expression of LINC00265 and CIRBP-AS1 had poor OS, in line with the result of Cox regression analysis (Figures 2D, S2C). In addition, we analyzed the expression of the four m5C-related lncRNAs in 33 tumors using a pan-cancer analysis. Compared to normal tissues, ZBTB20-AS4, CIRBP-AS1, GDNF-AS1, and LINC00265 were differentially expressed in a variety of tumors, including BRCA, COAD, KIRC, KIRP, LUAD, LUSC, and UCEC, suggesting that these m5C-related lncRNAs may be conserved in the tumor progression (Figure 2E). In summary, these results suggested that these four m5C-related lncRNAs were potential candidates for the development of a reliable prognostic model.

Next, we constructed the m5CrLS using the four m5c-related lncRNAs and their LASSO regression coefficients to quantify individual differences among LGG samples (Figure S3A). The distribution of m5CrLS score with survival time and status showed that the mortality increased with the score (Figures 3A, S3B). K-M curves suggested that the low m5CrLS score group had a significantly improved OS (p < 0.001) (Figure 3B). Notably, the m5CrLS group predicted prognosis was independent of age, WHO tumor grade, MGMT promoter status, and IDH status (Figures 3C, S3C, D). The ROC curves showed a robust time-dependent predictive power of the m5CrLS, with the 1-year AUC values of 0.854 and 0.740, the 3-year AUCs of 0.799 and 0.724, and the 5-year AUCs of 0.730 and 0.685, in the TCGA and CGGA datasets, respectively which was an improved predictor than clinical characteristics including age, gender, grade, IDH status and 1p/19q status (Figures 3D, S3E). The results of univariate [TCGA hazard ratio (HR) = 3.864, 95% confidence interval (CI) =2.825 - 5.286, p < 0.001; CGGA HR = 4.151, 95% CI = 2.881 - 5.980, p < 0.001] and multivariate Cox analysis (TCGA HR = 1.885, 95% CI = 1.178 – 3.015, p = 0.008; CGGA HR = 2.316, 95% CI = 1.395 - 3.844, p = 0.001) showed that the m5CrLS score was an independent risk factor (Figures 4A, S4A). The accuracy of the prognostic model can be determined by the calibration degree curve (Figure 4B). Moreover, to render the m5CrLS score a clinically applicable quantitative standard, a nomogram was constructed (Figure 4C).

In addition, the differential expression of the four prognostic m5C-related lncRNAs and their associated clinical features were exhibited. And samples of WHO grade II, 1p/19q co-deletion, IDH mutation, and MGMT promoter methylation were enriched in the low m5CrLS score group (p < 0.05) (Figures S4B, C). Particularly, WHO grade II, 1p/19q co-deletion, IDH mutation, MGMT promoter methylation, and oligodendroglioma had decreased m5CrLS scores (Figures 4D, S5A). The results of principal component analysis (PCA) in high and low m5CrLS score groups showed that those had different distributions based on the expression of prognostic biomarkers, all m5C-related lncRNAs, 13 m5C regulators, and all genes, suggesting significant differences in molecular and m5C methylation characteristics between patients in the two groups (Figure S5B).

To further uncover the molecular underpinning of m5CrLS stratification, functional enrichment analysis was performed. Firstly, 649 differentially expressed genes were identified and 520 genes were significantly up-regulated in high m5CrLS score LGG (Table S5). GO and KEGG enrichment analysis revealed that 520 up-regulated genes were mainly enriched in the extracellular matrix organization, extracellular structural organization, ECM-receptor interaction, and PI3K-AKT signaling pathway. Notably, high and low m5CrLS score patients have multiple immune response processes in different states of activation, such as humoral immune response, MHC class II protein complex (Figures 5A, B). In addition, we selected 20 classical tumor-associated pathways whose activation tended to correlate with the degree of tumor malignancy (24, 25). As a result, the high m5CrLS score group showed an increased ssGSEA score of angiogenesis, cell cycle, epithelial-mesenchymal transition (EMT), pan-fibroblast transforming growth factor-beta (Pan-F TBRS). In contrast, the low m5CrLS score group had increased ssGSEA score in the DNA damage repair, WNT target pathway (Figures 5C, D). The above results may suggest that tumor cells in the high m5CrLS score group were more malignant and aggressive.

Given that the differences in immune-related biological processes and signaling pathways between high and low m5CrLS score groups, the immune microenvironment was further explored. The stromal score, immune score, estimate score, and tumor purity were calculated based on the R package ‘ESTIMATE’. As a result, the high m5CrLS score group consistently scored higher, indicating a possible increase in the immune and stromal infiltration (Figures 6A, S6A). Thus, the fraction of six immune cells, including B cells, CD4⁺ T cells, CD8⁺ T cells, neutrophils, macrophages, and myeloid dendritic cells were inferred by TIMER2.0. Consistently, the infiltration of these immune cells was significantly increased in the high m5CrLS score group (Figure 6B). Spearman correlation analysis suggested that myeloid dendritic cells had the highest positive correlation with m5CrLS score (cor = 0.573, p < 0.001). Furthermore, myeloid dendritic cells showed a strong positive correlation with immune infiltration of neutrophils (cor = 0.86) (Figures 6C, S6B, C). Chemokines CCR1, CCR2, CCR3, CCR6, CCR5, CCL5, CXCL10, and XCL1 showed high expression in the high score group of m5CrLS to promote the migration of dendritic cells. In addition, the high expression of CCR7, CCL21 inhibited the function of dendritic cells and hindered the killing effect of T cells on tumor cells (31–33). The cytokine IL10, which showed high expression in the high m5CrLS score, could inhibit the release of IL12A from dendritic cells, thus causing anti-tumor immunosuppressive effects (34). In addition, IL6 and IL15 showed high expression in the m5CrLS high subgroup (Figures 6D, S6D). For validation, immune infiltration was also estimated using the ssGSEA algorithm and similar results were yielded that 24 out of 28 immune cells scored higher in high m5CrLS score patients (Figures 6E, S6E). To further explore the ability of anti-tumor immune response in high and low m5CrLS score groups, we quantified the activity of 7 steps of LGG anti-tumor immune response based on TIP. The results showed that the low m5CrLS score group scored higher in step 3, step 6, and step 7, while the high m5CrLS score group scored higher in the release of step 1, and step 5 (Figure 6F). The above results suggest that the stratification analysis based on m5CrLS can identify the antitumor immune activity of LGG to some extent.

Given the impact of immune checkpoints (ICPs) on anti-tumor immunity, we investigated the expression of those ICPs. Tumor necrosis factor superfamily ligands (TNFSF) and receptors (TNFRSF), including TNFSF9, TNFSF4, TNFSF15, TNFSF14, CD40, CD70, CD40LG, TNFRSF4, TNFRSF8, TNFRSF9, TNFRSF14, TNFRSF18, and TNFRSF25 were up-regulated in the high m5CrLS score group in both datasets. Programmed cell death protein-1 (PD-1) and two binding ligands, PD-L1 and PD-L2, were highly expressed in the high m5CrLS score group, indicating T-cell depletion and decreased effector function (35). CTLA4 was highly expressed in the high m5CrLS score group and could bind to equally highly expressed CD86 and CD80 to transduce T-cell suppressor signals and promote the suppressive function of Treg cells (36). TIM-3 (HAVCR2) and LGALS9 were significantly associated with T-cell depletion and impaired function, and low expression in the low m5CrLS score group suggested a stronger autoimmune and antitumor immune response than in the high m5CrLS score group (37). In addition, high expression of ICOS, ICOSLG, and IDO1 in high m5CrLS score patients contributed to Treg cell aggregation (38, 39). high expression of LAIR1 suppressed CD8⁺ T cell activity in high m5CrLS score patients (40). High expression of CD48, CD244 (41), CD200R1 (42), and BTLA (43) also had a suppressive effect on anti-tumor immunity (Figures 6G, S6F).

The tumor immune microenvironment was also associated with somatic mutation rates (44). The analysis of LGG mutation information allowed us to better understand patients with different m5CrLS. The mutation rate in the low m5CrLS score group was as high as 100%, while the mutation rate in the high m5CrLS score group was 90.8%. This implies a high release of tumor neoantigens and more effective anti-cancer immune activity in the low m5CrLS score group with higher somatic mutation rates (45–47). In addition, the two groups differed in terms of high IDH1 mutations (94% in the low m5CrLS score group and 60% in the high m5CrLS score group). It has been noted that IDH1 mutations suppress the infiltration of immune cells such as macrophages, microglia, and neutrophils (48). This further explains the lower number of cancer-promoting immune cells in patients with low m5CrLS scores. Interestingly, PTEN was mutated at a frequency of 8% in the high m5CrLS score group but did not appear in the top 15 mutated genes of the low m5CrLS score group in terms of mutation frequency. Mutations in PTEN result in its loss of tumor suppressive function (49). Notably, increased PTEN mutations and excessive activation of PI3K-AKT were associated with ICB treatment resistance mechanisms (50) (Figure 7A).

ICB has achieved surprising results in the treatment of tumors. Previously available studies have highlighted the impact of RNA modification on ICB treatment. With this in mind, we employed the TIDE algorithm based on the stratification of m5C-related lncRNAs to assess the efficacy of LGG on ICB treatment. The results showed that there was a significant difference in the response to ICB treatment between the two groups (p < 0.05), and that low m5CrLS score patients were more sensitive to ICB treatment (Figure 7B).

Next, the association of TMZ efficacy with m5CrLS score groups was examined. K-M analysis showed that in LGG of grade II who received TMZ, OS was significantly higher in the low m5CrLS score group. In contrast, there was no statistical difference between the high and low m5CrLS score groups that did not receive TMZ (Figures 7C, E). Similar results were also yielded in the samples of grade III that the low m5CrLS score group predicted improved OS in patients who received TMZ (Figures 7D, F).

In addition, the association between the stratification analysis and the efficacy of radiotherapy was evaluated. In the TCGA dataset, low m5CrLS score patients of grade II had a significant survival benefit after receiving radiotherapy, and such result can almost be validated by the CGGA cohort (Figures S7A, C). There was no significant difference in OS for patients in the high and low m5CrLS score groups of grade II who did not receive radiotherapy. In grade III, both datasets showed improved OS for the low m5CrLS score group regardless of whether they received radiotherapy or not (Figures S7B, D).