Gene expression, copy number variants, tumor somatic mutations, and matched clinical information of ccRCC from The Cancer Genome Atlas (TCGA) cohort were obtained. Gene expression data of 93 ccRCC tumors from the Clinical Proteomic Tumor Analysis Consortium (CPTAC) were obtained at https://proteomics.cancer.gov/programs/cptac. In addition, RNA-seq and clinical data of 91 ccRCC patients from the RECA-EU cohort were available from the International Cancer Genome Consortium (ICGC, https://dcc.icgc.org/) database and included in this study. Patients without overall survival information were removed from further analysis. In addition, 232 ccRCC samples with proteomics sequencing data with available clinical and pathologic electronic records were enrolled from our institute, Fudan University Shanghai Cancer Center (FUSCC, Shanghai, China). In total, 860 ccRCC tumor samples were included for analysis. Batch effects from nonbiological technical biases were corrected using the “ComBat” algorithm of sva package and the fragments per kilobase of transcript per million values were transformed into transcripts per kilobase million values.

A total of 17 m⁵C regulators were extracted from the integrated gene expression profiles to identify different m⁵C modification patterns. The 17 m⁵C regulators included three writers (DNMT1, DNMT3A, and DNMT3B), three erasers (TET1, TET2, and TET3), one regulating factor (DNMT3L), and 10 readers (MECP2, MBD1, MBD2, MBD3, MBD4, UHRF1, UHRF2, ZBTB4, ZBTB38, and ZFP57). The “ConsensusClusterPlus” R package was used to classify patients for further analysis, and 1000 times repetitions were conducted to ensure the stability of the classification (Wilkerson and Hayes, 2010). Overall survival was compared between patients with different m⁵C modification patterns.

GSVA, a commonly employed method for estimating the variation in pathways (Hänzelmann et al., 2013), was used to evaluate the potential biological differences between the m⁵C modification patterns using the “GSVA” R package. To estimate the TME cell infiltration, we applied the single-sample gene set enrichment analysis algorithm to evaluate the relative abundance of immune cells in the ccRCC TME. The reference gene sets for quantifying the immune cells were obtained from a previous study (Charoentong et al., 2017), and the examined immune cells included mast cells, monocyte, macrophage, activated CD4⁺ T cells, and other types of immune cells. Immune cell abundance was compared between m⁵C modification patterns, and the prognostic significance of each immune cell was also evaluated based on the overall survival information.

To explore the potential biological differences between m⁵C modification patterns, the limma package was used to identify differentially expressed genes (DEGs), and the threshold value was set as p < .05, |logFC|≥3 (Ritchie et al., 2015). Functional enrichment analyses were carried out to explore the potential functions of the DEGs. The expression profiles of DEGs were extracted, and unsupervised clustering was applied again to identify the subgroups; the subgroups were defined as m⁵C gene clusters.

A scoring system was constructed to evaluate the m⁵C modification patterns, and we termed it as m⁵C score. Univariate Cox regression was used to evaluate the prognostic value for each gene, and the genes with prognostic significance were extracted for further analysis. Random forest analysis and principal component analysis were used to construct the m⁵C relevant gene signature. Both principal components 1 and 2 were enrolled to calculate the signature scores, and the m⁵C score was defined as follows: m 5 C ⋅ score  =   ∑ ( PC 2 i + PC 2 i ) .

To explore potential associations between copy number variants and m⁵C score, Genomic Identification of Significant Targets in Cancer (version 2.0) was used to identify significantly amplified or deleted regions using TCGA copy number data (Beroukhim et al., 2007; Mermel et al., 2011). Q ≤ 0.05 was defined as significant, and the confidence interval was set to 0.95. Tumor immune dysfunction and exclusion (TIDE) was used to estimate the immunotherapy response based on the expression profiles (Jiang et al., 2018). Thus, the associations between m⁵C score and immunotherapy response were evaluated. The pRRophetic package was used to predict the half-maximal inhibitory concentration (IC50) of chemotherapy drugs in the high and low m⁵C score groups.

IHC was performed to evaluate the expression levels of Ki-67 (ab15580; Abcam), CD4 (RMA-0620, Maxim, China), CD8 (RMA-0514, Maxim, China), Glut-1 (ab115730; Abcam), PD-L1 (ab205921; Abcam), CXCL13 (ab246518; Abcam), TGF-β (ab189778; Abcam), FASN (ab99359; Abcam), CK (Kit-0009, Maxim, China), and FoxP3 (98,377, CST) following previously described procedures (Xu et al., 2021a; Xu et al., 2021b). Opal multispectral was implemented to identify differential immune cell infiltration and PD-L1 expression in different groups on a multispectral imaging system (Vectra^® Polaris™, Shanghai, China).

