Molecular Characterization and Clinical Relevance of m6A RNA Methylation Regulators in Cervical Cancer


 Background: RNA modification, such as methylation of N6 adenosine (m6A), plays a critical role in many biological processes. However, the role of m6A RNA modification in cervical cancer (CC) remains largely unknown. Methods: The present study systematically investigated the molecular signatures and clinical relevance of 20 m6A RNA methylation regulators (writers, erasers, readers) in CC. The mRNA expression and clinical significance of m6A-related genes were investigated using data from The Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression (GTEx) cervical cancer cohort. Mutations, copy number variation (CNV), differential expression, gene ontology analysis and the construction of a mRNA-microRNA regulatory network were performed to investigate the underlying mechanisms involved in the abnormal expression of m6A-related genes. Results: We found inclusive genetic information alterations among the m6A regulators and that their transcript expression levels were significantly associated with cancer hallmark-related pathways activity, such as the PI3K-AKT signaling pathway, microRNAs in cancer and the focal adhesion pathway, which were significantly enriched. Moreover, m6A regulators were found to be potentially useful for prognostic stratification and we identified FMR1 and ZC3H13 as potential prognostic risk oncogenes by LASSO regression. The ROC curves of 3, 5 and 10 years were 0.685, 0.726 and 0.741, respectively. The specificity for 3, 5 and 10 years were 0.598, 0.631 and 0.833, the sensitivity were 0.707, 0.752 and 0.811, respectively. Conclusions: Multivariable Cox regression analysis revealed that the risk score is an independent prognostic marker and can be used to predict the clinical and pathological features of CC.


Introduction
Cervical cancer (CC) is ranked fourth worldwide for causing morbidity and mortality in females [1].
Genetic and epigenetic alterations in CC are associated with complex factors, such as human papilloma virus (HPV) infection [2], tobacco abuse and oral contraceptives. On diagnosis, most CC patients have reached a late stage of cancer development when metastasis will most likely occur. Although surgical resection and radiotherapy techniques have been enhanced in recent years, the clinical prognosis and postoperative quality of life of patients with CC remain poor [3]. Thus, early diagnosis is critically important and therefore it is imperative to identify new potential biomarkers for CC. m6A RNA methylation is the most common internal modi cation of mRNA. Many proteins involved in m6A have been identi ed, including methyltransferase ('writer'), RNA binding protein ('reader') and demethylase ('eraser') [4]. These proteins and other markers have been shown to play important roles in the mRNA life cycle, as well as in cellular development and disease progression. Once the mechanisms involved in m6A modi cation are disrupted, they can induce a series of conditions including tumors, neurological diseases and a delay in embryonic developmental [5].
Recent advances in the study of m6A RNA methylation regulators has recognized their roles in the biological progression of various cancers. Lin et al. found that the m6A methyltransferases METTL3 can promote translation in lung cancer [6] and METTL14 has been shown to inhibit metastasis of hepatocellular carcinoma by regulating the level of miRNA modi cation [7]. In breast cancer, hypoxia stimulation promotes the expression of ALKBH5 and reduces the level of m6A modi cations, thereby improving the stability of mRNA and ultimately increasing the proportion of breast cancer stem cells [8].
However, little is known about the role of m6A RNA methylation regulators in CC. In our study, we employed the mRNA expression and clinical information from The Cancer Genome Atlas (TCGA) and the Genotype-Tissue Expression (GTEx) CC cohort [9] with the aim of establishing the relevance of m6A RNA methylation regulators dysregulation and their clinical prognosis value for CC.

