The gene mutation files and RNA sequencing data (count value) with the corresponding clinical information from ccRCC patients were downloaded from TCGA database (https://portal.gdc.cancer.gov/repository). Based on the clinical data, the 192 somatic mutation data and 261 RNA transcriptomes with early-stage ccRCC (pT1) were selected for inclusion. Data with VarScan2 analysis of somatic mutation and tumor mutation burden were downloaded from TCGA database. The result was displayed with the R packages ‘maftools’ (12). RNA sequencing counts in tumor tissues and adjacent relatively normal tissues were standardized and analyzed using the DESeq2 package (13).

To assess the reliability of the m6A methylation regulator expression across protein levels, 27 paired early-stage ccRCC tissues were then collected and analyzed by means of HPLC-MS. The present study was approved (approval no. KS2021034) by the Institutional Review Board of Peking Union Medical College Hospital, Chinese Academy of Medical Science and Peking Union Medical College (Beijing, China). Formal written consent was requested and signed by all patients to signify approval to participate. All patients received the partial or radical nephrectomy during the period of June 2019 to September 2019 and the ccRCC diagnosis was determined by at least two independent specialists for pathological analysis in our center. The HPLC-MS process was divided into the following steps:

i) Sample collection and preparation: A total of 27 paired ccRCC tissues, including 27 paired tumor and non-cancerous normal tissues (NATs), were collected for further analysis. The ccRCC tissues were resected from the center of the tumor without necrotic regions. NATs were defined as tissues at least 5 mm away from the cancer capsule and no tumor invasion was identified under microscopy. All tissue samples were collected from the operation room and transferred directly into the −80°C refrigerator for storage until final analysis took place.

A total of ~25–120 mg of each cryo-pulverized renal tumor tissue or NATs were homogenized separately in an appropriate volume of lysis buffer [i.e. 2% SDS, 20 mM Tris, Cocktail (1:100), DNAse (1:100), RNAse (1:1,000)] by repeated vortexing. Protein concentration was determined using Pierce™ BCA assay protein assay kit (Thermo Fisher Scientific, Inc.). A total of 100 mg of protein was reduced with 20 mM dithiothreitol (DTT) for 5 min at 95°C and subsequently alkylated with 50 mM iodoacetamide for 45 min at room temperature in the dark. Protein digestion was carried out using a filter-aided sample preparation technique method. Proteins were loaded onto 30-kDa filter devices (Pall Life Sciences). Trypsin (Trypsin Gold, mass spec grade, Promega Corporation) was added at an enzyme to protein ratio of 1:50, and samples were incubated at 37°C overnight.

ii) ESI-LC-MS/MS for proteome library generation: The pooled peptide samples (30 µl) of each group were separated by high-pH RPLC columns (4.6×250 mm, C18, 3 µm; Waters Corporation). Each pooled sample was loaded onto a column in buffer A1 (‘H₂O’, pH 10). The elution gradient was between 5–30% buffer B1 (90% ACN, pH 10; flow rate, 1 ml/min) for 30 min. Eluted peptides were collected at one fraction per min. After lyophilization, 30 fractions were resuspended in 0.1% formic acid, before being concatenated into 10 fractions by combining fractions 1, 11, 21 and so forth.

In order to generate a spectral library, fractions from RPLC were analyzed in the DDA mode. Parameters were set as follows: the MS was recorded at 350-1,500 m/z at a resolution of 60,000 m/z; the maximum injection time was 50 ms, the auto gain control (AGC) was 1e6, and the cycle time was 3 sec. MS/MS scans were performed at a resolution of 15,000 with an isolation window of 1.6 Da and a collision energy at 32% (HCD); the AGC target was 50,000, and the maximum injection time was 30 ms.

iii) ESI-LC-MS/MS for proteome data-independent acquisition analysis: Digested peptides were dissolved in 0.1% formic acid and separated on an RP C18 self-packing capillary LC column (75 µm × 150 mm, 3 µm). The eluted gradient was 5–30% buffer B2 (0.1% formic acid, 99.9% ACN; flow rate, 0.3 µl/min) for 60 min. For MS acquisition, the variable isolation window DIA method with 38 windows was developed. The specific window lists were constructed based on the DDA experiment of the pooled sample. The full scan was set at a resolution of 120,000 over the m/z range of 400 to 900, followed by DIA scans with a resolution of 30,000; the HCD collision energy was 32%, the AGC target was 1E6 and the maximal injection time was 50 ms.

