Decreased expression of METTL14 predicts poor prognosis and construction of a prognostic signature for clear cell renal cell carcinoma

Background METTL14, as one of N6-methyladenosine (m6A) related genes, has been found to be associated with promoting tumorigenesis in different types of cancers. This study was aimed to investigate the prognostic value of METTL14 in clear cell renal cell carcinoma (ccRCC). Methods We collected ccRCC patients’ clinicopathological parameters information and 13 m6A related genes expression from The Cancer Genome Atlas (TCGA) database. Univariate and multivariate Cox regression analyses were conducted to investigate whether METTL14 could serve as an independent factor correlated with overall survival (OS). Gene Set Enrichment Analysis (GSEA) was carried out to identify METTL14-related signaling pathways. Moreover, a risk score (RS) was calculated to predict the prognosis of ccRCC. Quantitative real-time PCR (qRT-PCR) was also utilized to verify the expression of METTL14 in clinical specimens. Results Differently expressed m6A related genes were identified between ccRCC tissues and normal tissues. Therein, METTL14 was lowly expressed in ccRCC tissues and verified by qRT-PCR (all p < 0.01). Survival analysis indicated that high expression of METTL14 was associated with better OS (p = 1e−05). GSEA results revealed that high METTL14 expression was enriched in ERBB pathway, MAPK pathway, mTOR pathway, TGF-β pathway and Wnt pathway. Moreover, METTL14 was proved to be an independent prognostic factor by means of univariate and multivariate Cox regression analyses. Nomogram integrating both the METTL14 expression and clinicopathologic variables was also established to provide clinicians with a quantitative approach for predicting survival probabilities of ccRCC. Furthermore, a METTL14-based riskscore (RS) was developed with significant OS (p = 6.661e−16) and increased AUC of 0.856. Besides, significant correlated genes with METTL14 were also provided. Conclusions Our results indicated that METTL14 could serve as a favorable prognostic factor for ccRCC. Moreover, this study also provided a prognostic signature to predict prognosis of ccRCC and identified METTL14-related signaling pathways.

rate of RCC is 74% and 16% of patients have been diagnosed with local invasion or distant metastasis, resulting in a poorer prognosis [3]. Hence, the early diagnosis and treatment of ccRCC are imperative and challenging. There is an urgent need to clarify the possible mechanisms underlying the development of ccRCC, to discover novel biomarkers for early screening, and to develop therapeutic targets.
Methylation of the N6 position of adenosine (m6A) has been found to be widespread throughout nuclear mRNA, rRNA, tRNA, and some snRNA of eukaryotes [4,5]. The modified RNAs, especially mRNA, play critical roles in the post-transcriptional regulation of gene expression. As the most common type of mRNA modification in eukaryote, m6A was proved to be highly conserved, as well as in human and mouse. The m6A sites were found mostly located in 3' untranslated regions (UTR) or near the stop codon [6,7]. RNA m6A modification has been shown to play vital for RNA splicing, export, stability, regulation of RNA-protein interaction and gene expression [4,[8][9][10]. m6A is a reversible modification of mRNA mediated by the dynamic activities of m6A "writers" and "erasers" [11]. The genes of "writers" compose a mRNA methyltransferase enzyme complex, including METTLE3, METTL14, WTAP, KIAA1429, ZC3H13 and RBM15, while the reversible process is catalysed by m6A "erasers", FTO and ALKBH5, conducting demethylation of m6A. In order to exert functions of m6A modifications, mRNA m6A sites can be recognized by RNA-binding proteins "readers", including YTHDF1, YTHDF2, YTHDC1, YTHDC2 and HNRNPC.
Recently, m6A modification has been found to regulate embryonic stem cells and myeloid differentiation, spermatogenesis, obesity, sex determination, neuronal development [12][13][14][15][16][17]. The m6A regulatory genes have been shown to plays an oncogenic role in both breast cancer and acute myeloid leukemia by means of cancer stem cell maintenance and dedifferentiation of cancer cells [4,18]. Previous study has suggested that METTL3 might be essential in lung cancer cells by directly promoting translation independently of methyltransferase activity and downstream m6A reader proteins [19]. According to another study, it has been proved that METTL14, a major RNA N6-adenosine methyltransferase, suppresses the metastatic potential of hepatocellular carcinoma by modulating m6A-dependent primary microRNA processing [20]. Therein, METTL14, as one of N6-methyladenosine (m6A) related genes, has been found to be associated with promoting tumorigenesis in different cancers. Recent studies reported that bioinformatics analysis could provide possible functions of many different novel genes in cancer research [21]. Hence, this study was aimed to investigate the prognostic value of METTL14 in ccRCC by bioinformatics analysis. Furthermore, we also provided a prognostic signature to predict prognosis of ccRCC patients and related signal pathways regulated by METTL14 by means of gene set enrichment analysis (GSEA).

