Establishment of a Risk Signature Based on m6A RNA Methylation Regulators That Predicts Poor Prognosis in Renal Cell Carcinoma

Introduction

Renal carcinoma, also known as renal cell carcinoma (RCC), is one of the most common malignant tumors of the urinary system, and it accounts for 2–3% of all malignant tumors.¹ Each year worldwide, ~295,000 individuals are diagnosed with RCC, resulting in 134,000 deaths.² The pathological types of RCC comprise clear cell, papillary, and chromophobe RCC, with >75% of cases being clear cell RCC.³ Although >50% of RCC patients are diagnosed with localized RCC, 30% have metastasis at initial diagnosis.^4,5 RCC is highly angiogenic, with a hypoxic tumor environment, and it is not sensitive to chemotherapy, so surgical resection is the only effective treatment.⁶ However, 20% of patients who have surgery experience local recurrence or distant metastases.⁷ Due to the lack of obvious symptoms in early RCC, early accurate diagnosis is a challenge, with one-third of patients displaying distant metastases at diagnosis, along with local recurrence being common after early surgical resection.^8,9 Therefore, identifying a risk signature to evaluate the prognosis of RCC patients at the diagnosis stage is of significance.

Epigenetic modification is highly associated with the occurrence and progression of tumors. It involves the regulation of gene function and expression levels via DNA methylation, thus affecting tumor progression. N6-methyladenosine (m6A) modification of RNA is the most common post-transcriptional modification in eukaryotic cells, accounting for 80% of RNA methylation modifications.¹⁰ m6A methylation is involved in the intricate regulation of important genes in many cellular activities. Additionally, several studies have demonstrated that m6A methylation has an important role during the occurrence and progression of tumors.^11–13

It has been reported that m6A modification is a dynamic reversible process.¹⁴ m6A-related proteins can be divided into three categories based on function: methyltransferases (writers), demethylases (erasers), and m6A-binding proteins/effectors (readers). Writers mediate RNA methylation, and they include METTL3, METTL14, WTAP, ZC3H13, RBM15, and KIAA1429 (also known as VIRMA). The first two are the most common writers, and they can catalyze m6A methylation of mRNA in vitro and in vivo¹⁵ Some writers, such as WTAP, do not display methyltransferase activity. However, WTAP interacts with the METTL3/METTL14 complex to influence m6A methyltransferase activity,¹⁶ ZC3H13 promotes m6A methylation by anchoring WTAP,¹⁷ and RBM15 and KIAA1429 facilitate the recruitment of the METTL3/METTL14 complex.¹⁸ Moreover, erasers are able to “erase” RNA methylation modification signals, mediating RNA demethylation. FTO, an eraser, uses ferrous iron as a cofactor along with α-ketoglutarate to oxidize the N-methyl at the m6A site to hydroxymethyl, thereby altering the cellular m6A level.¹⁴ Lastly, readers are responsible for “reading” the methylation information. They can be involved in downstream RNA translation and degradation. There are two “reading” modes: direct reading, where the readers (such as the YTH protein family, YTHDC1, YTHDC2, YTHDF1, and YTHDF2) selectively interact with the m6A site of RNA,¹⁹ and indirect reading. This involves the m6A methylation changing the secondary structure of the RNA and thereby rendering it accessible to certain “readers” such as heterogeneous nuclear ribonucleoprotein C (HNRNPC), which can bind to other proteins and promote mRNA maturation.²⁰

The interaction between m6A mRNA methylation and tumors has received increasing attention. Increasing evidence indicates that m6A mRNA methylation is closely associated with tumorigenesis and development, and m6A-related proteins are important regulators of these processes.¹⁰ In this study, RNA sequencing data from The Cancer Genome Atlas (TCGA) database were used to analyze the expression of 13 major m6A RNA methylation regulators in 893 RCC patients. In addition, a risk signature, based on two of the m6A methylation regulators, was established to predict the prognosis of the RCC patients. This risk signature was used to divide the RCC cases into the high- and low-risk groups (based on the median risk score). Thereafter, the associations between the high/low-risk groups and both overall survival and various clinicopathological characteristics were assessed.

