The methylation modification of m6A regulators contributes to the prognosis of ovarian cancer.

Background: Ovarian cancer (OV) is the leading cause of death in gynecological cancer. The dysregulation of N6-methyladenosine (m6A) modification is commonly found in cancers. However, there is a lack of research into m6A RNA methylation regulators in OV.
Methods: The RNA-Seq of 379 OV tissues and 88 healthy ovarian tissues was downloaded from The Cancer Genome Atlas (TCGA) and the Genotype-Tissue Expression (GTEx) databases, respectively. A Gene Ontology (GO) functional analysis was performed to verify the function of m6A RNA methylation regulators. Kaplan-Meier (K-M) curves and the log-rank (Mantel-Cox) test were used for the survival analysis. A Cox regression analysis was used to identify the genes related to overall survival (OS) and build the prediction model.
Results: m6A RNA methylation regulators were dysregulated in OV tissues compared with normal tissues (P<0.05), and patients with a high expression of KIAA1429 and YTHDC2 had a poor prognosis (P<0.05). A prognostic model was constructed based on the m6A RNA methylation regulators. Based on the risk signature, the patients were classified into high- and low-risk groups. The low-risk group's OS rate was significantly better than that of the high-risk group. The validity and accuracy of the prognostic model were verified by using TCGA and Gene Expression Omnibus (GEO) datasets, and the risk score from the prognostic model acted as an independent prognostic indicator in predicting the survival of OV patients. Conclusions: m6A RNA methylation regulators were dysregulated in OV tissues. More importantly, the prognostic model comprising the five selected m6A RNA methylation regulators could be a valuable tool for predicting the prognosis of OV patients. 2022 Annals of Translational Medicine. All rights reserved.

Introduction

Ovarian cancer (OV) is the leading cause of death in gynecological cancer since patients with early-stage OV are asymptomatic, and 75% will have reached an advanced stage (stage III or IV) by the time they are diagnosed (1). Although age-standardized rates are stable or falling in most high-income countries, they are rising in many low- and middle-income countries (2). According to the American Cancer Society, it was estimated that there were approximately 295,414 new cases and 184,799 deaths attributed to OV worldwide in 2018 (3).

Post-transcriptional modification of RNA is an important means of post-transcriptional processing of RNA, and it is an essential mechanism for post-transcriptional regulation of RNA (4).

N6-methyladenosine (m6A) modification is the most prevalent internal methylation in mRNA, and it is deposited by m6A methyltransferases, removed by m6A demethylases (5). m6A modification play important role in gene expression regulation of basic cellular processes and multiple physiological functions, and among of them, the most important role is modulating the eukaryotic transcriptome to influence mRNA splicing, export, localization, translation, and stability (6,7). Aberrant m6A modification can cause cancer, but it is also have the potential to yield new therapies (8).

New research has found that methylated nucleosides, especially m6A modification, which determine the fate of mRNAs by the alternatives of cell differentiation and embryonic development, exist in the mRNA of Novikoff hepatoma cells (9,10). Moreover, evidence is emerging that m6A modification and its regulatory proteins also play critical roles in various cancers including urological tumors, clear cell renal cell carcinoma, glioma, gastric cancer, etc. In urological tumors, many m6A modification-related proteins are dysregulated and promote cancer development (11). In clear cell renal cell carcinoma, alterations of m6A RNA methylation regulators were found to be related to the clinicopathological characteristics and survival of patients (12), and m6A RNA methylation regulators are also critical regulators of glioma progression (13). Further, m6A modification was shown to influence the malignant phenotypes of gastric cancer by regulating the Wnt/PI3K-Akt signaling pathway (14).

In ovarian cancer, m6A regulators are promising as potential molecular therapeutic targets for ovarian cancer (15). For example, YTHDF1 can promote tumorigenesis of ovarian cancer cells by regulating the translation of eIF3C via an m6A-dependent way, thus affects the global protein translation in ovarian cancer (16). Moreover, Han et al. showed that WTAP was significantly increased in ovarian tissues, and the highly expression is associated with the regulation of the cell cycle and the target of MYC (17).