Data acquisition
All of these normalized RNA-seq data and clinical data of kidney renal clear cell carcinoma (KIRC) patients were obtained from TCGA website (http://cance rgeno me.nih. gov/) by TCGA-assembler. Finally, we identified the gene expression profiles and clinical data of 539 ccRCC cases and 72 cases of normal control for further analysis. In this current study, overall survival (OS) was the primary outcome. A total of 13 m6A-related genes contained writers (KIAA1429, METTL3, METTL14, RBM15, WTAP and ZC3H13), erasers (FTO and ALKBH5) and readers (YTHDF1, YTHDF2, HNRNPC, YTHDC1 and YTHDC2). By utilizing the R programming language, we did an overlap by comparing the standardized RNA-seq data with the m6A-related genes.

Identification of differently expressed m6A-related genes
Differently expressed m6A-related genes between ccRCC and normal tissues were calculated by utilizing the R "Limma" package and adjusted P-value (FDR) less than 0.05 and at least twofold changes (FCs) were set as the threshold. Heatmap and boxplot were used to show the distribution and expression of differently expressed m6A-related genes respectively. Besides, the R "corrplot" package was used to calculate the pearson correlations between 13 m6A-related genes based on gene expression levels at the transcriptional level.

Kaplan-meier survival analysis and ROC curve of METTL14 in ccRCC
Based on the median expression of METTL14, ccRCC patients were divided into high-risk groups and lowrisk groups. By means of Kaplan-Meier (K-M) method and the log-rank test, we evaluated the survival differences between these two groups. The receiver operating characteristic curves (ROC) analyses were performed by using the R "survivalROC" package and the area under the curve (AUC) values were calculated to assess the specificity and sensitivity of METTL14.

Quantitative real-time PCR (qRT-PCR)
The total RNA from ccRCC tissues and adjacent normal tissues was acquired using TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) on the basis of a standard extraction protocol. The cDNA was synthesized by HiScript II (Vazyme, China), and qPCR analysis was performed on StepOne Plus Real-Time PCR system (Applied Biosystems, USA). The initial reaction was incubated at 95 °C for 10 min, followed by 40 cycles of 95 °C for 15 s and 60 °C for 1 min in accordance with the SYBR green method. The relative expression levels were analyzed by comparing 2 −ΔΔCT method. All reactions were carried out in triplicate. Actin was used for the internal reference. Primers were synthesized by TSINGKE (Beijing, China), including Actin (F:5' ATG ACT TAG TTG CGT TAC ACC 3' , R:5' GAC TTC CTG TAA CAA CGC ATC 3') and METTL14 (F:5' TTT CTC TGG TGT GGT TCT GG 3' , R:5' AAG TCT TAG TCT TCC CAG GATTG 3'). Six paired tumor tissues and their adjacent normal tissues were obtained from ccRCC patients from Affiliated Hospital of Nantong University. Ethical approval was obtained from the Institutional Research Ethics Committees of Affiliated Hospital of Nantong University and informed written consent was obtained from all of subjects.

