The prognostic values and clinical implications of m6A methylation regulators in epithelial ovarian cancer

Background: N6-methyladenosine (m6A) exists in both DNA and RNA modication. RNA m6A modication drives tumor initiation and metastasis through regulating cancer stem cells. However, the detailed mechanisms and the distinct m6A regulatory gene type underlying ovarian cancer mRNA modication remain unclear. Results: We analyzed copy number variation (CNVs) and mRNA expression of ovarian cancer cases in TCGA dataset to determine the copy number variation patterns of m6A regulatory genes, and the associations between m6A dysregulation or certain regulatory gene and overall survival. We showed the KIAA1429, as the writer gene, had highest amplication percentage and were associated with overall survival, or disease-free survival, whereas the associations with prognostic survival were independent of other prognostic factors including stage, grade, and debulking status of the tumour. Besides, METTL14 and YTHDC2, one as the writer gene and the other as reader gene, was also related with clinical outcome. Furthermore, subgroups analysis addressed that m6A upregulation especially writer gain contributed to prognosis in epithelial ovarian cancer. Conclusions: Collectively, our data addressed that m6A upregulation is likely to be critical to the clinical outcome, and KIAA1429 showed the highest correlation with clinical outcome in ovarian cancer among m6A regulatory genes.


Background
Epithelial ovarian cancer (EOC) is the most prevalent cause of death in gynecologic malignancies. The patients are usually diagnosed at advanced stages for the lacking of clinically meaningful biomarkers for early screening and understanding of molecular pathogenesis. Debulking surgery followed by chemotherapy can cause 50-80% complete clinical response in advanced stage cases. However, most cases will relapse and eventually die for recurrence or metastasis [1,2]. Recent studies focused on target therapy or immune therapy which is proven to prolong survival time in advanced-stage ovarian cancer patients. But the molecular pathogenesis and mechanism is underestimated. Hence, nding the functionally relevant molecular biomarkers and new therapeutic targets are still challenging.
M6A dysregulation exists in various process including stem cell fates [5], somatic cell reprogramming [6], and causes impaired cell proliferation, self-renewal capacity, developmental defects and cell death [7]. Articles showed that RNA m6A modi cation drives tumor initiation and metastasis through regulating cancer stem cells [8]. It is found to be associated with tumorigenesis in different cancers types, including breast cancer [9], hematologic malignancies [10,11].
Specially, the overall level of m6A regulators is found to be higher in ovaries compared to other organs and tissues [12], suggesting that m6A modi cation plays a critical role in reproduction system. One nding suggested that METTL3 has an oncogenic role in ovarian cancer development and aggressiveness [13]. The other nding showed that ALKBH5 was a candidate oncogene and a potential target for ovarian cancer therapy [14]. However, due to the heterogeneity and bias, comprehensive analysis is needed to analysis the expression of m6A regulators, the roles of m6A regulators in cancer progression and metastasis, the function in clinical survival, the correlation between m6A regulators and clinicopathological characteristics in ovarian cancer. Hence, we systematically analyzed the ovarian cancer cases from the TCGA database including the clinical and sequencing data of 579 cases, evaluated the expression spectrum of speci c m6A gene, and addressed the association between the genetic alteration or mRNA expression and clinicopathological features or clinical outcomes. In this study, we found that the alteration and mRNA expression level of m6A RNA methylation regulators play a critical role in EOC progression and metastasis.