Although m6A modification has been reported to play an important role in tumors, the relationship between m6A RNA methylation modulator signatures and OV prognosis remains to be verified. To confirm the precise m6A regulation pattern in OV, we systematically analyzed the expression of m6A RNA methylation regulators and constructed a prognostic model based on The Cancer Genome Atlas (TCGA) and the Genotype-Tissue Expression (GTEx) databases.

We present the following article in accordance with the TRIPOD reporting checklist (available at https://atm.amegroups.com/article/view/10.21037/atm-21-6462/rc).

Methods

Patient information

The clinical data of OV patients were downloaded from TCGA database (https://portal.gdc.cancer.gov/). There are 379 OV patients with RNA-Seq profiles and clinical data in TCGA database. The RNA-Seq profiles of 379 OV tissues and 88 healthy ovarian tissues were extracted from the University of California Santa Cruz (UCSC) database (http://xena.ucsc.edu/). Gene Expression Omnibus (GEO) datasets, including GSE66957 and GSE63885, were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The data from TCGA and GEO were downloaded on September 27, 2021. The study setting, eligibility criteria for participants, and treatments received were all obtained from TCGA. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

m6A RNA methylation modulator selection

A total of 13 genes are recognized as vital m6A methylation modulators, including METTL3, METTL14, WTAP, KIAA1429, RBM15, ZC3H13, YTHDC1, YTHDC2, YTHDF1, YTHDF2, HNRNPC, ALKBH5, and FTO (18).

RNA isolation and quantitative real‑time PCR (qRT‑PCR)

For tissue RNA isolation, 1 mL AG RNAex Pro Reagent (Accurate Biotechnology Co.) was added to 50 mg of tissue and total RNA samples were extracted according to the manufacturer’s instructions. Purified RNA was quantified resort to NanoVue (GE Healthcare Life Sciences). cDNAs were synthesized from total RNAs by using RT reagent Kit (Takara Co., LTD, Japan) and ReverTra Ace qPCR RT Kit (Toyobo Co., LTD, Japan) (19).

qRT-PCR of GAPDH, KIAA1429, HNRNPC, YTHDC2, ZC3H13 and WTAP was performed with the SYBR qPCR Mix (Toyobo Co., LTD, Japan). 10 µL reaction system was adopted according to the manufacturer’s instructions and amplified for 40 cycles. The expression levels were normalized by GAPDH. Relative expression was calculated using the method of 2−∆∆Ct (20). Primer names and primer sequences are listed in the following tables ().

Table 1

The forward and reverse primer sequence of the genes

Primer name	Primer sequence
GAPDH	Forward: CAACGTGTCAGTGGTGGACCTG
GAPDH	Reverse: GTGTCGCTGTTGAAGTCAGAGGAG
KIAA1429	Forward: GGAATGGACACGTTTATTCGAG
KIAA1429	Reverse: GATAGAGCACAGGAGCATATGT
ZC3H13	Forward: GATCAGTTAAAGCGTGGAGAAC
ZC3H13	Reverse: CTCTCTGTCGTGTTCATATCGA
WTAP	Forward: CTGACAAACGGACCAAGTAATG
WTAP	Reverse: AAAGTCATCTTCGGTTGTGTTG
YTHDC2	Forward: GAGAATTGGGCTGTCGTTAAAG
YTHDC2	Reverse: TGAAGCAGGATGAAATCGTACT
HNRNPC	Forward: CGTGTACCTCCTCCTCCTCCTATTG
HNRNPC	Reverse: CCCCGCTGTCCACTCTTAGAATTG

Bioinformatics analysis

A Gene Ontology (GO) functional enrichment analysis was performed to verify the function of the 13 m6A methylation modulators. Next, heatmap and violin plots were used to visualize the expression of the 13 regulators in 379 tumor tissues and 88 normal tissues. Spearman’s analysis was conducted to explore the correlations among these regulator genes, and a Kaplan-Meier (K-M) analysis was adopted to confirm the prognostic regulators. We performed a univariate Cox analysis to explore the prognostic role of m6A methylation regulators in OV patients. Then, a lasso regression analysis was used to construct the prognostic model. Five genes were identified as powerful prognostic factors. The risk score of each patient was calculated using the following formula: Risk score = 0.15987 × expression of KIAA1429 − 0.11282 × expression of HNRNPC + 0.069938 × expression of YTHDC2 + 0.109001 × expression of ZC3H13 + 0.123254 × expression of WTAP.

Statistical analysis

Unpaired t-tests were used to confirm the differential expression of the m6A RNA methylation modulators. The median risk score was set as the cut-off value to divide patients into high- or low-risk groups. The log-rank (Mantel-Cox) test was used for the survival analysis. Differences were considered statistically significant when the P value was <0.05.

Results

Confirmation of the function of m6A RNA methylation regulators

The level of m6A methylation is determined by the dynamic interplay between methyl transferases (‘writers’), binding proteins (‘readers’), and demethylases (‘erasers’). Writers of m6A include methyl transferase Like 3 (METTL3) (21), METTL14 (22,23), Wilms Tumor 1 Associated Protein (WTAP) (24), KIAA1429 (25), RNA Binding Motif Protein 15 (RBM15) (22), and ZC3H13 (26). Readers of m6A contain YTHDC1 (13,27), YTHDC2 (28), YTHDF1 (13), YTHDF2 (29), and HNRNPC (30). Erasers of m6A include ALKBH5 (31) and FTO (32).

We conducted a GO functional analysis to verify the function of the m6A RNA methylation regulators. Our results showed that the primary biological processes involved in these genes are RNA modification, methylation, metabolic process, and splicing. Moreover, these genes are involved in the formation of nuclear speck and methyltransferase complex (, ), and their molecular functions are catalytic activity and mRNA binding (, ).

Figure 1

GO functional analysis of m6A RNA methylation regulators. (A) Bar plot; (B) bubble plot; (C,D) circle plots. GO, Gene Ontology; m6A, N6-methyladenosine. BP, biological process; CC, cellular component; MF, molecular function.

Table 2

GO functional analysis of m6A RNA methylation regulators

Ontology	ID	Description	P value	Count	Gene ID
BP	GO:1903311	Regulation of mRNA metabolic process	6.82E-09	7	HNRNPC, RBM15, METTL14, METTL3, WTAP, YTHDC1, YTHDF2
BP	GO:0009451	RNA modification	5.22E-10	7	RBM15, METTL14, METTL3, ZC3H13, WTAP, FTO, ALKBH5
BP	GO:0008380	RNA splicing	1.72E-07	7	HNRNPC, RBM15, METTL14, METTL3, ZC3H13, WTAP, YTHDC1
BP	GO:0001510	RNA methylation	4.46E-08	5	RBM15, METTL14, METTL3, ZC3H13, WTAP
BP	GO:0080009	mRNA methylation	6.82E-09	4	METTL14, METTL3, ZC3H13, WTAP
BP	GO:0016556	mRNA modification	3.72E-08	4	METTL14, METTL3, ZC3H13, WTAP
CC	GO:0016607	Nuclear speck	8.66E-09	7	RBM15, METTL3, ZC3H13, WTAP, YTHDC1, FTO, ALKBH5
CC	GO:0034708	Methyltransferase complex	3.02E-08	5	RBM15, METTL14, METTL3, ZC3H13, WTAP
MF	GO:0140098	Catalytic activity, acting on RNA	8.92E-05	5	METTL14, YTHDC2, METTL3, FTO, ALKBH5
MF	GO:0003729	mRNA binding	0.002685	4	HNRNPC, RBM15, METTL14, METTL3

GO, Gene Ontology; BP, biological process; CC, cellular component; MF, molecular function.

The most significant biological processes in which m6A RNA methylation regulators are involved are RNA modification, regulation of mRNA metabolic process, and RNA splicing (). Among the genes related to m6A methylation, METTL3, METTL14, WTAP, and ZC3H13 participate in multiple RNA biological processes. However, YTHDF2, HNRNPC, YTHDC1, and FTO are only involved in RNA modification and the mRNA metabolic process ().

Dysregulation of m6A RNA methylation regulators in OV tissues

There are 379 OV tissues with RNA-Seq in TCGA database and 88 normal ovarian tissues with RNA-Seq in the GTEx database. The expression of m6A RNA methylation regulators was detected in OV tissues and healthy ovarian tissues, and all were dysregulated in OV tissues compared with normal ovarian tissues (, P<0.05). ZC3H13, ALKBH5, RBM15, YTHDF1, and YTHDF2 were upregulated in OV tissues, and the remainder were downregulated (, P<0.05).

Figure 2

The m6A RNA methylation regulators dysregulated in OV tissues. (A,B) The expression of m6A RNA methylation regulators in OV tissues versus healthy ovarian tissues; (C) the correlation between the m6A RNA methylation regulators; (D,E) the K-M plot curves of KIAA1429 and YTHDC2. ***P<0.0001. m6A, N6-methyladenosine; OV, ovarian cancer.

A co-expression analysis was performed to confirm if co-expression existed in the m6A RNA methylation regulators in OV tissues. The results showed the existence of positive and negative co-expressions. METTL3 (‘writers’) and HNRNPC proteins (‘readers’) have the most positive co-expression with an R value of 0.84, and YTHDF1 (‘readers’) and FTO (‘erasers’) have the strongest negative co-expression with an R value of −0.82 ().

Further, K-M analysis was applied to confirm the relationship between these genes’ expression and patients’ survival. Patients were divided into two groups according to the median value of the gene’s expression, and K-M curves were drawn to compare the OS time of the two patient groups. We found that a high expression of KIAA1429 and YTHDC2 predicted a poor prognosis in OV patients (P<0.05, ).

Construction of the prognostic model based on m6A RNA methylation regulators

A univariate Cox survival analysis was conducted to confirm the prognostic m6A RNA methylation regulators. We found five genes (KIAA1429, YTHDC2, ZC3H13, WTAP, and FTO) were highly expressed and were identified as high-risk factors for the poor prognosis of OV patients ().

Figure 3

Construction of the prognostic model based on m6A RNA methylation regulators. (A) The results and forest map of the univariate Cox survival analysis; (B,C) lasso regression analysis is used to screen the m6A RNA methylation regulators and confirm their regression coefficients in the prognostic model; (D) OV patients are divided into two groups according to the median value of the risk score; (E) the expression of the five m6A RNA methylation regulators in different clinicopathological parameters. **P<0.001. m6A, N6-methyladenosine; OV, ovarian cancer.

We also used a lasso regression algorithm to select genes and construct the prognostic model. According to the lasso regression analysis results, five genes (KIAA1429, YTHDC2, ZC3H13, WTAP, and FTO) were selected to construct the prognostic model for predicting the prognosis of OV patients (, ).

Table 3

The detailed information of the m6A RNA methylation regulators in the prognostic model

Gene	HR	HR.95L	HR.95H	P value	Coef
KIAA1429	1.386533	0.986996	1.947802	0.059499	0.15987
HNRNPC	0.762738	0.495315	1.174543	0.218851	−0.11282
YTHDC2	1.281947	0.871357	1.886009	0.207342	0.069938
ZC3H13	1.234432	0.961737	1.584448	0.098201	0.109001
WTAP	1.228301	0.923129	1.634359	0.158219	0.123254

m6A, N6-methyladenosine; HR, hazard rate.

A prognostic risk score formula was constructed according to the prognostic model. It was based on the linear combination of each gene’s expression multiplied by the coefficients from the results of the regression analysis as follows: Risk score = 0.15987 × expression of KIAA1429 − 0.11282 × expression of HNRNPC + 0.069938 × expression of YTHDC2 + 0.109001 × expression of ZC3H13 + 0.123254 × expression of WTAP. According to the risk score formula, each OV patient had a specific risk score value, and the patients were divided into two groups according to the median value of the risk score ().

Further, we compared the expression of the five m6A RNA methylation regulators in different clinicopathological parameters, including the two risk score groups. We found a significant difference between the high- and low-risk groups in OV patients’ survival status (P<0.05, ).

The risk score was closely related to the prognosis of OV patients

A K-M analysis was conducted to validate the risk score derived from the m6A RNA methylation regulators. Patients with a high-risk score had a poorer prognosis than patients with a low-risk score (P<0.05, ). Moreover, a three-dimensional heat map of survival time, risk score, and vital status showed more deaths among patients with a high-risk score ().

Figure 4

Validation of the risk score derived from the five m6A RNA methylation regulators. (A) Kaplan-Meier plot of patients with low- and high-risk scores; (B) a three-dimensional heat map of survival time, risk score, and vital status; (C-G) Kaplan-Meier plots of patients with different clinicopathological parameters, including age, grade (G3), tumor size, and White race; (H) the ROC curve of the risk score derived from the five m6A RNA methylation regulators; (I,J) the expression of m6A RNA methylation regulators in different age and risk score groups. m6A, N6-methyladenosine; OS, overall survival; ROC, receiver operating characteristic.

Next, we analyzed the effect of the risk score on the prognosis of patients with different clinicopathological parameters. As shown in the K-M plot curves, the risk score was a stable prognostic marker for OV patients of different ages (<65 or ≥65 years, ). Similarly, when we used the risk score to predict patients’ outcomes in other clinicopathological variables, we found that patients with a high-risk score who had grade 3 (G3) tumors, with tumor, and were White had a shorter survival time (). ROC curve of the risk model derived from the five m6A RNA methylation regulators showed that the risk model has good accuracy and specificity. (). Further, we found that the m6A RNA methylation regulators were dysregulated in different stratified clinicopathological parameters. Results showed that WTAP was upregulated in the higher age group (P<0.05, ). Additionally, all five m6A RNA methylation regulators showed significant differences between the low- and high-risk groups (P<0.05, ). HNRNPC was downregulated in the high-risk group, whereas the other four genes were upregulated in the high-risk group (P<0.05, ).

GEO datasets were used to verify the results from TCGA dataset and for further data mining. In dataset GSE66957, all five m6A RNA methylation regulators were downregulated in tumor tissues except ZC3H13, and the results are consistent with the results from TCGA dataset (). Moreover, the expression of m6A RNA methylation regulators were verified in clinical ovarian cancer tissues versus healthy ovarian tissues ().

Figure 5

The GEO datasets used to verify the results from TCGA dataset and for further data mining. (A) Verification of the expression of m6A RNA methylation regulators in GEO dataset GSE66957; (B) verification of the expression of m6A RNA methylation regulators in clinical ovarian cancer tissues versus healthy ovarian tissues; (C) four m6A RNA methylation regulators from our prognostic model show different expressions in different drug sensitivity groups; (D,E) YTHDC2 is upregulated in patients with larger tumors and those in an advanced stage; (F,G) univariate and multivariate Cox survival analyses reveal the prognostic factors for OV patients. GEO, Gene Expression Omnibus; TGCA, The Cancer Genome Atlas; m6A, N6-methyladenosine; OV, ovarian cancer.

Interestingly, in the stratified analysis in dataset GSE63885, we found that four of the genes in our prognostic model were expressed differently in the stratified clinicopathological parameters. In the stratification analysis of chemotherapy resistance, HNRNPC, KIAA1429, YTHDC2, and WTAP were upregulated in the moderately sensitive and resistant groups (). In the stratification analysis of tumor size and stage, YTHDC2 was more highly expressed in larger tumor sizes and stage 4, implying that YTHDC2 may be a key m6A methylation gene in OV patients although this will require further verification ().

Prognostic factors for OV patients

Univariate and multivariate Cox survival analyses were conducted to investigate the prognostic factors for OV patients. The univariate Cox survival analysis revealed that age, cancer status, and risk score were high-risk factors for the poor prognosis of OV patients (P<0.05, , ). Moreover, they were all independent risk factors affecting the prognosis of OV patients (, ).

Table 4

Univariate and multivariate analyses predicting the prognostic factors for ovarian cancer patients

Clinical feature	Univariate analysis					Multivariate analysis
	HR	HR.95L	HR.95H	P value		HR	HR.95L	HR.95H	P value
Age	1.352072	1.007566	1.81437	0.04441		1.353225	1.005279	1.821601	0.046077
Grade	1.396651	0.884412	2.205573	0.15184		1.302327	0.821149	2.065467	0.261623
Race	0.760836	0.447244	1.294305	0.313301		0.798667	0.468054	1.362808	0.409609
Cancer status	9.795922	4.335205	22.13507	4.10E-08		9.289214	4.107879	21.00585	8.61E-08
Risk score	3.143275	1.516093	6.516869	0.00208		2.3302	1.116067	4.865149	0.024301

HR, hazard rate.

Discussion

OV is one of the most lethal tumors for women, with a 5-year survival rate below 45% (2). It is difficult to diagnose OV in the early stages because the patient is generally asymptomatic. Therefore, more than half of patients are in an advanced stage (stage 3 or 4) by the time of diagnosis (1). Moreover, the treatments for OV are limited, and the death rates are higher than the incidence rates because of the resistance to radiotherapy and chemotherapy, especially in advanced stages (33). Thus, new, effective therapies urgently need to be identified.

m6A, as the most abundant post-transcriptional modification of RNA, is reportedly associated with different metabolic processes in various RNAs (5,10). Recently, several studies have revealed that m6A is involved in the development of many cancers (34,35). In this study, we analyzed the expression of the genes involved in modifying m6A in OV and found that all are dysregulated in OV tissues compared with normal ovarian tissues (). Moreover, positive and negative co-expression were found among the genes related to m6A methylation modification ().

KIAA1429, one of the writers of m6A, has been reported to facilitate the migration and invasion of hepatocellular carcinoma by inhibiting ID2 via upregulation of the m6A modification of ID2 mRNA (36). Additionally, YTHDC2, the reader of m6A, is dysregulated in many kinds of cancers and is implicated in the metastasis of colon tumors (37), the growth of tumor cells (38), and the susceptibility of pancreatic cancer (39). However, there is a lack of research on KIAA1429 and YTHDC2 in OV. In our research, we identified for the first time that a high expression of KIAA1429 and YTHDC2 is associated with the poor prognosis of OV patients ().

Multiple dysregulated gene signatures from patients’ tumor tissues can provide better accuracy in predicting cancer prognosis than one specific gene and enable more effective and individualized treatment plans (40,41). In this study, we found that all of the genes related to m6A methylation were dysregulated in OV patients. Therefore, we were interested in knowing whether multigene signatures of m6A methylation could provide a more accurate prognosis prediction for OV patients. A Cox survival analysis and lasso regression analysis were used to construct the prognostic model for predicting OV patients’ prognosis (). Finally, a risk score formula was constructed, where the risk score = 0.15987 × expression of KIAA1429 − 0.11282 × expression of HNRNPC + 0.069938 × expression of YTHDC2 + 0.109001 × expression of ZC3H13 + 0.123254 × expression of WTAP.

Moreover, K-M survival curves and a three-dimensional plot including patients’ risk score, survival time, and survival status were used to verify the accuracy and specificity of the prognostic model we constructed (). There was a greater number of deaths in the high-risk score group, and multivariate Cox survival analysis revealed that the risk score acted as an independent risk factor in the prognosis of OV patients (, ). The analysis from the GEO datasets was consistent with the results from TCGA dataset (). Of the identified genes, YTHDC2 was dysregulated in OV tissues and demonstrated a higher expression in larger tumor sizes and stage 4. This implies that YTHDC2 may be a key m6A methylation gene in OV, although this finding will require further verification (). The accuracy and specificity of the model constructed in this study will be improved by a larger sample of normal ovarian tissues, and the specific molecular mechanism of m6A needs to be further studied.

In the future, m6A modification could be used as diagnostic or prognostic targets for ovarian cancer. The upstream regulators of m6A modification, or their downstream targets, will provide new ideas and targets for the treatment of ovarian cancer. In general, the prospect of m6A modification research is helpful for the diagnosis, treatment and prognosis of ovarian cancer.

Conclusions

Collectively, this is the first study to combine TCGA and GTEx databases to explore the effect of m6A methylation genes on the prognosis of OV patients, which remedies the shortcoming of the small sample size of normal ovarian tissues in TCGA database. Moreover, we report for the first time that YTHDC2 (readers of m6A) and KIAA1429 (writers of m6A) are associated with the prognosis of OV patients. The prognostic model we have constructed could be a valuable clinical tool to indicate the prognosis for OV patients, and further research on the genes in our prognostic model will provide a future direction for the theoretical research and clinical treatment of OV.

The article’s supplementary files as