Univariate and multivariate cox hazard regression analyses
In order to identify independent prognostic factors, univariate and multivariate Cox regression analyses were utilized to exclude clinical characteristics with little prognostic values with OS from age, gender, race, grade, stage, staged T, staged N, staged M and METTL14 expression level by R package. Moreover, ROC curves of these nine factors and their AUC values were also calculated and compared by the R "survivalROC" package.

The establishment of a nomogram for prognosis prediction
In order to predict the 1-, 3-, and 5-year survival probabilities, a nomogram was conducted to visualize the associations between survival rates of OS and individual predictors (age, gender, race, grade, stage, T, M, N and METTL14) by means of the R "rms" package. By using the point scale in the nomogram, every parameter was assigned a point and the whole points were calculated by summing up the points of all factors.

Gene set enrichment analysis (GSEA)
GSEA is a systematic method used to figure out whether the hallmark gene sets predicted statistically significant differences between two different groups [22]. Here, GSEA was performed to analyze the significant survival differences between two groups divided by the expression of METTL14. The permutation tests were carried out 1000 times for discovering significant critical biological pathway. It had been considered to be significantly enriched that nominal p value less than 0.05 and FDR less than 25%.

Prognostic model (riskscore, RS) construction and evaluation
By means of multivariate COX regression analysis, we excluded some clinical characteristics (gender, race, T, M, N) with little prognostic value. Moreover, based on the correlation coefficient of each significant clinical characteristic (age, grade, and stage) and METTL14, we established an independent prognostic index (RS) for predicting the OS of ccRCC patients. It was calculated using the following formula β 1 × clinical characteristics 1 + β 2 × clinical characteristics 2 + ····· + β n × METTL14, where β was the correlation coefficient.
As for evaluation, ccRCC patients will be classified into low-risk and high-risk groups, according to the median risk score of our established prognostic model (RS). By means of Kaplan-Meier (K-M) method and the logrank test, we evaluated the survival differences between these two groups. The ROC analyses were performed by using the R "survivalROC" package based on each prognostic model. We also displayed the risk score distribution of different risk groups and the number of censored patients.