Alteration of m6A regulatory genes in EOC patients
With the analysis from sequencing data in TCGA, only 20 independent samples were identi ed with m6A genes mutations (Table 1). TP53 was shown with high mutation (88%) in this cohort, which is in line with published literatures [15]. However, copy number variations (CNVs) of m6A regulatory genes that included deep deletion/shallow deletion/ copy number gain/ ampli cation, had high frequency in EOC cases (Fig. 1A), with loss of copy number (3584/5733) in most of the CNV events (Fig. 1B, Table 2). Moreover, for alteration only including ampli cation and deep deletion as shown in cBioportal database, m6A regulatory genes showed higher frequency of ampli cation (Fig. 1C, Fig S1A). Among them, KIAA1429 or VIRMA (writer) showed the highest frequency of ampli cation (85/576) in this cohort (Fig. 1D, Fig S1B), WTAP (14/576) and YTHDC2 (14/576) had relative high frequency of deep deletion (Fig. 1E, Table 2), which are m6A "writer" and "reader" genes, showing the important of writers in m6A regulation. Correlation between M6A regulatory genes alterations and clinicopathological factors To address the correlation between alteration (deep deletion or ampli cation) of m6A regulatory genes and the clinicopathological and molecular characteristics, we assessed the related prognostic factors including age, stage, grade, debulking status, and chemotherapy sensitivity. The results depicts that the ampli cation or deep deletion of m6A regulatory genes were associated with grade (Table 3, p = 0.016), but not with age, stage, debulking status and chemotherapy response. It is reported that missense mutations of TP53 were found to be the most frequent, while mutations of BRCA1 and BRCA2 showed the critical prognostic value in ovarian cancers [16]. TP53 and BRCA mutation play a critical role in the progression and pathogenesis of ovarian cancer. Here we evaluated the association between the variation of m6A regulatory genes and the mutation of these genes. In fact, the alterations of m6A regulatory genes were not signi cantly correlated with BRCA and TP53 mutations (Table S1), which indicated that m6A regulation affects ovarian cancer through other critical genes and molecular pathways. Association between CNVs of m6A regulatory genes and survival of EOC patients To identify the association between m6A regulatory genes alteration and ovarian cancer patients' survival, we evaluated the clinical value of CNVs on the overall survival (OS) and progression-free survival (PFS) among all ovarian cancer patients. Figure  The results also showed that writer genes, a group of methyltransferase enzymes, might be important for patient survival.
Next, we conducted survival analysis with subgroups of CNVs (copy number gain of writer genes and deletion of eraser genes) to test the prognostic role of m6A regulatory genes especially the writers. Figure 3A&3B depicts that in eraser deletion patients, the group combination with writer genes gain had worse OS and DFS than those without writer genes gain. This evidence supported the correlation between up-regulated m6A level (especially writer genes gain) and poor survival. By univariable Cox regression analysis adjusted for prognostic factors including stage, grade, BRCA mutation, age, debulking status, and subgroups, we addressed that alteration of m6A regulator genes was an important risk factor independent of other factors (Fig. 3C, Table S2, HR = 1.51, 95% CI (1.12 − 2.1), p = 0.008). Furthermore, KIAA1429 was independent risk factor for overall survival when conducting multivariate analysis adjusted for prognostic factors (Fig. 3D, Table S3, HR = 1.30, 95% CI (1.13 − 1.51), p < 0.001).

Association between m6A genes level and survival in EOC patients
The effects of alterations in m6A regulatory genes on the mRNA expression were evaluated. Figure 4A-4C addressed that the mRNA expression were signi cantly correlated with the diverse CNV patterns. Higher mRNA level was shown to be related with ampli cation or copy number gains for all the m6A regulator genes, while decreased mRNA expression was shown in shallow deletion or deep deletions for the m6A regulator genes.
Considering the important role of m6A regulator genes in pathological feature and survival above with the CNV data, we investigated the relationship between individual m6A regulator genes mRNA level and clinical survival. K-M analysis was conducted with the expression data in the TCGA dataset. Figure 5A and Fig. 5B depicts the unadjusted HRs and 95% con dence intervals for medium of gene expression levels with OS and PFS. 376 patients with restriction to samples with mRNA data were included for analysis. KIAA1429 (HR = 1.27, p = 0.063, Fig. 5C), YTHDC2 (HR = 1.31, p = 0.039, Fig. 5D) were most signi cantly associated with poor overall survival, and YTHDC2 (HR = 1.39, p = 0.011, Fig. 5E) showed signi cant correlation with poor PFS. However, KIAA1429 showed no signi cant association with the poor PFS but with the tendency (Fig. 5F). In general, the results were in consistent with the prognostic value of m6A genes alteration above that copy number gain or ampli cation of writers was associated with poor prognosis, illustrating that high m6A expression correlated the poorer outcome.