Dataset
The CC mRNA expression pro les and clinical information were obtained from TCGA databases (http: cancergenome.nih.gov/) via the GDC-client, a total of 307 samples including 306 CC tissues and 3 paired adjacent non-cancerous tissues were included. Detailed clinical pathological of those CC patients are shown in Table 1. In order to avoid the ine ciency in various differential analyses caused by imbalance between the tumor and normal data, we acquired 10 cervical non-cancerous samples from the GTEx project, which produced RNA-seq data for over 8,000 normal samples, albeit from unrelated donors. We downloaded the GTEx project dataset from the UCXC Xena project (https://xena.ucsc.edu/public/). Both RNA-seq data in the TCGA project and GETx project were processed into the Fragments Per Kilobase of transcript per million mapped, reads (FPKM)-based gene expression, clinical information concluded the survival status, grades and survival time data. The two datasets were integrated via the R-package 'limma'.  [5]. We analyzed the integrated mRNA expression pro le of CC (including 306 tumors and 13 normal specimens) to explore the difference between normal and tumor samples of m6A target gene expression via the R-package 'limma'.
Constructing the microRNA-mRNA network The microRNAs, which interacted with m6A RNA methylation target genes, were predicted by Funrich software (version 3.0) (http://www.funrich.org/) [12]. A total of 433 microRNAs were found. We selected the top 100 closely interrelated microRNAs to each m6A RNA methylation target gene based on the score of interaction. Finally, the microRNA-mRNA network was visualized using Cytoscape software (ver. 3.6.1) [13].

Bioinformatic analysis
To investigate the molecular mechanisms of m6A RNA methylation regulators in CC, we divided 306 CC patients into 2 groups (cluster1 and cluster2) based on the characteristic similarity of m6A regulator expression via the 'ConsensusClsterPlus' R-package; we used 50 iterations and an 80% resample rate. Gene ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed to investigate the potential function of differentially expressed genes (DEGs) in the Database for Annotation, Visualization and Integrated Discovery (David)(https://david.ncifcrf.gov/). The association between m6A target gene expression and clinical pathological was assessed using a chisquared test.
Identi cation of prognostic m6A methylation regulators Univariate Cox regression analysis was employed to investigate the relationship between m6A target gene expression levels and the patients overall survival times (OS). We selected 2 genes with a P-value < 0.05 and constructed a prognostic risk score model using the LASSO Cox regression method [14]. The risk score model calculation formula was: Risk score = . Where the coe cient Coef i was determined by minimum criteria and X i is the expression level of each selected gene. According to the risk score calculation formula, we calculate the risk value of each CC patient, and all the patients were divided in the dataset into high-risk and low risk groups. The receiver operating characteristic (ROC) and Kaplan-Meier (KM) curve were plotted to evaluate the prognostic ability of the risk model.

Nomogram Analysis
We used rms R package to perform nomogram analysis by including those factors that signi cantly associated with OS of CRC patients in multivariate analysis. Calibration plot was applied to estimate the discrimination between actual and nomogram predicted OS probability.

Results
Expression of m6A regulators and clinical pathological correlation analysis Figure 1 shows a schematic diagram of the work ow analysis strategy used in the present study.
The aim was to examine the m6A molecular regulatory mechanisms involved in CC biological procession; we rst checked the correlation and interaction between the expression of individual m6A regulators. Spearman correlation analysis revealed that m6A regulators do not function independently; there was a high connection among writers, erasers and readers. For example, writer METTL3 was signi cantly positively correlated with reader HNRNPC and eraser FTO (r > 0.5), while eraser FTO was negatively correlated with writer RAM15 (r = -0.46). Taken together, these results indicated that m6A regulators play a critical role in the progression of CC, based on cross talk among writers, erasers, and readers. The Spearman correlation network results are presented in Fig. 3.

m6A target gene expression and genetic alterations in CC patients
Considering the importance of genetic alterations in various cancer types, we investigated m6A target gene mutations, copy number variations (CNVs) and alterations in mRNA expression in CC. The aim was to con rm whether genetic alteration could affect mRNA expression levels. Among the 278 CC cases, we determined that the m6A regulators were commonly dysregulated about 30% of cases (82) (Fig. 4A, B).
The m6A regulators overall mutation frequency in cervical squamous cell carcinoma was signi cantly higher than for normal cells (NC). Alterations in fusion only affected carcinoma samples. LRPPCR, ZC3H3 and VIRMA regulators displayed higher mutation frequencies among the m6A regulators. CNVs re ected 3 types including ampli cation, deep deletion and multiple alterations. We found that FMR1, METTL4 and VIRMA showed widespread CNVs ampli cation in CC. In contrast, ZC3H3 and RBM15B exhibited extensive deep deletions. According to the mRNA expression z-scores relative diploid samples results, we found that a m6A target gene with CNV ampli cation showed higher expression in cervical carcinoma cells (e.g. ELAVL1 and METTL4), while the genes with a CNV deletion had lower expression levels (e.g. ZC3H13 and METTL16) (Fig. 4C). All the data indicated that dysregulation of m6A RNA methylation regulators might be a prominent mechanism involved in cervical tumor progression.
Constructing a microRNA-mRNA regulatory network MicroRNAs are small, non-coding RNAs that function as regulators of gene expression [15]. Dysregulation of the expression of microRNA or microRNA mutations results in a gain or loss of microRNA functions, which leads to downregulation or upregulation of target proteins [16]. To understand the potential capability of m6A regulators, we constructed a microRNA-mRNA regulatory network [17]. We identi ed microRNAs that were related to m6A target genes by utilizing Funrich software and choose the top 100 microRNAs which were closely interrelated to each m6A target gene. Finally, the microRNA-mRNA regulatory network consisted of 100 microRNA and 17 mRNA nodes and was visualized using Cytoscape. As shown in Fig. 5, the key microRNA and mRNA molecules in the network had the highest degree of expression.

Analysis of consensus clustering of m6A regulators
We established consensus clustering of m6A regulators analysis and k = 2 seems to be the most proper parameter according to conjoining analysis consensus CDF and delta area results. We noticed that 151 out of 306 tumor samples clustered into the rst subgroup, 155 samples clustered into the second subgroup, namely cluster1 and cluster2 datasets. We set |logFC > 2| and P-value < 0.05 as the cut off values and gained 326 DEGs in cluster1 and 325 DEGs in cluster2. Then, we annotated the function of those DEGs by GO analysis using the DAVID database [18]. The results showed that the DEGs in the 2 subgroups were enriched into different biological processes, namely cellular and molecular components. Compared to the cluster1 dataset, the DEGs in cluster2 were more likely correlated with transcription, a DNA-template and metal ion binding GO terms [19]. A similar difference in the corresponding signaling pathways was also revealed by KEGG pathway analysis. The KEGG results of cluster1 found signi cantly enriched terms such as purine metabolism, endocytosis and mRNA surveillance pathway. In the cluster2 subtype, the PI3K-AKT signaling pathway [20] involved in cancer, microRNAs in cancer and focal adhesion pathways were signi cantly enriched. The results are shown in Fig. 6.

Prognostic value of m6A RNA methylation regulators
To investigate the prognostic role of m6A RNA methylation regulators in CC, we applied multivariate Cox proportional hazard model and identi ed that FMR1,and ZC3H13 were signi cantly associated with a patient's OS (P < 0.05; Fig. 7). FMR1 was risk genes for CC with HR < 1 and ZC3H13 was a protective gene with HR > 1. We picked out FMR1 and ZC3H13 to constructed a CC risk prognostic prediction model by using LASSO Cox regression, the risk score was calculated by the following formula: risk score = -0.0841 × FMR1 expression value + 0.1322 × ZC3H13 expression value. We calculated the risk score for all CC patients, and by using the median risk score as the cutoff value; we divided all patients into a high-risk group (n = 145) or a low-risk group (n = 146). The K-M survival curve showed that the risk score was a signi cant prognostic factor for OS (P < 0.05; Fig. 8A). Patients with low-risk scores were signi cantly longer lived than those with a high score. Meanwhile, the area under the ROC curve of this model was AUC = 0.685 (3 years), AUC = 0.726 (5 years), AUC = 0.741 (10 years) (Fig. 8B) and the speci city and sensitivity values ( Table 2) validated that our two-gene prognostic risk score model showed a good performance in predicting overall survival, and that the risk score had prognostic value for CC. We also explore if CC patient's risk score were associated with the clinical-pathological features which mainly included age, grade and stage. While the risk score was not signi cant different between elder and younger patients, lower grade and higher grade, as well as among stage1 to stage4. Besides, multivariate Cox regression analysis shown that the m6A regulators reliably predict CC patients prognosis independent of stage ( Fig. 9).

Construction and Validation of Nomogram
By four signi cant factors including age, grade, stage and risk score in multivariate Cox-regression analysis, we constructed a prognostic nomogram to predict the OS probability at 1, 3, and 5yeas. The calibration plots displayed that the risk score model have a good performance in predicting CC prognosis. (Fig. 10)

Discussion
A number of studies have revealed that elements of m6A RNA methylation involve complicated mechanisms that participate in the induction and development of multiple cancer types, such as hepatocellular carcinoma, breast and pancreatic cancer, and so on [21,22]. Therefore, we hypothesized that changes in the abundance of expression of these elements might be potential prognostic biomarkers for CC.
To explore the best candidates among all THE m6A-related players, the 2 CC datasets were surveyed and ZC3H13 and FMR1 were found to be the most commonly altered genes involved in m6A dynamics in CC.
Remarkably, CC samples displayed higher ZC3H13 and FMR1 transcript expression levels compared to NCs in both cohorts. The 2 m6A related genes have already been documented, but their roles in tumor inhibition and progression remain unclear. Zhu et al. found that ZC3H13 was associated with cell proliferation and invasion in colorectal cancer [23]. In contrast, the another m6A regulators have not been proven to be associated with cancer (FMR1) according to our extensive PubMed search. However, both were signi cantly associated with a patient's prognosis in CC. ZC3H13 (Zinc Finger CCCH-Type Containing 13) is a protein encoding gene and related component of the WMM complex, which is important for m6A methylation of mRNAs and can promote the e ciency of mRNA splicing and RNA processing at 3'-UTR [23]. The expression of ZC3H13 was upregulated in CC samples and was found to be an independent prognostic factor for OS. In contrast, FMR1 was were downregulated in CC samples and upregulated expression of them may be positively associated with OS. Furthermore, multivariable Cox regression analysis veri ed their prognostic role in CC. FMR1 was shown to regulate alternative splicing that speci cally recognizes and binds to m6A-containing RNAs. GO and KEGG analysis showed that the 2 subgroups not only were signi cantly correlated with clinical pathologicy but also with biological processes and pathways, such as cluster1 which was mostly enriched during endocytosis and mRNA surveillance, while cluster 2 was correlated with many cancer-related pathways, such as PI3K-Akt signaling and microRNAs in cancer pathways. The PI3K-Akt signaling pathway has been reported to be associated with different types of cancer by Deng et al. It is known that activation of the PI3K-Akt signaling pathway causes DNA damage, epithelial-mesenchymal transition and cancer stem cell promotion in ovarian cancer (OV) [25; 26]. In CC, inhibition of the PI3K/AKT signaling pathway can suppress tumor progression [27].
A 2-gene risk prognostic model based on m6A target genes for predicting OS of CC patients was developed. Patients were divided into high-risk and low-risk groups based on the risk score calculated by the 2-m6A related gene risk model. Survival, ROC curves and nomogram analysis indicated that the 2-m6A related gene risk prognostic model had high predictive accuracy. ZC3H13 and FMR1 transcript levels accurately discriminated between NCs and CCs, and may well serve as novel biomarkers for patient management. Multivariable Cox regression analysis proved that the risk score is an independent prognostic marker that predict the clinical pathological features of CC.
To explore the m6A regulatory mechanism in biological procession, we predicted the related microRNAs of the 20 m6A target genes in the Funrich database and constructed a complex miRNA-mRNA regulatory network. We found that these m6A target genes expression levels were directly regulated by many microRNAs, including mir-23a/b, mir-4262 and mir-30c/e. A number of studies have indicated that these microRNAs were signi cantly interrelated with cervical cancer procession. Several dysregulated microRNAs have been previously reported to exert their functions via miRNA-mRNA regulatory mechanisms, by dysregulation of ZC3H13, FMR1 activity [28]. Here, we built a miRNA-mRNA regulatory network that may be involved in the progression of CC, but further research is required to elucidate the molecular mechanisms involved.
There were a number of limitations to the present study. First, the clinical data of our samples were not complete. Second, the interaction of microRNA and mRNA was predicted by software, and requires further studies. Third, there was a lack of veri cation of animal experiments, thus further research is needed to support our study conclusions.
In conclusion, our work has proven that m6A RNA modi cation is signi cant for tumor initiation and progression of CC and its related genes. In particular, FMR1 and ZC3H13 may be potential prognostic biomarkers and be useful for the development of an effective treatment strategy for CC.

DATA AVAILABILITY STATEMENT
Our data was obtained from the publicly available datasets TCGA and UCXC Xena.

AUTHOR CONTRIBUTIONS
This study was conceived and designed by JC, XG, JH and XM. JC was responsible for data downloaded from the online database. JC and XM performed the experiments and analyzed data using the R program. XG and JH critically read and commented on the manuscript.