Zhijing Na
Ling Fan2*- 1Center of Reproductive Medicine, Shengjing Hospital of China Medical University, Shenyang, China
- 2Nursing Department, Shengjing Hospital of China Medical University, Shenyang, China
N6-Methyladenosine (m6A) is one of the most prominent modification regulating RNA processing and metabolism. Increasing studies have illuminated the vital role of m6A methylation in carcinogenesis. However, little is known about the interaction between m6A-related genes and survival of ovarian cancer (OC) patients. The purpose of this study was to obtain more reliable m6A-related genes that could be used as prognostic markers of OC using bioinformatics analysis performed on the RNA-seq data of OC. Gene expression datasets of all m6A-related genes as well as corresponding clinical data were obtained from the International Cancer Genome Consortium (ICGC) and The Cancer Genome Atlas (TCGA) databases. We detected differential expressed m6A-related candidate genes as well as their relationship and interaction. m6A RNA methylation regulator ALKBH5 and 35 m6A-related genes are dysregulated in OC. A gene set that could be used as a potential independent prognostic risk feature was further screened including NEBL, PDGFRA, WDR91, and ZBTB4. The results of mRNA expression analysis by PCR were consistent with those of bioinformatics analysis. We applied consensus clustering analysis on the expression of the four prognostic genes and obtained four OC subgroups TM1-TM4. There were significant differences in age, stage and grade among the subgroups, and the overall survival (OS) as well as Disease-free survival (DFS) of TM2 group were shorter than those of the other three groups. Further GO and KEGG enrichment analysis indicated that these differential genes were closely related to biological processes and key signaling pathways involved in OC. In summary, our study has indicated that m6A-related genes are key factors in the progression of OC and have potential effects on the prognostic stratification of OC and the development of treatment strategies.
Introduction
Ovarian cancer (OC) ranks the seventh most common cancer worldwide, with a total incidence of 239,000 each year (Bray et al., 2018). It is the leading cause of gynecologic cancer-related deaths among women, causing 152,000 deaths yearly (Torre et al., 2018). Due to the missing early symptoms and the absence of effective early detection strategies, approximately 70% of OC patients were diagnosis at advanced stage presenting with metastases (Oldak et al., 2019). Although the treatment for OC has been greatly improved in recent decades, the recurrence is frequent and the 5-year overall survival rate remains poor (Coleman et al., 2019). Hence, it’s urgent to explore specific early diagnostic and prognostic biomarkers and improve the unfavorable prognosis of patients suffering from OC.
N6-Methyladenosine (m6A) modification, the most prominent chemical modifications in eukaryotic mRNA and noncoding RNA (lncRNA), is a dynamic reversible process that regulates RNA processing and metabolism (Deng et al., 2018; Yu J. et al., 2018). m6A modification on RNAs was labeled by methyltransferases (writers), preferentially recognized and transmitted by binding proteins (readers), and erased by demethylases (erasers) (Zhao et al., 2020; Jiang et al., 2021). METTL3, METTL14, WTAP, KIAA1429, ZC3H13, and RBM15 are writers that catalyze the methylation of m6A into RNAs (Wang T. Y. et al., 2020). Readers are capable of selectively recognizing m6A-modified RNAs to mediate the translation and degradation of RNAs, including HNRNPC, YTHDC1, YTHDC2, YTHDF1, and YTHDF2 (Shi et al., 2019). Finally, FTO and ALKBH5, considered as erasers, can remove the methyl group from target RNAs to achieve the dynamics and reversibility of the m6A modification process (Livneh et al., 2020).
Currently, increasing studies have demonstrated the roles of m6A methylation in various carcinogenesis processes including cell self-renewal (Dai et al., 2018), differentiation (Sun et al., 2019), cell proliferation (Liu et al., 2020), migration (Cheng et al., 2019; Liu et al., 2020), invasion (Cheng et al., 2019), autophagy (Wang et al., 2019), apoptosis (Ianniello and Fatica, 2018), and metastasis (Yue et al., 2019). For instance, overexpression of METTL3 in AML significantly inhibits cancer cell differentiation and apoptosis via activating PI3K/AKT pathway (Ianniello and Fatica, 2018). ALKBH5 erases m6A modification of tumor suppressor gene FOXM1 to promote cancer cell maturation and tumorigenicity in glioblastoma (Zhang et al., 2017). The aberrant m6A modifications were reported to contribute to multiple cancers, including lung cancer (Li et al., 2019), breast cancer (Niu et al., 2019), glioblastoma (Zhang et al., 2017), acute myeloid leukemia (AML) (Ianniello and Fatica, 2018), liver cancer (Cheng et al., 2019), gastrointestinal cancer (Ni et al., 2019; Sun et al., 2019), endometrial cancer (Liu et al., 2018), et al. Thus, m6A has great potential as promising markers in the diagnosis, prognosis and personalized targeted therapies of cancers. Since the effects of m6A-related genes in OC have not been mentioned yet, we conducted our study to shed light on the expression pattern, prognostic value and potential mechanisms of m6A-related genes in OC patients.
Materials and Methods
The Expression Pattern of m6A-Related Genes in OC Patients
Data Source
We downloaded the RNA-seq transcriptome data and corresponding clinical information of 308 OC samples from the TCGA database from the UCSC’s xena database1. At the same time, RNA-seq data from 200 OC samples from the ICGC database independent of the samples used in the TCGA dataset were also obtained. The original RNA-seq data was standardized data, and it was uploaded online (DOI: 10.6084/m9.figshare.14399900), and the all TCGA codes of the patients used in this study was displayed in Supplementary Table 1. The corresponding clinical information of the validated ovarian cancer samples were downloaded, as shown in Table 1.
Table 1. The clinical information of the validated ovarian cancer samples.
Identification of Functionally m6A-Modified Genes
We first collected some m6A RNA methylation related genes from the known literature, which could be divided into three types according to the role they played in the methylation process: methyltransferases, binding proteins (readers) and erasers. A total of 21 methylation related genes were obtained, they are: METTL3, WTAP, ZC3H13, RBM15, METTL14, YTHDC1, YTHDC2, YTHDF2, YTHDF3, HNRNPA2B1, HNRNPC, HNRNPG, FTO, ALKBH5, IGF2BP2, IGF2, IGF2, ZNF217, elF3H, elF3J (Yang et al., 2018; Lence et al., 2019). Then, in the existing m6A database m6Avar2, there were 296 m6A-related genes related to OC disease. The two data sets of m6A-related genes were integrated together, after removing duplicates and genes that had no expression value in the sample or expression values that were less than 80% of all samples. The resulting m6A-related gene set contained a total of 267 m6A-related genes, including 18 m6A RNA methylation regulators and 249 m6A-related genes.
Expression Data
In order to determine the m6A RNA methylation regulatory factors that were differentially expressed according to different stages of OC, we compared the expression values of the TCGA data using the one-way ANOVA (p ≤ 0.05) with R language, version 3.6.0. UQ-FPKM (Upper-Quartile Fragments Per Kilobase Million reads) normalization allows for cross-sample comparison, thus we conducted a one-way ANOVA of the previously normalized absolute expression values in this study. The analysis workflow of our study is shown in Figure 1. First, we preprocessed the TCGA data and selected pathological features of stage (Stage I, Stage II, Stage III, Stage IV) from RNA-seq expression data set of m6A RNA methylation regulators and m6A-related genes related to OC, a data matrix containing 267 genes and 306 samples was obtained.
Figure 1. Analysis workflow of this study.
The Prognostic Value of m6A-Related Genes in OC Patients
Consensus Clustering
Consensus clustering analysis was conducted to classify OC samples, utilizing the m6A-related candidate gene set related to the prognosis of OC. Patients were then divided into subgroups according to age, stage, grade, new neoplasm event type, BRCA1 mutation status, BRCA2 mutation status, 19q13.2 CNV mutation and 19q13.42 mutation. Then, we conducted comparative analysis of subgroup survival differences.
Survival Analysis
Univariate Cox regression analysis was performed for each functional m6A modifier using R “Survival” package to screen gene sets related to prognosis, and genes with p ≤ 0.05 were candidate genes with potential independent prognostic efficacy. We applied LASSO regression to further screen for potential prognostic risk characteristics. Multivariate Cox regression analysis of risk characteristics was conducted to further screen for genes that were significantly related to survival. Meanwhile, risk scores were calculated based on the constructed risk characteristics (Supplementary Table 2). Each OC Risk score is calculated as follows: Risk Score = (Coef: regression coefficient; Exp: expression value; n: total number of samples; i:the identifier of the ith selected sample). Next, we split the patients with OC into low- and high-risk groups in terms of median risk score, and modeled these two categories as continuous variables to obtain the Hazard Ratio. Then, Kaplan-Meier test was adopted to test the significance of survival curves, and survival curves were drawn. Similarly, we analyzed the survival of OC samples in different clusters and compared the differences in survival between different subgroups. The predictive power of risk characteristics for 1-5-years survival estimates was achieved using nomogram. Nomograms utilize biological and clinical variables (such as tumor grade and patient age) to graphically depict statistical prognostic models that generate the likelihood of clinical events (such as cancer recurrence or death) for a given individual (Balachandran et al., 2015). The nomogram can visually display the results of Cox regression analysis.
The Potential Mechanisms of m6A-Related Genes in OC Patients
Construction of Interaction Networks
The interactions between the m6A-related genes were analyzed and proved with STRING database3. Meanwhile, correlations with a correlation coefficient threshold |r| > = 0.3 and a rank-sum test p ≤ 0.01 were selected using Spearman analysis.
Functional and Pathway Enrichment Analysis
To analyze potential functions and pathways of m6A RNA methylation related genes in OC, Enrichr4 was used for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of m6A-related genes related to OC.
Validation of Clinical Samples
Tissue Samples
For tissue microarray (TMA) cohort, tumor tissues including 8 OC tissue specimens and 8 normal ovarian tissue specimens, were obtained from May 2020 to July 2020 at the Shengjing Affiliated Hospital of China Medical University, China. None of the patients was administered any chemotherapy, immunotherapy, or radiotherapy prior to surgery. Patients with other types of malignant tumors, cardiovascular and cerebrovascular diseases, and mental illness were also excluded. Our study was approved by the Health Research Ethics Board of the Shengjing Affiliated Hospital of China Medical University, and all the cases were pathologically confirmed. The volume of a single tissue sample was about 0.5 cm3. Tissue samples were sharp dissected during surgery and quickly frozen with liquid nitrogen within 15 min after restriction, and then stored at −80∘C until RNA extraction.
Real-Time PCR Array
Trizol reagent (Invitrogen, CA, United States) was utilized to extract total RNA according to the manufacturer’s protocol. For the quantification of 4 m6A prognostic risk model genes, total RNA was then reverse-transcribed into the complementary DNAs (cDNAs) with the PrimeScriptTM RT Reagent Kit with gDNA Eraser (TaKaRa, Dalian, China) and amplified by GoTaq® qPCR Master Mix (Promega, Madison, WI, United States) using the ABI ViiA 7 Real-time PCR system (Applied Biosystems, United States). The specific primer sequences are listed in Supplementary Table 3. GAPDH was used as an internal control for the normalization of Gene expression. 2–ΔΔCT method was used to calculate the relative fold change of expression for samples.
Statistical Analysis
Statistical Analysis
We performed statistical analyses using SPSS 23.0 software (IBM Corp., Armonk, NY, United States) and GraphPad Prism 7 (San Diego, CA, United States). When the p-value (two-sided) ≤ 0.05, the difference was considered statistically significant.
Results
Expression Analysis of m6A RNA Methylation Regulators and m6A-Related Genes
Through univariate analysis of m6A-related genes related to OC in samples with different stages, a total of 36 candidate genes with significant differences in expression between different stages were selected from 267 m6A-related genes related to OC (p ≤ 0.05), of which 1 was an m6A methylation regulator ALKBH5 (Supplementary Figure 1).
Risk Model and Prognostic Analysis of m6A-Related Genes
In order to study the prognostic role of genes related to m6A RNA methylation regulation in OC, we used 306 samples with survival information in the TCGA dataset as the training data set, and all 36 related candidate genes were preliminarily screened for prognostic risk characteristics using Cox univariate regression analysis. We have initially obtained 5 genes that had potential effects on the survival of the samples. The expression level of WDR91 was positively correlated with the survival of OC patients, while the expression levels of NEBL, PDGFRA, ZBTB4, and FAM190A were negatively correlated with the survival of OC patients. The corresponding p-value and Hazard Ratio values of these 5 genes are shown in Supplementary Table 4 and Figure 2A, and the regression coefficients are shown in Figure 2B. LASSO regression analysis further proved that these five genes could constitute risk characteristics and the regression coefficients were obtained (Figure 2C). According to the median risk score calculated, samples were segregated into low- and high-risk groups, and survival differences between these two groups were statistically significant (p = 0.0002; Figure 2D). In order to eliminate the interference of other factors, we performed multivariate Cox regression analysis of candidate risk characteristics on these 5 genes using clinicopathological features such as age, stage, grade, etc. And 4 genes with significant impact on sample survival were obtained (p ≤ 0.05), respectively, NEBL, PDGFRA, WDR91, and ZBTB4. And these 4 genes could be used as independent prognostic markers. Subsequently, samples were split into low- and high-risk groups, and meaningful differences in survival were observed between the two groups (p = 0.0009; Figure 2E).
Figure 2. Construction of the prognostic signature with four m6A-related genes. (A) p-value, Hazard Ratio value and confidence interval for all risk genes with potential independent prognostic efficacy obtained by univariate Cox regression analysis. (B) Regression coefficients of 5 m6A related genes obtained by univariate Cox regression analysis. (C) LASSO regression map of 5 m6A related genes. (D) Survival curve drawn by dividing the TCGA training datasets into high and low risk groups based on the risk score calculated by the feature matrix constructed from the 5 candidate m6A-related genes. (E) Survival curve drawn by dividing the TCGA training datasets into high and low risk groups based on the risk score calculated by the feature matrix constructed from the 4 candidate m6A-related genes. (F) Survival curve drawn by dividing the ICGC validation datasets into high and low risk groups based on the risk score calculated by the feature matrix constructed from the 4 candidate m6A-related genes. (G) Survival time and risk store in TCGA training datasets. (H) Survival time and risk store in ICGC validation datasets.
Moreover, we used 200 samples from the ICGC test dataset to further validate the stability of the risk model and the potential independent prognostic efficacy of m6A-related genes. The results demonstrated that these genes could effectively distinguish the survival of low- and high-risk groups (p = 0.0072; Figure 2F). The sample survival and the model risk score in the TCGA training set and the ICGC validation set are shown in Figures 2G,H.
Validation of TCGA Expression Results Using Clinical Specimens
We examined the expression of 4 potential independent prognostic genes by qRT-PCR in 8 OC tissues and 8 normal ovarian tissues. We applied the unpaired t test to assess the differences between the two groups. The results showed that NEBL, PDGFRA and ZBTB4 were upregulated in OC tissues compared to in normal ovarian tissues, whereas WDR91 were downregulated in tumor tissues (Figure 3A). The mRNA expression results of qRT-PCR validation in 8 patients with OC were supporting effect on the establishment of the four-gene prognostic risk signature in OC.
Figure 3. Relationship between prognostic risk score and clinicopathological characteristics of ovarian cancer. (A) Expression of 4 prognostic predictors in ovarian cancer tissues and normal tissues. (B) The expression profile of clinical characteristics of high and low risk groups. The horizontal axis is four candidate prognostic genes. Red and blue indicate the expression value of each gene corresponding to each sample, and the vertical axis is the sample. Each top is shown from top to bottom. Age, stage, and risk score information for each sample, where all samples are ranked according to risk score from low to high. (C) Differences in risk scores between patients of different ages. (D) Differences in risk scores between patients of different stages. (E) Predicting patients 1–5-years survival with risk scores using nomograms. (F) Survival curves between different new tumor event samples in TCGA datasets.
Prognostic Risk Score and Clinicopathological Features in OC
The heat map depicted the expression of four candidate m6A-related genes in high-risk patients (Figure 3B). We studied the clinicopathological characteristics of OC in the low- and high-risk groups, including age, stage, grade, new neoplasm event type, BRCA1 mutation status, BRCA2 mutation status, 19q13.2 CNV status, and 19q13.42 CNV status. The results illuminated that no meaningful difference was observed in the clinicopathological characteristics between the two groups containing stage, grade, new neoplasm event type, 19q13.2 CNV status, and mutation status of BRCA1 and BRCA2. However, there was significant difference found in 19q13.42 CNV status (p = 0.0390) and different ages (p = 0.0023).
Additionally, we also detected the relationship between risk score and each clinicopathological feature, and found that the risk score was significantly different among patients of different ages (p = 0.0043; Figure 3C), and the difference was also significant among patients at different stages (p = 0.0342; Figure 3D).
We used nomograms to further demonstrate the 1–5 years survival rate predicted by the risk score, including the risk of sample illness, age, stage, and other factors, visualizing the results of Cox regression analysis. A more accurate understanding of survival by looking at the total number of points corresponding to a sample of a certain condition, and the predicted survival rate of the sample in 1–5 years is shown in Figure 3E. Survival differences between different new tumor states are illuminated in Figure 3F, including one case of distant metastasis, four cases of local regional disease, 12 cases of disease progression, and 146 tumor recurrences.
Consensus Clustering of m6A-Related Genes and Related Clinicopathological Characteristics and Survival Outcomes
To study the function of candidate m6A-related genes in OC, we separated the 308 TCGA OC samples into several subgroups using the expression similarity of 36 candidate m6A-related genes through the R “ConsensusClusterPlus” package. Based on the similarity of their expressions, k = 4 was the best k value for relatively stable clustering in a clustering range from 2 to 10 (Supplementary Figures 2A–C). All subgroups were named TM1, TM2, TM3, and TM4, respectively. After using Chi-Square test to analyze the clinicopathological characteristics of samples from the four subgroups of TM1 to TM4, it was found that the four subgroups had significant differences in age, stage and grade (p ≤ 0.05). But there was no notable difference in new neoplasm event type, BRCA1 mutation status, BRCA2 mutation status, 19q13.2 CNV status, and 19q13.42 CNV status. We further investigated the survival status between the four subgroups and uncovered that the difference in survival rates between these subgroups was not significant (p = 0.0998). What’s more, we separated the samples into several subgroups using the expression similarity of 4 prognostic genes. Utilizing consistent cluster analysis, the samples could be clearly divided into four categories with four prognostic genes, as shown in Figures 4A,B. Figure 4C demonstrated that the inflection point was larger when k = 4 and could be divided into four categories. What’s more, it was found that the new four subgroups had significant differences in age, stage and grade (p ≤ 0.05) (Figure 4D). We also investigated the survival status between the four subgroups and uncovered that the OS and DFS of TM2 group were significantly shorter than those of the other three groups (p = 0.0003, p = 0.002; Figures 4E,F).
Figure 4. Identification of consensus clusters by 4 prognostic genes. (A) Relative change in area under CDF curve for k = 2–10 classified by the prognostic risk model. (B) Consensus clustering cumulative distribution function (CDF) for k = 2–10 classified by the prognostic risk model. (C) Consensus clustering matrix for k = 4. (D) Heat maps and significant clinical features of the four subgroups. (E) Kaplan–Meier overall survival (OS) curves for 308 samples in the four subgroups classified by 4 prognostic genes. (F) Kaplan–Meier disease-free survival (DFS) curves for 308 samples in the four subgroups classified by four prognostic genes.
Interaction and Correlation Analysis of m6A Candidate Gene set
To further understand the interactions between the 36 m6A-related genes, we analyzed the interactions and correlations between these genes. The interactions between the 36 m6A-related genes are demonstrated in Figure 5A. When the correlation coefficient threshold of the expression amount was set to | r | ≥ 0.3, and the p-value of the rank sum test was set to p ≤ 0.01, a total of 67 pairs of significantly correlated interaction factors were obtained (Supplementary Table 5). The results of the correlation analysis of all 36 candidate genes are shown in Figure 5B, in which the blue ones are positive and the red ones are negative.
Figure 5. The interaction and correlation among m6A-related genes. (A) PPI network was constructed to evaluate the interaction among m6A-related genes. (B) The Pearson correlation analysis was used to determine the correlation among m6A-related genes.
GO and KEGG Pathway Analyses of Differentially Expressed m6A-Related Genes
We performed differential gene analysis on candidate genes between the four subgroups of TM1 to TM4 to find m6A-related genes that were differentially expressed in the four subgroups. We identified 9 genes with significant differential expression among the 36 m6A-related genes using one-way ANOVA (p ≤ 0.05; Supplementary Table 6). We performed GO and KEGG functional annotation of differential genes using functional enrichment analysis tool (Supplementary Tables 7, 8). The results demonstrated that these differentially expressed genes were significantly related to biological functions including regulation of actin filament organization, sphingomyelin metabolism, transcription of DNA templates, and apoptotic signaling pathways (Figure 6A). Meanwhile, the differential genes were closely related to sphingolipid metabolism pathways, which were vital pathways involved in OC initiation (Figure 6B).
Figure 6. Functional enrichment analysis of differentially expressed genes between subgroups. (A) GO function analysis of differential genes among four subgroups. (B) KEGG pathway analysis of differential genes among four subgroups.
Among these differential genes, NEBL is a candidate gene with potential independent prognostic efficacy. The design of subsequent experiments for these genes can be considered to further investigate the potential pathological value of these genes and the design of new drug targets.
Discussion
OC is one of the three most prevalent cancer types of female genitalia with the highest mortality rate, yet its pathogenesis is still unclear (Bamford and Webster, 2017). Thus, uncovering the intrinsic molecular mechanisms of OC tumorigenesis is of great significance. Aberrant expression of m6A-related genes has been identified in multiple cancers, which is proved to be closely related to pathological processes in cancer, including tumorigenesis, metastasis, and drug resistance (Ma et al., 2019). However, m6A-related genes have rarely been studied in OC. Thus, we conducted this study to investigate their roles in OC. Differential expression patterns for most m6A-related genes were identified in OC patients compared with controls based on RNA sequencing data from TCGA and ICGC. Significant difference was found in survival between the risk groups of samples. What’s more, a four-gene risk signature including NEBL, PDGFRA, WDR91 and ZBTB4 was confirmed, and this risk signature was considered as an independent predictor for prognosis of OC patients. Subsequently, we identified consensus clustering into four subgroups based on the expression of these four prognostic genes, and uncovered statistically significant differences of age, stage and grade among the subgroups. Moreover, the OS and DFS of TM2 group were shorter than those of the other three groups.
In terms of data from TCGA, m6A RNA methylation regulator ALKBH5 and other 35 m6A-related genes were dysregulated in OC patients, which indicated that m6A RNA methylation related genes might play a role in OC. Similar results were previously reported. Zhu et al. suggested that ALKBH5 regulated the proliferation, invasion and autophagy of OC cells via EGFR-PIK3CA-AKT-mTOR pathway and Bcl-2 (Zhu et al., 2019). In addition, METTL3 stimulates epithelial-mesenchymal transition (EMT) of cancer cells by stimulating receptor tyrosine kinase AXL, thereby enhancing the invasion and metastasis of OC (Hua et al., 2018). Han et al. (2020) studied the prognostic value of high-frequency genetic alterations of m6A RNA methylation regulators in OC. While we mainly focused on the prognostic value of m6A-related genes themselves for OC, and we identified a four-gene risk signature (Han et al., 2020). We also demonstrated distinct relationship between the expression of m6A-related genes and clinicopathological features in OC, such as 19q13.2_CNV, age and stage. 19q13.2_CNV was reported associated with decreased OC risk (Walker et al., 2017). What’s more, there were significant difference of age, stage and grade among the four subgroups using consensus clustering analysis, and the OS as well as DFS of TM2 group were shorter than those of the other three groups. The results above illustrate that m6A-related genes are likely to serve a vital function in the development/progression of OC, and further in-depth research is warranted to determine the underlying molecular mechanisms.
We then evaluated the effects of m6A-related gene alterations on survival of OC patients. Abundant studies have reported that the dysregulation of m6A-related genes was related to the prognosis of various cancers, such as breast cancer, gastric cancer, renal cell carcinoma, etc. (Liu et al., 2019; Su et al., 2019; Wang J. et al., 2020) In our study, we obtained 4 m6A-related genes that could be used as a feature of potential independent prognostic risk from the 36 candidate genes. Our four−gene prognostic signature implied that the expression level of NEBL, PDGFRA, ZBTB4 was negatively correlated with the survival of OC patients. Previous studies indicated that ZBTB4 and PDGFRA played roles in tumor progression, same as our prediction (Kim et al., 2012; Kang et al., 2015; Roussel-Gervais et al., 2017; Yu Y. et al., 2018; Ye et al., 2019). ZBTB4 is a mammalian DNA-binding protein acting as a transcriptional repressor (Yu Y. et al., 2018). ZBTB4 overexpression inhibits cancer cell proliferation and induces cell cycle arrest at G1 phase as well as apoptosis in Ewing’s sarcoma (Yu Y. et al., 2018). Roussel-Gervais et al. (2017) reported that decreased ZBTB4 expression correlated with the high genome instability among many frequent human cancers, which altered mitotic checkpoint, increased aneuploidy and promoted tumorigenesis. Similarly, underexpression of ZBTB4 is correlated with poor survival of breast cancer patients (Kim et al., 2012). PDGFRA is a gene encoding cell surface tyrosine kinase receptor (Ye et al., 2019). For PDGFRA, PDGFRA can promote downstream activation of the Notch1 pathway as well as the angiogenesis, proliferation and invasion of OC cells (Ye et al., 2019). And mutations in PDGFRA have been related to a variety of other cancers, including glioblastoma, melanoma, neuroendocrine carcinoma, etc. (Kang et al., 2015). However, the role of NEBL in cancer is complex. NEBL is a member of nebulin family of actin binding proteins (Zhang and Zhang, 2019). Zhang et.al found that NEBL promoted cancer cell proliferation, migration and invasion in cervical cancer by regulating PI3K/Akt pathway (Zhang and Zhang, 2019). On the contrary, upregulation of NEBL inhibited cancer cell migration and invasion and reversed TGF-β-induced EMT in prostate cancer (Wang et al., 2017). These indicate that m6A RNA methylation related genes might play different roles as tumor promoter or suppressor agents in different types of cancer, which need to be further elaborated.
Our GO analysis demonstrated that differentially expressed genes were significantly related to biological functions including regulation of actin filament organization and reorganization, sphingomyelin metabolic process, cardiac myofibril assembly, etc. Meanwhile, our KEGG analysis indicated that the differential genes were closely related to sphingolipid metabolism pathways. Previous studies have indicated that sphingolipid metabolic pathway participated in the regulation of vital cancer cellular processes, such as cell proliferation, migration, invasion, apoptosis and autophagy, playing an important role in the occurrence and development of OC (Hannun and Obeid, 2018; Jacob et al., 2018; Ogretmen, 2018). Meanwhile, key enzymes of sphingolipid metabolism were thought to be directly related to drug resistance in OC (Huang et al., 2016). These suggest that sphingolipid metabolism pathways may be a possible vital link in m6A-related genes involved in regulating OC.
However, there are also some potential limitations in the current study. First, there were no available datasets of tissues adjacent to OC from TCGA, and we used OC and normal ovarian tissue samples to verify the expression. Second, p ≤ 0.05 was considered statistically significant in our study which might have influence on the reliability and accuracy of the results (Colquhoun, 2017). Thus p ≤ 0.005 should be applied for our future research. Third, the m6A-related genes related to OC that we included were directly generated from the m6Avar database, and their regulatory functions might not have been verified enough. Subsequent experiments are needed to further verify the specific functioning mechanisms of these m6A-related genes, especially the 4 prognostic biomarkers we analyzed in this study.
Conclusion
In conclusion, our study for the first time analyzed the expression of m6A-related genes in OC, and discovered that m6A-related genes were tightly relevant to the prognosis of OC patients, highlighting their roles as prognostic biomarkers in OC patients, as well as their potential functions in the occurrence and progression of OC. Further research is required to investigate the regulatory mechanisms of m6A modification in OC, which will help develop m6A RNA methylation related genes as valuable therapeutic targets.
Data Availability Statement
Publicly available datasets were analyzed in this study, these can be found in The Cancer Genome Atlas (https://portal.gdc.cancer.gov/) the UCSC Xena Browser database (http://xena.ucsc.edu/).
Ethics Statement
The studies involving human participants were reviewed and approved by the Health Research Ethics Board of the Shengjing Affiliated Hospital of China Medical University. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.
Author Contributions
ZN and XW designed the study. ZN performed the data collection and the data analysis. ZN and LF drafted the manuscript. All authors read and approved the final version of the manuscript.
Funding
Work toward this research was supported by the National Natural Science Foundation of China (No. 81671423), the National Key Research and Development Program of China (No. 2018YFC1003100), and Freedom Researcher Fund of Shengjing Hospital (No. M0476). ZN was supported by the Ph.D. Scholarship from the Shengjing Hospital of China Medical University.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.542457/full#supplementary-material
Supplementary Figure 1 | The expression levels of m6A methylation regulators of samples at different stages in TCGA ovarian cancer datasets (For stages: dark blue represents stage I, yellow represents stage II, red represents stage III, and light blue represents stage IV. For expression levels: red represents high expression, and blue represents low expression).
Supplementary Figure 2 | Identification of consensus clusters by m6A-related genes. (A) Relative change in area under CDF curve for k = 2–10 classified by m6A-related genes. (B) Consensus clustering cumulative distribution function (CDF) for k = 2–10 classified by m6A-related genes. (C) Consensus clustering matrix for k = 4.
Footnotes
- ^ http://xena.ucsc.edu/
- ^ http://m6avar.renlab.org/
- ^ http://www.string-db.org/
- ^ http://amp.pharm.mssm.edu/Enrichr/
References
Balachandran, V. P., Gonen, M., Smith, J. J., and DeMatteo, R. P. (2015). Nomograms in oncology: more than meets the eye. Lancet Oncol. 16, e173–e180. doi: 10.1016/S1470-2045(14)71116-7
Bamford, J., and Webster, R. M. (2017). The ovarian cancer drug market. Nat. Rev. Drug Discov. 16, 451–452. doi: 10.1038/nrd.2017.92
Bray, F., Ferlay, J., Soerjomataram, I., Siegel, R. L., Torre, L. A., and Jemal, A. (2018). Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. Cancer J. Clin. 144, 1941–1953. doi: 10.1002/ijc.31937
Cheng, X., Li, M., Rao, X., Zhang, W., Li, X., Wang, L., et al. (2019). KIAA1429 regulates the migration and invasion of hepatocellular carcinoma by altering m6A modification of ID2 mRNA. Onco. Targets Ther. 12, 3421–3428. doi: 10.2147/OTT.S180954
Coleman, R. L., Fleming, G. F., Brady, M. F., Swisher, E. M., Steffensen, K. D., Friedlander, M., et al. (2019). Veliparib with first-line chemotherapy and as maintenance therapy in ovarian cancer. N. Engl. J. Med. 381, 2403–2415. doi: 10.1056/NEJMoa1909707
Colquhoun, D. (2017). The reproducibility of research and the misinterpretation of p-values. R. Soc. Open Sci. 4:171085. doi: 10.1098/rsos.171085
Dai, D., Wang, H., Zhu, L., Jin, H., and Wang, X. (2018). N6-methyladenosine links RNA metabolism to cancer progression. Cell Death Dis. 9:124. doi: 10.1038/s41419-017-0129-x
Deng, X., Su, R., Weng, H., Huang, H., Li, Z., and Chen, J. (2018). RNA N6-methyladenosine modification in cancers: current status and perspectives. Cell Res. 28, 507–517. doi: 10.1038/s41422-018-0034-6
Han, X., Liu, J., Cheng, G., and Cui, S. (2020). Gene signatures and prognostic values of m6a rna methylation regulators in ovarian cancer. Cancer Control. 27:1073274820960460. doi: 10.1177/1073274820960460
Hannun, Y. A., and Obeid, L. M. (2018). Sphingolipids and their metabolism in physiology and disease. Nat. Rev. Mol. Cell Biol. 19, 175–191. doi: 10.1038/nrm.2017.107
Hua, W., Zhao, Y., Jin, X., Yu, D., He, J., Xie, D., et al. (2018). METTL3 promotes ovarian carcinoma growth and invasion through the regulation of AXL translation and epithelial to mesenchymal transition. Gynecol. Oncol. 151, 356–365. doi: 10.1016/j.ygyno.2018.09.015
Huang, H., Tong, T. T., Yau, L. F., Chen, C. Y., Mi, J. N., Wang, J. R., et al. (2016). LC-MS based sphingolipidomic study on A2780 human ovarian cancer cell line and its taxol-resistant strain. Sci. Rep. 6:34684. doi: 10.1038/srep34684
Ianniello, Z., and Fatica, A. (2018). N6-Methyladenosine role in acute Myeloid Leukaemia. Int. J. Mol. Sci. 19:E2345. doi: 10.3390/ijms19082345
Jacob, F., Alam, S., Konantz, M., Liang, C. Y., Kohler, R. S., Everest-Dass, A. V., et al. (2018). Transition of mesenchymal and epithelial cancer cells depends on α1-4 galactosyltransferase-mediated glycosphingolipids. Cancer Res. 78, 2952–2965. doi: 10.1158/0008-5472.CAN-17-2223
Jiang, X. L., Liu, B. Y., Nie, Z., Duan, L. C., Xiong, Q. X., Jin, Z. X., et al. (2021). The role of m6A modification in the biological functions and diseases. Signal Transduct. Target Ther. 6:74. doi: 10.1038/s41392-020-00450-x
Kang, G., Bae, B. N., Sohn, B. S., Pyo, J. S., Kang, G. H., and Kim, K. M. (2015). Detection of KIT and PDGFRA mutations in the plasma of patients with gastrointestinal stromal tumor. Target Oncol. 10, 597–601. doi: 10.1007/s11523-015-0361-1
Kim, K., Chadalapaka, G., Lee, S. O., Yamada, D., Sastre-Garau, X., Defossez, P. A., et al. (2012). Identification of oncogenic microRNA-17-92/ZBTB4/specificity protein axis in breast cancer. Oncogene 31, 1034–1044. doi: 10.1038/onc.2011.296
Lence, T., Paolantoni, C., Worpenberg, L., and Roignant, J. Y. (2019). Mechanistic insights into m6A RNA enzymes. Biochim. Biophys. Acta Gene Regul. Mech. 1862, 222–229. doi: 10.1016/j.bbagrm.2018.10.014
Li, J., Han, Y., Zhang, H., Qian, Z., Jia, W., Gao, Y., et al. (2019). The m6A demethylase FTO promotes the growth of lung cancer cells by regulating the m6A level of USP7 mRNA. Biochem. Biophys. Res. Commun. 512, 479–485. doi: 10.1016/j.bbrc.2019.03.093
Liu, J., Eckert, M. A., Harada, B. T., Liu, S. M., Lu, Z., Yu, K., et al. (2018). m6A mRNA methylation regulates AKT activity to promote the proliferation and tumorigenicity of endometrial cancer. Nat Cell Biol. 20, 1074–1083. doi: 10.1038/s41556-018-0174-4
Liu, L., Liu, X., Dong, Z., Li, J., Yu, Y., Chen, X., et al. (2019). N6-methyladenosine-related genomic targets are altered in breast cancer tissue and associated with poor survival. J. Cancer 10, 5447–5459. doi: 10.7150/jca.35053
Liu, T., Yang, S., Sui, J., Xu, S. Y., Cheng, Y. P., Shen, B., et al. (2020). Dysregulated N6-methyladenosine methylation writer METTL3 contributes to the proliferation and migration of gastric cancer. J. Cell Physiol. 235, 548–562. doi: 10.1002/jcp.28994
Livneh, I., Moshkovitz, S. M., Amariglio, N., Rechavi, G., and Dominissini, D. (2020). The m6A epitranscriptome: transcriptome plasticity in brain development and function. Nat. Rev. Neurosci. 21, 36–51. doi: 10.1038/s41583-019-0244-z
Ma, S., Chen, C., Ji, X., Liu, J., Zhou, Q., Wang, G., et al. (2019). The interplay between m6A RNA methylation and noncoding RNA in cancer. J. Hematol. Onco. 12:121. doi: 10.1186/s13045-019-0805-7
Ni, W., Yao, S., Zhou, Y., Liu, Y., Huang, P., Zhou, A., et al. (2019). Long noncoding RNA GAS5 inhibits progression of colorectal cancer by interacting with and triggering YAP phosphorylation and degradation and is negatively regulated by the m6A reader YTHDF3. Mol. Cancer 18:143. doi: 10.1186/s12943-019-1079-y
Niu, Y., Lin, Z., Wan, A., Chen, H., Liang, H., Sun, L., et al. (2019). RNA N6-methyladenosine demethylase FTO promotes breast tumor progression through inhibiting BNIP3. Mol. Cancer 18, 46. doi: 10.1186/s12943-019-1004-4
Ogretmen, B. (2018). Sphingolipid metabolism in cancer signalling and therapy. Nat. Rev. Cancer 18, 33–50. doi: 10.1038/nrc.2017.96
Oldak, S., Ioannou, S., Kamath, P., Huang, M., George, S., Slomovitz, B., et al. (2019). Polypharmacy in patients with ovarian cancer. Oncologist 24, 1201–1208. doi: 10.1634/theoncologist.2018-0807
Roussel-Gervais, A., Naciri, I., Kirsh, O., Kasprzyk, L., Velasco, G., Grillo, G., et al. (2017). Loss of the Methyl-CpG-binding protein ZBTB4 alters mitotic checkpoint, increases aneuploidy, and promotes tumorigenesis. Cancer Res. 77, 62–73. doi: 10.1158/0008-5472.CAN-16-1181
Shi, H., Wei, J., and He, C. (2019). Where, when, and how: context-dependent functions of RNA methylation writers, readers, and erasers. Mol. Cell. 74, 640–650. doi: 10.1016/j.molcel.2019.04.025
Su, Y., Huang, J., and Hu, J. (2019). m6A RNA methylation regulators contribute to malignant progression and have clinical prognostic impact in gastric cancer. Front. Oncol. 9:1038. doi: 10.3389/fonc.2019.01038
Sun, T., Wu, R., and Ming, L. (2019). The role of m6A RNA methylation in cancer. Biomed. Pharmacother. 112:108613. doi: 10.1016/j.biopha.2019.108613
Torre, L. A., Trabert, B., DeSantis, C. E., Miller, K. D., Samimi, G., Runowicz, C. D., et al. (2018). Ovarian cancer statistics, 2018. Cancer J. Clin. 68, 284–296. doi: 10.3322/caac.21456
Walker, L. C., Marquart, L., Pearson, J. F., Wiggins, G. A., O’Mara, T. A., Parsons, M. T., et al. (2017). Evaluation of copy-number variants as modifiers of breast and ovarian cancer risk for BRCA1 pathogenic variant carriers. Eur. J. Hum. Genet. 25, 432–438. doi: 10.1038/ejhg.2016.203
Wang, B., Zhang, L., Zhao, L., Zhou, R., Ding, Y., Li, G., et al. (2017). LASP2 suppresses colorectal cancer progression through JNK/p38 MAPK pathway meditated epithelial-mesenchymal transition. Cell Commun. Signal 15:21. doi: 10.1186/s12964-017-0179-9
Wang, J., Zhang, C., He, W., and Gou, X. (2020). Effect of m6A RNA methylation regulators on malignant progression and prognosis in renal clear cell carcinoma. Front. Oncol. 10:3. doi: 10.3389/fonc.2020.00003
Wang, T. Y., Kong, S., Tao, M., and Ju, S. Q. (2020). The potential role of RNA N6-methyladenosine in cancer progression. Mol. Cancer 19:88. doi: 10.1186/s12943-020-01204-7
Wang, X., Wu, R., Liu, Y., Zhao, Y., Bi, Z., Yao, Y., et al. (2019). m6A mRNA methylation controls autophagy and adipogenesis by targeting Atg5 and Atg7. Autophagy 26, 1–15. doi: 10.1080/15548627.2019.1659617
Yang, Y., Hsu, P. J., Chen, Y. S., and Yang, Y. G. (2018). Dynamic transcriptomic m6A decoration: writers, erasers, readers and functions in RNA metabolism. Cell Res. 28, 616–624. doi: 10.1038/s41422-018-0040-8
Ye, W., Ni, Z., Yicheng, S., Pan, H., Huang, Y., and Xiong, Y. (2019). Anisomycin inhibits angiogenesis in ovarian cancer by attenuating the molecular sponge effect of the lncRNA Meg3/miR 421/PDGFRA axis. Int. J. Oncol. 55, 1296–1312. doi: 10.3892/ijo.2019.4887
Yu, J., Chen, M., Huang, H., Zhu, J., Song, H., Zhu, J., et al. (2018). Dynamic m6A modification regulates local translation of mRNA in axons. Nucleic Acids Res. 46, 1412–1423. doi: 10.1093/nar/gkx1182
Yu, Y., Shang, R., Chen, Y., Li, J., Liang, Z., Hu, J., et al. (2018). Tumor suppressive ZBTB4 inhibits cell growth by regulating cell cycle progression and apoptosis in Ewing sarcoma. Biomed. Pharmacother. 100, 108–115. doi: 10.1016/j.biopha.2018.01.132
Yue, B., Song, C., Yang, L., Cui, R., Cheng, X., Zhang, Z., et al. (2019). METTL3-mediated N6-methyladenosine modification is critical for epithelial-mesenchymal transition and metastasis of gastric cancer. Mol. Cancer 18:142. doi: 10.1186/s12943-019-1065-4
Zhang, S., Zhao, B. S., Zhou, A., Lin, K., Zheng, S., Lu, Z., et al. (2017). m6A Demethylase ALKBH5 maintains tumorigenicity of glioblastoma stem-like cells by sustaining FOXM1 expression and cell proliferation program. Cancer Cell 31, 591–606.e6. doi: 10.1016/j.ccell.2017.02.013
Zhang, Y., and Zhang, L. (2019). Knockdown of LASP2 inhibits the proliferation, migration, and invasion of cervical cancer cells. J. Cell Biochem. 120, 15389–15396. doi: 10.1002/jcb.28806
Zhao, Y. C., Shi, Y. F., Shen, H. F., and Xie, W. Z. (2020). m6A-binding proteins: the emerging crucial performers in epigenetics. J. Hematol. Oncol. 13:35. doi: 10.1186/s13045-020-00872-8
Keywords: subgroup, prognosis, m6A-related genes, m6A RNA methylation regulators, ovarian cancer
Citation: Na Z, Fan L and Wang X (2021) Gene Signatures and Prognostic Values of N6-Methyladenosine Related Genes in Ovarian Cancer. Front. Genet. 12:542457. doi: 10.3389/fgene.2021.542457
Received: 14 March 2020; Accepted: 19 April 2021;
Published: 18 August 2021.
Edited by:
Shankar Subramaniam, University of California, San Diego, United StatesReviewed by:
Dechao Bu, Institute of Computing Technology, Chinese Academy of Sciences (CAS), ChinaGonzalo Riadi, University of Talca, Chile
Copyright © 2021 Na, Fan and Wang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Ling Fan, fanl@sj-hospital.org; Xiuxia Wang, wangxxsj@sina.cn