Methods

Datasets

Transcriptome (RNA sequencing) data and relevant clinical information on RCC patients (including those with clear cell, papillary, and chromophobe RCC) were acquired from the TCGA database (cancergenome.nih.gov). The TCGA dataset included information on 893 RCC patients and 128 normal controls. Additionally, a Gene Expression Omnibus (GEO) dataset (GSE29609; ncbi.nlm.nih.gov/geo), involving 39 RCC patients, was used to validate the prognostic value of the risk model constructed based on TCGA data.

Selection of m6A RNA Methylation Regulators

Thirteen major m6A regulators were analyzed: ZC3H13, WTAP, METTL3, METTL14, RBM15, KIAA1429, HNRNPC, YTHDC1, YTHDC2, YTHDF1, YTHDF2, ALKBH5, and FTO. An expression matrix of the 13 genes was established, and data on the clinicopathological characteristics were also obtained and utilized in the bioinformatics analyses.

Cell Culture

RCC cell lines (786–0 and 769-P) were cultured in Roswell Park Memorial Institute (RPMI)‐1640 medium (Gibco, CA, USA) with 10% fetal bovine serum (FBS; Life Technologies, Carlsbad, CA, USA), 100 U/mL penicillin, and 100 mg/mL streptomycin. Another RCC cell line (Caki-1) was cultured in McCoy’s 5A (Gibco, USA) with 10% FBS (Life Technologies), 100 U/mL penicillin, and 100 mg/mL streptomycin. The 293T cell line (the control kidney cell line) was cultured in Dulbecco’s Modified Eagle Medium (DMEM; Gibco) with 10% FBS (Life Technologies), 100 U/mL penicillin, and 100 mg/mL streptomycin. The cells were cultured at 37°C in 5% CO₂. All cells were acquired from the American Type Culture Collection (ATCC; Rockville, MD, USA).

Real-Time PCR

The RNA was extracted from the cell lines using TRIzol Reagent (Invitrogen, CA, USA) and reverse transcribed (PrimeScript™ RT Master Mix; Takara, Kyoto, Japan) to acquire cDNAs. Real-time PCR was performed using SYBR^® Premix Ex Taq™ GC (Takara). β-actin was used as the internal control. The primers for METTL14 and WTAP (m6A RNA methylation regulators) are listed in Table S1. The 2^−ΔΔCt method was used to calculate the relative expression of these genes.

Bioinformatics Analyses

The limma package (version 3.22.7) in R was used to analyze the expression of the 13 m6A regulators in the TCGA dataset.²¹ Subsequently, the vioplot package (version 3.4) was used to visualize the expression of these genes in the RCC and normal tissues. A consensus clustering algorithm in the ConsensusClusterPlus package (version 1.52.0) inR²² was used to evaluate the stability of clustering and determine the optimum number of clusters. Least absolute shrinkage and selection operator (LASSO) Cox regression was used to establish a risk signature for predicting overall survival based on the m6A regulators. The coefficients of the selected genes were determined, with penalty parameter (λ) tuning based minimum criteria. The glmnet package (version 3.0) and the survival package were used to execute this program.²³ The following formula was used to calculate the risk score for each patient: Risk score = β1×ExprGene1 + β2×ExprGene2 + … + βn×ExprGenen, where β represents the coefficients (1 to n). The median risk score was used as the cut-off value to determine the high- and low-risk groups. The predictive performance of the risk model was assessed using a receiver operating characteristic (ROC) curve analysis.²⁴ Furthermore, Kaplan–Meier analysis was used to compare overall survival between the high- and low-risk groups, and univariate and multivariate Cox regression analyses of the association between the risk score and overall survival were conducted. The clinicopathological characteristics of the high/low-risk groups (based on the median risk score) were also investigated.

Statistical Analysis