Validation of prognostic role of KIAA1429 in EOC patients
To validate the prognostic signi cance of KIAA1429 expression in ovarian cancer patients, we integrated the expression pro ling in Kaplan-Meier plotter online database (K-M plotter). We evaluated the prognostic value of KIAA1429 at mRNA level by Kaplan-Meier plotter analysis with cases enrolled from multiple GEO (Gene Expression Omnibus) datasets. Figure  Gene set enrichment analysis of KIAA1429 expression level As mentioned above, m6A upregulation, especially KIAA1429 gain, is important to promote cancer pathogenesis and progression. Herein we evaluated the molecular mechanism of m6A upregulation in EOC by GSEA analysis (Gene set enrichment analysis). We explored the gene set enrichment between samples with different KIAA1429 mRNA expression. The GSEA results showed that high KIAA1429 expression was associated with some critical pathway. Nuclear hormone receptor binding, small conjugating protein ligase activity, acid amino acid ligase activity, positive regulation of transcription from RNA polymerase II promoter, and positive regulation of RNA metabolic process were related with KIAA1429 upregulation in EOC patients, as Table 4 and Fig. 7A-E shown.
However, the further mechanism and regulation are needed to be illustrated. Discussion N6-methyladenosine (m6A) plays important roles in cancer progression and metastasis by controlling cell differentiation and cell pluripotency especially cancer stem cells [17]. It has been shown that m6A modi cation affect the breast cancer stem cells enrichment through changing the NANOG mRNA stability, relating with a shorter survival of breast cancer cases [18,19]. In human lung cancer, METTL3 can promote oncogene translation by decreasing epidermal growth factor receptor (EGFR) protein expression [20]. Besides, m6A regulatory genes also show an critical oncogenic role in acute myeloid leukemia [21], glioblastoma [8,22], and hepatocellular carcinoma [18].
Especially, Xiaoyao Lin et al. revealed that m6A modi cation of mRNAs regulates epithelial-mesenchymal transition (EMT) which is an important process for cancer cell metastasis [23]. Upregulation of METTL3, one of the writers in m6A regulatory genes, was proved to be shown in ovarian carcinoma, and signi cantly related with some prognostic clinicopathological factors including grade, lymph-node metastasis, and stage. Interestingly, METTL3 can promote EMT by upregulating the receptor tyrosine kinase AXL to promote ovarian cancer progression [13]. Moreover, other m6A modi cation genes, such as m6A demethylases FTO and ALKBH5, can regulate the Wnt/β-catenin pathway, contributing to PARP inhibitor resistance in BRCA-mutated EOC cells [24]. The published data showed that m6A regulation play an important role in ovarian cancer progression. In this article, we integratively analyzed the CNVs and mRNA expression data from m6A regulator genes in TCGA dataset. We explored the prognostic value of distinct writers, readers, and erasers, the association of copy number alteration and prognosis, and the speci c gene regulator for the ovarian cancer prognosis.
Surprisely, we found that KIAA1429, METTL14 and YTHDC2 were associated with cliniopathological features and clinical outcomes. Firstly, we addressed that ampli cation of KIAA1429 which showed highest among the thirteen genes' ampli cation, and METTL14 had correlation with clinicopathological grade and poor prognosis, which is consistent with the data shown in previous published studies. Besides, we elucidated the KIAA1429 mRNA expression and its potential clinical signi cance in epithelial ovarian cancer. The results depicted that KIAA1429 upregulation was related with poor prognosis, which is in accordance with the relationship between CNV and clinical prognosis. This was proved in breast cancer that KIAA1429 was highly expressed in breast cancer tissues and down-regulated in non-cancerous breast tissues. The prognosis, especially overall survival of breast cancer was associated with KIAA1429 expression. METTL14, one of the writers in m6A regulatory genes and one important component of the methyltransferase complex, plays a critical oncogenic role in human AMLs and liver cancer by positioning RNA substrates for methylation medicated by a complex comprising RBM15-WTAP-METTL3-METTL14 [8,18,[25][26][27]. Further-more, our ndings also found that ampli cation of YTHDC2, as an m6A modi cation reader, was correlated with poor clinical outcome in ovarian cancer, which is shown in studies that YTHDC2 promotes cancer metastasis [28] We identi ed the mechanism of KIAA1429 mRNA expression on clinical prognosis. Finally, 19624 genes were enrolled into the GSEA process. We addressed that positive regulation of transcription from RNA polymerase II promoter and positive regulation of RNA metabolic process were critical pathways, while articles showed that N6-methyladenosine promote cancer progression through regulating RNA metabolism [17]. As shown by the subgroups analysis, writers gain plus eraser loss showed association with poor prognosis in EOC. As methyltransferase enzymes, writer genes play an important role and the most part in m6A regulation process. The results implied that the regulation of m6A level might be associated with writers' expression and function.
Our articles showed high ALKBH5 loss in EOC but lower alteration including deep deletion and ampli cation. It is ambiguous about the role of ALKBH5 or the m6A regulatory genes loss in ovarian cancer. ALKBH5 can induce the stem cell phenotype in breast by mediating the m6A-demethylation of NANOG mRNA [19], and further studies proved that ALKBH5 regulated the malignant behavior of glioblastoma [29]. In ovarian cancer, ALKBH5 may be a potential oncogene in EOC by inhibiting autophagy of epithelial ovarian cancer through regulating miR-7 and BCL-2 [14]. Our article showed no relationship between ALKBH5 CNVs or mRNA expression level with prognostic value including overall survival and prognostic free survival, which is addressed with TCGA database and validated with K-M plotter database. The reason for the inconsistent may be for the limited clinical samples and cell lines which can be solved. In the current study, we used a vast spectrum of tumour samples from the publicity online studies and datasets. This indicated that our conclusions are general.
But there were still some shortcomings in this article. We used the CNVs and expression datasets from TCGA dataset in which the sample size may be limited.
Then we used GEO as the validation set, which includes a large sample size to improve repeatability. This approach may increase statistical power, but limitation exists insofar as some patients' information will be lost and data integrity will be impaired. Also, GEO samples have the possibility that patients are likely to have differences in clinical characteristics, and treatment. These differences affected the association between the m6A regulator genes and clinical outcome, with the possibility that we increased the false negative ratio.
In summary, this study addressed the critical role of m6A upregulation, especially writers' gain or ampli cation or mRNA upregulation, in clinical outcome of ovarian cancer, providing the prognostic role and possible targets for tumor diagnosis and prognosis. Analysis with large-scale genomic data from databases can reveal the role of m6A methylation regulator genes that is an important part of RNA stability, methylation and functions.