iv) Spectral library generation: To generate a comprehensive spectral library, a pooled sample from each group was processed. DDA data were processed using Proteome Discoverer software (Thermo Fisher Scientific, Inc.) and searched against the human UniProt database (https://sparql.uniprot.org/) appended with the iRT fusion protein sequence (Biognosys AG).

A maximum of two missed cleavages for trypsin were used, cysteine carbamidomethylation was set as a fixed modification, and methionine oxidation deamination and +43 on Kn (Carbamyl) were used as variable modifications. Parent and fragment ion mass tolerances were set to 10 ppm and 0.02 Da, respectively. The applied false discovery rate (FDR) cutoff was 0.01 at the protein level. The results were then imported to Spectronaut Pulsar software (Biognosys AG) to generate the spectral library.

Additionally, DIA data were imported into Spectronaut Pulsar software and searched against the human UniProt database to generate DIA library. The final library was generated by combining DDA and DIA libraries of ccRCC and controls.

v) Data analysis: DIA-MS data were analyzed using Spectronaut Pulsar (Biognosys AG) with default settings. All results were filtered with a Q-value cutoff of 0.01 which corresponded to an FDR of 1%. Proteins identified in more than 50% of the samples in each group were retained for further analysis. Missing values were imputed based on the k-nearest neighbor method. Raw proteomics data were transformed using log2 and then centralized. Wilcoxon signed rank test was implemented with the software R version 4.1.1 (R-project.org).

Selection of m6A RNA methylation regulators: After integrating RNA sequencing and proteomic data, 19 m6A RNA methylation regulators were included for further analysis and included: ALKBH5, CBLL1, ELAVL1, FMR1, FTO, FXR2, HNRNPA2B1, HNRNPC, METTL 14, METTL 3, RBM15, RBM15B, RBMX, RBMXL1, WTAP, YTHDC1, YTHDC2, YTHDF2 and YTHDF3.

To assess the survival impact of every single m6A methylation regulators on early stage ccRCC, univariate survival analysis was applied. Firstly, pROC package (14) was applied to determine the optimal threshold of every single gene expression and then all the ccRCC patients were divided into two subgroups, namely high and low subgroups. Then, Kaplan-Meier survival analysis with log rank test was calculated and displayed with package ‘survminer’ (15).

Aiming to construct a comprehensive prognostic model, all 19 m6A methylation regulators were included into multivariate Cox's regression analysis with risk scores for each sample being calculated. Prognostic risk scores for ccRCC patients were calculated using the following formula: Risk score=∑ðβi × Exp iÞ, with i representing the number of prognostic m6A expression levels.

According to the median of risk score, all samples were divided into two subgroups: high risk and low risk, and survival curves and time-dependent ROC curves were generated. The Human Protein Atlas (16) was searched to assess the protein expression level of m6A methylation regulators within the prognostic model.

In order to assess the independence of the diagnostic model, certain critical clinical parameters of the patients were included into the multivariate analysis, including sex, age, tumor side, pathological T stage, N stage and tissue grade. Furthermore, these factors were applied using the ‘rms’ package (17) to develop a survival time nomogram.

