A novel m a gene signature associated with regulatory immune function for prognosis prediction in clear-cell renal cell carcinoma

Introduction

It is estimated that there will be 73, 750 new cases of renal malignant tumors in the United States in 2020 ( Siegel et al., 2020 ). As one of the most aggressive malignancies, renal cell carcinoma (RCC) accounts for 2–3% of the malignancies ( Ljungberg et al., 2015 ). Clinically, clear-cell RCC (ccRCC) is the most common type of renal carcinoma, representing ~80% of RCC. Although the therapy methods of ccRCC, including the surgical and targeted therapy, have been improved, the poor survival prognosis is also noticed in ccRCC, which is obligated to most metastatic cases of RCC ( Reuter, 2006 ). Even for the localized ccRCC, the tumor recurrence and progression can also be observed in ~25% of patients after primary treatment ( De et al., 2014 ). Therefore, there is still an urgent need to explore effective prognostic markers and fully clarify the molecular mechanism underlying the tumorigenesis of ccRCC.

Currently, the risk stratification and prognostic prediction for the ccRCC patients are mainly conducted using the Fuhrman grade and tumor–node–metastasis stage system. Even patients with similar clinicopathologic characteristics could also suffer from variable survival outcomes. Taking into account the heterogeneity in ccRCC, several potential biomarkers for diagnosis and survival prediction have been reported in recent years, such as carbonic anhydrase 9 ( Genega et al., 2010 ), PBRM1, and BAP1 mutation ( Varela et al., 2011 ). However, few of these potential biomarkers could finally be transformed into clinical practice due to the unsatisfactory accuracy and sensibility.

The modification and epigenetic alteration of RNA have recently been primarily identified, including N6-methyladenosine (m ⁶ A) ( Roundtree et al., 2017 ; Boccaletto et al., 2018 ). M ⁶ A RNA methylation acts as one kind of conserved internal modifications in eukaryotic nuclear RNAs ( Dubin and Taylor, 1975 ), which is one of the most common RNA modifications and plays diverse roles in various biological processes ( Bokar et al., 1997 ; Dominissini et al., 2012 ). The m ⁶ A status in a cell is regulated by three types of m ⁶ A RNA methylation genes, including methyltransferases called writers (RBM15, METTL14, KIAA1429, and others), readers (FMR1, IGF2BP1, YTHDC1, and others), and demethylases called erasers (ALKBH5 and FTO) ( Li A. et al., 2017 ). At present, the low expression levels of FTO and ALKBH5 have been proved to be associated with worse cancer-specific survival of RCC patients after nephrectomy ( Strick et al., 2020 ). However, the specific role of m ⁶ A regulators in carcinogenesis and progression of ccRCC has not been fully understood by now.

In this study, we carried out a comprehensive analysis based on the expressions of m ⁶ A RNA methylation regulator to identify their important roles in the diagnosis and survival prediction for the ccRCC patients.

Materials and Methods

Datasets Source and Screen

Data sets from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) were strictly screened for this study. The inclusion criteria were as follows: studies with more than 60 ccRCC patients and studies with complete clinical data and gene expression data of m ⁶ A RNA methylation genes. Five GEO datasets, including GSE46699 ( Eckel-Passow et al., 2014 ), GSE53757 ( von Roemeling et al., 2014 ), GSE22541 ( Wuttig et al., 2012 ), GSE17895 ( Dalgliesh et al., 2010 ), and GSE40435 ( Wozniak et al., 2013 ), were finally selected and downloaded with gene expression data and clinical data. Four datasets (GSE46699, GSE53757, GSE22541, and GSE17895) with microarray data were detected based on the Affymetrix Human Genome U133 Plus 2. 0 Array. We firstly performed background adjustment and quantile normalization from the raw CEL data of the four datasets by the Robust Multi-array Average in the R environment ( Irizarry et al., 2003 ). And then, we annotated the probe of each gene by using the hgu133plus2. db in R. Finally, the processed microarray data from the four datasets were merged as the Meta-cohort after batch normalization. The processed expression data of GSE40435 based on the Illumina HumanHT-12 V4. 0 was directly downloaded from the GEO database, as it had been normalized in a previous study ( Wozniak et al., 2013 ), which was defined as the GSE-cohort. TCGA-KIRC data set, including normalized RNA-sequencing data and clinical data, was directly downloaded from the TCGA database. Expression data files of 21 m ⁶ A RNA methylation genes were further extracted from the processed data sets mentioned earlier. Basic clinical characteristics of data sets in this study are shown in Supplementary Table 1 .