A Wilcox test was used to compare differences between two groups. The Kaplan–Meier method was used to conduct survival analysis, and the cutoff value was defined via the survminer package. A log-rank test was used to detect the significance. The receiver operating characteristic curve (ROC) was drawn to evaluate the predictive ability for immunotherapy response.

A total of 17 m⁵C regulators including three writers, three erasers, one regulating factor, and 10 readers were manually identified in this study. The dynamic reversible process of m⁵C DNA methylation mediated by regulators as well as their potential biological functions for ccRCC are summarized in Figure 1A. We detected significant differences in the expressions of m⁵C regulators between ccRCC and para-cancer tissues (p < .05) (Figure 1B). Analysis of CNV frequency indicated that CNV alterations were prevalent in the 17 m⁵C regulators, and half of the m⁵C regulators more frequently showed copy number amplification compared with copy number loss (Figure 1C). Besides this, in DNA variation profiles, we found 19 experienced samples of m⁵C regulators with a frequency of 5.65% among 336 ccRCC samples from TCGA (Figure 1D). The location of CNV alterations of the m⁵C regulators on chromosomes is shown in Figure 1E. Notably, ccRCC samples could be distinguished from normal samples completely based on the expression pattern of these m⁵C regulators (Figure 1F). These findings suggest a high degree of m⁵C modification–mediated intertumoral heterogeneity of genetic and expressional alteration landscape between ccRCC and adjacent normal samples, suggesting that the aberrant expression of m⁵C regulators may play an essential role in ccRCC malignancy.

Three data sets (both proteome and transcriptome) with available survival and clinicopathological information (TCGA, CPTAC, and RECA-EU) were included in one meta-cohort. Figure 2A shows the comprehensive landscape of the interaction of the 17 m⁵C regulators, interaction network, and the prognostic implications for ccRCC patients. The results identified TET2, MBD1, MBD2, MECP2, ZBTB4, ZBTB38, and UHRF2 as protumorigenesis indicators for ccRCC, and MBD3, UHRF1, and DNMT3B were identified as significant favorable factors for ccRCC. We also found that expression of the m⁵C regulators in the same functional category exhibited remarkable correlations, and a marked association was displayed among writers, regulators, erasers, and readers. For instance, whether ccRCC tumors with a high writer gene expression exhibit a high eraser gene expression normally depended on the different writer and eraser genes. However, we found that tumors with high expression of the m⁵C reader gene ZFP57 showed low expression of some reader genes (ZBTB4, UHRF2, MBD3, and MBD2) although the high expression of other reader genes was not affected. These results imply that a cross-talk among the genes encoding m⁵C writers, readers, regulators, and erasers could play essential roles in the malignancy of different m⁵C modification patterns and tumor immune microenvironment characterization among individual ccRCC samples.

We next used the ConsensusClusterPlus R package to identify ccRCC patients with qualitatively different m⁵C modification patterns based on the transcriptional expression of 17 m⁵C regulators. Two distinct modification patterns were classified using unsupervised clustering: m⁵C cluster A (including 271 cases) and m⁵C cluster B (354 cases) (Figure 2B). Survival analysis of patients in the two clusters revealed that patients with the m⁵C cluster A modification pattern showed improved survival compared with patients with the m⁵C cluster B pattern (Figure 2C).

To investigate the clinical differences and biological processes between the two distinct m⁵C modification patterns, we constructed a clustering heat map showing differentially expressed m⁵C regulators and clinical information, including age, sex, stage, and survival status in the two m⁵C modification patterns (Supplementary Figure S1). GSVA enrichment analysis indicated that m⁵C clusters mainly differ in heterochromatin, peptidyl modification pathways, and microRNA post-transcriptional regulation (Figure 3A). The ccRCC samples in m⁵C cluster A showed prominent upregulation in E2F1, miR-147B, miR-3910, miR-4261, miR-3689-3p, miR-4719, PBXIP1, and ZNF184 targeted regulation and downregulation in peptidyl modification processes, such as histone binding, peptide amino acid modification, protein autoubiquitination, ubiquitin-like protein ligase, and transferase activities (Supplementary Figure S2A).

We next examined immune cell infiltration to assess differences in the immune context of the TME between m⁵C modification patterns. m⁵C cluster A was remarkably rich in innate immune cell infiltration and the active immune response process with a high abundance of activated CD4 T cells, immature B cells, regulatory T cells, Tfh cells, dendritic cells, eosinophils, macrophages, mast cells, natural killer cells, and neutrophils (Figure 3B). The results from GSVA analyses demonstrate that the m⁵C cluster A modification pattern, which predicts favorable clinical outcome, was significantly associated with antitumor immune responses. Therefore, we hypothesized that the peptidyl modification inactivation in m⁵C cluster A may be involved in the antitumor effects of immune cell infiltration related to this cluster.