One-way analysis of variance (ANOVA) was used to compare the expression levels of the 13 m6A regulators in the TCGA dataset. Student’s t-test was used to explore the associations between the high/low-risk groups of RCC patients (based on the median risk score) and clinicopathological characteristics, such as age, gender, tumor stage, and tumor, node, and metastasis (TNM) status. Kaplan–Meier analysis and the log rank test were used to analyze the overall survival of RCC patients in the two clusters in the TCGA dataset and in the high- and low-risk groups (based on the median risk score) in the TCGA and GEO datasets. Furthermore, univariate and multivariate Cox regression analyses were used to assess the associations of the m6A RNA methylation regulator expression levels, the risk score, and other clinicopathological characteristics with overall survival.

Cases with missing data in the TCGA dataset were excluded from the Kaplan–Meier and Cox regression analyses of overall survival, leaving 634 cases for these analyses. P<0.05 was considered to indicate statistical significance. The statistical analyses were performed using R software (version 3.6; www.r-project.org).

Results

Expression of m6A RNA Methylation Regulators in RCC

The m6A RNA methylation regulators play key roles in the formation and development of various tumors. First, the expression levels of these regulators (Table 1) in 893 RCC and 128 normal tissues obtained from the TCGA database were compared. Of the 13 genes, ten genes displayed significant differences. Six genes (METTL3, RBM15, YTHDF1, YTHDC2, ALKBH5, and FTO) were significantly upregulated and four genes (METTL14, ZC3H13, HNRNPC, and YTHDF2) were significantly downregulated in RCC patients compared to normal controls (Figure 1A and B). Furthermore, the expression levels of the 13 regulators were positively correlated (Figure 1C).

Table 1 Expression of the 13 N6-Methyladenosine (m6A) RNA Methylation Regulators in Renal Cell Carcinoma (RCC)

Figure 1 Expression levels of m6A RNA methylation regulators in RCC. (A) Heatmap of the expression of 13 m6A RNA methylation regulators in RCC tissues compared to normal tissues. (B) Violin plots of the expression of the 13 regulators in RCC. (C) Spearman correlation analysis of the 13 regulators in RCC. *P<0.05, **P<0.01, and ***P<0.001. Abbreviations: RCC, renal cell carcinoma; m6A, N6-methyladenosine.

Identification of Two Clusters of RCC Cases by Consensus Clustering Analysis

The ConsensusClusterPlus package was used to perform a consensus clustering analysis based on the expression levels of the 13 m6A RNA methylation regulators in the 893 RCC tissues in the TCGA database. The results indicated that the relative change in the area under the cumulative distribution function (CDF) curve obviously decreased when k=3 (Figure 2A and B). However, the consensus matrix showed that there was considerable classification overlap among the clusters when the cases were divided into three clusters. Therefore, the cases were divided into two clusters (Figure 2C). For k=2–9, the distributions of the RCC cases are presented in Figure 2D. Principal component analysis of the two clusters was used to assess the performance of the classification. The results indicated that the cluster 1 cases clustered together, as did the cluster 2 cases (Figure 2E).

Figure 2 Consensus clustering analysis of RCC cases. (A) CDF curves for k=2–9. (B) Relative change in the area under the CDF curves for k=2–9. (C) Consensus matrix for k=2. (D). Tracking plot of cases for k=2–9. (E). Principal component analysis plot based on the RNA expression profiles of the 13 major m6A regulators, which grouped the cases into two clusters. Red: cluster 1; blue: cluster 2. Abbreviations: CDF, cumulative distribution function; RCC, renal cell carcinoma.

Cluster 1 and 2 are Associated with a Clinicopathological Characteristics (T Status) in RCC

To further understand the clustering results, a Kaplan–Meier analysis of the overall survival of 634 RCC patients in the two clusters was conducted. The results indicated that cluster 2 had a worse prognosis compared to cluster 1, but the difference was not significant (P=0.084; Figure 3A). Moreover, the clusters were significantly associated with T status (Figure 3B). Collectively, the results indicated that the clusters were closely associated with RCC.

Figure 3 Overall survival and clinicopathological characteristics of RCC patients in the two clusters. (A) Kaplan–Meier curves of overall survival of RCC patients in clusters 1 (red) and 2 (blue). (B) Heatmap showing the associations between clinicopathological characteristics and the expression of each of the 13 m6A RNA methylation regulators in clusters 1 and 2. Red: upregulated; green: downregulated. **P<0.01. Abbreviations: m6A, N6-methyladenosine; RCC, renal cell carcinoma.