Construction of the Diagnostic and Prognosis Model Based on the m ⁶ A Regulators for Clear-Cell Renal Cell Carcinoma

Differential expression analysis of the 21 m ⁶ A RNA methylation genes between normal and tumor samples was performed in the Meta-cohort and the GSE-cohort. Only m ⁶ A RNA methylation genes differentially expressing in both the two cohorts were further used for the construction of the diagnostic model. We carried out the least absolute shrinkage and selection operator (LASSO) analysis via the glmnet package to identify ccRCC-related m ⁶ A regulators and calculate coefficients of each gene in the Meta-cohort. The number of lambda values in LASSO was set as 1, 000 to ensure the robustness of our diagnostic models. The m ⁶ A diagnostic score was then calculated as follow:

$\begin{array}{ }\text{m}6\text{A diagnostic score}={\sum }_{i=1}^n\left(\text{Coefi}*\text{Di}\right)\end{array}$

Coefi means the coefficient of each ccRCC-related m ⁶ A RNA methylation gene, whereas Di represents the related messenger RNA (mRNA) expression. The diagnostic model was further verified in two independent cohorts (GSE-cohort and TCGA-cohort).

To construct the prognosis model based on the m ⁶ A regulators for ccRCC, 531 ccRCC patients from the TCGA cohort were randomly assigned to a training cohort and a testing cohort by the random number method. We performed a LASSO Cox regression analysis to identified survival-related m ⁶ A RNA methylation genes and calculated their coefficients in the TCGA training cohort. The m ⁶ A risk score was then calculated as follow:

$\begin{array}{ }\text{m}6\text{A risk score}={\sum }_{i=1}^n\left(\text{Coefi}*\text{Ri}\right)\end{array}$

Coefi means the coefficient of each survival-related m ⁶ A RNA methylation gene, whereas Ri represents the related mRNA expression. The prognosis model was further validated in the TCGA testing cohort and the independent validation cohort (GSE22541). The cutoff value of high and low m ⁶ A risk score was defined as the median score of respective cohorts.

Development and Evaluation of the Prognostic Nomogram Model

We developed a prognostic nomogram model based on m ⁶ A risk score, age, tumor stage, and tumor grade using the rms and nomogramEx packages in the R environment. The calibration with bootstrapping was conducted to verify the nomogram-predicted probabilities of the 1-, 3-, and 5-year overall survival (OS) via plotting on the x-axis, with actual OS probabilities plotting on the y-axis. The time-dependent receiver operating characteristic (ROC) curve analysis and the decision curve analysis were also performed to identify the specificity, sensitivity, and clinical utility of the prognostic model.

Weighted Gene Co-expression Network and Functional Enrichment Analysis

The significantly differentially expressed genes (DEGs) were firstly identified between the ccRCC tissue and the normal renal tissue from GSE40435, GSE53757, and GSE46699. The weighted gene co-expression network analysis (WGCNA) was then performed based on the DEGs by the WGCNA package in R ( Langfelder and Horvath, 2008 ). When the soft-thresholding power of β value was defined as 8, the recruited DEGs were hierarchically clustered into seven gene modules. The correlation between gene module and clinical characteristic was further analyzed to identify the optimal gene module with the highest correlation. We finally conducted the Kyoto Encyclopedia of Genes and Genomes pathway and the Genetic Ontology analysis by the Metascape ( Zhou et al., 2019 ) to explore the potential biological mechanisms in which the m ⁶ A risk score might be involved in.

Statistical Analysis