We further assessed the prognostic implications of immune cell infiltration in ccRCC (Figure 3C). Univariate Cox analysis indicated that T follicular helper cells (p = .022), immature B cells (p = .013), mast cells (p = .006), type 17 T helper cells (p = .036), and activated CD8 T cells (p = .036) could serve as independent prognostic protective factors in ccRCC, and MDSC (p < .001) was a remarkable risk indicator for 616 ccRCC patients from the TCGA and CPTAC cohorts (Figure 3C). When clinicopathological factors were analyzed, we found no significant differences in the pathology types and genetic variations between the patients in the two m⁵C modification pattern groups, which suggests that DNA m⁵C methylation modification did not influence clinical and pathologic features of tumors (Supplementary Figure S2B–G).

To further explore the biological consequences of the distinct m⁵C modification patterns, we then investigated the genetic constitution of individual m⁵C clusters patterns and identified 180 m⁵C phenotype–related DEGs using the Limma package of R software. Random forest analysis was performed to determine the most important m⁵C gene signatures in identifying m⁵C modification patterns (Figure 4A). To investigate the regulation mechanism of DNA m⁵C modifications on ccRCC, we then performed unsupervised clustering analyses based on the obtained 180 m⁵C phenotype–related signatures to classify patients into different genotypes. Consistent with the clustering of m⁵C modification patterns, the unsupervised clustering algorithm also revealed two distinct m⁵C modification genomic subtypes, named as m⁵C gene clusters A and B (Figure 4B). Kaplan–Meier analysis of ccRCC cases in the combined discovery TCGA and test CPTAC cohorts revealed that patients in the m⁵C gene cluster B group (n = 247) showed significantly poor survival compared with cases in m⁵C gene cluster A (n = 369) (Figure 4C). Prominent differences in the expression of m⁵C regulators between the distinct m⁵C gene clusters were confirmed using unpaired t test, and the results were in accordance with the results of DNA m⁵C methylation modification patterns (Figure 4D). These results revealed the presence of distinct m⁵C methylation modification patterns in ccRCC and showed that these patterns could distinguish aggressiveness in ccRCC.

Next, the clusterProfiler package was used to perform GO and KEGG functional enrichment analysis for the m⁵C DEGs. The biological processes, cellular components, and molecular functions with significant enrichment are summarized in Figure 5A. Enriched terms in biological processes were related to m⁵C modification, neutrophil activation–related immune response, and response to hypoxia, which provided a basis that m⁵C modification may play an important role in the immune regulation of the ccRCC TME (Figure 5A). We further found that ccRCC samples in m⁵C gene cluster B showed advanced clinical stages and exhibited higher mortality (Figure 5B). Older patients were concentrated in the m⁵C gene cluster B, and the distinct genotype clusters were characterized by different m⁵C signature genes.

The above findings demonstrate that the m⁵C methylation modification plays a key regulatory role in reshaping different TME landscapes. Nevertheless, these results were determined on the patient population and might not provide accurate information on survivorship based on m⁵C modification patterns in individual ccRCC patients. Considering the individual intratumor heterogeneity of m⁵C methylation and using the phenotype-related genes, we establish a scoring system for easy quantification of the m⁵C modification patterns for individual ccRCC patients and named this system m⁵C score. The alluvial diagram was applied to visualize the alterations of individual patients (Figure 6A). The m⁵C score clusters prominently classified the patients into two prognostic groups (good and poor) and enabled stratification of patients in both the discovery TCGA and validation real-world FUSCC cohorts. Survival analysis indicated that high m⁵C score was significantly correlated with poor overall survival (HR = 0.3 with 95% CI from 0.22 to 0.41, p < .0001) in 516 patients with ccRCC from TCGA (Figure 6B) and correlated with worse overall survival in 266 patients with ccRCC from FUSCC (Figure 6C).

We next investigated the relationship of m⁵C score with clinical and pathologic characteristics in ccRCC patients from the training, testing, and validation cohorts. Consistent with its prognostic value, the m⁵C score significantly increased with advancing clinical stages and aggressive ISUP grade and reached the highest level at stage IV or grade 4 (Figures 7A,B). There was no difference in age between the two clusters. The proportion of males in the high m⁵C score group was markedly higher than that of females, which is consistent with the result that male patients have a worse prognosis than female patients with ccRCC (Figures 7C,D).