To assess the relationship between prognostic model and tumor purity, R package ‘estimate’ (18) was used to calculate tumor purity, the scores of stromal cells and the infiltration level of immune cells. To further explore the immune cell composition within the tumor microenvironment, Tumor Immune Estimation Resource (TIMER, http://timer.cistrome.org/) was used to compute the proportion of six kinds of immune cells. Spearman's correlation analysis was used to investigate correlations between risk score and immune infiltrating cells. In addition, associations between clinical characteristics, tumor purity and risk scores were analyzed using the standard Kruskal-Wallis test for multi-group (Dunn's test was used for the post hoc analysis) and Wilcoxon signed rank test for two-group comparisons, respectively.

All the data analysis and the figures were completed by R (version 4.1.1) as well as corresponding R packages illustrated in the Methods. The non-parametric Spearman's rank correlation analysis was used to calculate the correlation coefficient. All the tests were two-sided, and P<0.05 was considered to indicate a statistically significant difference.

Clinical characteristics for 261 TCGA participants and 27 paired proteomic samples have been provided in Table I. At the end of the follow-up period, 44 of 261 patients had succumbed while none of the 27 patients had experienced tumor recurrence or succumbed. The number of men in this cohort was slightly greater than women with a ratio of 158:103 for TCGA and 20:7 in the proteomic sample. In addition, the average age within the two datasets was similar at 59 years for the TCGA group and 57 years for the proteomic sample.

Firstly, the gene mutation pattern of the 19 m6A methylation regulator genes was examined (Fig. 1A). It was identified that the mutation frequency was relatively low with a rate of 4.69%. Among the 9 mutated samples, YTHDC2 demonstrated the highest mutation frequency at 2% and was followed by WTAP with 1%. Furthermore, missense mutation was the most common mutation type with a rate of 3.6%.

Expression levels for RNA sequencing and proteomic data were presented in Fig. 1B and C, respectively. It was revealed that there existed 11 differentially expressed genes, including 4 upregulated and 7 downregulated m6A regulators (Fig. 1B). According to the level of proteins, 15 aberrantly expressed m6A regulators were detected, except for FMR1, FTO, METTLE14 and METTLE3 (Fig. 1C).

In order to identify survival-related m6A methylation regulators, single-variate Cox analysis was performed and the results were provided in a forest plot (Fig. 2A). It could be clearly observed that the expression of FMR 1, METTL14, RBMXL1 and YTHDF 2 were significantly associated with early-stage ccRCC survival. Notably, all 4 regulators were negatively associated with ccRCC overall survival having hazard ratios (HR) of less than 1. After calculating the best cut-off value of each gene with ROC curves, the survival curves for m6A regulators were plotted. The results suggested that 8 gene expression levels, namely RBMXL1, YTHDF3, FMR1, METTLE14, YTHDF2, ELAVL1, FXR2 and RBM 15 were associated with ccRCC patient survival and the top three are revealed in Fig. 2B-D.

In order to comprehensively describe the predictability of m6A methylation regulators, all 19 m6A regulators were included to construct a prognostic model. After applying multivariate Cox regression analysis, a risk score algorithm was established: Risk score=(−0.54) × expr (YTHDF2) + 0.67 × expr (YTHDC2) + (−1.05) × expr (FMR1) + 1.43 × expr (ELAVL1). Then, all 261 patients were classified into two subgroups according to the median risk score, namely high risk and low risk. Through survival analysis plots, it could be clearly observed that a higher risk score was significantly associated with poorer ccRCC survival (Fig. 3A). To assess the YTHDF2, YTHDC2, FMR1 and ELAVL1 expression, it was found that the immunohistochemistry results from the Human Protein Atlas were in accordance with the risk score. Compared with normal kidney tissue, the expression levels of YTHDF2 and FMR1 within tumor cells were low. In contrast to this situation, the YTHDC2 and ELAVL1 expression levels were upregulated, suggesting the unfavorable role in predicting ccRCC survival (Fig. S1).

To further assess the accuracy of the prognostic model, a ROC curve was generated with an AUC of 65.46% (Fig. 3B). To examine the independency of the model, multivariate Cox analysis with several confirmed critical clinical characteristics was performed and the results showed that the model remained an independent indicator of the survival risk (Fig. 3C).

Furthermore, the association between risk scores and clinical characteristics were also assessed (Fig. 4A-F). The results suggested that risk scores are significantly associated with survival status. A cluster heat map was generated to show RNA expression level for all regulators included in the model (Fig. 5).

To predict survival status more comprehensively and easily, a nomogram was constructed after integrating prognostic risk score and clinical information such as sex, age, laterality, pathological grade, T stage and N stage (Fig. 6). It was revealed that the nomogram is suitable for predicting 1-year, 2-year, 3-year and 5-year survival for ccRCC patients.

m6A methylation regulators participate in the regulation of tumor microenvironment immune status during tumorigenesis (6). After comparing the tumor purity within different risk subgroups, it was found that the higher risk group was significantly associated with higher immune score (P<0.001) and lower tumor purity (P=0.0066). However, the stromal score within high and low risk groups did not show significant difference (P=0.68) (Fig. 7A-C).

To further confirm which kind of immune cells maybe associated with the prognostic model, Spearman's correlation analysis was implemented to assess the association between risk scores and immune cell composition (Fig. 7D-I). The results revealed that risk scores were statistically associated with various immune cells, including dendritic cells, CD4⁺ T cells, CD8⁺ T cells and macrophages as opposed to B cells and neutrophils.