The integrated dataset containing seven independent PAAD datasets has been described in our previous studies (9). We collected the PAAD datasets from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/), The Cancer Genome Atlas (TCGA) (https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga), International Cancer Genome Consortium (ICGC) (https://dcc.icgc.org/), and ArrayExpress databases (https://www.ebi.ac.uk/arrayexpress). Patients without survival information were eliminated from the integrated datasets. Seven independent datasets, namely, the GSE28735, GSE57495, GSE62452, MTAB-6134, TCGA-PAAD, ICGC-AU and ICGC-CA datasets, were collected as the training cohort, and the datasets included 930 patients with PAAD. The details of these seven datasets were summarized in sheet 1, 3, 8 and 9 of Supplementary Table 1 . XLX (In the supplementary data). For it’s difficult to find more PAAD datasets with clinical data, we validated our results in three datasets of the seven datasets. The ICGC-AU and ICGC-CA datasets, including 295 patients, and the TCGA-PAAD dataset were used as the validation cohort for further analysis. The ComBat method from the “SVA” package was used to correct the batch effects from nonbiological technical biases (21). In this process, some information was lost, and we could not detect the expression of YBX2 in the integrated dataset with 930 PAAD patients.

Twelve m⁵C regulators were determined based on previous studies (14, 22). For GSE62165 dataset has no clinical data, it’s excluded from the integrated dataset. For GSE62165 has both PAAD and corresponding normal tissues, it was used to analyze the mRNA expression of these twelve regulators. Additionally, the protein expression of these regulators was acquired from the Human Protein Atlas (HPA) database (http://www.proteinatlas.org/). Genetic alteration was thought to modulate gene expression. To uncover the underlying mechanism of m⁵C regulators expression pattern, the cBioPortal database (http://www.cbioportal.org/) was used to examine m⁵C regulator genetic alterations, including amplification, deep deletion, and missense mutations.

Based on m⁵C regulator expression levels, the patients with PAAD were divided into different clusters using the unsupervised clustering package “ConsensusClusterPlus” (23). To validate the classification stability, this process was repeated 1000 times. The principal component analysis was performed using the “PCA” package in R software.

The ssGSEA algorithm was used to analyze the infiltration level of 24 immune cell types that play critical roles in cancer immunity (24). The marker genes for each type of immune cell were determined according to previous studies (25–27). Based on the immune infiltration level, the patients in the dataset were divided into high-, moderate- and low-infiltration groups using unsupervised clustering.

The ESTIMATE algorithm was used to calculate the immunoscore and stromal score (28). By analyzing specific mRNA expression, we also determined tumor cellularity and tumor purity (28). The TIDE score was determined using the TIDE algorithm (29). A high TIDE score indicates that cancer cells are more likely to resist immunotherapy and promote immune evasion.

Firstly, we extracted differently expressed genes (DEGs) between normal and tumor tissues in PAAD patients from GEPIA database (http://gepia2.cancer-pku.cn/#index) (30). Then, we obtained 2439 differentially expressed genes among three m⁵C clusters by using limma packages in R software and an adjusted p<0.05 was set as the cutoff value (31). Lastly, we took intersection from DEGs (from GEPIA) and 2439 DEGs (from m⁵C clusters) to increase the specificity of the target genes. Generally, we obtained 1720 merged differentially expressed m⁵C-related genes (m⁵C DEGs). We performed unsupervised clustering of these m⁵C DEGs in training cohort and divided the cohort into three clusters, named gene cluster A, B, and C. The gene signatures that were positively correlated with gene cluster were named as m⁵C gene signature A, and the rest m⁵C DEGs were termed as m⁵C gene signature B. We performed functional enrichment analysis of m⁵C gene signature A and B using over representation analysis. In order to reduce the noise or redundant genes, we used the Boruta algorithm to perform dimension reduction in the m⁵C gene signature A and B (32).

To construct the m⁵C score, principal component analysis (PCA) was utilized to construct m⁵C score from the m⁵C gene signature A and B. (1) the PC1A represent the first component of m⁵C gene signature A, and (2) the PC1B represent the first component of m⁵C gene signature B. m⁵C score = ∑PC1A-∑PC1B. 930 patients in training cohort were stratified into two subgroups as high- and low- m⁵C score by using surv_cutpoint function in the survminer package in R software.

GSE91061 was downloaded to explore the ability of the m⁵C score to predict immunotherapy response. The GSE91061 dataset contains patients with melanoma receiving PD-1 treatment (33).

Gene set variation analysis (GSVA) was used to explore the underlying characteristic of the high and low m⁵C score groups of patients with PAAD (34). The hallmark gene set (h.all.v6.2.entrez.gmt) was obtained from the Molecular Signatures Database (MSigDB). The cancer immunity cycle and immunotherapy-related pathways were detected according to previous studies (35, 36). Additionally, we obtained cancer driver genes from the Integrative OncoGenomics database (IntOGen, https://www.intogen.org/) (37).

The "pRRophetic" package in R was used to analyze the chemotherapy response (IC50) based on the Genomics of Drug Sensitivity in Cancer (GDSC) database (38). IC50 indicates the effects of the tested drugs on suppressing cell growth. If p value < 0.05, it indicates the sensitivity of the drug is different between high m⁵C and low m⁵C groups

PAAD tissues and control tissues were obtained from patients who underwent surgical resection at Tangdu Hospital. Written informed consent was obtained from every patient, and all procedures performed in patients were approved by the Research Ethics Committee of Tangdu Hospital. The tissues were cut into small pieces using scissors and lysed with lysis buffer (C1053-100, Applygen, Beijing, China) supplemented with protease and phosphatase inhibitors (04906837001 and 04693159001, Roche, Basel, Switzerland). After incubation on ice for 0.5 h, the samples were centrifuged, and the supernatant was collected. The protein concentration was analyzed using a BCA kit (23225, Thermo, MA, US). The samples were analyzed using 10% SDS gel electrophoresis and transferred to polyvinylidene fluoride (PVDF) membranes. The blots were incubated with primary antibodies at 4°C overnight. As YBX2 was not detectable in the integrated dataset, we only examined 11 m5C modification regulators in PAAD and corresponding tissues. Anti-DNMT1 (A19679), anti-DNMT3A (A19659), anti-DNMT3B (A11079), anti-NSUN2 (A3443), anti-NSUN3 (A12892), anti-NSUN4 (A14983), anti-NSUN5 (A5992), anti-NSUN6 (A7205), anti-TRDMT1 (A10535), anti-YBX1 (A7704) and anti-GAPDH (A19056) antibodies were purchased from Abclonal (Wuhan, China). The anti-NSUN7 (17546-1-AP) antibody and secondary HRP antibodies (SA00001-1 and SA00001-2) were purchased from Proteintech (Wuhan, China). The blots were washed with TBST buffer three times and incubated with the corresponding secondary HRP antibodies for 2 h at room temperature. Finally, the blots were assessed using a GelDoc XR+ gel imager (Bio-Rad, USA).

R-4.0.4 was used to perform the statistical analyses in this study. For quantitative data, Student’s t tests were used for normally distributed variables, and the Wilcoxon rank-sum test was used for nonnormally distributed variables. To compare more than two groups, Kruskal–Wallis tests and one-way ANOVA were used as nonparametric and parametric methods, respectively. The Kaplan–Meier method and the Cox proportional hazards model were used to analyze survival. The “survival” and “Survminer” packages were used to examine the association between the factors and prognosis. The surv-cutpoint function from the “survival” package was applied to divide the included PAAD samples into different subtypes. Correlation analysis in the present study was performed with Pearson or Spearman correlation coefficients and using the “corrplot” package in R software. All comparisons were two-sided with an alpha level of 0.05, and the Bonferroni method was used to control the false discovery rate (FDR) for multiple hypothesis testing.

m⁵C regulators play a critical role in maintaining target RNA stability, translation, or translocation (39). However, the landscape of expression and genetic alterations of m⁵C regulators has not been elucidated in PAAD. We analyzed the expression pattern of m⁵C regulators in GSE62165, which contains 131 PAAD tissues and adjacent normal tissues. The present results indicated that the levels of DNMT3B, NSUN6, NSUN7 and YBX2 mRNA were significantly decreased while the levels of DNMT1, DNMT3A, NSUN2, NSUN3, NSUN5 and YBX1 mRNA were significantly increased in PAAD patients ( Figures 2A, B ). We also analyzed the protein expression of m⁵C regulators using the Human Protein Atlas (HPA) database ( Supplementary Figure 1 ). Our results showed that NSUN6 levels were decreased in PAAD, in line with previously established expression patterns (20). Genetic alterations are critical factors that influence disease progression (40). We examined the genetic alterations of these regulators in cBioPortal database. These results indicated that 51 of the 259 (19.69%) PAAD patients in 2 studies experienced genetic alterations of 12 m⁵C regulators ( Figure 2C ). The results confirmed that m5C regulators are differentially expressed in PAAD patients, which may play a critical role in PAAD progression.

As shown in our previous studies, we conducted further analysis in the integrated dataset with 930 PAAD patients, which was composed of seven dependent PAAD dataset with clinical data (9). According to the similarity displayed by the m⁵C regulator expression levels, k = 3 was identified, with optimal clustering stability from k = 2 to 9 ( Supplementary Figure 2 ). The integrated dataset was divided into three independent m⁵C clusters ( Figure 3A ). As shown in Figure 3B , m⁵C clusters were significantly related to PAAD patient survival status. m⁵C cluster 2 had the best overall survival (OS) when compared to clusters 1 and 3 ( Figure 3B ). This result confirms that m⁵C regulators are critical factors influencing the prognosis of PAAD. It has been reported m⁵C modification is involved in immune microenvironment regulation (22). Then we analyzed the correlation between m⁵C regulators and immune microenvironment of PAAD. The result indicated immune cell infiltration was differentially distributed between the m⁵C clusters. Specifically, m⁵C cluster 2 had the highest proportions of B cells, NK CD56 bright cells, Tfh cells, Th17 cells and the lowest proportion of Treg cells when compared to clusters 1 and 3 ( Figure 3C ). These cells have been reported to activate the immune response and kill cancer cells (41). Additionally, we analyzed the relationships between the m⁵C clusters and immune-related genes. In the present study, CD274, CTLA4, HAVCR2, IDO1, LAG3, and PDCD1 were included in the immune checkpoint-related signature, and CD8A, CXCL10, CXCL9, GZMA, GZMB, PRF1, TBX2, and TNF were included in the immune activity-related signature (42). The results showed that cluster 2 had the lowest levels of CTLA4 and IDO1, which play an inhibitory role in immune cell activation ( Figure 3D ). These results showed that m⁵C cluster 2 had the most potent immune response when compared with cluster 1 and 3, which partly account for the best prognosis of patients in the m⁵C cluster2. These results confirmed that m⁵C regulators are related to the PAAD immune microenvironment. The tumor immune microenvironment plays a role in the effect of m5C regulators on patient prognosis.

To further probe into the biological behaviors accounting for the important role of m⁵C regulators in the PAAD TME and PAAD development, we tried to identified the m⁵C DEGs. As described in the method section, DEGs from GEPIA dataset was intersected with DEGs from three m⁵C clusters. Finally, we identified a total of 1720 m⁵C DEGs ( Figure 4A ). Based on the unsupervised clustering of these m⁵C DEGs, 930 PAAD patients were divided into three m⁵C-related gene clusters, namely, m⁵C gene clusters A, B, and C ( Figure 4B and Supplementary Figure 3 ). The results showed that the m⁵C gene clusters were strongly correlated with survival status: m⁵C gene cluster C had the best OS when compared to clusters A and B ( Figure 4C ). We also analyzed the correlations with immune cell infiltration. m⁵C gene cluster C had the highest proportion of NK CD56bright cells and the lowest proportion of Treg cells ( Figure 4D ). Additionally, these two kinds of immune cells were overlapped with that in m⁵C cluster2. Functionally, NK CD56bright cells are mainly responsible for innate immune response and Treg cells are one of the most canonical immune-suppressive cells in the body (43, 44). Moreover, m⁵C gene cluster C had the lowest levels of CD274, CTLA4 and PDCD1 ( Figure 4E ). These results confirmed m⁵C gene cluster C had the most potent immune response when compared with cluster A and B. Collectively, these results further confirmed that m⁵C regulators are critical factors influencing the PAAD TME and PAAD development. The coherence between the prognostic and TME features in the three gene clusters indicated that this classification was reliable and reasonable.

To acquire quantitative indicators of m⁵C regulators, a heatmap was generated to visualize the correlation among the m⁵C clusters, m⁵C gene clusters, patient status, and immune infiltration. As shown in the left axis of the heatmap, m⁵C gene signatures that were positively correlated with m⁵C gene cluster were named as m⁵C gene signature A, and the rest m⁵C DEGs were termed as m⁵C gene signature B ( Figure 5A ). This method has been verified in previous study (42). Then the dimension reduction was presented by the Boruta algorithm based on the m⁵C gene signatures A and B and principal component analysis (PCA) was utilized to extract the principal component from the m5C gene signature A and B. (1) the PC1A represent the first component of m5C gene signature A, and (2) the PC1B represent the first component of m5C gene signature B. m5C score = ∑PC1A-∑PC1B. 930 patients in training cohort were stratified into two subgroups as high- and low- m⁵C score by using surv_cutpoint function in the survminer package in R software. The results indicated that higher m⁵C score was related to longer OS than lower score ( Figure 5B ). Consistently, cluster C had the highest m⁵C score and the best OS rate ( Figure 5C ). These results suggested m⁵C score was a potential indicator to predict PAAD prognosis. It has been validated m⁵C modification was related to immune microenvironment of PAAD, we also analyzed the correlation between the m⁵C score, infiltrating immune cells, and immune-related genes. The high m⁵C score group had a higher proportion of NK CD56 bright cells and a lower proportion of Treg cells, which was consistent with the result of m⁵C cluster 2 and gene cluster C ( Figure 5D ). Moreover, the high m⁵C score group had lower levels of CD274, CTLA4, IDO1, LAG3, and PDCD1 expression despite the fact that the differences in some factors were not statistically significant ( Figure 5E ). Additionally, we determined the relationship between the m⁵C score and cancer-associated fibroblasts (CAFs), T cell dysfunction, microsatellite instability (MSI), and the TIDE score. It has been validated an increased number of CAFs promotes cancer development by hindering immune cell infiltration into cancer tissues, while MSI and the TIDE score are precise predictors of chemotherapy and immunotherapy response (45–47). High m⁵C score was found to be related to lower CAF levels, T cell dysfunction, and lower TIDE scores but were also related to higher MSI ( Figures 6A–H ). Moreover, the high m⁵C score group had a lower stromal score, immune score, and ESTIMATE score but higher tumor purity ( Figure 6I ). These results further confirmed m⁵C score was related with immune microenvironment of PAAD. High m⁵C score indicated a better immune response.

The cancer immunity cycle is useful for determining which steps are involved in the immune response (48). Next we aimed to explore the immune steps regulated by the m⁵C score via using the cancer immunity cycle tool. The results suggested that the high m⁵C score group had high proportions of basophils and neutrophils ( Figure 6J ). Moreover, we aim to explore the mechanism underlying m⁵C score. The correlation of m⁵C score and immunotherapy-related pathways was analyzed. The m⁵C score was related to base-excision repair, the cell cycle, DNA replication, the p53 signaling pathway and other pathways ( Figures 6K, L ). These results assisted to elucidate the mechanisms of m⁵C score in regulating immune response. To comprehensively validate the utility of the m⁵C score, a heatmap was constructed to delineate the correlation between the m⁵C score and the TME. The m⁵C score was related to the proportions of CD8⁺ T cells, NK cells and B cells, the expression of immune checkpoints and MHC-II, the EMT process, the cell cycle and mismatch repair ( Supplementary Figure 4 ). m⁵C modification were closely related with immune microenvironment and participated in immune response modulation via multiple signaling pathways.

The m⁵C score reflects regulation and modulation of the immune microenvironment. However, whether the m⁵C score could predict immunotherapy response had not been explored. We examined the utility of the m⁵C score in predicting immunotherapy response. The results showed that patients with high m⁵C score had better sensitivity and response to immunotherapy. A high m⁵C score was related to a higher PR (partial response)/CR (complete response) rate and better OS when compared with low m⁵C score ( Figures 7A–C ). It has been validated chemotherapy agents activated immune responses in patients with tumors, which assisted to improve immunotherapy response (49, 50). It propelled us to examine the chemotherapy agents which were responded differently in high and low m⁵C score group. The results indicated IC50 values of 5-Fluorouracil, Erlotinib, Mitomycin C, Trametinib, SB52334 and AS605240 were decreased in high m5C score group ( Figures 7D–I ). It has been proved some of these regents exerted an inhibitory role on PAAD development (51–53). Combination with these agents may greatly enhance the immunotherapy responsiveness. Our present results confirmed that the m⁵C score is a precise indicator for the prediction of immunotherapy benefit.

The m⁵C score was verified to predict PAAD patient prognosis in the training cohort. Next, we validated the utility of the m⁵C score using the TCGA PAAD dataset, the ICGC-AU and the ICGC-CA datasets. Consistently, high score was related to better OS, disease-specific survival (DSS), disease-free interval (DFI), and progression-free interval (PFI) in the TCGA PAAD dataset ( Figures 8A–D ). OS was also demonstrated to be related to the m⁵C score in the ICGC-AU and ICGC-CA datasets ( Figure 8E ). We further examined the relationship between the m⁵C score and tumor mutation burden (TMB). TMB has been reported to be related to cancer development and act as an indicator of the immunotherapy response and prognosis (54). In the present study, there was no significant relationship between these two factors, indicating that the m⁵C score was an independent immunotherapy response indicator ( Figures 8F, G ).

To further uncover the underlying molecular mechanism of the m⁵C score, we analyzed the correlation of cancer-related driver genes with immune infiltration, the m⁵C clusters, the m⁵C gene clusters and the m⁵C score ( Figure 9A ). The results indicated that ACVR2A, ARID1A, ARID2, BAP1, BCL11B, BCORL1, BIRC6, CDH10, FAT1, FAT4, FBXW7, GNAS, KDM6A, KRAS, MARK2, NRAS, PCBP1, RBM10, RNF43, SETD1B, SETD2, SF3B1, SMAD3, SMARCA4, SPEN, WT1 and ZNF521 were differentially expressed between the high and low m⁵C score groups ( Figure 9B ). Gene set variation analysis (GSVA) was also conducted in the high and low m⁵C score groups. The terms reactive oxygen species pathway, DNA repair, unfolded protein response, p53 pathway and PI3K-AKT-MTOR pathway were differentially activated. Moreover, terms related to metabolic processes such as fatty acid metabolism, cholesterol homeostasis, glycolysis, and adipogenesis were significantly altered between the high and low m⁵C score groups. Additionally, inflammation-related pathways were also altered between the high and low m⁵C score groups, including the interferon alpha response, interferon gamma response, TNFα signaling via NFκB, IL6-JAK-STAT3 signaling and IL2-STAT5 signaling pathways ( Figure 9C ). Finally, we verified the expression pattern of m⁵C regulators in PAAD patient tissues and corresponding normal tissues ( Figure 10 ). Consistently, DNMT1, DNMT3A, NSUN2, NSUN3, NSUN5 and YBX1were increased while DNMT3B, NSUN6 and NSUN7 were decreased in PAAD samples when compared with mRNA and immunohistochemical results. Consistent with previous studies, DNMT1, DNMT3A and YBX1 were increased while NSUN6 was decreased in PAAD tissues when compared with adjacent normal tissues (20, 55–57).