To reveal the role of the m⁵C score phenotype in the comprehensive molecular landscape of ccRCC, we examined tumor somatic mutation and evaluated DNA variation in the m⁵C score clusters. Patients with mutation in PBRM1, a gene frequently mutated in ccRCC, showed a prominently lower m⁵C score compared with patients with wild-type PBRM1 (Figure 7E). The m⁵C score did not show a significant association with tumor mutation burden in patients with ccRCC (Figure 7F).

We next evaluated the differences in the DNA variation landscape in the two m⁵C score clusters. The top 20 frequently mutated genes in the m⁵C score clusters are shown in Figure 7G,H. VHL (mutation frequency, 40%), BAP1 (13%), SETD2 (13%), TTN (13%), and MTOR (11%) were the five most frequently mutated genes in the m⁵C score^high group (Figure 7G), whereas VHL (24%), PBRM1 (19%), TTN (15%), SETD2 (7%), and MTOR (6%) were the five most frequently mutated genes in the m⁵C score^low group (Figure 7H). Thus, we speculate that the significantly higher mutation frequency of BAP1 in the high m5C score group may contribute to the poor prognosis for ccRCC patients and the low mutation frequency of PBRM1 may reduce immunotherapy efficiency for ccRCC patients. Copy number variant features are depicted in Figure 7I,J. In addition to the common mutation site located in 5q35.3, copy number variant in m5C score^high samples were generally located in 3q25.33, 2q10.53, and 9p12.3 loci.

To define the role of m⁵C-related phenotypes in regulation of the TME, we first investigated cancer-related pathways characterizing m⁵C gene clusters based on training and testing cohorts. As shown in Figure 8A, TGF-β signaling, oxidative phosphorylation, and fatty acid metabolism were significantly downregulated in ccRCC samples in the m⁵C score^high group compared with the m⁵C score^low group, whereas pathways involved in protumorigenesis responses of the TME, such as hypoxia, glycolysis, epithelial-mesenchymal translation, and IL6-JAK/STAT3 signaling, were markedly upregulated in the m⁵C score^high group. We next evaluated the immune cell infiltration in the TME in m⁵C-related phenotype clusters. The results indicated that CD4⁺ T cell memory resting, mast cell resting, and monocyte and NK cell infiltration significantly correlated with a high m⁵C score, whereas plasma cell, M0 macrophage, Treg cell, and neutrophil infiltration were significantly associated with low m⁵C score in ccRCC patients (Figure 8B). To evaluate the regulatory role of m⁵C score in TME, we explored the expression of chemokine, cytokine, and immune checkpoints in m⁵C score clusters (Figure 8C). We found the expression levels of immune checkpoint factors were significantly different in the m⁵C score^high group, suggesting that the high m⁵C score cluster may indicate an immune-suppressive microenvironment.

Immunotherapies, including anti-immune checkpoints, are revolutionizing the field of cancer therapy. RCC is resistant to traditional cytotoxic chemotherapy but can be responsive to immunotherapy. Therefore, we investigated whether the m⁵C modification signature could predict the responses to chemotherapy and ICTs in the combined ccRCC cohorts (n = 860, TCGA, CPTAC, and FUSCC). Evaluation of the ICC50 of cisplatin showed that the low m⁵C score group was significantly correlated with a higher IC50 value, which indicates that the low m⁵C score group may be less sensitive to cisplatin (Figure 9A). However, no significant differences were observed in predicting IC50 values of gemcitabine between the m⁵C modification groups (Figure 9B). The TIDE algorithm was used to predict intratumoral heterogeneity and responsiveness to immunotherapy. The findings indicate that a higher m⁵C score was significantly correlated with an elevated TIDE score, suggesting that the high m⁵C score group may show a reduced response to immunotherapy, such as PD-1 and PD-L1 blockade (Figure 9C). The ROC curve showed a relatively stable ability for predicting the immunotherapy response of m⁵C score with an AUC of 0.676 (Figure 9D).

To further test the stability of m⁵C score model, we applied the m⁵C score signature established in the real-world FUSCC proteomics cohort and evaluated TME characteristics by IHC staining analysis of 30 consecutive ccRCC tissue sections. IHC staining revealed significantly decreased CD8, PD-L1, and GLUT-1 expression and elevated FoxP3, CXCL13, and FASN expression and Ki-67 staining in tumors from the FUSCC cohort (p < .05) in the m⁵C score^high group (Figure 10), suggesting immune-suppressive characteristics of the TME. Furthermore, we found a significantly decreased number of infiltrated CD4⁺ T cells and CD8⁺FoxP3⁺ Treg cells and downregulated PD-L1 expression in the immune-cold m⁵C score^high group using opal multimarker IHC staining (Figure 10). In general, the data from multiomics bioinformatics to the real world demonstrate that lower m⁵C score predicts better responses to immunotherapy for ccRCC patients.