The continuous variable with abnormal distribution was analyzed by non-parametric tests (Mann–Whitney U test for comparison between two groups and Kruskal–Wallis test for comparisons among more than two groups). The log-rank test was used for analyzing the Kaplan–Meier curves of OS and disease-free survival. The Cox regression analyses were carried out to identify the m ⁶ A risk score as a prognostic factor of OS. R (3. 6. 2) and Statistical Package for Social Sciences 24. 0 software (SPSS Inc., Chicago, IL, USA) were used for statistical analysis and data visualization. A two-tailed P -value of <0. 05 was considered significant.

Results

Important Roles of m ⁶ A Regulators in Clear-Cell Renal Cell Carcinoma

We explored and verified 21 m ⁶ A regulators with differential mRNA expressions between the tumor samples and normal samples in the Meta-cohort and the GSE-cohort ( Figure 1A ). Only 11 m ⁶ A regulators were significantly differentially expressed in both the two cohorts, including five writers (RBM15, RBM15B, WTAP, KIAA1429, and CBLL1), one eraser (FTO), and five readers (YTHDC1, YTHDC2, YTHDF2, IGF2BP1, and HNRNPA2B1). The 11 m ⁶ A regulators were further analyzed by the LASSO analysis in the Meta-cohort ( Figure 1B ). Finally, eight ccRCC-related m ⁶ A RNA methylation genes, including RBM15, WTAP, CBLL1, FTO, YTHDC1, YTHDC2, YTHDF2, and HNRNPA2B1, were selected for the construction of the m ⁶ A-related diagnostic model ( Figure 1C ). The ROC analysis revealed that the area under the curve (AUC) of the diagnostic model reached 0. 924 [95% confidence interval (CI): 0. 897–0. 946] in Meta-cohort ( Figure 1D ). This diagnostic model was further validated in two independent cohorts, including GSE-cohort and TCGA-cohort. AUC values reached 0. 867 (95% CI: 0. 812–0. 910) in GSE-cohort ( Figure 1E ) and 0. 795 (95% CI: 0. 761–0. 827) in TCGA-cohort ( Figure 1F ), which proved the stability of the diagnostic efficiency of our m ⁶ A-related diagnostic model for the ccRCC patients.

FIGURE 1

Important role of m ⁶ A regulators in ccRCC.(A)Exploring (Meta-cohort) and verifying (GSE-cohort) different mRNA expression of 21 m ⁶ A regulators between normal and tumor samples in ccRCC.(B)10-fold cross-validated error for LASSO analysis. First vertical line equals minimum error, whereas second vertical line shows cross-validated error within one standard error of minimum.(C)Profile of coefficients in model at varying levels of penalization plotted against log (lambda) sequence.(D)ROC curve of diagnostic score in training cohort.(E, F)ROC curves of diagnostic score in validation cohorts (GSE-cohort and TCGA-cohort). m ⁶ A, N ⁶ -methyladenosine; ccRCC, clear-cell renal cell carcinoma; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

Effective Survival Prediction of Clear-Cell Renal Cell Carcinoma Patients by the m ⁶ A Risk Score

As shown in Figure 2A , the first vertical line pointed at 7 in the TCGA training cohort, which equaled the minimum 10-fold cross-validated error, indicating that seven optional prognosis-related m ⁶ A regulators, including WTAP, ZC7H13, FTO, ELAVL1, HNRNPA2B1, HNRNPC, and IGF2BP1, were selected by the LASSO Cox regression analysis ( Figure 2B ). A significant difference was observed in the OS [hazard ratio (HR) = 3. 474, 95% CI: 2. 260–5. 342, p < 0. 0001] between the ccRCC patients with high and low m ⁶ A risk score in the TCGA training cohort ( Figure 2C ), which was further verified in the TCGA testing cohort ( Figure 2D ) with HR value of 1. 679 (95% CI: 1. 113–2. 532, p = 0. 015) for the OS and HR value of 2. 101 (95% CI: 1. 186–3. 724, p = 0. 001) for the disease-free survival in the independent GSE22541 cohort ( Figure 2E ). A higher m ⁶ A risk score was correlated with the higher tumor stage, higher tumor grade, and death ( Figure 2F ).