Statistical analysis
All statistical data and figures were analyzed by using SPSS 24.0 (IBM, Chicago, USA), R 3.3.1 (https ://www.rproje ct.org/) and GraphPad Prism 6.0 (San Diego, CA, USA). The correlations between two different genes were analyzed by Pearson's correlation method. The association between clinicopathologic parameters and METTL14 was estimated by the Wilcoxon signed rank test and logistic regression. The Kaplan-Meier curve and log-rank test were carried out to estimate survival predictive performance of METTL14 and RS. Univariate and multivariate Cox regression analyses were used to evaluate the relationship between variables and OS. A nomogram was also created using the rms package of R software. ROC and its AUC values were calculated to evaluate the prognostic ability of METTL14 or RS by using the package of "survivalROC" in R. All statistical results with p < 0.05 were considered statistically significant.

Relative expression of m6A-related genes in ccRCC
We analyzed the mRNA expression levels of 13 m6Arelated genes (METTL3, METTL14, WTAP, RBM15, ZC3H13, KIAA1429, FTO, ALKBH5, YTHDC1, YTHDC2, YTHDF1, YTHDF2, HNRNPC) in TCGA database to identify the different expressed genes between ccRCC tissues and normal tissues. Figure 1a is the heat map which indicated the expression levels of m6A-related genes. The expression of each m6A-related gene in ccRCC tissues and adjacent normal tissues is displayed in Fig. 1b. We further investigated the correlation between each two types of m6A-related genes. As displayed in Fig. 1c, positive correlations were shown as red bubbles and negative correlations were shown as blue bubbles. The size of each bubble represented the P value. The result revealed that METTL14 and YTHDC1 were the two most relevant m6A-related genes (Pearson's r = 0.66). METTL14 also had a positive correlation with KIAA1429 (Pearson's r = 0.53) and ZC3H13 (Pearson's r = 0.53), respectively.

The expression of METTL14 in ccRCC and qRT-PCR verification
As presented in 72 pairs of ccRCC and matched normal tissues, METTL14 had a low expression in tumor tissues in TCGA dataset (p = 0.004, Fig. 2a). Moreover, the mRNA expression level of METTL14 was significantly decreased in ccRCC tissues (n = 539) compared with that in normal tissues (n = 72) in Fig. 2b (p < 0.001). Receiver operating characteristic (ROC) curves for 5-year survival were carried out to demonstrate the performance of the METTL14 expression level for survival prediction (Fig. 2c), and its 5-year AUC was 0.669. The K-M curve indicated that patients in the low-METTL14 group had a shorter OS time than patients in the high-METTL14 group (p < 0.001, Fig. 2d). As showed in Fig. 2e and Additional file 1. Table S1, the qRT-PCR results verified that the expression of METTL14 was significantly down-regulated in the ccRCC tissues (n = 6) compared with adjacent normal tissues (n = 6; p < 0.001).

Association between METTL14 expression and clinicopathologic parameters
The results of independent sample t-tests showed that the METTL14 expression was lower in high grade (p < 0.001; Fig. 3a), pathologic stage (p < 0.001; Fig. 3b), T stage (p < 0.001; Fig. 3c), and M stage (p < 0.001; Fig. 3d). In addition, there is no difference of METTL14 expression identified between male and female, Asian and Caucasian, or N0 stage and N I-III stage.
Based on logistic regression analysis, METTL14 expression was revealed to be associated with prognostic clinicopathologic characteristics (Table 1). High METTL14 expression in ccRCC patients is significantly These analyses revealed that ccRCC patients with low METTL14 expression tend to progress to a more advanced grade and distant metastasis than those with high METTL14 expression.

METTL14 could serve as an independent prognostic factor
Univariate Cox and multivariate Cox regression analyses were conducted to investigate whether the METTL14 expression was an independent factor correlated with Fig. 2 The expression of METTL14 in ccRCC tissues; a pairwise boxplot of the METTL14 expression between the ccRCC (n = 72) and matched normal tissues (n = 72) in TCGA dataset; b relative expression levels of the METTL14 expression between the ccRCC (n = 539) and normal tissues (n = 72) in TCGA dataset; c ROC curve performed to assess the OS predictive performance of METTL14 expression in TCGA dataset; d Kaplan-Meier curves for OS in the high-risk and low-risk groups when stratified by METTL14 expression; e qRT-PCR verification of METTL14 expression between ccRCC tissues (n = 6) and adjacent normal tissues (n = 6). The bar graphs represent means ± standard deviation. ***p < 0.001 OS ( Table 2). Univariate Cox analysis revealed that the METTL14 expression, age, grade, pathological stage, T stage and M stage were all associated significantly with the overall survival of ccRCC patients (all p < 0.05; Fig. 4a). Multivariate Cox regression analysis showed high METTL14 expression correlated significantly with a better OS (HR = 0.808; p = 0.008). Other clinicopathological parameters correlated with prior overall survival consist of old age, low grade and advanced stage (Fig. 4b).
Hence, our results indicated that the METTL14 expression might be an independent prognostic factor of OS when adjusted by those variables. Moreover, ROC curves for OS were also carried out to demonstrate the predictive ability of the METTL14 expression and clinicopathological variables (Fig. 4c). The AUC of the METTL14 expression (AUC = 0.667) was obviously higher than that of age (AUC = 0.

Establishment of ccRCC prognostic prediction nomogram
In order to predict the 1-, 3-, and 5-year survival probabilities, a nomogram was established based on eight clinicopathological parameters (age, gender, race, grade, stage, T, M, N) and METTL14 gene (Fig. 4d). By using the point scale in the nomogram, every parameter was

GSEA identified METTL14 -related signaling pathways
In consideration of METTL14 as an independent prognostic factor of OS for ccRCC patients, we were eager to explore how METTL14 was involved in ccRCC pathogenesis. We carried out Gene Set Enrichment Analysis (GSEA) between tissues with different METTL14 expression levels. Based on the normalized enrichment score (NES) and FDR q-val (FDR < 0.01), we selected the most significantly enriched biological pathways. The results indicated that high METTL3 expression was associated with some essential signaling pathways including ERBB pathway, MAPK pathway, mTOR pathway, renal cell carcinoma, pathway in cancer, TGF-β pathway and Wnt pathway ( Fig. 5 and Table 3), giving a clue of the underlying mechanism in the pathogenesis of ccRCC.

Construction and evaluation of the RS
Based on the multivariate Cox regression analysis with TCGA dataset, METTL14 expression and three clinicopathological variables consisting of age, pathologic stage, and histological grade were indicated as independent prognostic parameters (p < 0.05) ( Table 2). Thus, METTL14 expression and three clinicopathological variables including age, pathologic stage, and histological grade were integrated to develop the riskscore (RS) that could predict the prognosis of ccRCC by using a multivariate Cox regression model ( Table 4). The formula of RS was as follows: RS=(0.036132 × age) + (0.310219 × grade) + (0.509547 × stage) + (− 0.20258 × METTL14 expression).The ccRCC patients in TCGA dataset were divided into high-risk and low-risk groups by the median expression value of RS. In order to verify the performance of RS in predicting the OS of ccRCC, we constructed a Kaplan-Meier curve to analyze the different survival time between the highand low-risk groups. Patients in the low-risk group had much better OS compared with that in the high-risk group (P < 0.001, Fig. 6a). Figure 6b, c indicated that the higher the risk scores, the higher the patients in high-risk subgroups, and the higher the numbers of dead persons. In order to validate the prognostic predictive performance of RS for ccRCC, we also performed a ROC curve analysis based on data in the TCGA dataset. The AUC for the RS was 0.856 (Fig. 6d), indicating the satisfactory performance of the RS for prognostic prediction.

Discussion
Recently, more and more research has focused on the field of cancer [23]. As the most common histological type of RCC, clear cell renal cell carcinoma accounts for 90% of kidney neoplasms [24]. Although the development of operational methods and medical oncology, it remained eager to explore new reliable prognostic targets for survival prediction and treatment of ccRCC. Recently, some studies have reported the pivotal roles that m6A methylation played in tumorigenesis and development [4,25]. However, few researches focused on the association between m6A methylation regulatory genes and the prognostic prediction in ccRCC. Hence, our present study focused on the prognostic role of METTL14 in ccRCC.
In the current study, we explored the association among METTL14 expression, clinicopathologic parameters and patient survival based on TCGA database. The results revealed that METTL14 showed obviously lower expression in ccRCC tissues than that in adjacent normal tissues. Additionally, high METTL14 expression was significantly associated with better pathologic stage, histological grade, T stage, M stage, and satisfactory survival time. The METTL14 expression level was proved to be an independent predictive factor of overall survival for ccRCC patients. In order to explore how METTL14 was involved in ccRCC pathogenesis, we carried out GSEA between tissues with different METTL14 expression levels and found that high METTL14 expression was associated with some essential signaling pathways including ERBB pathway, MAPK pathway, mTOR pathway, renal cell carcinoma, pathway in cancer, TGF-β pathway and Wnt pathway. We further constructed a risk score that could predict the prognosis of ccRCC by combining the METTL14 expression with three clinicopathological variables including age, pathologic stage, and histological grade. Then, we confirmed the performance of RS in the TCGA database and the AUC for the RS was 0.856, which indicated the satisfactory ability of the RS for prognostic prediction. Finally, we investigated the relationship of METTL14 and related genes based on TCGA dataset and discovered that the five most negatively relevant genes are NUTF2, C19orf53, POLR2J, PGLS, and ADRM1. In contrast, the five most negatively relevant genes are ZNF24, UBE3A, IRE82, AQR, and UBR1.
Given the results of correlation analysis based on TCGA cohort, METTL14 and YTHDC1 were the two most relevant m6A regulatory genes (Pearson's r = 0.66), which indicated that the process of RNA m6A modification in ccRCC may be the result of combined action by m6A writer gene (METTL14) and reader gene (YTHDC1). Recently, a study suggested that numerous concurrence of CNVs in two m6a regulatory genes was observed in ccRCC, hence, m6A writer gene and reader gene may have synergistic effects on tumorigenesis and progression [25]. Unlike METTL3, METTL14 have been demonstrated to show tumorsuppressive functions in most types of cancer. In hepatocellular carcinoma, reduced METTL14 expression, but not METTL3 was associated with adverse patient prognosis and tumor metastasis by regulating the process of pri-miR126 [20]. A similar finding is indicated for glioblastoma stem cells, METTL14 knockdown could enhance glioblastoma stem cell growth, self-renewal and tumorigenesis by enhancing m6A modification of ADAM19 and decreasing its expression, suggesting a tumor suppressor role of METTL14 [26]. Recently, it is found that mutation of METTL14 in endometrial cancer cells could down regulate m6A mRNA methylation and enhance proliferation and tumorigenicity by regulating AKT activity [27]. In contrast, one recent study demonstrated that high METTL14 expression is present in normal hematopoietic stem/progenitor cells (HSPCs) and acute myelogenous leukemia (AML) cells with t(11q23), t(15;17) or t(8;21) translocation, and suppresses terminal myeloid differentiation and promotes leukemogenesis by positively regulating expression of MYB and MYC through m6A-based post-transcriptional regulation [28].
In our study, the results demonstrated that low METTL14 expression was significantly associated with poor clinical pathological variables and adverse survival time. Hence, METTL14 expression may be a potential prognostic marker of favorable survival and a tumor suppressor in ccRCC. Similar as our results, a previous study revealed that METTL14 could suppress P2RX6 mRNA and protein level. Moreover, the m6A-suppressed P2RX6 activation could promote renal cancer cells migration and invasion through ATP-induced Ca2 + influx modulating ERK1/2 phosphorylation and MMP9 signaling pathway [29].
Our results discovered that high METTL14 expression had association with ERBB pathway, MAPK pathway, mTOR pathway, renal cell carcinoma, TGF-β pathway and Wnt pathway by GSEA. These pathways  are crucial biological processes in tumorigenesis and development of ccRCC. As we all know, epidermal growth factor receptor (EGFR), as known as ErbB1, plays vital roles in enhancing tumor cell proliferation, suppressing apoptosis and contributing to the development and metastasis of renal cell carcinoma [30][31][32], and translation of EGFR has been recently found to be promoted by METTL3 in cancer cells [19]. A previous study reported an integrated molecular study of ccRCC (106 ccRCC cases). mTOR was found to be mutated in 6 cases [33]. Furthermore, 26 percent of ccRCC cases were discovered to have mutations involving in PI3K/ AKT/mTOR signaling. Nishida et al. [34] revealed that loss of TGFBR3, one of the TGF-β pathway molecular, could enhance metastatic abilities in ccRCC through TGF-β-dependent pathway. However, the relationship between METTL14 expression with these signaling pathways is the first to be reported, and the regulatory mechanism needs to be further investigated. In addition, Wnt pathway is another important pathway associated with high METTL14 expression in ccRCC by GSEA. A recent study proved knockdown of m6A writer gene METTL14 in gastric cancer cell enhanced proliferation and invasiveness by Wnt and PI3K-Akt signaling [35]. Therefore, further investigations need to be carried out to explore that whether progression in ccRCC could be affected by METTL14 through Wnt/ PI3K-Akt signaling pathway.