The work flow chart in the present work is presented in Supplementary Figure 1 . The Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO) databases were searched to obtain the clinical annotation and related gene expression data. Patients who had no survival data were eliminated from this study. The eligible PTC cohorts [including GSE33630 (32), GSE65144 (33), and GSE29265, together with TCGA-PTC (The Cancer Genome Atlas- papillary thyroid carcinoma)] were collected into the present work. With regard to Affymetrix^® microarray data, raw “CEL” files were downloaded to adjust the background and normalize the quantile using the multiarray averaging approach by affy and simpleaffy packages. In terms of microarray data of additional sources, matrix files after normalization were collected directly. For TCGA datasets, the RNA sequencing information (FPKM values) of gene levels was obtained based on the Genomic Data Commons (GDC, https://portal.gdc.cancer.gov/) by TCGAbiolinks of R package, a software designed to comprehensively analyze GDC data (34). Thereafter, the FPKM values were converted to the transcripts per kilobase million (TPM) values. At the same time, the GSE65144 (12 tumor and 13 normal samples), GSE33630 (60 tumor along with 45 normal samples), and GSE29265 (29 tumor and 20 normal samples) datasets were also downloaded. R package (version 3.6.1) was utilized for data analysis.

Altogether, 13 regulators were obtained based on TCGA datasets to identify the diverse m⁵C regulator-mediated m⁵C modification patterns. All the 13 genes, except for ALYREF and NSUN1, were with available expression profiles. The remaining 11 m⁵C regulators contained one eraser (TET2) and 10 writers (NSUN2-7, DNMT1-2, and DNMT3A-3B). Of our 493 patients from TCGA-PTC, 9 among those 11 genes were differentially expressed between tumor and normal tissues (with the exception for NSUN3 and DNMT3A) ( Supplementary Table 1 and Supplementary Figure 2 ). Later, consensus clustering was adopted for identifying the different m⁵C modification patterns according to nine m⁵C regulator expression levels, and then patients were classified accordingly. The above procedure was performed using the ConsensuClusterPlus package (35) for 1000 iterations to guarantee the classification stability.

To investigate the heterogeneities in biological process among the m⁵C modification patterns, GSVA was carried out by the “GSVA” R package. Notably, GSVA is the unsupervised, non-parametric approach usually used to estimate variations of pathways and biological processes within the expression dataset samples (36). The “c2.cp.kegg.v7.0.symbols” gene sets were extracted based on MSigDB database to conduct GSVA. The adjusted p < 0.05 indicated statistical significance. Meanwhile, functional annotation was performed using WebGestaltR package (37), and the threshold was FDR < 0.05.

Estimate R package was utilized to calculate immune and stromal scores for all samples to reflect the immune and stromal cell infiltration degrees on the whole. Besides, CIBERSORT algorithm (38) was adopted for quantifying cell infiltration relative abundance within the TME of PTC. Thereafter, the gene set used to mark the TME-infiltrating immunocyte type was acquired to score different human immunocyte subtypes, like the activated CD8 T cells, regulatory T cells, natural killer T cells, activated dendritic cells (DCs), and macrophages (39).

To identify the m⁵C-associated genes, the patients were divided into three different m⁵C modification patterns according to nine m⁵C regulator expression levels. DEGs were determined across the diverse modification patterns using the empirical Bayesian method in the limma R package (40). The adjusted p < 0.05 served as the significance criterion to determine DEGS.

To quantify m⁵C modification patterns among individual tumors, the scoring system, m⁵C-score, was built based on the m⁵C gene signature, as shown below:

First of all, DEGs obtained based on the diverse m⁵C-clusters were subjected to normalization across all samples, and then, those overlapped genes were selected. Afterwards, all cases were divided into different groups via the unsupervised clustering approach, so as to analyze the overlapped DEGs. In addition, the gene cluster number and the stability were defined using the consensus clustering algorithm. Later, prognostic analysis was carried out for all genes selected in our constructed signature by the use of the univariate Cox regression model. Later, those significant genes were obtained in subsequent analysis. In this study, p < 0.01 was selected as the criterion to screen 49 genes. Supplementary Table 2 shows the results of single factor survival analysis for the 49 genes. Then, principal component analysis (PCA) was utilized for constructing the m⁵C signature. PC1 and PC2 were adopted as the signature scores; as a result, the score was focused on the set that had the greatest number of well-correlated (or anticorrelated) genes, and the contributions of genes not tracking with other members in the set were down-weighted. Later, the m⁵C-score was defined by the GGI-like approach (41):

m 5 C − score  =   ∑ ( PC 1 i +  PC 2 i )

In the formula, i represents the 49 m⁵C phenotype-associated gene expression levels.

First of all, this study mined the signature genes significantly correlated with the m⁵C-score (correlation coefficient > 0.4), and the PTC diagnostic model was constructed by the SVM method. Thereafter, the accuracy of this model was verified using samples from TCGA and GEO databases.

The Spearman and distance correlation analysis was adopted to calculate the correlation coefficients of TME-infiltrating immunocytes with m⁵C regulator expression levels. Thereafter, Kruskal–Wallis test and one-way ANOVA were applied in comparing the heterogeneities among three groups. Based on correlations of the m⁵C-score with patient survival, the survminer R package was utilized to determine the threshold value for every dataset. In addition, the “surv-cutpoint” function was used for dichotomizing the m⁵C-score by testing the possible threshold values to find the maximal rank statistic. Later, all cases were classified as high or low m⁵C-score group according to the maximal log-rank statistics for reducing the calculation batch effect. In the meantime, the log-rank test and Kaplan–Meier approach were adopted for identifying the significance of differences, so as to generate survival curves. The univariate Cox regression model was used for calculating hazard ratios (HRs) for the m⁵C regulators as well as the m⁵C phenotype-associated genes. Meanwhile, the receiver operating characteristic (ROC) curve was plotted to assess the sensitivity and specificity of our diagnostic model and the m⁵C-score, and pROC R package was utilized to quantify the area under the curve (AUC). The two-sided statistical p < 0.05 indicated statistical significance. The R 3.6.1 software was employed for data processing.

According to nine m⁵C regulators with expression profiles in the TCGA-PTC dataset, PTC samples were identified from normal samples ( Figure 1A ). Afterwards, the expression profile data of these nine m⁵C regulators were carried out with z-score standardization using the scale function in mosaic package. Then, three different m⁵C modification patterns were discovered according to those nine m⁵C regulator expression patterns ( Figures 1B, C ). These three patterns were named m⁵C-clusters 1–3. It was observed from Figure 1D that the expression level of these nine m⁵C regulators showed significant differences among the three distinct subtype samples.

Prognostic analysis was also carried out for these three major m⁵C modification patterns, which suggested that the m⁵C-cluster 2 modification pattern showed survival advantage ( Figure 1E ). However, due to the speciality of PTC and the good overall prognosis, there was no significant statistical difference among these three subtypes. Besides, average survival time of samples in these three subtypes was also analyzed, which discovered that the average survival time of C2 subtype samples was 1307.657 days, that of C1 subtype samples was 1125.877 days, and that of C3 subtype samples was 1202.695, with that in C2 higher than those in C1 and C3.

To explore those biological behaviors in the different m⁵C modification patterns, GSVA was conducted. It was illustrated from Figure 2A that m⁵C-cluster 1 significantly associated with the amino acid metabolic pathways; m⁵C-cluster 2 was enriched to the endocrine system, lipid metabolism, and cancer, whereas m⁵C-cluster 3 was associated with cell cycle, DNA repair, and nucleic acid metabolism.

Furthermore, the distribution of clinical features of samples in the above three subtypes was statistically analyzed. The statistical results are displayed in Supplementary Table 3 and Figure 3 . It was found from the results that, multiple clinical features in the three subtype samples were randomly distributed, with no significant difference.

In addition, the ESTIMATE algorithm was applied in quantifying the differences in stromal cell infiltration among the three subtype samples. As shown in Figure 2B , the stromal score in m⁵C-cluster 2 was the highest, followed by m⁵C-cluster 3, while m⁵C-cluster 1 had the lowest score. In addition, there were significant differences among them. Thereafter, the CIBERSORT deconvolution algorithm was utilized for comparing the heterogeneities in immunocyte components of the three m⁵C modification patterns ( Figure 2C ). Meanwhile, the support vector regression was used to determine the immunocyte types in tumors. As a result, high levels of Tregs and monocytes were detected in m⁵C-cluster 1 and m⁵C-cluster 3, whereas excessive resting/activated DCs were found in m⁵C-cluster 2. Recently, research has particularly focused on the RNA modification mechanism in the regulation of DC activation. DCs function to present antigen and to activate the naive T cells, which connect the intrinsic immunity with the adaptive one (42).

Finally, this study analyzed the expression of the 34 known immune checkpoints in the three subtype samples. As found from Figure 2D , there were significant differences in the expression of these 34 immune checkpoints among the three subtypes. Most immune checkpoint genes were highly expressed in m⁵C-cluster 2, followed by m⁵C-cluster 3, while m⁵C-cluster 1 had the lowest expression, which was consistent with the average survival time of samples in the three subtypes.

To better investigate the possible biological behaviors of all the m⁵C modification patterns, the limma package was used to determine 690 m⁵C phenotype-associated DEGs ( Supplementary Figure 3 ). In addition, KEGG pathway enrichment analysis was carried out on DEGs using the WebGestaltR package. It was surprising that these genes were enriched to cell cycle, DNA repair, cell adhesion molecules, and immune inflammatory response-related pathways. These findings verified the important role of m⁵C modification in cancer cells themselves and in TME immunomodulation ( Figure 4A ). To better validate such regulatory mechanism, the unsupervised clustering analysis was performed using those 690 m⁵C phenotype-associated genes, for the sake of classifying cases to distinct genome subtypes. Similar to clustering analysis of m⁵C modification patterns, three different m⁵C modification genome phenotypes were found, which were referred to as m⁵C gene-clusters A–C, separately ( Supplementary Figure 4 ). According to such results, there were three different m⁵C methylation modification patterns in PTC. Besides, there were diverse signature genes in the three different gene clusters ( Supplementary Figure 4 ). The m⁵C regulator expression levels were significantly different among the three m⁵C gene-clusters ( Supplementary Figure 5 ), consistent with the results obtained for m⁵C methylation modification patterns. The expression quantities of these nine genes were the highest in gene-cluster B samples, followed by gene-cluster A samples, and were the lowest in the gene-cluster C samples.

First of all, we analyzed the stromal scores of three m⁵C gene-cluster subtypes. The results suggested that ( Figure 4B ) there were significant differences in the stromal score of three subtypes, among which, gene-cluster C had the highest score, followed by gene-cluster B, while gene-cluster A had the lowest score. Then, we analyzed the distribution of 22 immunocytes in the three m⁵C gene-cluster subtypes. As observed from Figure 4C , the distribution of 15 immunocytes in three subtypes showed statistically significant differences. These findings revealed the important role of m⁵C methylation modification in the formation of diverse TME landscapes and tumor-related immune regulation.

Nonetheless, the above results were obtained from patient population alone, which might not precisely estimate the m⁵C methylation modification patterns of individual cases. Due to the m⁵C modification complexity and heterogeneity in individual samples, this study established the scoring system (m⁵C-score) using the phenotype-associated genes for quantifying m⁵C modification patterns in individual PTC cases. Besides, those attribute alterations in individual patients were visualized by the alluvial diagram ( Figure 5A ). It was discovered from the figure that, among the 3 m⁵C-cluster subtypes, samples in m⁵C-cluster 2 and m⁵C-cluster 3 subtypes were mostly distributed in the low m⁵C-score group, while those in the high m⁵C-score group were basically derived from the m⁵C-cluster 1 subtype. In the three m⁵C gene-cluster subtypes, the m⁵C-score values of Cluster A and Cluster C samples were lower. Samples aged over 40 years were mostly classified into the low m⁵C-score group, while females mostly belonged to the high m⁵C-score group.

To further evaluate the differences between low and high score samples, the limma package was used to analyze the DEGs between the two groups. Using the thresholds of logFC > log₂(1.2) and p < 0.05, 67 DEGs were screened, including 58 upregulated and 9 downregulated ones ( Figure 5B ). Moreover, the WebGestaltR package was utilized for the GO and KEGG enrichment analyses of DEGs, with p < 0.05 as the threshold. A total of 62 biological processes (BP), 2 cellular components (CC), 6 molecular functions (MF), and 9 pathways were selected. As shown in Supplementary Figure 6 , these genes were mainly involved in tumor proliferation and immune response-related biological processes/molecular functions and signaling pathways, such as MAPK, TNF, and IL-17.

Subsequently, this study observed the correlation of the m⁵C-score with patient survival and analyzed the difference in prognosis between high and low m⁵C-score samples. The results suggested that, samples with low m⁵C-scores had better prognosis than those with a high score, regardless of DFS or OS ( Figures 5C, D ). In addition, it was also discovered that there was no difference in the clinical features (such as T, M, and stage) between high and low m⁵C-score samples ( Supplementary Figure 7 ). The expression levels of nine m⁵C regulators in the high m⁵C-score group were significantly higher than those in the low score group, and there was significant difference between the two groups ( Supplementary Figure 8 ).

Subsequently, this study observed the correlation of the m⁵C-score with TME. First of all, the CIBERSORT method was adopted to evaluate the infiltration level of each immunocyte type in the high and low m⁵C-score TCGA-TPC samples. The results are presented in Supplementary Figure 9A . There were significant differences in six cell types between high and low m⁵C-score groups. In addition, this study also calculated the stromal score, immune score, and ESTIMATE score in different samples. As presented in Supplementary Figure 9B , in the low m⁵C-score group, the immune score was significantly higher than that in the high m⁵C-score group, which was consistent with the previous results that show that the low m⁵C-score group had better prognosis than the high score group. Moreover, it was discovered through expression of immune checkpoint genes that there were significant differences in 16 immune checkpoint gene expression levels between high and low m⁵C-score groups ( Supplementary Figure S10 ). Based on these findings, low m⁵C-score showed close correlation with immune activation. Furthermore, the m⁵C-score helped to assess m⁵C modification patterns in individual tumors and better assess the TME cell infiltration features of tumors, thus contributing to distinguishing the true or false TME immune infiltration.

Last, this study integrated the influences of the m⁵C-score and various immunocyte infiltration levels on the prognosis for PTC patients. From Figure 6 , it was discovered that resting CD4⁺ memory T cells and CD8⁺ T cells were mainly enriched in low m⁵C-score samples, while activated NK cells and monocytes were mostly enriched in the high m⁵C-score group. Then, the median infiltration level of the above four cell types was used to divide all samples into high and low immunocyte infiltration level groups. It was discovered that samples with low m⁵C-score and low infiltration level of resting CD4⁺ memory T cells had the best prognosis, while those with high m⁵C-score and low infiltration level of resting CD4⁺ memory T cells had the poorest prognosis. In addition, samples with low m⁵C-score and high CD8⁺ T cell infiltration had the best prognosis, while those with high m⁵C-score and low CD8⁺ T cell infiltration had the poorest prognosis. Furthermore, it was found that samples with low m⁵C-score and high monocyte infiltration had the best prognosis, while those with high m⁵C-score and high monocyte infiltration had the poorest prognosis. According to the prognostic prediction model, we analyzed the correlation between the m5C-score and Treg expression in 24 PTC patients. The m5c-score showed a negative relationship with CD3+CD4+/CD3+CD8+ (r = −0.9543, p < 0.0001; Figure 5E ), but a positive relationship with CD4+CD25+ Tregs percentage (r = 0.4477, p = 0.015; Figure 5F ).

First of all, this study calculated the correlation of 49 m⁵C phenotype-related genes with the m⁵C-score. Then, 10 signature genes related to the m⁵C-score were screened by the threshold of correlation coefficient >0.4, which were used as the features to construct the SVM classification model.

In order to verify the classification efficiency and accuracy of the model, we used the expression profile data of TCGA tumor samples as the training set. The m⁵C-score was utilized to classify the samples into high and low groups. Then, the expression profile data of these 10 genes were used to construct the SVM classification model to classify the TCGA-TPC samples. It was discovered that, compared with the m⁵C-score classification results, the accuracy reached 98.3%, and the sensitivity was up to 88.9%. The 493 samples were accurately classified, with an area under the ROC curve (AUC) of 0.936 ( Figure 7A ). The above results demonstrated that the classification model constructed based on these 10 signature genes well simulated the classification results of the m⁵C-score. The gene number was substantially reduced, which significantly improved the classification efficiency.

Thereafter, all the 551 TCGA samples (including 493 tumor samples and 58 normal samples) were used as the verification set 1. The abovementioned 10 genes were used as the features to construct the SVM classification model to classify the samples. Surprisingly, it was discovered that the model accurately classified TCGA-TPC samples into tumor samples and para-carcinoma tissue samples, with a classification accuracy of 89.7% and a sensitivity of 98.6%. Of the 551 samples in verification set 1, 538 were accurately classified, with an AUC of 94.2% ( Figure 7B ).

To further verify the model classification efficiency and accuracy, another three sets of microarray data were also downloaded, and the 10 signature genes were used for SVM verification. The GSE29265 dataset was utilized as verification set 2, which included 49 samples (20 normal samples and 29 tumor samples), with a model classification accuracy of 95%. Of the 49 samples, 48 were accurately classified, the model sensitivity to high and low scores was up to 100%, and the AUC was 97.5% ( Figure 7C ). Meanwhile, the GSE33630 dataset was used as verification set 3, which included 105 samples (45 normal samples and 60 tumor samples). The model classification accuracy reached up to 100%, all the 105 samples were accurately classified, the model sensitivity to high and low scores was 100%, and the AUC was 100 ( Figure 7D ). The GSE65144 dataset was used as verification set 4, which contained 25 samples (13 normal samples and 12 tumor samples). The model classification accuracy was 84.5%, all the 25 samples were accurately classified, the model sensitivity to high and low scores was 100%, and the AUC was 92.3% ( Figure 7E ).

We firstly used the L1000 fireworks display (L1000FWD) tool, and a reverse drug screening method for deferentially expressed genes in high- and low-risk groups of the m⁵C-score and obtained small molecules (drugs, Supplementary Table 4 ). In the interaction database between CMAP drug and gene expression, we analyzed 67 drugs that may interact with genes with different changes in the risk model constructed by the m⁵C-score, and selected 55 small molecules (drugs, Supplementary Table 5 ). We compared the potential drug overlap between L1000 and CMAP annotation, and found that there were five overlapping small molecules (S8), namely cephaeline, emetine, anisomycin, ouabain, and thapsigargin. CCK8 was used to detect the effect of five potential drugs on the growth and metabolic activity of PTC tumor cells. It was found that compared with the control group, the five drugs could inhibit the growth of thyroid cancer cells in different degrees ( Figure 8A ). Consistent with this, results of subcutaneous transplantation model also showed that intraperitoneal injection of these five drugs could significantly inhibit the growth of tumor, respectively ( Figure 8B ).