FIGURE 2

Prognosis model based on m ⁶ A regulators for ccRCC.(A)10-fold cross-validated error for LASSO analysis. First vertical line equals minimum error, whereas second vertical line shows cross-validated error within one standard error of minimum.(B)Profile of coefficients in model at varying levels of penalization plotted against log (lambda) sequence.(C–E)Kaplan–Meier survival analysis of overall survival stratified by m ⁶ A risk score for ccRCC patients in TCGA training cohort, TCGA validation cohort, and disease-free survival in independent validation cohort (GSE22541).(F)Heatmap summarized expression of seven m ⁶ A methylation genes selected by LASSO Cox regression analysis and distribution of clinicopathologic factors in whole TCGA cohort. Cutoff value of high and low m ⁶ A risk score was defined as median score of respective cohorts. m ⁶ A, N ⁶ -methyladenosine; ccRCC, clear-cell renal cell carcinoma; LASSO, least absolute shrinkage and selection operator; TCGA, The Cancer Genome Atlas.

The m ⁶ A risk score was identified as a risk factor for the OS in ccRCC patients by the univariate Cox regression analysis ( Figure 3A ), which was further verified in both the TCGA training cohort ( Figure 3B ) and the independent GSE22541 cohort ( Figure 3C ). There existed significant differences of the m ⁶ A risk score among different tumor grades ( P < 0. 0001, Figure 3D ), different tumor stages ( P < 0. 0001, Figure 3E ), different tumor lymph node metastasis status ( P = 0. 06, Figure 3F ), and different tumor distant metastasis status ( P < 0. 0001, Figure 3G ).

FIGURE 3

Evaluation of prognosis model based on m ⁶ A regulators.(A–C)Univariate Cox regression analyses of m ⁶ A risk score and clinicopathologic factors in TCGA training cohort, TCGA validation cohort, and independent validation cohort.(D–G)Difference in distribution of m ⁶ A risk score of different tumor stages, tumor grades, metastasis status, and lymph node metastasis status in whole TCGA cohort. m ⁶ A, N ⁶ -methyladenosine; TCGA, The Cancer Genome Atlas.

Prognostic Accuracy of the m ⁶ A Risk Score Integrated With the Clinicopathologic Factors

To improve the accuracy of our survival prognostic model, we developed an integrated nomogram by combining the m ⁶ A risk score and predictable clinicopathologic factors, including age, tumor grade, and tumor stage. The integration nomograms for the OS prediction are shown in Figure 4A . The calibration plots revealed that the 1-, 3-, and 5-year OS probabilities predicted by the integrated nomogram model had an excellent agreement with the actual observations ( Figure 4B ), indicating good performance in predicting the survival status of the ccRCC patients. The decision curve analysis illustrated that when the threshold probability was more than 0. 2, the integrated nomogram for the OS prediction could be more favorable than either the m ⁶ A risk score and the predictable clinicopathologic factors ( Figure 4C ). Further time-dependent ROC curve revealed that the AUC values of the integration nomogram for the OS prediction in 1, 3, and 5 years arrived at 85. 2, 82. 4, and 78. 3%, respectively ( Figure 4D ).

FIGURE 4

Construction and evaluation of nomogram prognosis model in TCGA cohort.(A)Nomogram based on m ⁶ A risk score and clinicopathologic factors for 1-, 3-, and 5-year overall survival prediction of ccRCC patients.(B)Evolution of prognostic nomogram model for 1-, 3-, and 5-year overall survival prediction.(C)Decision curve analysis compared overall survival benefits among nomogram, tumor stage, tumor grade, and m ⁶ A risk score. Black line: all victims dead. Gray line: none victims dead.(D)Time-dependent ROC curve of 1-, 3-, and 5-year overall survival prediction for prognostic nomogram model. m ⁶ A, N ⁶ -methyladenosine; TCGA, The Cancer Genome Atlas; ccRCC, clear-cell renal cell carcinoma; ROC, receiver operating characteristic.

