m6A RNA methylation regulators were associated with the malignancy and prognosis of ovarian cancer

ABSTRACT

N6-methyladenosine (m6A) RNA methylation regulators play a regulatory role in tumor pathogenesis and development. However, the role of m6A regulator genes in ovarian cancer (OC) has not been fully elucidated. This study aims to investigate the mRNA expressions, clinicopathological features, and prognostic values of m6A regulators in OC. Here, we demonstrate that the 17 m6A RNA methylation regulators are differentially expressed in ovarian cancer and normal tissues. By using consensus clustering, all ovarian cancer patients can be divided into two subgroups (cluster 1 and 2) based on the expression of 17 m6A RNA methylation regulators. Using Gene Set Enrichment Analysis, we identified that cluster 1 was most connected to oxidative phosphorylation pathways. Regression models identified that prognosis is associated with HNRNPA2B1, KIAA1429, and WTAP. qRT-PCR result show that the expression trends of HNRNPA2B1 and KIAA1429 are consistent with the predicted results. Multivariate Cox regression analysis results show that the risk score was an independent predictive factor in OV. The overall survival of high-risk patients was significantly shorter than that of low-risk patients. ROC curve analysis showed that the prognostic signature precisely predicted the 5-year survival of OV patients. A nomogram was developed to predict each patient’s survival probability and well calibrated and showed a satisfactory discrimination. Dendritic fraction, macrophage fraction, and neutrophil fraction showed higher fraction in high-risk patients. In conclusion, m6A RNA methylation regulators are vital participants in ovarian cancer pathology, and three-gene mRNA levels are valuable factors for prognosis predictions.

Graphical Abstract

KEYWORDS:

• Ovarian cancer
• m6A
• prognosis
• immune cells infiltration
• tcga

Background

Epigenetic modification is a change in the expression of a nucleotide sequence. Previous studies have shed light on epigenetic pathways such as histone modification, chromosome remodeling, DNA methylation, and non-coding RNA regulation. RNA-level modification can be accomplished by N7-methyladenosine, N1-methyladenosine, pseudouridine, 5-methylcytosine, N6, 2ʹ-O-methylation, and 2ʹ-O-dimethyladenosine (m6A). Among them, m6A is a form of RNA methylation discovered in the 1970s.

RNA methylation is a dynamic and reversible process involving methyl-transferases ‘writers’, binding proteins (‘readers’), and demethylases (‘erasers’), just like DNA methylation. The prominent m6A methylation regulators include ‘writers’ like METTL3, METTL14, WTAP, KIAA1429, RBM15, and ZC3H13; ‘readers’ like YTHDC1, YTHDC2, YTHDF1, YTHDF2, and HNRNPC and ‘erasers’ like FTO and ALKBH5. More and more studies have found that regulators in m6A RNA methylation are associated with tumorigenesis. For example, Chen M et al. found that METTL3 promoted liver cancer progression [1]. Li J et al. found that FTO promoted the growth of lung cancer cells by regulating the m6A level of USP7 mRNA [2]. Hua W et al. found that METTL3 promoted ovarian carcinoma growth and invasion [3]. Mei Chen et al. found that m6A RNA methylation regulators can promote malignant progression and affect the prognosis of bladder cancer [4]. Shuai Ma et al. found that there is a meaningful interaction between m6A RNA methylation and non-coding RNA in cancer [5]. The mechanism of action between m6A RNA methylation and specific tumors has increasingly become the focus of in-depth research.

Molecular technologies for the treatment of ovarian cancer are ongoing. Ovarian cancer poses a great threat to women’s life. Abnormal modification of RNA may act in tumor development. The study by Zhao Ma et al. found that METTl3 independently of METT114 and WTAP regulates m6A in endometrioid epithelial ovarian cancer [6]. Jie Li et al. found that YTHDF2 is a protein inhibited by miR-145, which can regulate the proliferation, apoptosis, and migration of ovarian cancer cells [7]. Takeshi Fukumoto et al. found that N 6-Methylation of Adenosine of FZD10 mRNA Contributes to PARP Inhibitor Resistance [8]. These research results all show that m6A RNA methylation plays a pivotal role in the development and treatment of ovarian cancer.

