The mRNA (RNA sequencing) Fragments Per Kilobase of transcript per Million Fragments standardized expression data and corresponding clinicopathological features were retrieved for 489 ccRCC tissues and 72 adjacent non-tumor tissues from TCGA (http://cancergenome.nih.gov/) and 91 ccRCC tissues from the International Cancer Genome Consortium (ICGC; http://icgc.org/). Patients without prognostic information were excluded from analysis. Overall survival (OS) and disease-free survival (DFS) were the primary end points of this study. OS was defined as the time between diagnosis and death or was censored at the last follow-up. DFS was defined as the time between diagnosis or surgery and the recurrence of disease or death (for any reason) at the last follow-up.

A total of 19 m6A RNA methylation regulators that were available in TCGA and ICGC cohorts were selected. The expression of 19 m6A RNA methylation regulators was systematically compared in 539 ccRCC tissues with different pathological features and with 72 adjacent non-tumor tissues.

To identify the optimal molecular subgroups, the present study divided patients with ccRCC into different clusters by applying the ‘ConsensusClusterPlus’ package in R v.1.50.0 (http://www.bioconductor.org/packages/release/bioc/html/ConsensusClusterPlus.html) based on the expression of the 19 m6A RNA methylation regulators. The resampling program was used to sample 80% of the samples 50 times; the similarity distance between samples was estimated by the Euclidean distance (25), and ‘km dist’ was used as the clustering algorithm to select the reliable and stable subgroup classification. The criteria to determine the optimal number of clusters were a relatively high consistency between clusters and no appreciable rise in the area under the cumulative distribution function curve. The ‘proportion of ambiguous clustering’ was used to select an optimal value of clusters (k value) (26). Venn online software (http://bioinformatics.psb.ugent.be/webtools/Venn/) was used to identify the overlapping differentially expressed genes (DEGs) between clusters. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment and Gene Ontology (GO) analyses were performed and visualized using the Database for Annotation, Visualization and Integrated Discovery software (DAVID; version 6.7; http://david.abcc.ncifcrf.gov/home.jsp) to provide comprehensive pathway interpretations and functional annotation of DEGs (|log₂FC|>2 and adjusted P value <0.05) between different clusters. Search Tool for the Retrieval of Interacting Genes software (STRING; http://string-db.org/) was used to analyze the interactions and evaluate the level of interactions (including interactions determined by experiments or obtained from curated databases, and interactions determined by text mining or co-expression analyses) among the m6A RNA methylation regulators.

Univariate Cox regression analysis was used to identify the prognostic value of the expression of the m6A RNA methylation regulators in the training cohort (TCGA). Regulators associated with OS in univariate analyses were subsequently selected for least absolute shrinkage and selection operator (LASSO) Cox regression to construct a m6A-related risk signature for clinical prognosis (27). As a result, seven m6A RNA methylation regulators with their corresponding coefficients were determined by the minimum mean cross-validated error, choosing the optimal penalty parameter λ related to the minimum 10-fold cross validation within the training set. The risk score of each patient with ccRCC in the training and validation (ICGC) cohorts was calculated using the following formula:

Where x_i is the standardized expression value of each selected m6A RNA methylation regulator, and Coef_i is the corresponding coefficient of the gene. All patients were divided into low- and high-risk groups based on the median value of the risk scores in the training and validation cohorts.

The human RCC cell lines SW839, SN12C, 786-O and OSRC-2 and human normal renal tubular epithelial cell line HK-2 were obtained from Cell Bank of the Chinese Academy of Sciences. The human RCC cell lines were cultured in RPMI-1640 (Gibco; Thermo Fisher Scientific, Inc.) with 10% FBS (HyClone; GE Healthcare Life Sciences) and 1% penicillin/streptomycin (P/S; Gibco; Thermo Fisher Scientific, Inc.). HK-2 cells were cultured in Keratinocyte Medium (ScienCell Research Laboratories, Inc.) with 1% Keratinocyte Growth Supplement (ScienCell Research Laboratories, Inc.) and 1% P/S (ScienCell Research Laboratories, Inc.). All cells were cultured at 37°C in 5% CO₂.

Total RNA was isolated from cells using TRIzol^® reagent (Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions. Reverse transcription was performed using SuperScript III Reverse Transcriptase (Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions, and qPCR was performed using a SYBR^® Green PCR Master Mix (Applied Biosystems; Thermo Fisher Scientific, Inc.) and a 7900HT Fast Real-Time PCR System (Applied Biosystems; Thermo Fisher Scientific, Inc.). The thermocycling conditions were as follows: 95°C for 10 min, followed by 40 cycles of 95°C for 10 sec and at 60°C for 1 min. Relative mRNA levels were normalized against β-actin. Data were analyzed using the 2^−ΔΔCq method (28). Primer sequences are presented in Table SI.

Total m6A RNA modification was detected in 200 ng aliquots of total RNA extracted from cells using the EpiQuik m6A RNA Methylation Quantification Kit (EpiGentek Group, Inc.) according to the manufacturer's instructions. Briefly, total RNA was bound to wells using the RNA Binding Solution. m6A RNA modification was detected using capture and detection antibodies. The detected signal was enhanced and quantified by reading absorbance at 450 nm using a microplate reader (BioTek Instruments, Inc.). The amount of m6A RNA modification was proportional to the OD intensity. The experiments were performed in triplicate.

The differences in the expression levels of m6A RNA methylation regulators between ccRCC tissues and adjacent non-tumor tissues were compared using Wilcoxon test. Differences in the expression levels of m6A RNA methylation regulators and m6A RNA modification levels between HK-2 and human RCC cell lines were analyzed using one-way ANOVA followed by Dunnett's post-hoc test. The distributions of age, sex, histological grade and TNM stage between clusters and between risk subgroups were analyzed using the Kruskal-Wallis test and the Chi-square test, respectively. One-way ANOVA was performed to compare the risk scores in ccRCC with different T stages. Student's t-test was employed to compare the risk scores in patients grouped by binary clinical variables. Univariate and multivariate Cox regression analyses were performed to identify the prognostic value of the risk score and other clinicopathological features. Comparisons of survival between the clusters or risk groups were performed using Kaplan-Meier curves with the log-rank test, followed by pairwise comparisons between clusters using the ‘survminer’ R package v.0.4.6 with Bonferroni correction. Receiver operating characteristic curves (ROCs) were used to test the prediction accuracy of the risk score and clinicopathological features for 5-year survival and recurrence. Variables significant in the multivariate Cox regression analyses were used to construct prognostic nomograms validated by the concordance index (c-index) and calibration plots using the ‘Rms’ package (https://cran.r-project.org/web/packages/rms/index.html) in R. SPSS 22.0 (IBM Corp.), GraphPad Prism 6 (GraphPad Software, Inc.) and R v3.6.1 (https://www.r-project.org/) were used for statistical analysis. P<0.05 was considered to indicate a statistically significant difference, with a Bonferroni correction (p<0.05/3) applied for pairwise comparisons.

Among 19 m6A RNA methylation regulators, 15 DEGs were identified, including 9 upregulated and 6 downregulated genes in 539 ccRCC tissues compared with 72 adjacent non-tumor tissues in TCGA cohort (Fig. 1A), which indicated that m6A RNA methylation regulators exerted important biological functions in the tumorigenesis of ccRCC. The main clinicopathological characteristics of patients with ccRCC are presented in Table I. In the present study, G1-G2 was defined as low histological grade, G3-G4 was defined as high histological grade, stage I–II was defined as low pathological stage, and stage III–IV was defined as high pathological stage. The associations between the expression levels of m6A regulatory regulators and the clinicopathological features of ccRCC, including histological grade and pathological stage, were analyzed in TCGA cohort. The results demonstrated that the expression of most m6A RNA methylation regulators was significantly associated with histological grade (Fig. 1B) and pathological stage (Fig. 1C). In addition, different expression levels of FTO, IGF2BP2, IGF2BP3, KIAA1429, YTHDC1 and ZC3H13 were identified by quantitative analyses according to histological grade or pathological stage (Fig. 1D and E). These results suggested that the m6A RNA methylation regulators contributed to the malignant progression of ccRCC.

The network and correlation analyses among 19 m6A RNA methylation regulators were conducted to determine their interactions (Fig. 2A and B). ‘Writers’ had a wide range of interactions, including interactions determined by experiments or obtained from curated databases, with the m6A RNA methylation regulators. In addition, the ‘writers’ METTL14, KIAA1429, WTAP, RBM15 and ZC3H13 were co-expressed with more than half of the 19 m6A RNA methylation regulators. Based on these results, ‘writers’ were considered to be the hub gene set in the network.

WTAP was the only ‘writer’ with various known interactions with the other five ‘writers’, revealing that it may be a hub gene of the ‘writers’. WTAP expression was also positively associated with the ‘writers’ KIAA1429, ZC3H13 and RBM15 in ccRCC (Fig. 2B). In addition, the expressions levels of KIAA1429, YTHDF3, RBM15, YTHDC2, FTO, ZC3H13, YTHDC1, METTL14 and YTHDF2 were positively associated with each other.

By contrast, the interactions between ‘erasers’ and other regulators were identified by text mining or co-expression analyses, but lacked known interactions. With the exception of FTO, the ‘erasers’ ALKBH3 and ALHBH5 lacked co-expression with other m6A RNA methylation regulators. Independent interaction groups were identified within the ‘readers’, indicating the diverse functions of the ‘readers’; YTHDC1, YTHDC2, YTHDF1, IGF2BP2 and HNRNPA2B1 exhibited co-expression with each other in ccRCC.

According to the ‘proportion of ambiguous clustering’ method and the criteria for selecting the number of clusters, k=3 was selected as the optimal value in TCGA cohorts (Fig. 2C-E). Thus, the patients from TCGA cohort were divided into three clusters, namely cluster 1, 2 and 3. Differences in clinicopathological features and survival distributions between different clusters were identified, with the expression values of the m6A RNA methylation regulators screened out as heatmaps (Fig. 2F). Cluster 1/2/3 subgroups were significantly associated with sex (P<0.001), M stage (P=0.048), N stage (P=0.039), pathological stage (P=0.021) and survival outcome (P<0.001) (Table II). Additionally, patients in cluster 1 exhibited the shortest OS (P=0.02), followed by cluster 3 and cluster 2 (Fig. 2G). These results revealed that patients with ccRCC in different clusters exhibited significantly diverse clinicopathological features, and the clustering results were associated with survival.

To identify different biological processes among the three clusters, DEGs were identified among the clusters, and their pathway interpretations and functional annotations were visualized. As pathological features between cluster 2 and cluster 3 were similar, and patients with ccRCC in cluster 1 exhibited the worst survival outcomes among all clusters and largely different clinicopathological features from cluster 2 and cluster 3, DEGs between cluster 1 and cluster 2 and between cluster 1 and cluster 3 were investigated. A total of 1,386 DEGs (736 upregulated and 650 downregulated) were identified between cluster 1 and cluster 2, and 2,287 DEGs (1,756 upregulated and 531 downregulated) were identified between cluster 1 and cluster 3. Venn online software was used to identify the overlapping DEGs between the two groups of DEGs (Fig. 3A and B). As a result, 603 upregulated and 149 downregulated overlapping DEGs were selected and analyzed by DAVID software to perform GO enrichment and KEGG pathway analyses.

The results of GO biological processes analysis demonstrated that upregulated and downregulated overlapping DEGs were mainly enriched in RNA metabolism, including ‘RNA export from the nucleus’, ‘mRNA export from the nucleus’, ‘RNA processing’, ‘regulation of RNA splicing’, ‘mRNA splice site selection’ and ‘translation’ (Fig. 3C and E). In KEGG pathway analysis, similar changes were observed in the corresponding signaling pathways, including ‘spliceosome’, ‘ribosome’, ‘transcriptional misregulation in cancer’, ‘mRNA surveillance pathway’ and other malignancy-related pathways including ‘primary immunodeficiency’, ‘regulation of autophagy’ and ‘response to oxidative stress’ (Fig. 3D and F). These results suggested that cluster 1/2/3 subgroups were associated not only with the clinicopathological features and survival, but also with malignancy-related biological processes.

The prognostic significance of the m6A RNA methylation regulators in patients with ccRCC was further analyzed. In TCGA cohort, univariate Cox regression analysis identified 14 and 11 m6A RNA methylation regulators significantly associated with OS and DFS, respectively (Fig. 4A and B). Among the 14 m6A RNA methylation regulators in OS, six regulators were associated with poor OS (IGF2BP1, IGF2BP2, IGF2BP3, HNRNPA2B1, METTL3 and ALKBH3), and eight were associated with favorable OS (YTHDF2, YTHDF3, FTO, YTHDC1, YTHDC2, ZC3H13, KIAA1429 and METTL14). For DFS, two regulators were associated with poor DFS (IGF2BP2 and IGF2BP3), and nine regulators were associated with favorable DFS (YTHDF1, YTHDF2, YTHDF3, YTHDC1, YTHDC2, ZC3H13, KIAA1429, METTL14 and RBM15). The prognostic values of the m6A RNA methylation regulators in patients with different histological grades and pathological stages in TCGA cohort were subsequently analyzed by univariate Cox regression. In both low and high histological grade ccRCC, IGF2BP3 (P=0.030 and P<0.001, respectively) and KIAA1429 (P=0.045 and P=0.031, respectively) were significantly associated with OS; in both low and high pathological stage, METTL14 was significantly associated with OS (P=0.035 and P=0.003, respectively) (Fig. S1).

Since the OS data was available in both TCGA and ICGC cohorts, the OS-associated m6A RNA methylation regulators in TCGA cohort were used to construct a m6A-related risk signature, and its accuracy was validated in the ICGC cohort. As a result, seven genes (IGF2BP2, IGF2BP3, METTL3, METTL14, HNRNPA2B1, KIAA1429 and ALKBH3) were identified by LASSO Cox regression in the training cohort (TCGA) to build the m6A-related risk signature (Fig. 4C and D). Subsequently, the risk score of each patient with ccRCC in the two cohorts was calculated using these seven genes and the formula presented in Materials and methods. Patients with ccRCC from the two cohorts were divided into low- and high-risk subgroups based on the median risk score. In TCGA cohort, the Kaplan-Meier method revealed that patients with high risk scores exhibited shorter OS (P<0.001; Fig. 4F) and DFS (P<0.001, Fig. 4G) compared with those with low risk scores. In the ICGC cohort, patients with high risk scores also exhibited shorter OS (P=0.005; Fig. 4H) compared with those with low risk scores.

Survival analyses were performed to compare the two risk groups of patients with different histological grades and pathological stages. In TCGA cohort, patients in the high-risk group exhibited poor OS and DFS compared with patients in the low-risk group regardless of histological grade and pathological stage (Fig. S2A-H). In the ICGC cohort, patients with high risk scores exhibited shorter OS compared with patients with low risk scores with a low pathological stage, but not a high pathological stage (P=0.001 and P=0.549, respectively; Fig. S2I and J). ICGC data of histological grade and DFS were unavailable.

Based on these results, the risk signature derived from the seven m6A regulators was identified to be a powerful prognostic tool for patients with ccRCC. Further research is needed to evaluated whether the prognostic value of the risk signature was also independent of the known prognostic factors, such as histological grade and pathological stage.

The heat map of the expression levels of the seven selected m6A RNA methylation regulators in the high- and low-risk subgroup patients in TCGA cohort is presented in Fig. 5A and B. Patients in the high-risk group had a higher proportion of males (P=0.015), higher pathological stage (P<0.001), higher histological grade (P<0.001), higher T stage (P<0.001), higher M stage (P<0.001), higher N stage (P=0.012) and poorer OS compared with patients in the low-risk group (Table III). The levels of risk score according to different clinicopathological features were compared; significantly different risk scores were determined between patients stratified by pathological stage, histological grade, T stage, N stage and M stage in TCGA cohort (all P<0.001; Fig. 5C-H).

ROCs were used to determine the prognostic accuracy of the risk score, histological grade and pathological stage in patients with ccRCC. The results revealed that in TCGA cohort, the accuracy of the risk score [area under the curve (AUC), 0.736] was superior compared with those of histological grade (AUC, 0.681) and pathological stage (AUC, 0.720) in predicting 5-year survival (Fig. 5I). In predicting 5-year recurrence, the risk score (AUC, 0.728) was superior compared with histological grade (AUC, 0.722) but inferior to pathological stage (AUC, 0.820) (Fig. 5J). The combination of these three characteristics improved the accuracy of predicting the 5-year survival and recurrence of the same cohort (AUC, 0.782 and 0.859, respectively). Therefore, the risk score was able to correctly predict the prognosis for patients with ccRCC.

In TCGA cohort, univariate Cox regression analyses demonstrated that the risk score, pathological stage, histological grade, T, M and N stage had prognostic value in OS and DFS, whereas age was associated exclusively with OS (Fig. 6A and B). These factors were used for multivariate Cox regression analysis. The risk score, histological grade and pathological stage remained significantly associated with OS and DFS. Additionally, age, M and N stage were independent prognostic factors of OS, whereas T stage was an independent prognostic factor of DFS (Fig. 6A and B). Similar results were observed in the validation cohort, where the risk score (P=0.014) and pathological stage (P<0.001) were associated with OS in multivariate analyses (Fig. 6C).

Nomograms for 3- and 5-year OS and DFS were constructed based on the m6A-related signature and other independent prognostic factors in multivariate Cox regression analyses (Fig. 7A and B). The c-indices of nomograms were 0.783±0.018 (mean ± SEM) and 0.819±0.018 for OS and DFS, respectively, indicating that the prognostic prediction of nomograms was largely consistent with the actual OS and DFS in patients with ccRCC. Additionally, calibration plots demonstrated that nomogram prediction exhibited an agreement with actual 3-year OS and DFS (Fig. 7C and E) and a relative agreement with actual 5-year OS and DFS (Fig. 7D and F).

RT-qPCR was used to validate the expression levels of the seven selected m6A RNA methylation regulators in four human RCC cell lines (SW839, SN12C, 786-O and OSRC-2) and one human normal renal tubular epithelial cell line (HK-2). The results demonstrated significant differences in the expression levels of five m6A RNA methylation regulators (KIAA1429, ALKBH3, HNRNPA2B1, IGF2BP2 and METTL3) between RCC and normal renal cells (Fig. 8A). In addition, m6A RNA modification levels in the cell lines were detected. The m6A RNA modification levels in RCC cell lines SW839 and 786-O were significantly higher compared with that in HK-2 cells (Fig. 8B).