Immune-Related Pathways Were Associated With the m ⁶ A RNA Methylation Risk Signature

A total of 1, 152 DEGs were identified in the GSE40435 cohort, GSE53757 cohort, and GSE46699 cohort ( Figure 5A ), including 577 upregulated genes and 570 downregulated genes ( Figure 5B ). The DEGs were hierarchically clustered into seven gene modules by the WGCNA method ( Figure 5C ), and the brown model was found to perform the highest correlation to the m ⁶ A RNA methylation risk signature ( Figure 5D ). Significantly different expressions of 128 genes in the brown model were observed between the patients with high and low m ⁶ A risk scores ( Figure 6A ). As shown in Figure 6B , some immune-related pathways, including the immune response regulating pathway, the cytokine-mediated pathway, and the lymphocyte activation-associated pathway, were enriched in genes from the brown model, suggesting that m ⁶ A RNA methylation might affect clinical prognosis through regulating immune functions of ccRCC patients. Further correlation analysis revealed that the m ⁶ A risk score was significantly correlated with the abundances of T regulatory cells (Tregs) and T follicular helper cells (Tfhs) ( Figure 6C ). Higher abundances of the Tregs and the Tfhs were found in the patients with high m ⁶ A risk scores ( Figures 7A, B ).

FIGURE 5

WGCNA for differentially expressed genes in ccRCC.(A)Volcano plots showing significantly upregulated and downregulated genes in GSE40435, GSE53757, and GSE46699.(B)Venn diagrams presenting intersections of significantly upregulated and downregulated genes among GSE40435, GSE53757, and GSE46699.(C)Gene dendrogram with different colors showing modules identified by WGCNA.(D)Relationship between gene modules and clinical characteristic. WGCNA, weighted gene co-expression network analysis; ccRCC, clear-cell renal cell carcinoma.

FIGURE 6

Potential mechanisms of co-expressed genes associated with m ⁶ A risk score.(A)Heatmap visualizing expressions of co-expressed genes in brown module.(B)Potentially enriched pathways of co-expressed genes in brown module.(C)Correlation analysis of m ⁶ A risk score and abundance of 22 types of immune cells. m ⁶ A, N ⁶ -methyladenosine.

FIGURE 7

Comparing abundance of T regulatory cells(A)and T follicular helper cells(B)between patients with high-risk scores and low-risk scores. Cutoff value of high- and low-risk scores was defined as median score.

Discussion

Due to the high incidence and mortality of ccRCC, the accurate diagnosis and survival prediction of ccRCC patients were urgently needed. Thanks to the innovative developments in high-throughput genetic diagnostic techniques for oncology, there has been a revolutionary improvement in the efficient diagnosis of a malignant tumor at the molecular level ( Cancer Genome Atlas Research Network, 2013 ; Linehan et al., 2016 ; Zehir et al., 2017 ).

In this study, we designed a diagnostic model for ccRCC based on the m ⁶ A regulators. The ccRCC and normal renal tissue could be accurately told out by the m ⁶ A diagnostic score in Meta-cohort, and the result was further verified in the independent GSE-cohort and TCGA-cohort, with the AUC values of 0. 924, 0. 867, and 0. 795, respectively. Effective survival prediction of ccRCC by the m ⁶ A risk score was also identified in the TCGA training cohort and validated in the testing cohort and the independent GSE22541 cohort, with the HR values of 3. 474, 1. 679, and 2. 101 for clinical survival, respectively. It was proved that the m ⁶ A risk score in ccRCC was associated with higher tumor stage, higher tumor grade, and tumor metastasis. All of these results suggested the crucial roles of m ⁶ A regulators in regulating tumorigenesis and tumor progression, as previously reported ( He et al., 2019 ; Wang et al., 2019 ; Tian et al., 2020 ).

An integrated nomogram by combining m ⁶ A risk score and predictable clinicopathologic factors, including age, tumor grade, and tumor stage, was constructed in this study to improve the accuracy of our survival prognostic model. The time-dependent ROC curve revealed that the AUC values of the integration nomogram for OS prediction in 1, 3, and 5 years arrived at 85. 2, 82. 4, and 78. 3%, respectively, indicating advantageous usability of our survival prediction model in clinical practice.

