m6A RNA methylation regulators can contribute to malignant progression and impact the prognosis of bladder cancer

Abstract N6-methyladenosine (m6A) is the most common form of messenger RNA (mRNA) modification. An increasing number of studies have proven that m6A RNA methylation regulators are overexpressed in many cancers and participate in the development of cancer through the dynamic regulation of m6A RNA methylation regulators. However, the prognostic role of m6A RNA methylation regulators in bladder cancer (BC) is poorly understood. In the present study, we downloaded the mRNA expression data from The Cancer Genome Atlas (TCGA) database and the corresponding clinical and prognostic information. The relationship between m6A RNA methylation regulators and clinicopathological variables of BC patients was assessed by the Kolmogorov–Smirnov test. The expression of the m6A RNA methylation regulators was differentially associated with different clinicopathological variables of BC patients. The least absolute shrinkage and selection operator (LASSO) Cox regression model was then applied to identify three m6A RNA methylation regulators. The risk signature was constructed as follows: 0.164FTO − (0.081YTHDC1+0.032WTAP). Based on the risk signature, the risk score of each patient was calculated, and the patients were divided into a high-risk group and a low-risk group. The overall survival (OS) rate of the high-risk group was significantly lower than that of the low-risk group. The risk signature was not only an independent prognostic marker for BC patients but also a predictor of clinicopathological variables. In conclusion, m6A RNA methylation regulators can participate in the malignant progression of BC, and a risk signature with three selected m6A RNA methylation regulators may be a promising prognostic biomarker to guide personalized treatment for BC patients.


Introduction
Bladder cancer (BC) is the most common malignant tumor of the urinary system. In recent years, both morbidity and mortality have increased [1]. BC is classified as muscle-invasive BC (MIBC) and non-MIBC (NMIBC) according to the degree of tumor invasion. Approximately 30% of patients experience invasion of the muscular layer at the time of diagnosis [2]. The incidence of BC in males is higher than that in females [3]. BC is formed by the interaction of various internal and external factors [4]. So far, the regulatory mechanism of BC occurrence and development has not been fully elucidated.
Recently, reversible N6-methyladenosine (m6A) RNA methylation regulators have been found to provide a new method of post-transcriptional regulation [5]. Geneticists found the methylation of m6A in eukaryotic messenger RNA (mRNA) [6]. RNA methylation modification accounts for more than 60% of all RNA modifications, while m6A RNA methylation is the most common type of RNA methylation modification on higher biological mRNAs. m6A RNA methylation regulators play an important role in regulating mRNA splicing, translation and stability [7]. m6A RNA methylation regulators are widely distributed in various types of RNA, such as mRNA, transport RNA (tRNA) and ribosomal RNA (rRNA) [8]. The dysregulation of m6A regulators leads to decreased cell proliferation, a loss of self-renewal capacity, developmental defects and cell death [9]. m6A RNA methylation regulators are involved in the occurrence and development of many cancers, such as liver cancer [10,11], glioblastoma [12], osteosarcoma [13] and colorectal cancer [14]. The level of modification of transcript m6A is regulated by methyltransferases, binding proteins and demethylases [15]. The methyltransferases, known as the 'writers' , are responsible for the addition of the methyl group to the nitrogen on the sixth carbon of the aromatic ring of an adenosine residue and mainly include KIAA1429, METTL3, METTL14, RBM15, WTAP and ZC3H13 [16]. The m6A binding proteins, known as the 'readers' , can preferentially bind the RNA to confer its fate and regulate downstream functions and include HNRNPC, YTHDC1, YTHDC2, YTHDF1 and YTHDF2 [17]. The demethylases, known as the 'erasers' , can specifically target RNA m6A and mainly include ALKBH5 and FTO [18,19]. METTL3 and METTL14 can form a stable heterodimer core complex that functions in cellular m6A deposition on mammalian nuclear RNAs [20]. METTL14 expression is positively correlated with METTL3 expression and is highly expressed in breast cancer [21]. KIAA1429 is an interacting partner of the methyltransferase complex components [22]. WTAP is a regulatory subunit that is important for the localization of METTL3 and METTL14 into nuclear speckles [23]. FTO is associated with human obesity [24]. FTO is highly expressed in acute myeloid leukemia and can inhibit the differentiation of acute myeloid leukemia cells induced by all-trans retinoic acid [25]. ALKBH5 is highly expressed in glioblastoma, and ALKBH5 silencing inhibits the proliferation of glioblastoma [26]. METTL3 is highly expressed in BC and promotes BC proliferation in an m6A-dependent manner [27]. METTL3 can promote BC progression through the AFF4/NF-κB/MYC signaling network [28]. In vivo and in vitro experiments showed that the METTL3-mA-CDCP1 axis and chemical carcinogens could promote the tumorigenesis of BC [29]. Although an increasing number of studies have proven that m6A RNA methylation regulators play an important role in the occurrence and development of cancer, the relationship between m6A RNA methylation regulators and BC is still not fully clear. Many studies have predicted the prognosis of BC by constructing a prognostic signature based on miRNAs [30] and long noncoding RNAs (lncRNAs) [31], but BC prognosis has not been predicted by constructing a prognostic signature based on m6A RNA methylation regulators.
In the present study, we used data from The Cancer Genome Atlas (TCGA) database to analyze the expression of 13 m6A RNA methylation regulators in BC and their relationship with the clinicopathological features of BC. After Cox univariate analysis and least absolute shrinkage and selection operator (LASSO) Cox regression analysis, we selected three m6A methylation regulators to construct a risk signature. Then, we analyzed the prognostic role of the risk signature in BC. We found that m6A RNA methylation regulators play an important role in the malignant progression of BC and that the risk signature can predict the prognosis of BC patients.