Prognostic Risk Model Based on the m6A RNA Methylation Regulators

Subsequently, the associations between the expression levels of the m6A RNA methylation regulators and prognosis in RCC patients were assessed. Univariate Cox regression analyses were conducted according to the expression levels of the regulators in the TCGA database. The results indicated that WTAP (hazard ratio [HR], 1.051; 95% confidence interval [CI], 1.021–1.082), HNRNPC (HR, 1.016; 95% CI, 1.003–1.028), METTL3 (HR, 1.086; 95% CI, 1.017–1.161), YTHDC1 (HR, 0.955; 95% CI, 0.90–0.98), and METTL14 (HR, 0.730; 95% CI, 0.638–0.836) were significantly associated with overall survival (Figure 4A).

Figure 4 Identification of risk signature based on m6A RNA methylation regulators. (A) Univariate Cox regression results (hazard ratio, 95% confidence interval) for the 13 m6A RNA methylation regulators. (B and C) Coefficients from the multivariate LASSO Cox regression. LASSO: least absolute shrinkage and selection operator. Abbreviation: m6A, N6-methyladenosine.

Thereafter, two of the genes with P<0.001 in the univariate analyses of TCGA data (METTL14 and WTAP) were selected, and the LASSO Cox regression algorithm was subsequently applied. These two genes were selected for the risk signature via the minimum criteria. The risk scores were then determined using the coefficients acquired from the LASSO algorithm (Figure 4B and C). The patients were then divided into high- and low-risk groups based on the median risk score.

Risk Scores are Associated with Clinicopathological Characteristics in RCC

To further analyze the associations between the risk score and the clinical traits of RCC patients, the associations between the high/low-risk groups and both overall survival and various clinicopathological characteristics were assessed. First, the results indicated that the high-risk group had a worse overall survival compared to the low-risk group (Figure 5A). Second, the associations between the high/low-risk groups and clinicopathological characteristics, including TNM status, stage and gender, were determined (Table 2), and strong associations between the high/low-risk groups and M status, T status, stage, and gender were identified (Figure 5C). The results also indicated that WTAP was downregulated while METTL14 was upregulated in the high-risk group compared to the low-risk group.

Table 2 Clinical Characteristics of the Patients in High- and Low-Risk Groups (Based on the Median Risk Score)

Figure 5 Overall survival and clinicopathological characteristics of RCC patients in the high- and low-risk groups. (A). Kaplan–Meier curves of overall survival of patients in the high- and low-risk groups (based on the median risk score involving two m6A RNA methylation regulators). (B) Receiver operating characteristic curve showing the predictive performance of the risk model. (C) Heatmap showing the associations between clinicopathological characteristics and the expression of two m6A RNA methylation regulators (METTL14 and WTAP) in the high- and low-risk groups. Red: upregulated; green: downregulated. *P<0.05 and ***P<0.001. Abbreviations: RCC, renal cell carcinoma; m6A, N6-methyladenosine.

A ROC curve analysis was used to further verify the predictive performance of the risk signature. The area under the ROC curve was 0.701 (Figure 5B), which indicated that the model had high predictive performance.

Validation in the Gene Expression Omnibus (GEO) Dataset

The GSE29609 GEO dataset (involving 39 RCC patients) was used to further validate the prognostic value of the risk signature in RCC patients. Based on the median risk score, the 39 patients were divided into high- and low-risk groups. The results indicated that the patients in the high-risk group had a significantly worse prognosis than the patients in the low-risk group (P<0.05; Figure 6). The result was consistent with the TCGA results, indicating that the risk model was accurate and effective.

Figure 6 Kaplan–Meier curves of overall survival of RCC patients in the high- and low-risk groups (based on the median risk score involving two m6A RNA methylation regulators) in the Gene Expression Omnibus dataset. Abbreviations: m6A, N6-methyladenosine; RCC, renal cell carcinoma.

Cox Regression Analysis Indicates That the Risk Signature is an Independent Prognostic Indicator