Different types of m ⁶ A RNA methylation regulators work differently in tumorigenesis. For example, FTO could promote the progression of lung carcinoma by releasing the m ⁶ A modification in MZF1 mRNA and strengthening its stability ( Liu et al., 2018 ). However, METTL14 was found to suppress RCC by downregulating P2RX6 protein translation ( Gong et al., 2019 ). Generally, writers could irritate m ⁶ A modifications in the mRNA of tumor suppressor genes or oncogenes. And then, these modifications could be recognized by readers, resulting in downregulating tumor suppressor or upregulating oncogene expression ( Wang et al., 2020 ).

A total of 128 DEGs associated with the m ⁶ A RNA methylation risk signature in the brown model were enriched in immune regulated pathways, including the immune response regulating pathway, the cytokine-mediated pathway, and the lymphocyte activation-related pathway. In addition, the m6A risk score was significantly correlated with the abundance of Tregs and Tfhs. The important role of m ⁶ A modification in immune cell-related pathogenesis through the m ⁶ A mediated degradation of Naïve T cell had been recently identified ( Li H. B. et al., 2017 ). Our enrichment results indicated that the m ⁶ A RNA methylation in ccRCC might affect the prognosis through regulating immune function, especially the functions of Tregs and Tfhs.

There are also several defects in this study. Firstly, our main findings were based on integrated bioinformatics analyses, without experimental verification. The experiment verification at the cellular level, including regulation mechanism research and function verification, are still needed. Secondly, all of the data sets analyzed in this study were acquired from the public database, resulting in a potential bias in genetic and clinical data. However, cross-validation among independent datasets has been performed as much as possible to reduce the potential bias. Thirdly, the cutoff value of the high and low m ⁶ A risk score group was defined as the median score of respective cohorts. Actually, the optional cutoff value of m ⁶ A risk score to distinguish the patients with high survival risk is supposed to be identified and verified in larger clinical patient cohorts. Finally, our work still requires further clinical validation to verify its application to the clinic in ccRCC patients.

Conclusions

In summary, we identified the important role of the m ⁶ A regulator in ccRCC. Immune-related pathways might be involved in the regulation of ccRCC through the m ⁶ A RNA methylation. The novel m ⁶ A gene signature might act as an effective biomarker for prognosis prediction of the ccRCC patients, which still requires further clinical validation and experimental verification.

Data Availability Statement

The original contributions generated for the study are included in the article/ Supplementary Material , further inquiries can be directed to the corresponding authors.

Author Contributions

JZ, XW, and SC conceived the study, designed the research, and wrote the paper. SC, EZ, NZ, and TW conducted and analyzed experiments. LJ and XW provided samples. JZ supervised the research. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the National Natural Science Foundation of China (81772705, 81972393, and 82002665).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank the TCGA, GEO database for gene expression and clinical data.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www. frontiersin. org/articles/10. 3389/fcell. 2020. 616972/full#supplementary-material

References

Boccaletto, P., Machnicka, M. A., Purta, E., Piatkowski, P., Baginski, B., Wirecki, T. K., et al. (2018). MODOMICS: a database of RNA modification pathways. 2017 update. Nucleic Acids Res . 46, 303–307. doi: 10. 1093/nar/gkx1030

Bokar, J. A., Shambaugh, M. E., Polayes, D., Matera, A. G., and Rottman, F. M. (1997). Purification and cDNA cloning of the AdoMet-binding subunit of the human mRNA (N6-adenosine)-methyltransferase. RNA 3, 1233–1247.

Cancer Genome Atlas Research Network (2013). Comprehensive molecular characterization of clear cell renal cell carcinoma. Nature 499, 43–49. doi: 10. 1038/nature12222

Dalgliesh, G. L., Furge, K., Greenman, C., Chen, L., Bignell, G., Butler, A., et al. (2010). Systematic sequencing of renal carcinoma reveals inactivation of histone modifying genes. Nature 463, 360–363. doi: 10. 1038/nature08672

