The mRNA expression dataset and clinicopathological characteristic of OSCC patients were extracted from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) and the NCBI GEO database (https://www.ncbi.nlm.nih.gov/geo/). In total, 558 OSCC samples, 32 normal samples, and 101 oral leukoplakia samples were screened out for further evaluation. 322 OSCC samples and 32 normal control samples were extracted from TCGA-HNSCC cohort, the cancerous sites included alveolar ridge, hard palate, floor of the mouth, the base of the tongue, oral cavity, buccal mucosa, and oral tongue. GSE65858 includes 83 OSCC patients; GSE41613 contains 97 OSCC samples; GSE26549 includes 86 oral leukoplakia patients; GSE31056 includes 22 OSCC patients; and GSE85195 contains 34 OSCC patients and 15 oral leukoplakia patients. Detailed information on these datasets is provided in Supplementary Table S1. The clinicopathological information of OSCC patients is provided in Supplementary Table S2. The data processing were performed as we previously described (Gao et al., 2021). To remove the batch effects, the RNA sequencing data (FPKM format) data of OSCC patients in TCGA HNSCC was transformed into transcripts per kilobase million (TPM) format and then combined with microarray data in GSE65858. Copy number variations (CNVs) were obtained from the UCSC Public Hub and analyzed using the RCircos package in R. Somatic mutation expression data were downloaded from TCGA and analyzed using the maftools package in R.

We conducted consensus clustering analysis to classify OSCC patients into distinct m5C-modification clusters using the ConsensusClusterPlus R package. In total, 16 m5C methylation regulators, including “writers” (DNMT1, DNMT3A, DNMT3B, NOP2, NSUN2-7, and TRDMT1), “erasers” (ALYREF and YBX1), and “readers” (TET2, TET3 and ALKBH1) retrieved from literature were analyzed in this study (Chen et al., 2021c) (Trixl and Lusser, 2019) (Bohnsack et al., 2019). Clustering heatmaps were plotted using the pheatmap package.

We conducted single-sample gene set enrichment analysis (ssGSEA) to quantify the relative amount of 27 immune cell types in the TIME of OSCC using the GSVA and GSEABase packages. These gene sets acquired from recent studies were provided in Supplementary Table S3 (Charoentong et al., 2017) (He et al., 2018). ESTIMATE method was used to estimate the proportions of stromal cells (SCs), immune cells (ICs) in the TIME of OSCC. Tumor purity was predicted based on the SCs, ICs, and ESTIMATE scores. The deconvolution approach CIBERSORT was used to investigate the differences in 22 immune cell subtypes among the m5C-modification clusters. The “signature matrix” of 547 genes was acquired from CIBERSORT (https://cibersort.stanford.edu/index.php), namely, LM22.txt. The differences in physiological processes in different m5C-modification clusters were evaluated by gene set variation analysis (GSVA). Gene sets “c2.cp.kegg.v7.4.symbols” acquired from the MSigDB (http://www.gsea-msigdb.org/gsea/msigdb) were used for the GSVA analysis.

We used WGCNA to explore m5C-modification cluster-related modules. Based on the top 50% of the most variable genes (7,949 genes), the gene co-expression network was constructed. The correlation of the expression of these filtered genes was calculated, and β = 6 was selected as optimal soft threshold for building a scale-free topological network. Topological overlap matrix (TOM) analysis was performed to cluster the adjacency matrix of gene expression. The modules were identified with the criterion set as the mini-size of module gene numbers was 40, and a cut height of 0.9. The module with the strongest relationship with m5C-modification clusters was selected for further analysis.

M5Cscore was adopted for the quantitative analysis of m5C-modification patterns in an individual patient with OSCC. This method was established as we previously described (Gao et al., 2021). Briefly, the prognostic value of each gene in the selected module was evaluated by Cox regression analysis. Subsequently, principal component analysis (PCA) was conducted to calculate the principal components (PCs) 1 and 2, which were used for m5Cscore calculations. m5Cscore = ∑ ( PC 1 i + P C 2 i ) (1) where _i is the expression of the selected genes.

All statistical analyses were implemented using R software version3.5.3 (https://www.r-project.org/). The Kaplan–Meier analysis and log-rank analysis were performed to evaluate the progression-free survival (PFS) and overall survival (OS) between clusters. We conducted Pearman correlation analysis to explore the correlation between m5Cscore and known gene signatures and the correlation coefficient of gene expression. Biological gene signatures were retrieved from previous studies (Couzin-Frankel, 2013) (Melero et al., 2021). Student’s t-test was used to evaluate the statistical differences between two groups. Kruskal-Wallis tests and one-way ANOVA tests were performed when there were more than two groups. In all cases, the p-value was two-sided, and a p-value < of 0.05, was considered statistically significant.

In this study, sixteen m5C RNA methylation regulators, including “writers” (DNMT1, DNMT3A, DNMT3B, NOP2, NSUN2-7, and TRDMT1), “erasers” (ALYREF and YBX1), and “readers” (TET2, TET3 and ALKBH1) were analyzed to analyze the expression and prognosis significance of m5C RNA methylation regulators in OSCC. Between OSCC samples and normal controls, the expression levels of m5C regulators were significantly different, except for TET2 (Figure 1A). The expression of NSUN7 was markedly downregulated in OSCC samples, while that of other regulators was significantly upregulated. Among these m5C regulators, CNV alterations were relatively prevalent. NSUN3, NSUN2, DNMT3B, NOP2, YBX1, ALKBH1, NSUN4, and TET3 were predominantly amplified, while ALYREF had a higher frequency of CNV deletion (Figures 1B,C). The frequency of somatic mutations in DNMT3B, DNMT1, DNMT3A, TET3, NSUN7, NOP2, NSUN6, TET2, and ALYREF was >1%. However, among the 316 OSCC samples, only 27 samples (8.54%) showed m5C regulators mutation (Figure 1D). The aforementioned findings suggest that there is a marked difference between OSCC and normal controls in the expression and genetic variations in m5C regulators.

We also evaluated the prognostic value of m5C regulators in OSCC (Figures 1E,F). Most m5C regulators showed tumor-promoting effects, except for TET2 and NSUN7. The high expression of nine m5C regulators suggested poorer survival outcomes in OSCC patients: NOP2 (p = 0.034), NSUN3 (p = 0.007), NSUN4 (p < 0.001), NSUN5 (p = 0.019), NSUN6 (p = 0.002), TRDMT1 (p = 0.013), ALYREF (p = 0.010), ALKBH1 (p = 0.004) and YBX1 (p = 0.010) (Figures 1G–O). In the validation cohort GSE41613, the high expression of seven m5C regulators suggested worse outcomes, namely NSUN2 (p = 0.010), NSUN3 (p = 0.046), TRDMT1 (p = 0.032), DNMT3A (p = 0.046), DNMT3B (p < 0.001), ALKBH1 (p < 0.001), and ALYREF (p < 0.001); however, high expression of TET2 (p = 0.011) was related to better prognosis (Supplementary Figure S1). These findings indicate that crosstalk among m5C regulators might affect the formation and effect of distinct m5C modification clusters and might further affect the oncogenesis, development, and prognosis of OSCC.

To comprehensively evaluate the association among these m5C regulators, we conducted consensus clustering analysis to stratify 405 OSCC patients into two independent clusters: clusters A and B (Figure 2A; Supplementary Figures S2 A–C). Figure 2B shows that the expression of m5C regulators in the two clusters was significantly distinct. The expression of m5C regulators was noticeably higher in cluster B. The OS analysis indicated that compared with cluster A, prognosis of patients in cluster B was markedly poorer (p = 0.020) (Figure 2C). What’s more, in GSE41613, OSCC patients were also classified into two clusters (Supplementary Figures S2D–F). Most m5C regulators were highly expressed in cluster A (Supplementary Figure S2G). The OS outcomes in cluster A were worse than those in cluster B (p = 0.047) (Figure 2D).

We further evaluated the effect of m5C modification patterns on the TIME in OSCC. CIBERSORT analysis results indicated that compared with cluster A, cluster B had markedly lower fractions of resting dendritic cells (DCs), activated natural killer (NK) cells, neutrophils, and CD8⁺ T cells, and markedly higher fractions of resting memory CD4⁺ T cells and M0 macrophages (Figure 2E). Furthermore, ssGSEA analysis showed compared with cluster A, cluster B had lower immune activity in most immune cell types, included CD8⁺ T cells, tumor-infiltrating lymphocytes (TILs), and B cells, which indicated a deeper immunosuppressed state in cluster B (Figure 2F). GSVA enrichment analysis revealed that cluster B was significantly enriched in ubiquitin-mediated proteolysis, mismatch repair-related pathways, progesterone-mediated oocyte maturation, and cell-cycle. Cluster A was markedly enriched in cytokine-cytokine receptor interaction, neuroactive ligand receptor interaction, and PPAR signaling pathway-related pathways (Figure 2G). These findings imply that m5C-related modification patterns might participate in TIME regulation in OSCC.

To identify the m5C-modification cluster-related modules, we constructed a gene co-expression network (Figure 3A). β = 6 was selected as the optimal soft threshold for building a scale-free topological network (Supplementary Figure S3A). As illustrated in Figure 3B; 11 modules were identified. The yellow module demonstrated the highest correlation with the m5C cluster (r = 0.64, p = 3e-48). In addition, genes in the yellow module were markedly co-expressed (weighted correlation = 0.83, P < 1e-200) (Figure 3C). 53 genes in this module were screened out for further analysis, with the screening criteria of |MM| ≥ 0.8 and |GS| ≥ 0.2 (Supplementary Table S4).

Next, we conducted a consensus clustering analysis based on m5C cluster-related genes. OSCC patients were divided into two m5C gene clusters: gene clusters A and B (Figure 3D; Supplementary Figures S4B–D). Cluster B had the highest expression of most m5C regulators and poorer prognosis (p = 0.032) (Figures 3E,F). We verified and obtained the same outcomes in the validation cohort (p < 0.001) (Figure 3G; Supplementary Figures S3E–H).

We calculated the m5Cscore for individual OSCC patients using PCA. Patients with high-m5Cscore showed poorer OS than those with low-m5Cscore (Figures 4A,B). A consistent prognosis result was verified using GSE41613 (Figure 4C). The PFS of patients in high-m5cscore group was also poorer than those in low-m5Cscore group in GSE65858 (Figure 4D). We further evaluated the m5Cscore between patients in different tumor stages, nodal stages, and clinical stages. As shown in Figures 4E–G, m5Cscore progressively increased with the progression of the T stage (Figure 4E), N stage (Figure 4F), and clinical stage (Figure 4G). At different T and N stages, OSCC patients with a high m5Cscore showed poorer OS (Supplementary Figures S4 A–D). m5Cscore and m5C clusters and m5C gene clusters had a positive correlation, i.e., m5Cscore were higher in m5C cluster B and m5C gene cluster B (Supplementary Figures S4 E–F). In addition, compared with low-m5Cscore group, the SCs, ICs, and ESTIMATE scores in the high-m5Cscore group were noticeably lower, while the tumor purity score was noticeably higher (Figures 4H–K). What’s more, the immune activity of most types of immune cells in the high-m5Cscore group were significantly lower, compared with low-m5Cscore group (Figure 4L), which further suggests an association between a suppressive TIME and high m5Cscore in OSCC patients.

To verify the correlation between tumor mutational burden (TMB) and m5Cscores, we performed TMB quantification analysis in the high- and low-m5Cscore groups. In both m5Cscore groups, the somatic mutation rates of TP53, TTN, CDKN2A, NOTCH1, and CASP8 were relatively higher. Generally, the proportion and extent of somatic mutations in the high-m5Cscore group were higher and more extensive (Figures 5A,B). In addition, there was a high correlation between high TMB and high m5Cscore (p = 0.0028) (Figures 5C,D), as well as a worse prognosis (Figure 5E). The combination of a high m5Cscore and high TMB suggested an extremely poor OS for patients with OSCC (Figure 5F). Besides, m5Cscore was positively correlated with antigen processing machinery, CD8⁺ T effector, and immune checkpoint, but negatively correlated with cell cycle, Wnt target, mismatch repair, and cell cycle regulator by Spearman correlation analysis (Figure 5G).

We further investigated other predictive roles of m5Cscore in OSCC progression. Tumor recurrence is a crucial prognostic factor for OSCC. In this study, we compared the m5C regulators expression and m5Cscores between OSCC patients with and without recurrence in GSE31056. As illustrated in Figure 6A, there were five regulators—NOP2, DNMT1, DNMT3A, DNMT3B, and TET3—that were upregulated in OSCC patients with recurrence. Surprisingly, for patients with recurrence, m5Cscore was markedly higher than those without recurrence (Figure 6B). Moreover, 100% of patients in the high-m5Cscore group experienced recurrence in the course of OSCC (Figure 6C). We then analyzed the m5Cscore in patients with oral leukoplakia (Figures 6D–F). The m5Cscore was also markedly higher in patients with oral leukoplakia who developed OSCC. Three m5C regulators—NOP2, DNMT3B, and ALYREF—were upregulated (Figure 6D). We also verified the expression levels of m5C regulators in the microarray dataset GSE85195, which contained 34 OSCC patients and 15 oral leukoplakia patients. The expression of eight m5C regulators was upregulated in OSCC patients compared with that in oral leukoplakia patients (Figure 6G). These findings suggest that m5C methylation regulators might affect OSCC progression and that m5Cscore might has a predictive effect on the recurrence of OSCC and the deterioration of oral leukoplakia.