Univariate (Figure 7A) and multivariate (Figure 7B) Cox regression analyses indicating the associations between the clinicopathological characteristics (including risk score) and overall survival were performed. The results indicated that age, stage, and risk score were significantly associated with overall survival in RCC (Figure 7A and B). Collectively, these results indicated that the risk score obtained from the risk signature may serve as an independent prognostic factor for overall survival in RCC patients.

Figure 7 Associations between the clinicopathological characteristics and overall survival of RCC patients. (A) Univariate and (B) multivariate Cox regression results (hazard ratio, 95% confidence interval) indicating the associations between the clinicopathological characteristics and the overall survival of RCC patients. Abbreviation: RCC, renal cell carcinoma.

Validation of mRNA and Protein Expression Levels of the Two m6A RNA Methylation Regulators

To further verify our results, information from the Human Protein Atlas (HPA) database was used to analyze the protein expression of the two m6A RNA methylation regulators: METTL14 and WTAP. The protein expression of METTL14 was higher in RCC tissues than normal tissues (Figure 8A), while the protein expression of WTAP was higher in normal tissues than RCC tissues (Figure 8B). These results were consistent with the expression of METTL14 and WTAP at the mRNA level (Figure 8C and D); our results showed that compared to the control kidney cell line (293T), METTL14 mRNA was significantly decreased in the RCC cell lines (786–0 and 769-P), while WTAP mRNA was significantly increased in RCC cell lines (786–0 and Caki-1). These results demonstrate that METTL14 and WTAP showed significant changes in RCC at both the mRNA and protein levels, indicating that they have important roles in RCC.

Figure 8 METTL14 and WTAP levels. (A) METTL14 and (B) WTAP protein expression in normal and RCC tissues from the Human Protein Atlas, as detected by immunohistochemical analysis. (C) METTL14 and (D) WTAP mRNA expression in normal (293T) and RCC cell lines. **P<0.01, ***P<0.001, and ****P<0.0001. Abbreviation: RCC, renal cell carcinoma.

Discussion

The majority of RCC cases are clear cell RCC.³ Patients with early RCC have no obvious symptoms, so some of the patients already have metastases at pathological diagnosis. Molecularly targeted therapies are prone to drug resistance. RCC patients typically have a poor prognosis and the 5-year survival rate is <40%.²⁵ Interferon and interleukin-2-based immunotherapy, as the first-line treatment, displays 15% efficacy.²⁶ The occurrence and development of RCC is a complicated process involving numerous factors, such as smoking, obesity, genetics, and cell cycle changes.²⁷ Thus, the precise mechanisms underlying RCC are complex. Developing novel therapeutic strategies based on these mechanisms is critical.

m6A modification of RNA refers to methylation of the nitrogen atom at position 6 of adenine (A), catalyzed by methyltransferase.²⁸ Each mRNA has 3–5 m6A methylation sites that regulate RNA stability, localization, transport, splicing, and translation at the post-transcriptional level.²⁹ Increasing evidence has demonstrated the role of m6A modification in acute myeloid leukemia (AML), hepatocellular carcinoma (HCC), breast cancer, and other types of cancer. In malignant tumors, m6A modification can induce a carcinogenic or anticancer effect. However, little is known about the roles of m6A methylation regulators in RCC. The risk model involving m6A RNA methylation regulators in RCC in this study provides new insights for further research into RCC.

According to the 893 RCC and 128 normal tissues analyzed in this study, the expression of the 13 m6A RNA methylation genes was significantly altered in RCC patients compared to normal controls, which indicated that m6A RNA modification may serve a role in RCC. Among the 13 regulators, the change in FTO expression was the most significant. FTO is a member of the non-heme Fe (II) and α-ketoglutarate-dependent dioxygenase ALKB protein family.³⁰ Recent studies have suggested that FTO expression may be associated with several tumors, such as breast, thyroid, endometrial, and gastric cancer.^31,32 Li et al reported that FTO affects the occurrence and resistance of leukemia by reducing m6A levels in the ASB2 and RARA target genes.³³ Zhou et al demonstrated that FTO expression was significantly increased in patients with cervical squamous cell carcinoma (CSCC), and FTO and β-catenin expression have prognostic value in CSCC.³⁴ These studies above have investigated the biomolecular mechanism underlying FTO-induced m6A modification during cancer development, but further investigation is required.