Data processing
There are 579 ovarian cancer patients with CNV data from the TCGA database. Segmentation analysis and GISTIC algorithm were used to identify the copynumber alteration [30]. There are 316 samples with mutation data and CNVs information. 379 cases were included with mRNA expression data and 3 replicates were removed by patient id. Clinical data and prognostic factors included survival data (overall survival and progression free survival), age, stage, grade, debulking status (optimal and suboptimal), and chemotherapy response (sensitivity and resistant).

Gene set enrichment analysis (GSEA)
Gene set enrichment analysis (http://software.broadinstitute.org/-gsea/index.jsp) is a method to detect the critical biological process in ovarian cancer by genes enrichment analysis comparing two distinct biological status. In our study, two groups were included according the quartile of KIAA1429 mRNA expression level, with one group are cases with lowest quartile expression and the other group with highest quartile expression. Process with normalized pvalue < 0.05 and FDR < 0.25 were identi ed to be signi cantly enriched.

Kaplan-Meier (K-M) plotter database
The Kaplan Meier plotter (http://kmplot.com/analysis) is based on an online database [31] and is capable to assess the association of genes on survival in ovarian cancer. K-M plotter was mainly used in validation for the prognosis value of KIAA1429 mRNA expression. The correlation of individual KIAA1429 mRNA expression and survival in ovarian cancer was analyzed online and presented with the hazard ratio, 95% con dence intervals and computed log rank pvalue.

Statistical analysis
All the statistical analysis was conducted by R language 3.6.1 (https://www.r-project.org/). We analyzed the association between clinicopathological features or prognostic factors with the CNV or mRNA expression level of m6A regulatory genes with chi-square test. Kaplan-Meier curve and Cox regression analysis were used for univariate or multivariate analysis and for evaluation of the prognosis value of m6A regulatory gene. Cox proportional hazard regression model was performed using R language. Results with a p-value < 0.05 were considered to be signi cant. Availability of data and material The datasets generated and analysed during the current study are available from TCGA, GEO, and Kaplan-Meier plotter that provide free online tools and resources Competing interest: The authors declare that they have no competing interests.  Evaluation of the effects of alterations in m6A regulatory genes on the mRNA expression  The prognostic value of KIAA1429 at mRNA level by Kaplan-Meier plotter analysis with cases enrolled from multiple GEO (Gene Expression Omnibus) datasets.