De, P., Otterstatter, M. C., Semenciw, R., Ellison, L. F., Marrett, L. D., and Dryer, D. (2014). Trends in incidence, mortality, and survival for kidney cancer in Canada, 1986-2007. Cancer Causes Control 25, 1271–1281. doi: 10. 1007/s10552-014-0427-x

Dominissini, D., Moshitch-Moshkovitz, S., Schwartz, S., Salmon-Divon, M., Ungar, L., Osenberg, S., et al. (2012). Topology of the human and mouse m6A RNA methylomes revealed by m6A-seq. Nature 485, 201–206. doi: 10. 1038/nature11112

Dubin, D. T., and Taylor, R. H. (1975). The methylation state of poly A-containing messenger RNA from cultured hamster cells. Nucleic Acids Res . 2, 1653–1668. doi: 10. 1093/nar/2. 10. 1653

Eckel-Passow, J. E., Serie, D. J., Bot, B. M., Joseph, R. W., Cheville, J. C., and Parker, A. S. (2014). ANKS1B is a smoking-related molecular alteration in clear cell renal cell carcinoma. BMC Urol . 14: 14. doi: 10. 1186/1471-2490-14-14

Genega, E. M., Ghebremichael, M., Najarian, R., Fu, Y., Wang, Y., Argani, P., et al. (2010). Carbonic anhydrase IX expression in renal neoplasms: correlation with tumor type and grade. Am. J. Clin. Pathol . 134, 873–879. doi: 10. 1309/AJCPPPR57HNJMSLZ

Gong, D., Zhang, J., Chen, Y., Xu, Y., Ma, J., Hu, G., et al. (2019). The m(6)A-suppressed P2RX6 activation promotes renal cancer cells migration and invasion through ATP-induced Ca(2+) influx modulating ERK1/2 phosphorylation and MMP9 signaling pathway. J. Exp. Clin. Cancer Res . 38: 233. doi: 10. 1186/s13046-019-1223-y

He, L., Li, H., Wu, A., Peng, Y., Shu, G., and Yin, G. (2019). Functions of N6-methyladenosine and its role in cancer. Mol. Cancer 18: 176. doi: 10. 1186/s12943-019-1109-9

Irizarry, R. A., Bolstad, B. M., Collin, F., Cope, L. M., Hobbs, B., and Speed, T. P. (2003). Summaries of Affymetrix GeneChip probe level data. Nucleic Acids Res . 31: e15. doi: 10. 1093/nar/gng015

Langfelder, P., and Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinform . 9: 559. doi: 10. 1186/1471-2105-9-559

Li, A., Chen, Y. S., Ping, X. L., Yang, X., Xiao, W., Yang, Y., et al. (2017). Cytoplasmic m(6)A reader YTHDF3 promotes mRNA translation. Cell Res . 27, 444–447. doi: 10. 1038/cr. 2017. 10

Li, H. B., Tong, J., Zhu, S., Batista, P. J., Duffy, E. E., Zhao, J., et al. (2017). m(6)A mRNA methylation controls T cell homeostasis by targeting the IL-7/STAT5/SOCS pathways. Nature 548, 338–342. doi: 10. 1038/nature23450

Linehan, W. M., Spellman, P. T., Ricketts, C. J., Creighton, C. J., Fei, S. S., Davis, C., et al. (2016). Comprehensive molecular characterization of papillary renal-cell carcinoma. N. Engl. J. Med . 374, 135–145. doi: 10. 1056/NEJMoa1505917

Liu, J., Ren, D., Du, Z., Wang, H., Zhang, H., and Jin, Y. (2018). m(6)A demethylase FTO facilitates tumor progression in lung squamous cell carcinoma by regulating MZF1 expression. Biochem. Biophys. Res. Commun . 502, 456–464. doi: 10. 1016/j. bbrc. 2018. 05. 175

Ljungberg, B., Bensalah, K., Canfield, S., Dabestani, S., Hofmann, F., Hora, M., et al. (2015). EAU guidelines on renal cell carcinoma: 2014 update. Eur. Urol . 67, 913–924. doi: 10. 1016/j. eururo. 2015. 01. 005