To further study the effect of the 13 m6A regulators on the clinical characteristics and prognosis of RCC, clinical data were downloaded from the TCGA database for analysis. Consensus clustering was performed to divide the RCC cases into two clusters. The results indicated that the clusters were associated with a clinicopathological characteristic (T status) in RCC. The results also indicated that the expression of the m6A RNA methylation regulators differed between the two clusters. Therefore, this analysis suggested that the m6A RNA methylation regulators were strongly associated with RCC.

The associations between the m6A RNA methylation regulators and the prognosis of RCC were further investigated by performing univariate Cox regression analyses. Based on the results, two m6A RNA methylation regulators with P<0.01 (METTL14 and WTAP) were selected for LASSO Cox regression analysis to establish a risk signature. Additionally, based on the heatmap, METTL14 was upregulated and WTAP was downregulated in the high-risk group compared to the low-risk group. Although METTL14 and WTAP are both writers, their expression levels and roles in tumors are not the same. WTAP is upregulated in AML, which promotes tumorigenesis and development. In contrast, METTL14 is downregulated in HCC, where it has a suppressive role.^35,36 Each m6A RNA methylation regulator works differently in different tumors. For example, METTL3 is downregulated in glioblastoma and has an inhibitory role, but it is upregulated in lung cancer and HCC.^37–39 The roles of METTL14 and WTAP in RCC are not completely understood and require further investigation.

Our results demonstrated that the high/low-risk groups (based on the median risk score) were significantly associated with RCC prognosis. Moreover, the area under the ROC curve for the risk signature was 0.701, which indicated that the predictive performance of the model was good. In addition, the high/low-risk groups were associated with clinicopathological characteristics (P<0.05; T status, M status, stage, and gender). Finally, univariate and multivariate Cox regression analyses were performed to assess the associations between the risk score (and other clinicopathological characteristics) and overall survival. The results indicated that the risk signature may serve as an independent prognostic factor in RCC. Therefore, WTAP and METTL14 may serve as independent prognostic factors in RCC.

WTAP is an essential protein in m6A methylation modification⁴⁰ as it coordinates METTL3-METTL14 heterodimer localization and promotes m6A methylation.¹⁶ This indicates that the carcinogenic effect of WTAP may be closely related to the m6A methyltransferase complex. In addition, WTAP is linked to the development of other cancers. WTAP expression was found to be significantly increased in pancreatic ductal adenocarcinoma (PDAC) and was an independent prognostic factor in PDAC.⁴¹ WTAP may affect RCC cell proliferation by regulating the stability of CDK2 mRNA, leading to the occurrence and development of cancer.⁴¹ Therefore, WTAP may serve as a novel target for the diagnosis and treatment of RCC.

METTL14 is an important component of the m6A methyltransferase complex.¹⁵ Recent studies have reported its role in the malignant phenotype of tumors. One study demonstrated that METTL14 knockout in AML cell lines effectively inhibited proliferation.⁴² Moreover, METTL14 and FTO were downregulated in HCC,³⁶ which indicated that m6A modification involved a complex feedback regulation mechanism involving both the writer and the eraser. To date, few studies have investigated the role of METTL14 in RCC. This study indicated that METTL14 was upregulated in the high-risk group compared to the low-risk group, which suggested that METTL14 upregulation is associated with poor RCC prognosis.

The specific regulation of the m6A methylation regulators and their potential biological functions require further investigation, especially regarding the two m6A genes (which encode writers) included in the risk signature in this study.

Conclusions

In conclusion, the results indicated that the risk score may serve as an independent prognostic factor and may also predict clinicopathological characteristics of RCC. However, the roles of m6A RNA methylation regulators in RCC are not completely understood, so further investigation of the precise underlying mechanisms is required. The results of this study may aid in investigating prognosis, diagnosis, and therapeutic strategies for RCC.