Data sources
RNA-seq transcriptome data and the corresponding clinicopathological and prognostic information for 408 BC patients were obtained from the TCGA database (https://cancergenome.nih.gov/).

Selection of m6A RNA methylation regulators
There are 13 m6A-related genes in BC and 13 genes are recognized as m6A methylation regulators. Although 13 m6A-related genes have been systematically analyzed in gliomas [32], they have not been systematically analyzed in BC. We obtained mRNA expression data for 13 m6A-related genes from the BC cohort of the TCGA database. We compared the relationship among the expression of 13 m6A-related genes and the clinicopathological variables of patients with BC.

Bioinformatics analysis
The interaction between m6A RNA methylation regulators was analyzed using the STRING database (http://www. string-db.org/). To determine the prognostic role of m6A RNA methylation regulators in BC patients, we performed Cox univariate analysis using data from the TCGA database. Three m6A RNA methylation regulators were selected for further analysis, and a risk signature was developed using the LASSO Cox regression algorithm [33,34]. The risk score was calculated as follows: where Coefi is the coefficient and x i is the expression value of each selected gene. This formula was used to calculate the risk score for each patient.

Statistical analysis
All statistical analyses were performed using R software (version 3.5.1). The Wilcoxon's test was used to compare the expression of m6A RNA methylation regulators between cancer and normal tissues. The relationship between m6A RNA methylation regulators and the clinicopathological characteristics of BC patients was analyzed using the Kolmogorov-Smirnov test. The median risk score was used as a cut-off value to classify patients into a high-risk group and a low-risk group, and the Kaplan-Meier method was used to analyze the difference in overall survival (OS) between the high-and low-risk groups. The chi-square test was used to compare the relationship between the risk score and clinicopathological variables. Cox univariate and multivariate analyses of the relationship between clinicopathological variables and risk score were performed.

Expression of m6A RNA methylation regulators in BC
We analyzed the expression of 13 m6A RNA methylation regulators in BC using data from the TCGA database. We found that six m6A RNA methylation regulators (METTL3, HNRNPC, YTHDF2, YTHDF1, ZC3H13 and FTO) were differentially expressed in BC tissues compared with normal tissues. Of these regulators, the expression of METTL3 (P<0.001), HNRNPC (P<0.01), YTHDF2 (P<0.01) and YTHDF1 (P<0.001) in BC tissues was significantly higher than that in normal tissues, and the expression of ZC3H13 (P<0.01) and FTO (P<0.01) in normal tissues was significantly higher than that in BC tissues (Figure 1).

Relationship between the expression of m6A RNA methylation regulators and the clinicopathological features of BC
Because m6A RNA methylation regulators play an important role in tumorigenesis, we further studied the relationship between m6A RNA methylation regulators and BC clinicopathological variables to explore the clinical significance of m6A RNA methylation regulator expression. The results showed that expression of ALKBH5 (P=0.029), FTO (P=0.027), KIAA1429 (P=0.001), RBM15 (P=0.000), YTHDF1 (P=0.002) and METTL3 (P=0.011) was significantly correlated with grade (Figure 2A-F). As the grade increased, the expression of ALKBH5, FTO, KIAA1429, RBM15 and YTHDF1 increased, while the expression of METTL3 decreased. ALKBH5 expression was significantly up-regulated in M1 compared with M0 (P=0.006, Figure 2G). The expression of WTAP was significantly correlated with N stage (P=0.001, Figure 2H  To further understand the relationship among 13 m6A RNA methylation regulators, we analyzed the interaction and correlation between these genes. In the interaction network, WTAP is at the center and interacts with other m6A RNA methylation regulators ( Figure 3A). The correlation between YTHDC1 and METTL14 is the most significant ( Figure 3B).

Prognostic role of m6A RNA methylation regulators in BC
To investigate the effect of m6A RNA methylation regulators on BC prognosis, we performed Cox univariate analysis ( Figure 4A). The risk signature was established by selecting WTAP, YTHDC1 and FTO, and the risk scores of BC patients were calculated through the LASSO Cox regression model, where the coefficient of YTHDC1 was −0.081, the coefficient of FTO was 0.164, and the coefficient of WTAP was −0.032 ( Figure 4C). According to the median risk score, patients were divided into low-and a high-risk groups. There was a significant difference in the OS rate between the two groups, and the OS rate of the high-risk group was significantly lower than that of the low-risk group ( Figure 4B, P=2.304e-04). The 5-year OS rate was 31.8% in the high-risk group and 54.0% in the low-risk group.

The role of the risk signature in BC patients prognosis
The heat map shows the expression of the three selected m6A RNA methylation regulators and clinicopathological variables in the high-and low-risk groups. We found significant differences between the two groups in stage (P<0.05) and fustat (P<0.01) ( Figure 5A). We used Cox univariate and multivariate analyses to determine whether the risk signature was an independent predictor. Univariate analyses showed that age, stage, T stage, N stage and risk score were significantly linked with OS (P<0.001, Figure 5B). Multivariate analyses showed that age and risk score were significantly linked with OS (P<0.001, Figure 5C). These results suggest that the risk signature is a risk factor for BC patients and can independently predict the prognosis of BC patients.
To detect whether the risk signature adds prognostic value to the clinical system, the patients were also grouped according to clinicopathological variables (gender, grade, M stage, N stage, T stage and stage). The high-and low-risk groups were distinguishable in both the N0 (P=2.399e-02, Figure 6A) and N1/2/3 (P=3.312e-02, Figure 6B) patient groups, and high-risk patients had a significantly lower OS rate. We found that female patients (P=2.15e-05) and those with a high-grade (P=4.889e-04), M0 stage (P=1.82e-03), T3/4 stage (P=3.574e-04) and stage III/IV (P=7.62e-04) in the high-risk group had a significantly lower OS rate than those in the low-risk group ( Figure 6C-G); however, for male patients and those with a low-grade, M1 stage, T1/2 stage and stage I/II, no significant difference was present (Supplementary Figure S1), which may be due to the small sample size.

Discussion
BC is one of the most common malignant tumors in the urinary system. The high recurrence rate and drug resistance make the treatment of BC difficult in clinical research. Studies have shown that multiple genes always interact with each other to regulate the development of tumors [35]. m6A RNA methylation regulators play an important role in the development of cancer [36]. Therefore, it is necessary to explore the influence of m6A RNA methylation regulators on BC.
In the present study, we first analyzed the expression of m6A RNA methylation regulators in BC and normal tissues and the relationship between their expression and different clinicopathological variables. The expression of the m6A RNA methylation regulators was differentially associated with different clinicopathological variables of BC patients. METTL3 is an effective therapeutic target for several cancers, such as pancreatic cancer [37], melanoma [38] and lung adenocarcinoma [39]. METTL3 is overexpressed in hepatocarcinoma cells, and its overexpression is associated with a poor prognosis of liver cancer patients [40]. METTL3 is highly expressed in ovarian cancer, and the expression of METTL3 is significantly associated with grade, pT status, pN/pM status and International Federation of Obstetrics and Gynecology (FIGO) stage [41]. All the above studies indicate that METTL3 is a potential oncogene, but in our study, we found that although METTL3 is highly expressed in cancer tissue, it is expressed at low levels in high grade. KIAA1429 is highly expressed in breast cancer and can promote breast cancer proliferation and metastasis in vivo and in vitro [42]. We found that the expression of KIAA1429 did not differ between BC and normal tissues but was meaningful in grade and highly expressed in high grade. YTHDF1 is overexpressed in colorectal cancer, and YTHDF1 knockdown significantly inhibits Wnt/β-catenin pathway activity in colorectal cancer cells [43]. We found that YTHDF1 is highly expressed in BC and high grade, suggesting a potential carcinogenic effect of YTHDF1 in cancer. ALKBH5 is highly expressed in ovarian cancer. Silencing ALKBH5 can enhance autophagy and inhibit the proliferation and invasion of ovarian cancer cells [44]. We found that ALKBH5 expression did not differ between BC and normal tissues, but the expression of ALKBH5 was positively correlated with grade and M1 stage. These results indicate that ALKBH5 is a potential oncogene that promotes tumor progression. HNRNPC knockdown inhibits the proliferation of breast cancer cells [45]. In our study, HNRNPC was highly expressed in BC. FTO is highly expressed and linked with a poor prognosis in breast cancer [46], but we found that FTO is expressed at low levels in BC. The carcinogenic effects of m6A regulatory genes appear to be controversial. This may be due to the heterogeneity of cancer, which results in the differential expression of m6A regulatory genes in different cancers. The same gene has different roles in different cancers.
Whether m6A RNA methylation regulators have prognostic value in cancer is of great significance [47]. After Cox univariate analysis and LASSO Cox regression analysis, we selected three m6A RNA methylation regulators to construct a risk signature that stratified BC patients into a low-risk group and a high-risk group. Univariate analyses showed that age, stage, T stage, N stage and risk score were significantly associated with OS. Multivariate analyses showed that age and risk score were significantly linked with OS. When patients were also grouped according to clinicopathological variables, the risk signature could also distinguish the OS outcomes of patients with N0 stage, patients with N1/2/3 stage, females, patients with high-grade, patients with M0 stage, patients with T3/4 stage and patients with stage III/IV. The study also has some limitations. We only use the data of TCGA database to verify Cox regression model, and we should also verify it in other public databases.

Conclusions
m6A RNA methylation regulators are closely related to malignant clinicopathological features, and a risk signature with three selected m6A RNA methylation regulators can independently predict the prognosis of BC patients. Our study provides important evidence for further study of m6A RNA methylation regulators in BC.