The aim of our study is to find a direct and fixed algorithm for batch prediction of m6A RNA regulators that have an impact on ovarian cancer. Our goal is to provide a new molecular mechanism for the development of ovarian cancer.

Materials and methods

Data Sources

RNA-seq data and corresponding clinicopathological data were obtained from TCGA (http://cancergenome.nih.gov/http://cancergenome.nih.gov/) including 379 OC patients [9]. The expression dataset (N = 379 for ovarian cancer from TCGA and N = 88 for normal ovarian tissues from GTEx) [10] was downloaded from the UCSC Xena project (http://xena.ucsc.edu/http://xena.ucsc.edu/) [11]. Data were renormalized based on total reads for each sample to generate RPKM (Reads Per Kilobase of transcripts per Million mapped reads) [12] and then the expressions of between normal ovarian and OC tissues were compared. ‘limma’ package [13] was performed to solve the imbalance between the tumor and normal data and then analyzed different expressions among normal ovarian and OC.

Selection of m6A RNA methylation regulators

17 m6A RNA methylation regulators from published papers have been identified for subsequent studies.

Bioinformatic analysis

We divided the samples into different groups using ‘ConsensusClusterPlus’ [14–16]. Enrich Database (https://amp.pharm.mssm.edu/Enrichr) was conducted for functional analysis. Interactions among m6A RNA methylation regulators were analyzed by using the Search Tool for the Retrieval of Interacting Genes Database (STRING) database (http://www.string-db.org/). Gene Set Enrichment Analysis (GSEA) [17] was based on different subgroups of OC for identifying the functions. Univariate Cox regression analysis [18] was used to determine the prognostic value of m6A RNA methylation regulators. From this, it has been proved that three genes were significantly associated with survival through this study (P < 0.05), which are selected for further functional analysis and development of a potential risk signature with the LASSO Cox regression algorithm [19]. The minimum criteria [20] was set to determine the 3 genes and their coefficients. The risk score [21] for the signature was calculated by using the formula: Risk score = ΣCoefi*xini = 1, where Coefi is the coefficient, and xi is the z-score-transformed relative expression value of each selected gene. This formula was used to calculate a risk score for each patient in TCGA datasets. To reveal potential Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways [22] of the high- and low-risk groups and three prognosis-related genes, GSEA [23] was utilized to find enriched terms in C2 and in TCGA‐OV database. And p < 0.05 were considered to be statistically significant.

Complex cancer genomic profiles are accessible from the cBioPortal tool (http://www.cbioportal.org/http://www.cbioportal.org/) [24,25,].

The abundant six subtypes of tumor-infiltrating immune cells including CD4 T cells, CD8 T cells, B cells, neutrophils, macrophages, and dendritic cells can be visualized by TIMER online database [26].

Construction and validation of the nomogram

A nomogram and calibration curve were constructed by using the ‘rms’ package [27] on R. Harrel’s concordance index (C-index) [28] was measured. The nomogram was then subjected to bootstrapping validation (1,000 bootstrap resamples) to calculate a relatively corrected C-index [29]. Finally, we used decision curve analysis (DCA) to determine the clinical usefulness [30].

Statistical analysis

Patients were divided into two clusters by consensus expression of m6A RNA methylation regulators. Chi-square tests were used to compare the distribution of grade, race, stage age, and tumor status between the two risk groups.

Patients were divided into two groups based on the median levels of risk score [31]. This prognostic model and patient survival information were merged. Kaplan–Meier survival curves were used to compare the prognostic ability of the prediction models [32]. Area under the curve (AUC) [33] value for the Receiver operating characteristic (ROC) curves of each prognostic model was calculated by survival.ROC package in R. Besides, univariate and multivariate Cox regression analyses were conducted to compare the hazard ratio (HRs) [34] of prognostic models and important clinical features for OC. Differences among clinical parameters (age, grade, stage, and tumor status) were tested by using independent t-tests and P < 0.05 were considered to be statistically significant.

qRT-PCR to verify the expression of hub genes in clinical samples

16 ovarian cancer tissues and 8 normal ovarian tissues were obtained. The acquisition process and sample storage and processing are the same as before [35]. The extraction of RNA and the implementation of qRT-PCR were also operated strictly in accordance with the protocols. TABLE S1 showed the primer sequences.

Results

The present study aims to identify m6A RNA modulators associated with prognosis in OC. We analyzed the association between each m6A RNA methylation regulator and the clinical features of OC. Survival analysis showed that ALKBH5, METTLE14, METTLE16, YTHDF1, YTHDF2, YTHDF3, and ZC3H13 were significantly associated. The expression profiles of METTLE3, YTHDC2, and YTHDF3 were all related to the clinicopathological features in OC. By using consensus clustering, all ovarian cancer patients can be divided into two subgroups (cluster 1 and 2) based on the expression of 17 m6A RNA methylation regulators. Using Gene Set Enrichment Analysis, we identified that cluster 1 was most connected to oxidative phosphorylation pathways. Regression models identified that prognosis is associated with HNRNPA2B1, KIAA1429, and WTAP. Combined with clinical characteristics, a nomogram for predicting prognosis was constructed and calibrated. The immune infiltration analysis was then performed to assess the impact of the three m6A RNA modulators on tumor behavior.

Expression of m6A RNA methylation regulators was correlated with OC clinicopathological features

The correlation between each m6A RNA methylation regulator and OC clinical features was analyzed. These clinical features include age, grade, stage, and tumor-status. FMR1, METTL16, RBM15, and WTAP had significant correlations with age (Supporting Figure 1). FTO and YTHDC2 had significant correlation with OC grade (Supporting Figure 2). ALKBH5, METTL3, METTL14, METTL16, RBM15, and YTHDF1 had significant correlation with OC stage (Supporting Figure 3). HNRNPA2B1 and METTL16 had significant correlation with OC status (Supporting Figure 4). For ovarian cancer, BRCA1 gene mutation is a common pathogenic factor and marker for ovarian cancer [36]. The expression of the m6A RNA methylation regulator in ovarian cancer patients with and without BRCA1 mutation was also analyzed, and it has been found that the expression of both FMR1 and YTHDF1 showed differences between BRCA1-mutation and non-BRCA1-mutation groups (Supporting Figure 5). Seventeen m6A RNA methylation regulators in OC tissue samples and normal tissue samples were analyzed later, and it has been shown that all the regulators are differentially expressed (Figure 1). In addition, survival analysis found that the overall survival (OS) and progression-free survival (PFS) of ALKBH5, METTLE14, METTLE16, YTHDF1, YTHDF2, YTHDF3, and ZC3H13 were meaningful (Supporting Figure 6 and Supporting Figure 7). In order to assess the diagnostic value of 17 m6A RNA regulators, an ROC curve was generated by using expression data from ovarian cancer patients and healthy participants (Supporting Figure 8). The area under the ROC curve (AUC) indicates a modest diagnostic value.

Figure 1. Expression of m6A RNA methylation regulator in OC samples and normal samples. (a) Heatmap showed that the 17 m6A RNA methylation regulators expressed differently between OC samples and normal samples. (b) Expression level of 17 m6A RNA methylation regulators in OC samples and normal samples

Figure 2. Differential clinicopathological features and overall survival of OC in the cluster ½ subgroups. (a) Consensus clustering matrix of 379 TGGA samples for k = 2. (b) Consensus clustering cumulative distribution function (CDF) for k = 2 to 10. (c) Relative change in area under CDF curve for k = 2 to 10. (d) Heatmap and clinicopathologic features of the two clusters defined by the m6A RNA methylation regulators consensus expression. (e) Kaplan–Meier overall survival (OS) curves for 315 out of 379 OC samples in the TCGA dataset

Figure 3. Interaction among m6A RNA methylation regulators and functional annotation of OC in cluster ½ subgroups. (a) The m6A modification-related interactions among the 17 m6A RNA methylation regulators. (b) Spearman correlation analysis of the 17 m6A modification regulators. (c) GO analysis by GSEA of cluster 1 and cluster 2

Figure 4. Gene set enrichment show genes with higher expression in cluster 1 were enriched for KEGG of malignant tumors

Figure 5. Risk signature with 17 m6A RNA methylation regulators. (a) The process of building the signature containing 17 m6A RNA methylation regulators. The hazard ratios (HR), 95% confidence intervals (CI) calculated by univariate Cox regression. (b) LASSO regression analysis was used to calculate the coefficient of interferon gamma response genes. (c) Three genes were selected as active covariates to determine the prognostic value after 10-fold cross-validation for the LASSO model. (d-e) The risk scores for all patients in TCGA cohort are plotted in ascending order and marked as low risk (blue) or high risk (red), as divided by the threshold (vertical black line). (f) The distribution of risk score, survival status, and the expression of 3 genes of each patient in TCGA cohort by z-score, with red indicating higher expression and light blue indicating lower expression

Figure 6. Screening of Prognosis-related m6A RNA methylation regulators. (a) Kaplan–Meier overall survival (OS) curves for patients in the TCGA datasets assigned to the low- and high-risk groups. (b) ROC curve for 5-year survival prediction and clinical characteristics, including age, stage, grade, and risk score. (c) Univariate Cox regression analysis of the associated between clinicopathological factors (including risk score) and overall survival of patients. (d) Multivariate Cox regression analysis of the associated between clinicopathological factors (including risk score) and overall survival of patients. (e) WTAP expression levels in different age groups. (f) The expression levels of HNRNPA2B1 in the TUMOR FREE group and the TUMOR group. (g) RISK SCORE in different age groups

Figure 7. GSEA results and KEGG enrichment. (a) Top enriched KEGG pathways in the high risk group are represented by the curves above the x-axis in the graph. Top enriched KEGG pathways in the low risk group are represented by the curves below the x-axis in the graph (p-value < 0.05) in TCGA dataset. The names of enriched KEGG pathways are listed on the right side. (b) GSEA plots of KEGG Pathways in which the WTAP, KIAA1429 and HNRNPA2B1 were co-enriched

Figure 8. The nomogram to predict 3‐ or 5‐year OS in the entire set. (a) The nomogram for predicting proportion of patients with 3‐ or 5‐year OS. (b-c) The calibration plots for predicting patient 3‐ or 5‐year OS. Nomogram‐predicted probability of survival is plotted on the x‐axis; actual survival is plotted on the y‐axis. (d) DCA for assessment of the clinical utility of the nomogram. The x‐axis represents the percentage of threshold probability, and the y‐axis represents the net benefit. DCA: decision curve analysis; OS: overall survival

Two clusters of m6A RNA methylation regulators were associated with distinct OC clinical outcomes and clinicopathological features

With clustering stability increasing from k = 2 to 10, k = 2 seemed to be an adequate selection based on the expression similarity of m6A RNA methylation regulators (Figure 2a-c). Then, 379 OC samples were clustered into two subgroups in the TCGA dataset (cluster 1:216; cluser2:163). The two subgroups were named as cluster 1 and cluster 2 and clinicopathological features of the 2 two clusters were compared by k = 2 (Figure 2d). A significantly shorter OS was displayed in Cluster 1 than in Cluster 2 (Figure 2e).

Categories identified by consensus clustering are closely associated with the progression of OC

In order to understand the interactions among the 17 m6A RNA methylation regulators, the interaction (Figure 3a) and correlation (Figure 3b) among these regulators were analyzed. METTLE3 lies at the core of the network of m6A RNA methylation regulators. Its interactions and co-expressions with KIAA1429, METTLE14, WTAP, YTHDF1, YTHDF3, YTHDF2, and YTHDC1 were constructed and displayed in the String database. METTLE3 was also significantly correlated with YTHDC2 and YTHDF3. Three genes might co-work to regulate the progression of OC. The functional analysis of the two clusters was further performed. GSEA was functioned to show that the most relative GO terms were chemokine activity, macrophage chemotaxis, macrophage migration, monocyte chemotaxis, and tumor necrosis factor biosynthetic process in cluster 1 and cluster 2 (Figure 3c). The most relative pathways were allograft rejection, asthma, graft-versus-host disease, oxidative phosphorylation, and ribosome (Figure 4). The above findings suggest that the two categories identified by consensus clustering are closely associated with the progression of OC.

m6A RNA methylation regulators had prognostic significance

The prognostic ability of m6A RNA methylation regulators in OC was investigated. Univariate Cox regression analysis was conducted based on the expression levels of m6A RNA methylation regulators (Figure 5a). The results indicated that the three genes were significantly correlated with OS (P < 0.05). WTAP and KIAA1429 were two risky genes with HR > 1, while HNRNPA2B1 was a risky gene with HR < 1. To evaluate the ability of m6A RNA methylation regulators in predicting the clinical outcomes of OC, the LASSO Cox regression algorithm on three prognosis-associated genes was performed (Figure 5b-c), which were selected to construct the risk signature based on the minimum criteria. Coefficients obtained from LASSO algorithm were used to calculate the risk score: HNRNPA2B1*-0.01+ KIAA1429*0.085+ WTAP*0.03. The OC samples (n = 374) into low- and high-risk groups were separated based on the median risk score. The distribution of risk score, survival status, and the expression of three genes from each patient were also displayed (Figure 5d-f). Significant difference was observed in OS between the two groups (Figure 6a). ROC curves for 5-year survival were used to reveal the predictive performance of the three gene risk signatures. The 5-year AUC of the signature was 0.649, which was obviously higher than that of stage (AUC = 0.512), grade (AUC = 0.525) and age (AUC = 0.539) (Figure 6b). The results showed the three gene risk signatures had a stronger ability to predict OC survival than clinical factors.

Prognostic value and clinical utility of three m6A regulators

the TGGA dataset, univariate and multivariate regression models were constructed to identify whether the risk signature was an independent prognostic factor. Univariate analysis showed that tumor status and risk score were both correlated with OS (Figure 6c). Having absorbed three genes into the multivariate regression analysis, tumor status and risk score remained significantly correlated with OS (Figure 6d). Furthermore, the clinical features were associated with the three genes and Risk score (Table 1). It has also been found that the WTAP expression level is significantly different in different age groups (Figure 6e). The expression levels of HNRNPA2B1 in the TUMOR FREE group and the TUMOR group were also significantly different (figure 6f). Risk scores of patients in different age groups were also significantly different (Figure 6g). The stratification analysis was then performed based on grade, age, stage, and tumor status. Patients were stratified into Grade I/II and Grade III/IV subgroups and Stage I/II and Stage III/IV subgroups. As shown in Supporting Figure 9a, the prognosis of high-risk patients was significantly worse than that of low-risk patients in the Stage III/IV subgroup, which was consistent with the results of Grade II/IV subgroup (Supporting Figure 9b). However, there is no statistical significance in Stage I/II subgroup and Grade I/II subgroup. The prognostic ability of the three-gene signature combined with age and tumor status was also assessed. The patients were also stratified into different subgroups, including subgroups which is above 60 years and below 60 years. Interestingly, it has been revealed that high-risk patients in two subgroups were inclined to unfavorable OS (Supporting Figure 9c-e). Most of the immunity-related pathways were enriched in high-risk group, like T cell receptor signaling pathway, cytokine receptor interaction, and TOLL-like receptor signaling pathway. Most of the immunity-unrelated pathways were enriched in the low-risk group, like DNA replication and linoleic acid metabolism (Figure 7a). The three genes from the risk score model were co-enriched in cell adhesion molecule cams and chemokine signaling pathway (Figure 7b). The expression levels of three hub genes in 16 ovarian cancer clinical tissues and 8 normal ovarian tissues were also verified by the research (Figure 9). The results showed that the expression of HNRNPA2B1 was higher in normal ovarian tissues, and the expression of KIAA1429 was higher in ovarian cancer tissues, and both of them were consistent with our predicted trends. The expression level of WTAP was higher in ovarian cancer tissues, which was contrary to our prediction.

Table 1. Clinical significance of three prognosis-related genes

Gene	Age (≥60/<60)		Stage (I–II/III–IV)		Grade (1–2/3-4)		Tumor Status (with tumor/tumor free)
	T	P	T	P	T	P	T	P
HNRNPA2B1	−0.165	0.869	1.49	0.154	−0.571	0.570	2.844	0.005
KIAA1429	−0.615	0.539	1.909	0.072	−1.268	0.209	0.082	0.935
WTAP	3.641	3.264e-04	1.337	0.198	1.218	0.230	−0.146	0.884

Bold values indicate P < 0.05.

Note: t: t value of student’s t test; P: P‐value of student’s t test.

Figure 9. Expression level of WTAP, KIAA1429 and HNRNPA2B1 in clinical group

Association between three m6A regulators and immune infiltration

TCGA dataset was used to search the most significant tumor-infiltrating immune cells. Risk score was calculated to indicate the association between immune infiltration and three m6A regulators. Dendritic fraction, Macrophage fraction, and Neutrophil fraction have been discovered to be mostly enriched in high-risk group (Supporting Figure 10A-C).

A nomogram based on three m6A regulators

Encompassing age, stage, grade, tumor status, and risk score, a nomogram was constructed to predict the three-year or five-year OS of OC (Figure 8a). The calibration curve form Figure 8b-c suggests that the nomogram exhibited a performance as good as that of the Kaplan–Meier estimates. The C-index for this nomogram was 0.789 and became 0.773 after bootstrapping validation, showing its good discriminating ability. Meanwhile, DCA was created to estimate the clinical utility of the nomogram. The results of DCA showed that the nomogram containing three mRNAs’ signature had better prediction ability, with a threshold ranging from 2% to 83% (Figure 8d).

Genetic information of the seventeen genes

The genetic alteration harbored in the 17 genes was analyzed with cBioPortal software. The network was exhibited, which is constructed by METTL3, HNRNPA2B1, HNRNPC, FMR1, and their 50 most associated neighbor genes (only four out of the 17 genes had a joint node, while the remaining three genes had no junctions and were not shown) (Supporting Figure 11A). Supporting Figure 11B-C illustrates that the 17 genes were altered in 471 (79%) of the 594 patients (606 in total); YTHDF1, WTAP, and ZC3H13 showed the most diverse alterations, including amplification, missense mutation etc.

Discussion

Ovarian cancer is the most common gynecological malignancy. Most of the women have developed advanced stage when diagnosed [37]. Early diagnosis is critical to improve OC prognosis because the 5-year relative survival rate at the local stage is 93%. Specific biological diagnostic markers have been defined [38]. m6A modification has been implicated in mRNA turnover, localization, or translation [39–41]. As mainstream RNA regulators, m6A modulators have been widely proven to coordinate proteins related to coding functions in mRNA, and ultimately affecting tumors such as colorectal cancer, melanoma, etc. [42–47]. We hope to find m6a regulators related to ovarian cancer in batches. Therefore, it was reasonably speculated that m6A RNA methylation regulators are associated with ovarian cancer. It was also identified that two OC subgroups (Cluster 1 and Cluster 2) based on the expression of m6A RNA methylation regulators during the research. Two clusters were not only related to OC prognosis and clinicopathological features, but also related to some functional pathways, including graft-versus-host disease and oxidative phosphorylation. Coincidentally, these functional pathways have been found to regulate the development of OC. For example, Bay JO et al. found that an OC patient developed acute graft-versus-host disease and from this time her tumor diminished progressively [48]. Hänel M et al. also discovered a graft-versus-tumor effect in refractory ovarian cancer [49]. Pastò A et al. found oxidative phosphorylation in the stem cells from epithelial ovarian cancer patients [50]. Oxidative phosphorylation has also been validated as a therapeutic target for ovarian cancer [51].

METTL3 is the most widely studied writer. The experiments of Hua W et al. demonstrated that METTL3 promoted ovarian carcinoma growth and invasion through regulating AXL translation and epithelial to mesenchymal transition [3]. Cai X et al. demonstrated that HBXIP-elevated METTL3 expression and promoted the progression of breast cancer via inhibiting tumor suppressor let-7 g [52]. Vu LP et al. found that METTL3 curbed myeloid differentiation of normal hematopoietic and leukemia cells [53]. It has been identified that METTLE3 co-worked with others in OC development as the central gene in the network of m6A RNA methylation regulators. This conclusion has been proven in previous studies [54,55]. Studies have shown that N6-methyladenosine (m6A) is installed by the METTL3-METTL14-WTAP methyltransferase complex [56]. It has been identified that METTLE3 co-worked with others in OC development as the central gene in the network of m6A RNA methylation regulators. This conclusion has been proven in previous studies [57]. Wang Qiang and others found that the expression of METL3 promotes tumor angiogenesis and glycolysis in gastric cancer, and Mettl3 may be a cancer-promoting factor for gastric cancer [58]. Mettl3 has also been proven to promote the progression of cervical cancer, which is also a gynecological tumor [59].

Another writer associated with ovarian cancer is WTAP [60]. Consistently, our results showed that WTAP had a mutation rate of 16% and was highly expressed in patients aged over 60. The older the age, the greater the cumulative mutational load. The prognosis of bladder cancer and malignant glioma is also affected by WTAP [61,62]. Jing Wang et al. found that WTAP functions as an oncogenic factor that promotes the progression of ovarian cancer in which WTAP-HBS1L/FAM76A axis may be involved [63].

METTL14 is the main factor involved in aberrant m6A modification. Ma JZ et al. found that METTL14 as a writer suppressed the metastatic potential of hepatocellular carcinoma by modulating N6-methyladenosine-dependent primary microRNA processing [64]. In the study, METTL14 was lowly expressed in ovarian cancer samples and highly expressed in patients of Stage 1/2 subgroup and correlated with OC prognosis; thus, it is clear that the same M6A regulator also exerted different effects in different tumors. Mettl14 has also been shown to be related to breast cancer, colorectal cancer, and pancreatic cancer [65–67]

KIAA1429, METTL16, RBM15, and ZC3H13 were less studied in tumors. Qian JY et al. found that KIAA1429 acted as an oncogenic factor in breast cancer by regulating CDK1 in an N6-methyladenosine-independent manner [68]. KIAA1429 also has been proven to regulate cell proliferation by targeting c-Jun messenger RNA directly in gastric cancer [69] and participate in the migration and invasion of hepatocellular carcinoma [70]. ZC3H13 was found to suppress colorectal cancer proliferation and invasion [71]. ZC3H13 was also found to be Tumor Suppressor Genes in Breast Cancer [72]. We found that ZC3H13 showed a mutation rate of 18% and highly expressed in OC samples, and its expression level was negatively correlated with OC prognosis. It has been discovered that METTL16 was lowly expressed in OC tissues and was positively correlated with the prognosis via this study. METTL16 was lowly expressed in the samples collected from patients younger than 60. The same trend was shown in patients of stage III–IV subgroup and with tumor subgroup. These results support our hypothesis that METTL16 suppressed the development of OC. These findings indicate that increasing the level of m6A enrichment which was conducted by the writer can indeed alter the development of the tumor. However, a specific mechanism still needs to be tapped.

As complementary factors of writers, erasers also exert effects on a variety of tumors. Obesity is a high-risk factor for many tumors [73], and fat mass and obesity (FTO) are associated with obesity [74]. Akbari ME et al. found that FTO gene affected obesity and breast cancer through similar mechanisms [75]. FTO is associated with the occurrence and prognosis of gastric cancer [76]. Alkylation repair homolog protein 5 (ALKBH5) is also associated with pancreatic cancer [77], gastric cancer [78] and breast cancer [79]. Zhu H et al. found that ALKBH5 inhibited autophagy of epithelial ovarian cancer through regulating miR-7 and BCL-2 [80].

Three readers (YTHDF1, YTHDF2, and YTHDF3) were all highly expressed in OC samples and negatively related to prognosis. YTHDF1, YTHDF2, and HNRNPC have been intensely studied. YTH domain family 1 (YTHDF1) has a mutation rate of 27% and high expression rate in OC samples, which is associated with the poor prognosis of OS and DFS. YTHDF1 has been shown to be involved in the regulation of colorectal and pancreatic cancer [81,82]. Overexpression of YTHDF1 is associated with the poor prognosis of hepatocellular carcinoma [83]. In pancreatic cancer cells, YTHDF2 orchestrates epithelial–mesenchymal transition/proliferation dichotomy [82]. Li J et al. found that downregulation of N6-methyladenosine binding YTHDF2 protein mediated by miR-493-3p suppressed prostate cancer [84]. YTHDF2 also has a certain regulatory effect in lung cancer and gastric cancer [85,86]. HNRNPC can serve as a candidate biomarker for chemoresistance in gastric cancer [87]. Kleemann M et al. demonstrated that MiR-744-5p could induce cell death by directly targeting HNRNPC and NFIX in ovarian cancer [88]. The BRCA gene mutation is a feature of hereditary ovarian cancer. The samples were classified according to the presence of BRCA gene mutation [89], and it showed that the expression of FMR1 was related to BRCA gene mutation. Gleicher N et al. also found that BRCA/FMR1 had a correlation with ovarian cancer [90]. The mechanism of BRCA/FMR1 mutation causing ovarian cancer deserves further study.

A prognostic regression analysis found that HNRNPA2B1, KIAA1429, and WTAP have the strongest correlation with OC. Subsequently, OC patients were stratified into two subgroups with statistically different survival outcomes. In addition, univariate and multivariate Cox analyses identified the prognostic signature as an independent factor. In this study, due to the lack of an external validation cohort, to validate the prognostic performance of the 3-mRNA signature cannot be achieved. Thus, bootstrapping with 1,000 resamples was applied to internally validate the performance of three mRNA signature. The C-index for the internal validation was 0.773, indicating its good performance in clinical use. Moreover, a nomogram containing the 3-mRNA signature and other clinical features of ovarian cancer was built. The nomogram showed a moderate performance in predicting the survival of OC patients. Meanwhile, the results of DCA suggested that the nomogram showed better prediction ability, with a threshold ranging from 2% to 83%. Except WTAP, the role of both HNRNPA2B1 and KIAA1429 in ovarian cancer has not been thoroughly studied and can be used to direct further research. GSEA showed that the samples from the high-risk group were mainly enriched in immune-related pathways. Interestingly, the study showed that Dendritic fraction, Macrophage fraction, and Neutrophil fraction were related to the three m6A regulators. Surprisingly, there are experiments demonstrating that dendritic cell (DC) immunotherapy can induce anti-tumor T cell immunity [91]. Macrophage can regulate the progress of OC through multiple mechanisms like CD47 [92] and NF-κB activation [93]. Meta-analysis by Chen S et al. showed neutrophil-to-lymphocyte ratio is a potential prognostic biomarker in patients with ovarian cancer [94]. It can be seen that distorted immune microenvironment may induce tumors to some extent.

This research has the following shortcomings: 1. The model is not validated with external data. 2. Lack of verification of in vitro and in vivo experiments. Prospective clinical trials are necessary in the future to reconfirm the findings.

Conclusion

This study analyzed the association between m6A regulators and clinical features of OC. The three selected m6A RNA methylation regulators (HNRNPA2B1, KIAA1429, and WTAP) showed high prognostic value for OC and were also enriched in the biological processes and signaling pathways that drive the malignant progression of OC. High-risk patients who had a dendritic fraction, macrophage fraction, and neutrophil fraction were also found in this study. Future clinical and experimental research is warranted to further verify the results of this study. In brief, this study provides novel markers for evaluating OC prognosis.

Abbreviations

Ovarian cancer (OC) Search Tool for the Retrieval of Interacting Genes Database (STRING)	=
Receiver operating characteristic (ROC)	=
Gene set enrichment analysis (GSEA)	=
Decision curve analysis (DCA)	=
Kyoto Encyclopedia of Genes and Genomes (KEGG)	=
Area under the curve (AUC)	=
Hazard ratio (HR)	=

Declaration

Ethics approval and consent to participate

The study was approved by the Ethics Committee of the Jiangsu Province Hospital of Chinese Medicine and samples were obtained with written informed consent from all patients, and all informed consent were conducted in accordance with the Declaration of Helsinki.

Author contribution

Hongyu Guo: Conceptualization, Data Curation

JinHui Liu: Data Curation

Dandan Hong: Writing – Original draft preparation

Jing Ji: Visualization

Qin Zhang: Supervision

Qun Guan: Supervision

Qingling Ren: Visualization

Cheng Zhang: Revision

Highlights

1.Three m6A regulatory genes were identified to predict OS of ovarian cancer patients.

2.The prognostic signature of three m6A regulatory genes are an independent prognostic factor.

3.The immune infiltrations are enriched in different risk groups.

Consent for publication

Written informed consent for publication was obtained from all participants.

Acknowledgements

Not applicable.

Supplementary material

Supplemental data for this article can be accessed here.

Disclosure statement

The authors declare that they have no competing interests.

Data availability statement

The data and materials are available in The Cancer Genome Atlas (https://cancergenome.nih.gov/https://cancergenome.nih.gov/).

Additional information

Funding

The authors have no funding to report.