Reuter, V. E. (2006). The pathology of renal epithelial neoplasms. Semin. Oncol . 33, 534–543. doi: 10. 1053/j. seminoncol. 2006. 06. 009

Roundtree, I. A., Evans, M. E., Pan, T., and He, C. (2017). Dynamic RNA modifications in gene expression regulation. Cell 169, 1187–1200. doi: 10. 1016/j. cell. 2017. 05. 045

Siegel, R. L., Miller, K. D., and Jemal, A. (2020). Cancer statistics, 2020. CA Cancer J. Clin . 70, 7–30. doi: 10. 3322/caac. 21590

Strick, A., von Hagen, F., Gundert, L., Klümper, N., Tolkach, Y., Schmidt, D., et al. (2020). The N(6) -methyladenosine (m(6) A) erasers alkylation repair homologue 5 (ALKBH5) and fat mass and obesity-associated protein (FTO) are prognostic biomarkers in patients with clear cell renal carcinoma. BJU Int . 125, 617–624. doi: 10. 1111/bju. 15019

Tian, J., Zhu, Y., Rao, M., Cai, Y., Lu, Z., Zou, D., et al. (2020). N(6)-methyladenosine mRNA methylation of PIK3CB regulates AKT signalling to promote PTEN-deficient pancreatic cancer progression. Gut. 69, 2180–2192. doi: 10. 1136/gutjnl-2019-320179

Varela, I., Tarpey, P., Raine, K., Huang, D., Ong, C. K., Stephens, P., et al. (2011). Exome sequencing identifies frequent mutation of the SWI/SNF complex gene PBRM1 in renal carcinoma. Nature 469, 539–542. doi: 10. 1038/nature09639

von Roemeling, C. A., Radisky, D. C., Marlow, L. A., Cooper, S. J., Grebe, S. K., Anastasiadis, P. Z., et al. (2014). Neuronal pentraxin 2 supports clear cell renal cell carcinoma by activating the AMPA-selective glutamate receptor-4. Cancer Res . 74, 4796–4810. doi: 10. 1158/0008-5472. CAN-14-0210

Wang, Q., Chen, C., Ding, Q., Zhao, Y., Wang, Z., Chen, J., et al. (2019). METTL3-mediated m(6)A modification of HDGF mRNA promotes gastric cancer progression and has prognostic significance. Gut 69, 1193–1205. doi: 10. 1136/gutjnl-2019-319639

Wang, T., Kong, S., Tao, M., and Ju, S. (2020). The potential role of RNA N6-methyladenosine in Cancer progression. Mol. Cancer 19: 88. doi: 10. 1186/s12943-020-01204-7

Wozniak, M. B., Le Calvez-Kelm, F., Abedi-Ardekani, B., Byrnes, G., Durand, G., Carreira, C., et al. (2013). Integrative genome-wide gene expression profiling of clear cell renal cell carcinoma in Czech Republic and in the United States. PLoS ONE 8: e57886. doi: 10. 1371/journal. pone. 0057886

Wuttig, D., Zastrow, S., Füssel, S., Toma, M. I., Meinhardt, M., Kalman, K., et al. (2012). CD31, EDNRB and TSPAN7 are promising prognostic markers in clear-cell renal cell carcinoma revealed by genome-wide expression analyses of primary tumors and metastases. Int. J. Cancer 131, E693–704. doi: 10. 1002/ijc. 27419

Zehir, A., Benayed, R., Shah, R. H., Syed, A., Middha, S., Kim, H. R., et al. (2017). Mutational landscape of metastatic cancer revealed from prospective clinical sequencing of 10, 000 patients. Nat. Med . 23, 703–713. doi: 10. 1038/nm. 4333

Zhou, Y., Zhou, B., Pache, L., Chang, M., Khodabakhshi, A. H., Tanaseichuk, O., et al. (2019). Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat. Commun . 10: 1523. doi: 10. 1038/s41467-019-09234-6