RNA 5-Methylcytosine Regulators Contribute to Metabolism Heterogeneity and Predict Prognosis in Ovarian Cancer
Jie Xu
Xiaoyi Liu
Yanjie Chen
Yuya Wang
Tao Liu 
Ping Yi- Department of Obstetrics and Gynecology, The Third Affiliated Hospital of Chongqing Medical University, Chongqing, China
5-Methylcytosine (m5C) is an abundant and highly conserved modification in RNAs. The dysregulation of RNA m5C methylation has been reported in cancers, but the regulatory network in ovarian cancer of RNA m5C methylation-related genes and its implication in metabolic regulation remain largely unexplored. In this study, RNA-sequencing data and clinical information of 374 ovarian cancer patients were downloaded from The Cancer Genome Atlas database, and a total of 14 RNA m5C regulators were included. Through unsupervised consensus clustering, two clusters with different m5C modification patterns were identified with distinct survivals. According to enrichment analyses, glycosaminoglycan and collagen metabolism–related pathways were specifically activated in cluster 1, whereas fatty acid metabolism–related pathways were enriched in cluster 2, which had better overall survival (OS). Besides the metabolism heterogeneity, the higher sensitivity to platinum and paclitaxel in cluster 2 can further explain the improved OS. Ultimately, a least absolute shrinkage and selection operator prediction model formed by ALYREF, NOP2, and TET2 toward OS was constructed. In conclusion, distinct m5C modification pattern exhibited metabolism heterogeneity, different chemotherapy sensitivity, and consequently survival difference, providing evidence for risk stratification.
Introduction
Ovarian cancer (OVC) is the most lethal gynecological cancer (Siegel et al., 2021). Because of asymptomatic onset and lack of efficient screening tests, more than 75% of patients are diagnosed at an advanced stage with a 5-year survival rate of 29%, in contrast to 92% for early stage (Singer et al., 2003). The standard frontline care is the debulking surgery to no tumor residual and platinum-based adjuvant chemotherapy, with antiangiogenic therapy applied in patients who have suboptimal tumor reduction and stage IV disease (Lheureux et al., 2019). The poly-ADP-ribose polymerase (PARP) inhibitors have been applied in frontline care for maintenance therapy and in patients with recurrence (Kristeleit et al., 2017; Pujade-Lauraine et al., 2017; Moore et al., 2018). However, the moderate activity of PARP inhibitors was found in patients with homologous recombination dysfunction, and a worse therapeutic effect was observed in homologous recombination–proficient patients (Kaufman et al., 2015). Despite initial response to the first-line treatment, 25% of patients have a relapse within 6 months (Jemal et al., 2003), and more than 80% of patients eventually have a recurrence (Kim et al., 2020). Immunotherapy has demonstrated modest response rates of 10% to 15%, despite a large proportion of OVCs with high expression of programmed death ligand 1 (Pujade-Lauraine, 2017). Facing those challenges in diagnosis and treatment, seeking predictive biomarkers could enable early diagnosis, survival prediction, and identification of patient subgroups who would maximally benefit from those treatments (Lheureux et al., 2019). Consequently, investigations are devoted to molecular and function profiling of OVC for optimal biomarkers.
Dysregulation of RNA expression profile is an important hallmark of tumors (Chai et al., 2019a). RNA 5-methylcytosine (m5C) modification, the methylation of the fifth carbon in cytosine base in RNA sequences, has emerged as one of the critical posttranscriptional regulators of gene expression and has been identified in tRNA, rRNA, and mRNA (Nombela et al., 2021). The distribution of m5C site in mRNA has been reported to be mainly deposited in the coding sequences and enriched around the translation initiation codon (Amort et al., 2017; Yang et al., 2017; Yang et al., 2019; Tang et al., 2020). RNA m5C modification is a reversible and dynamic process mediated by a group of proteins named “writers,” “erasers,” and “readers,” which work as methyltransferases (NSUN, DNMT, and TRDMT families), demethylases (TET family), and binding proteins (ALYREF and YBX1), respectively. RNA m5C modification has been involved in the regulation of gene expression (Roundtree et al., 2017; Chai et al., 2019b) and thus has participated in a series of physiological and pathological processes including cancers (Chen et al., 2021).
Cancer is considered as a disease characterized by the accumulation of genetic or epigenetic alterations of different oncogenes and tumor suppressors (Nombela et al., 2021). Metabolism reprogramming is another indispensable hallmark of cancer. Mutation of oncogene and tumor suppressors drives the reprogramming of metabolism and rewiring of epigenetic modification. Cancer cell fate can also be modified by epigenetic modification and alteration of metabolites. m5C regulator dysregulation has been reported in multiple cancers such as breast cancer, leukemia, bladder cancer, and skin squamous cell carcinoma (Freeman et al., 1991; Blanco et al., 2016; Cheng et al., 2018; Chen et al., 2019). It has also been demonstrated that RNA m5C modification could promote glucose metabolism through enhancing PKM2 mRNA stability in bladder cancer (Wang et al., 2021). However, the role of RNA m5C regulators-mediated m5C methylation modification, as well as its effect on metabolism reprogramming in OVC, remains unclear.
In this study, we revealed the landscape of genetic variation and gene expression of m5C regulators in OVC and established a prognostic prediction model formed by ALYRER, NOP2, and TET2 for overall survival (OS). We also dissected the potential roles of m5C modification in metabolism heterogeneity and altered chemotherapeutic drug sensitivity, which could result in survival differences of OVC patients.
Materials and Methods
Data Resources
The workflow of our study is shown in Supplementary Figure S1A . The fragments per kilobase of exon model per million mapped fragments (FPKM) files of RNA-seq transcriptome data, as well as clinical information of 374 cases of OVC, were downloaded from The Cancer Genome Atlas (TCGA) database. The SOFT formatted matrix files of three Gene Expression Omnibus (GEO) datasets (GSE27651, GSE52037, GSE54388, and GSE19829) (Clough and Barrett, 2016) were downloaded using R package getGEO. The Masked Copy Number Segment data of DNA copy number variation (CNV) data of OVC were downloaded from the Genomic Data Commons (https://portal.gdc.cancer.gov).
RNA m5C Regulators
Fourteen m5C regulators including eight writers (NOP2, NSUN2, NSUN3, NSUN4, NSUN5, NSUN6, NSUN7, TRDMT1); four erasers (TET1, TET2, TET3, ALKBH1); and two readers (YBX1, ALYREF) were enrolled in this study (Chen et al., 2021). DMNT3A and DMNT3B were excluded because they have only been reported in Arabidopsis thaliana for now.
m5C Regulators Mutation and CNV Analysis
The somatic mutation investigation of m5C regulators in pan-cancer and the CNV analysis in OVC were performed using cBioPortal website (www.cbioportal.org) (Gao et al., 2013). The Pan-Cancer Project of TCGA was enrolled for somatic mutation evaluation in pan-cancer. The Pan-Cancer Project of TCGA-OV with both somatic mutation and mRNA data was enrolled for CNV analysis in OVC.
Differentially Expressed Gene Analysis
Principal component analysis (PCA) using R package FactoMineR and differential gene expression analysis using R package limma (Ritchie et al., 2015) were conducted, in order to display the different profiles of m5C regulators between human ovarian surface epithelium (HOSE) and OVC. Differential analysis was also utilized in seeking DEGs that were specifically up-regulated in each cluster. DEGs were defined as genes with p < 0.05 and |fold change| >1.2.
Interaction Between 14 m5C Regulators
The protein–protein interaction (PPI) network plot was constructed using the STRING 11.0 b website (https://string-db.org/). The correlation analysis of the m5C regulators among mRNA expression and CNV level and between them was conducted by R package corrplot.
Clustering Analysis of 14 m5C Regulators
The ConsensusClusterPlus package (Wilkerson and Hayes, 2010) was performed to identify distinct m5C phenotype based on the expression of 14 m5C regulators, and 1,000 times repetitions were conducted to ensure the stability of the classification.
Cluster Function Annotation and Exploration of Cluster Metabolism Heterogeneity
The cluster function annotation was conducted using R package Gene Set Variation Analysis (GSVA) to explore the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment among different m5C clusters. GO hallmark and KEGG hallmark gene sets were downloaded from MsigDB dataset (http://www.gsea-msigdb.org/gsea/msigdb). Further functional annotation of each m5C cluster was performed by R package ClusterProfiler (Yu et al., 2012) using the top 500 expressed genes and DEGs in each cluster for GO and KEGG pathways. Then, Gene Set Enrichment Analysis (GSEA) was performed using cluster DEGs by ClusterProfiler package to further identify up-regulated pathways in the individual cluster. The up-regulated pathways identified in multiple methods were finally visualized in circos plot using R package circlize. Then, metabolic pathways were downloaded from KEGG database including 1,653 human genes assigned to 91 pathways. The GSVA scores of metabolic pathways were calculated using the GSVA package for further correlation with mRNA level of m5C regulators.
Prediction of Drug Sensitivity
The drug sensitivity was predicted using calcPhenotype function in R package oncoPredict (Maeser et al., 2021) based on Genomics of Drug Sensitivity in Cancer cell line dataset (https://www.cancerrxgene.org/). The prediction ability by drug sensitivity score calculation was validated in OVC clinical trial (GSE51373) with area under the curve (AUC) of 0.786. Imputed lower sensitivity score represents higher sensitivity of the drug.
Cell Culture and Cell Growth
OVCAR3 cells were cultured in Dulbecco modified eagle medium (GIBCO, United States) supplemented with 10% fetal bovine serum (GIBCO, United States), penicillin (GIBCO, United States), and streptomycin (GIBCO, United States) and maintained at 37°C in 5% CO2 cell culture incubator. The cells transfected with siRNA targeting TET2 or control (TsingkeBiotechnology, China) were seeded in 96-well plates for the cell viability test. CCK-8 reagent (DOJINDO, Japan) was added into the plate and incubated for 2 h. The cell absorbance at 450-nm wavelengths was measured by using the microplate reader (BioTek, United States) at 0, 1, 2, 3, and 4 days. All experiments were performed in triplicate.
Western Blot
Cells were collected and lysed with cell lysis buffer (Beyotime, China) on ice for 30 min, and the lysate was obtained by centrifugation at 12,000g for 10 min. Proteins were fractionated by sodium dodecyl sulfate–polyacrylamide gel electrophoresis and then transferred onto 0.22-μm NC membranes. The membranes were blocked with 5% nonfat milk in TBS/Tween-20 and blotted with the antibody (anti-TET2; ProteinTech) at 4°C overnight. Corresponding secondary antibodies (ZSGB-BIO, China) were added on the membrane at room temperature for 1.5 h. Immunoreactive bands were visualized using enhanced chemiluminescence detection reagent (Millipore, United States).
Statistical Analysis
Correlations among different m5C regulators were evaluated by Spearman correlation analyses using R package corrplot. Correlations between m5C regulators and clinicopathological parameters were evaluated by Spearman correlation analyses using SPSS 25.0. A χ2 test was conducted to compare the clinicopathological parameters of clusters. Kruskal–Wallis test was used to compare gene expression among different samples. R packages survival and survminer (Scrucca et al., 2007) were used to perform the univariate Cox proportional hazards analysis and Kaplan–Meier analysis for OS. R package forestplot and survminer were used for visualized the Cox analysis results and survival curves, respectively. Genes with p < 0.05 in univariate analysis were selected to the least absolute shrinkage and selection operator (LASSO) method regression analysis using R package glmnet. Patients with survival information were randomly divided into two groups (75% in the training group and 25% in the test group) by createDataPartition function from R package caret. Three gene signature and their corresponding coefficients were determined in the training group by glmnet package, and the risk score was calculated for each patient using the prediction function. The AUC of receiver operating characteristic (ROC) curve was calculated by R package survivalROC (Heagerty and Zheng, 2005). True positive (TP) and false positive (FP) of every patient in the training group were calculated through survivalROC function, and the minimum value of the formula (TP-1)*2 + FP*2 was determined as the best cutoff value. This cutoff value was used in the internal training set, internal testing set, and external testing set to divide the samples into the high-score group and the low-score group. R 4.0.3 was used for all the statistical analyses in this study. p < 0.05 is the significance threshold for all the data.
Results
Profiles of Genetic Variation and Gene Expression of RNA m5C Regulators in Ovarian Cancer
First, we comprehensively studied the profile of the genetic variation frequency of m5C regulators in the pan-cancer cohort. The amplification is the prevalent variation pattern of m5C regulator genes in OVC, and 13 of 14 (92.9%) regulators were amplified (Figures 1A,B). Among those regulators, YBX1, NOP2, and NSUN4 genes exhibited the highest amplification frequencies of 7%, 6%, and 5%, respectively (Figure 1B). The significant positive correlation of CNV among regulators was demonstrated especially between ALYREF and writers, as well as TET2 and writers (Figure 1C). Then, we explored the correlation between CNV and mRNA levels of each regulator and found a significant positive correlation in all 14 regulators (Figure 1D). Differential analysis was further performed to profile the expressions of 14 m5C regulators between HOSE and OVC. Consistent with the CNV pattern, most regulators were significantly up-regulated in OVC compared with HOSE tissues, whereas TET2 expression was decreased in two GEO cohorts (Figures 1F,G). Besides, a significant distinction of m5C regulators’ expression profiles among HOSE and OVC was illustrated by PCA (Figure 1E).
FIGURE 1. Genetic variation and gene expression of RNA m5C regulators in ovarian cancer. (A) Genetic alteration frequencies of m5C regulators in pan-cancer. (B) Copy number variation of m5C regulators in ovarian cancer. (C) Correlation of copy number variations among m5C regulators. (D) Principal component analysis for the expression profile of m5C regulators to distinguish OVC from HOSE samples in GSE54388 cohort. OVC, ovarian cancer, HOSE, human ovarian surface epithelium. (E) The differential expression analysis of 14 m5C regulators between OVC and HOSE samples in three independent GEO cohorts. Up-regulated in OVC: orange; up-regulated in HOSE samples: blue. (F) Correlation between copy number variation and mRNA level of m5C regulators. (G) The boxplot of expression of 14 m5C regulators in OVC and HOSE samples in GSE27651 cohort.*p < 0.05, **p < 0.01, ***p < 0.001.
The Interaction and Correlation Analysis of RNA m5C Regulators in Ovarian Cancer
To testify whether m5C regulators have correlations among each other, Spearman correlation analysis of mRNA levels of 14 regulators indicated significant positive correlation among most regulators (Figure 2A). The comprehensive landscape of m5C regulators was depicted with the PPI network according to the STRING 11.0 b website (Figure 2B). The writers and erasers had remarkable interactions within each other except the readers (YBX1 and ALYREF). To further investigate the correlations of writers and erasers who work as methyltransferases and demethyltransferases and affect the amount and distribution of m5C modification, comparisons of writer gene expression were performed in patients with high and low eraser gene expression (Figures 2C–J) (Zhang et al., 2020). The results showed that writer genes exhibited different correlations with eraser genes. NSUN6 expression was positively correlated with TET1/2/3, whereas NSUN7 expression was negatively correlated with eraser genes TET1/2 (Figures 2H,I). As writer genes NOP2, NSUN4, and NSUN2 have relatively higher amplification (Figure 1B), we analyzed whether the CNV of those writer genes is correlated with eraser genes. Of these, only TET3 was significantly down-regulated in patients with NSUN2 amplification compared with wild type (Figures 2K-M). These data indicate the complex cross-talk among m5C regulators in OVC.
FIGURE 2. The interaction and correlation analysis between RNA m5C regulators in ovarian cancer. (A) Correlation of mRNA levels among m5C regulators. (B) Protein–protein interaction plot among the 14 m5C regulators. (C–J) Correlations between “writers” (NUSN1, NSUN2, NSUN3, NSUN4, NSUN5, NSUN6, NSUN7, TRDMT1) and “erasers” (TET1, TET2, TET3, ALHBK1) at mRNA level (RPKM). (K–M) Differential mRNA levels of “erasers” between amplified type and wild type of “writers” with the highest CNV frequencies (NOP2, NSUN4, and NSUN2).
Consensus Clustering of m5C RNA Methylation Regulators Identifying Two Clusters With the Distinct OS
Based on the expression profile of m5C regulators in 374 OVC patients (TCGA), we used unsupervised consensus clustering analysis to distinguish different m5C modification patterns, and two clusters were identified (Figures 3A–C). k = 2 is the optimal stable clustering when compared with k = 3–8 (Supplementary Figure S1B). Distinct expression profiles of m5C RNA methylation regulators, clinicopathological parameters, and log2 (fold change) of regulators are illustrated in Figure 3D. According to χ2 test results, no statistical differences were found in lymphatic invasion (p = 1), venous invasion (p = 0.60), tumor residual (p = 0.28), International Federation of Gynecology and Obstetrics stage (p = 0.59), and grade (p = 0.13) between two clusters. TET2 was significantly up-regulated in cluster 1; ALYREF, NOP2, NSUN4, NSUN5, and YBX1 were substantially up-regulated in cluster 2. Despite similar clinicopathological parameters, the OS of patients in cluster 2 was better than cluster 1 (p = 0.015) (Figure 3E). We also examined the correlations of m5C regulators and clinicopathological parameters, and the results showed weak correlation between them (Supplementary Figure S2A).
FIGURE 3. Unsupervised consensus clustering analysis of RNA m5C regulators. (A) Consensus clustering matrix for the most suitable k (k = 2). (B) Relative changes in the area under the clustering cumulative function (CDF) curve at k = 2–10. (C) Consensus clustering CDF for k = 2–10. (D) The m5C regulator expression profiles and clinicopathological characteristics in two clusters. (E) Kaplan–Meier overall survival curves for ovarian cancer patients of two clusters in the TCGA cohort.
Functional Annotation Revealing the Metabolism Heterogeneity of the Two Subtypes
To further investigate the activated biological processes that may result in distinct survival in each cluster, functional annotations were conducted through three different methods. First, GSVA scores of GO and KEGG hallmark gene sets were calculated, and the differential analysis was performed to explore the discrepancies in pathways among two clusters, of which metabolism-related pathways occupied more than half of the top 30 enriched pathways (Figures 4A,B). Then, the top 500 expressed genes of each cluster and phenotype-related DEGs were selected for GO enrichment analysis. The results revealed that glycosaminoglycan and collagen metabolism–related pathways were up-regulated in cluster 1, whereas fatty acid metabolism–related pathways were up-regulated in cluster 2, which were further supported by the GSEA results (Figures 4C–F). To further explore the metabolism heterogeneity between the two clusters, the enrichment scores of 91 metabolic pathways downloaded from the KEGG database were calculated, and the varied metabolic pathways between the two clusters were analyzed. Differences were found in 27 of 91 (29.7%) metabolic pathways between two clusters. In addition, a correlation analysis was performed between those metabolic pathway enrichment scores and mRNA levels of m5C regulators. A significant negative correlation with these metabolism processes was found in eraser genes and positive correlation in writer genes except for NSUN6 and NOP2 (Figure 4G). Taken together, these results indicate that m5C modification might contribute to OVC survival difference through regulating metabolism heterogeneity.
FIGURE 4. Metabolic heterogeneity in ovarian cancers with distinct m5C modification patterns. (A,B) Bubble plot showing GSVA enrichment analysis of the top 30 changed GO and KEGG hallmark pathways between two m5C clusters. (C) The most enriched metabolism-related GO pathways using differentially expressed genes in each m5C cluster. (D) The most enriched metabolism-related GO pathways using the top 500 up-regulated genes in each m5C cluster. (E,F) GSEA analysis showing the significant metabolism-related pathways up-regulated in each m5C cluster. (G) Correlation analysis between GSVA enrichment scores of the changed metabolism-related pathways in each m5C cluster and mRNA levels of m5C regulators.
Representative Metabolism-Related Genes With the Prognostic Value in Different m5C Modification Patterns
An overlapping analysis was made to explore the leading metabolic genes that contributed to survival difference of two clusters (Figure 5A). 13 m5C cluster-specific metabolism genes with prognostic values were identified, as is illustrated in the Venn plot. According to the results of GO enrichment analysis using phenotype-related DEGs, GSEA, and GSVA, the collagen, glycosaminoglycan, and aminoglycan metabolism were up-regulated in cluster 1, whereas fatty acid metabolism was up-regulated in cluster 2. We visualized those genes in activated pathways of two m5C modification phenotype clusters (Figures 5B,C), and those important genes with prognostic value were marked in the circos plot. The fold change and p of 13 representative genes are shown in the volcano plot (Figure 5D). Kaplan–Meier analysis curves for OS in eight representative genes are presented (Figures 5E–L).
FIGURE 5. Selection of representative genes differentially expressed in m5C cluster, related to metabolism and with prognostic significance. (A) Venn plot showing the overlapping analysis of genes that were differentially expressed (fold change >1.5 and p < 0.05) in two clusters, enriched in up-regulated metabolic pathways in each cluster, and had the prognostic significance. (B,C) Circus plot exhibiting the metabolic pathways up-regulated in each m5C cluster and those 13 genes with prognostic value were marked. Blue, up-regulated in cluster 1. Orange, up-regulated in cluster 2. (D) Volcano plot distributes the relative expressions 13 representative genes with log2 (fold change) and log10 (FDR). The reference is gene expression of cluster 1. Blue, up-regulated in cluster 1. Orange, up-regulated in cluster 2. (E–L) Kaplan–Meier overall survival curves of eight representative genes.
Different Sensitivity to Chemotherapies Between Two Clusters
In order to detect whether the two clusters with m5C modification patterns had different drug sensitivity, we made the sensitivity prediction of 199 drugs for 374 OVC patients from TCGA dataset. Two clusters exhibited different sensitivities in a total of 44 drugs via t test for imputed sensitivity score, including Paclitaxel_1,080, Docetaxel_1,007, and Cisplatin_1,005, which were included in the standard chemotherapy of OVC. Higher sensitivities of these three chemotherapeutic drugs were observed in cluster 2 (Figure 6).
FIGURE 6. The boxplot of imputed sensitivity score of three chemotherapeutic drugs in two clusters.
Establishment and Evaluation of a Risk Score Signature With RNA m5C Methylation Regulators
As clusters with different m5C modification patterns exhibited altered metabolism pathways and chemotherapeutic drug sensitivity that were associated with the differences in prognosis, we further explored the predictive value of RNA m5C methylation regulators for OS. The univariate survival analyses using Cox proportional hazards models were performed, and three genes with p < 0.1 were selected for the LASSO Cox algorithm (Supplementary Figure S2B). Then, the LASSO algorithm with 10 folds of cross-validation was applied to establish a risk score prediction model for OS. Two hundred eighty of 374 OVC patients in TCGA cohorts formed the internal training set, and the rest of 94 patients formed the internal testing set. The GSE19829 formed the external testing set. ALYREF, NOP2, and TET2 finally entered in the model (Figures 7A,B), and the risk core model is (ALYREF*-0.00349425170618238) + (NOP2*0.0200758391547147) + (TET2*0.145822477408262). Patients were then divided into subgroups of low score and high score according to the cutoff value. The distribution of risk scores in the internal training set and internal testing set is demonstrated in Figure 7C and Supplementary Figure S2C, and the results suggest the risk score model could distinguish those patients with poor survivals. Patients with high score exhibit poor OS, which were verified in internal training set (p = 0.016, hazard ratio [HR] = 3.2), internal testing set (p = 0.012, HR = 10), and external testing set (p = 0.02, HR = 2.1) (Figures 7D–F). The AUCs of ROC curves of the prediction were approximately 0.6 (Supplementary Figures S2D,E). To confirm the role of key genes in the risk score signature, in vitro experiments were performed. The results showed that knockdown of TET2 in OVC cells could hinder the cell malignant growth, suggesting its potential in the prediction of OS prognosis (Supplementary Figure S2F). Furthermore, multivariate analysis of Cox proportional hazards was performed to further confirm the performance of m5C risk score prediction for OS, and the results showed that the m5C risk score was an independent prognostic factor (p = 0.014, HR = 38.34) (Figure 7G).
FIGURE 7. Construction and evaluation of prognostic prediction model with three RNA m5C regulators in TCGA cohort and GEO cohort (GSE19829). (A,B) LASSO Cox regression analysis results showing the identification of three prognostic risk signature genes, and the risk score model is (ALYREF*-0.00349425170618238) + (NOP2*0.0200758391547147) + (TET2*0.145822477408262). (C) The distribution of prognostic signature-based risk score in internal training set. (D–F) The Kaplan–Meier overall survival analysis for patients with high score and low score in internal training set, internal testing set, and external testing set. (G) The m5C risk score in Cox multivariate analysis for OS of ovarian cancer patients.
Discussion
In the present study, the metabolism heterogeneity in two clusters based on m5C expression profile was annotated through enrichment analyses and found to be significantly correlated with the mRNA levels of m5C regulators. Then, 13 metabolism-related DEGs of two m5C clusters related RNA m5C methylation were identified, which were also associated with OS. Besides metabolism reprogramming, RNA m5C regulators could also trigger altered chemotherapeutic drug sensitivity and consequently influence survival in OVC. Ultimately, a prognostic model comprising of ALYREF, NOP2, and TET2 for OS was developed for further verifying predictive value of m5C regulators for prognosis in OVC.
As the present diagnosis and screening program in OVC still have limitations, efforts have been made to seek predictive biomarkers related to cancer occurrence and progression. Through these biomarkers, the distinction of subtypes with different prognoses and molecular characteristics will make the identification of patient subgroups who could well respond to a certain treatment or who had worse survival. In the present study, incorporating mRNA data of ALYREF, NOP2, and TET2 could well stratify those patients with worse survivals, providing evidence for clinical practice. The predictive value of m5C regulators has also been confirmed in glioma (Wang et al., 2020), breast cancer (Huang et al., 2021), and head and neck squamous cell carcinoma (Xue et al., 2020), and our study made extending support for their predictive role in OVC.
Metabolism reprogramming in OVC has been implicated in the pathogenesis, progression, and target therapy for cancer. In our study, glycosaminoglycan and collagen metabolism–related pathways were activated in cluster 1 and fatty acid metabolism–related pathways were activated in cluster 2. Fatty acids are important components of lipids such as fats, sterol esters, and phospholipids. Lipid metabolism dysregulation has been verified to participate in cancer progression (Li et al., 2021), but it has a more special implication in OVC, as an almost symbiotic relationship exists between OVC and the fat-containing cells in the omentum. FABP4, ELOVL2, and ACSL5 as important genes involved in biosynthesis and transport of fatty acids were found to be up-regulated in cluster 2. FABP4 was highly expressed on the membrane of metastasis OVC cells at the adipocyte–cancer cell interface and mediated lipid accumulation and effect on invasion (Zhao et al., 2019). ELOVL2 has been reported to participate in the biosynthesis of polyunsaturated fatty acids and be involved in tumorigenicity in glioma cancer stem cells (Gimple et al., 2019). ACSL family is responsible for activating long-chain fatty acids, and family members have opposite functions toward carcinogenesis. ACSL5 is nuclear-coded and expressed in the mitochondria and physiologically participates in the proapoptotic sensing of cells acting as a tumor suppressor, which could possibly explain the relatively better prognosis of cluster 2 where ACSL5 was dominantly up-regulated (Quan et al., 2021). The collagen metabolism alteration influences the distribution of collagen and extracellular matrix (ECM) structure, thus affecting cancer progression (Xu et al., 2019). The glycosaminoglycans, another important component of ECM, were found to be involved in multiple signaling pathways related to angiogenesis, cancer invasion, and metastasis (Morla, 2019). VCAN, which is one member of glycosaminoglycan gene sets with leading expression in cluster 1, was previously reported to be up-regulated in ovarian stromal cells and associated with increased microvessel density and poorer survival (Ghosh et al., 2010). In this study, we found the heterogeneity of lipid, glycosaminoglycans, and collagen metabolism in two clusters with distinct m5C modification. The stepwise accumulation of altered metabolism at mRNA levels in different m5C clusters eventually resulted in distinct prognoses, indicating that the metabolism alteration has prognostic significance.
The role of epigenetic modifications in cancer metabolism reprogramming has been broadly reported, but there is still a lack of disclosure of how RNA m5C modification functions in cancer metabolism. It was demonstrated that ALYREF binds the 3′-UTR of PKM2 mRNA and promotes the glucose metabolism of bladder cancer in an m5C-dependent manner (Wang et al., 2021). In this study, we found 13 representative metabolic genes that were related to m5C RNA methylation in OVC. All of them have positive or negative correlations with RNA m5C regulators at the mRNA level (Supplementary Figure S2G). Consequently, experimental verification could be done in the future to verify the regulatory role of RNA m5C methylation of those metabolism genes in cancer.
In conclusion, our study depicted the landscape of genetic variation and gene expression of m5C regulators in OVC and established a prognostic prediction model formed by ALYRER, NOP2, and TET2 for OS. We also uncovered the indispensable roles of m5C modification in metabolism heterogeneity and altered sensitivity to chemotherapeutic drugs.
Data Availability Statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.
Author Contributions
PY put forward the ideas of this article. PY, TL, and JX drafted and reviewed the article. JX performed the bio-informatics analysis. YC, YW, and XL helped with acquisition of data and analysis and interpretation of data. All authors read and approved the final manuscript.
Funding
This work was sponsored by the National Natural Science Foundation of China (82072886), the Natural Science Foundation of Chongqing (cstc2020jcyj-msxmX0344) and the Natural Science Foundation of Chongqing (cstc2021jcyj-bsh0022).
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.
Acknowledgments
The authors acknowledge contributions from TCGA databases.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2022.807786/full#supplementary-material
References
Amort, T., Rieder, D., Wille, A., Khokhlova-Cubberley, D., Riml, C., Trixl, L., et al. (2017). Distinct 5-methylcytosine Profiles in Poly(A) RNA from Mouse Embryonic Stem Cells and Brain. Genome Biol. 18 (1), 1. doi:10.1186/s13059-016-1139-1
Blanco, S., Bandiera, R., Popis, M., Hussain, S., Lombard, P., Aleksic, J., et al. (2016). Stem Cell Function and Stress Response Are Controlled by Protein Synthesis. Nature 534 (7607), 335–340. doi:10.1038/nature18282
Chai, R. C., Li, Y. M., Zhang, K. N., Chang, Y. Z., Liu, Y. Q., Zhao, Z., et al. (2019a). RNA Processing Genes Characterize RNA Splicing and Further Stratify Lower-Grade Glioma. JCI Insight 5, e130591. doi:10.1172/jci.insight.130591
Chai, R. C., Wu, F., Wang, Q. X., Zhang, S., Zhang, K. N., Liu, Y. Q., et al. (2019b). m6A RNA Methylation Regulators Contribute to Malignant Progression and Have Clinical Prognostic Impact in Gliomas. Aging 11, 1204–1225. doi:10.18632/aging.101829
Chen, X., Li, A., Sun, B.-F., Yang, Y., Han, Y.-N., Yuan, X., et al. (2019). 5-methylcytosine Promotes Pathogenesis of Bladder Cancer through Stabilizing mRNAs. Nat. Cel Biol 21 (8), 978–990. doi:10.1038/s41556-019-0361-y
Chen, Y. S., Yang, W. L., Zhao, Y. L., and Yang, Y. G. (2021). Dynamic Transcriptomic M5 C and its Regulatory Role in RNA Processing. Wiley Interdiscip. Rev. RNA 12 (4), e1639. doi:10.1002/wrna.1639
Cheng, J. X., Chen, L., Li, Y., Cloe, A., Yue, M., Wei, J., et al. (2018). RNA Cytosine Methylation and Methyltransferases Mediate Chromatin Organization and 5-azacytidine Response and Resistance in Leukaemia. Nat. Commun. 9 (1), 1163. doi:10.1038/s41467-018-03513-4
Clough, E., and Barrett, T. (2016). The Gene Expression Omnibus Database. Methods Mol. Biol. 1418, 93–110. doi:10.1007/978-1-4939-3578-9_5
Freeman, J. W., McGrath, P., Bondada, V., Selliah, N., Ownby, H., Maloney, T., et al. (1991). Prognostic Significance of Proliferation Associated Nucleolar Antigen P120 in Human Breast Carcinoma. Cancer Res. 51 (8), 1973–1978.
Gao, J., Aksoy, B. A., Dogrusoz, U., Dresdner, G., Gross, B., Sumer, S. O., et al. (2013). Integrative Analysis of Complex Cancer Genomics and Clinical Profiles Using the cBioPortal. Sci. Signal. 6 (269), pl1. doi:10.1126/scisignal.2004088
Ghosh, S., Albitar, L., LeBaron, R., Welch, W. R., Samimi, G., Birrer, M. J., et al. (2010). Up-regulation of Stromal Versican Expression in Advanced Stage Serous Ovarian Cancer. Gynecol. Oncol. 119 (1), 114–120. doi:10.1016/j.ygyno.2010.05.029
Gimple, R. C., Kidwell, R. L., Kim, L. J. Y., Sun, T., Gromovsky, A. D., Wu, Q., et al. (2019). Glioma Stem Cell-specific Superenhancer Promotes Polyunsaturated Fatty-Acid Synthesis to Support EGFR Signaling. Cancer Discov. 9 (9), 1248–1267. doi:10.1158/2159-8290.cd-19-0061
Heagerty, P. J., and Zheng, Y. (2005). Survival Model Predictive Accuracy and ROC Curves. Biometrics 61 (1), 92–105. doi:10.1111/j.0006-341x.2005.030814.x
Huang, Z., Pan, J., Wang, H., Du, X., Xu, Y., Wang, Z., et al. (2021). Prognostic Significance and Tumor Immune Microenvironment Heterogenicity of m5C RNA Methylation Regulators in Triple-Negative Breast Cancer. Front. Cel Dev. Biol. 9, 657547. doi:10.3389/fcell.2021.657547
Jemal, A., Murray, T., Samuels, A., Ghafoor, A., Ward, E., and Thun, M. J. (2003). Cancer Statistics, 2003. CA Cancer J. Clin. 53, 5–26. doi:10.3322/canjclin.53.1.5
Kaufman, B., Shapira-Frommer, R., Schmutzler, R. K., Audeh, M. W., Friedlander, M., Balmaña, J., et al. (2015). Olaparib Monotherapy in Patients with Advanced Cancer and a Germline BRCA1/2 Mutation. J. Clin. Oncol. 33 (3), 244–250. doi:10.1200/JCO.2014.56.2728
Kim, H., Xu, H., George, E., Hallberg, D., Kumar, S., Jagannathan, V., et al. (2020). Combining PARP with ATR Inhibition Overcomes PARP Inhibitor and Platinum Resistance in Ovarian Cancer Models. Nat. Commun. 11 (1), 3726. doi:10.1038/s41467-020-17127-2
Kristeleit, R., Shapiro, G. I., Burris, H. A., Oza, A. M., LoRusso, P., Patel, M. R., et al. (2017). A Phase I-II Study of the Oral PARP Inhibitor Rucaparib in Patients with Germline BRCA1/2-Mutated Ovarian Carcinoma or Other Solid Tumors. Clin. Cancer Res. 23 (15), 4095–4106. doi:10.1158/1078-0432.ccr-16-2796
Lheureux, S., Braunstein, M., and Oza, A. M. (2019). Epithelial Ovarian Cancer: Evolution of Management in the Era of Precision Medicine. CA Cancer J Clin. 69, 280. doi:10.3322/caac.21559
Li, L-Y., Yang, Q., Jiang, Y-Y., Yang, W., Jiang, Y., Li, X., et al. (2021). Interplay and Cooperation between SREBF1 and Master Transcription Factors Regulate Lipid Metabolism and Tumor-Promoting Pathways in Squamous Cancer. Nat. Commun. 12 (1), 4362. doi:10.1038/s41467-021-24656-x
Maeser, D., Gruener, R. F., and Huang, R. S. (2021). oncoPredict: an R Package for Predicting In Vivo or Cancer Patient Drug Response and Biomarkers from Cell Line Screening Data. Brief Bioinform 22. doi:10.1093/bib/bbab260
Moore, K., Colombo, N., Scambia, G., Kim, B.-G., Oaknin, A., Friedlander, M., et al. (2018). Maintenance Olaparib in Patients with Newly Diagnosed Advanced Ovarian Cancer. N. Engl. J. Med. 379 (26), 2495–2505. doi:10.1056/nejmoa1810858
Morla, S. (2019). Glycosaminoglycans and Glycosaminoglycan Mimetics in Cancer and Inflammation. Int. J. Mol. Sci. 20 (8), 1963. doi:10.3390/ijms20081963
Nombela, P., Miguel-López, B., and Blanco, S. (2021). The Role of m6A, m5C and Ψ RNA Modifications in Cancer: Novel Therapeutic Opportunities. Mol. Cancer 20 (1), 18. doi:10.1186/s12943-020-01263-w
Pujade-Lauraine, E., Ledermann, J. A., Selle, F., Gebski, V., Penson, R. T., Oza, A. M., et al. (2017). Olaparib Tablets as Maintenance Therapy in Patients with Platinum-Sensitive, Relapsed Ovarian Cancer and a BRCA1/2 Mutation (SOLO2/ENGOT-Ov21): a Double-Blind, Randomised, Placebo-Controlled, Phase 3 Trial. Lancet Oncol. 18 (9), 1274–1284. doi:10.1016/S1470-2045(17)30469-2
Pujade-Lauraine, E. (2017). New Treatments in Ovarian Cancer. Ann. Oncol. 28 (Suppl. l_8), viii57–viii60. doi:10.1093/annonc/mdx442
Quan, J., Bode, A. M., and Luo, X. (2021). ACSL Family: The Regulatory Mechanisms and Therapeutic Implications in Cancer. Eur. J. Pharmacol. 909, 174397. doi:10.1016/j.ejphar.2021.174397
Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res. 43 (7), e47. doi:10.1093/nar/gkv007
Roundtree, I. A., Evans, M. E., Pan, T., and He, C. (2017). Dynamic RNA Modifications in Gene Expression Regulation. Cell 169 (7), 1187–1200. doi:10.1016/j.cell.2017.05.045
Scrucca, L., Santucci, A., and Aversa, F. (2007). Competing Risk Analysis Using R: an Easy Guide for Clinicians. Bone Marrow Transpl. 40 (4), 381–387. doi:10.1038/sj.bmt.1705727
Siegel, R. L., Miller, K. D., Fuchs, H. E., and Jemal, A. (2021). Cancer Statistics, 2021. CA Cancer J. Clinicians 71 (1), 7, doi:10.3322/caac.21654
Singer, G., Oldt, R., Cohen, Y., Wang, B. G., Sidransky, D., Kurman, R. J., et al. (2003). Mutations in BRAF and KRAS Characterize the Development of Low-Grade Ovarian Serous Carcinoma. JNCI J. Natl. Cancer Inst. 95 (6), 484–486. doi:10.1093/jnci/95.6.484
Tang, Y., Gao, C.-C., Gao, Y., Yang, Y., Shi, B., Yu, J.-L., et al. (2020). OsNSUN2-Mediated 5-Methylcytosine mRNA Modification Enhances Rice Adaptation to High Temperature. Dev. Cel 53 (3), 272–286. doi:10.1016/j.devcel.2020.03.009
Wang, J. Z., Zhu, W., Han, J., Yang, X., Zhou, R., Lu, H. C., et al. (2021). The Role of the HIF‐1α/ALYREF/PKM2 axis in Glycolysis and Tumorigenesis of Bladder Cancer. Cancer Commun. 41 (7), 560–575. doi:10.1002/cac2.12158
Wang, P., Wu, M., Tu, Z., Tao, C., Hu, Q., Li, K., et al. (2020). Identification of RNA: 5-Methylcytosine Methyltransferases-Related Signature for Predicting Prognosis in Glioma. Front. Oncol. 10, 1119. doi:10.3389/fonc.2020.01119
Wilkerson, M. D., and Hayes, D. N. (2010). ConsensusClusterPlus: a Class Discovery Tool with Confidence Assessments and Item Tracking. Bioinformatics 26 (12), 1572–1573. doi:10.1093/bioinformatics/btq170
Xue, M., Shi, Q., Zheng, L., Li, Q., Yang, L., and Zhang, Y. (2020). Gene Signatures of m5C Regulators May Predict Prognoses of Patients with Head and Neck Squamous Cell Carcinoma. Am. J. Transl Res. 12, 6841. ISSN:1943-8141/AJTR0116260.
Xu, S., Xu, H., Wang, W., Li, S., Li, H., Li, T., et al. (2019). The Role of Collagen in Cancer: from Bench to Bedside. J. Transl Med. 17 (1), 309. doi:10.1186/s12967-019-2058-1
Yang, L., Perrera, V., Saplaoura, E., Apelt, F., Bahin, M., Kramdi, A., et al. (2019). m5C Methylation Guides Systemic Transport of Messenger RNA over Graft Junctions in PlantsC Methylation Guides Systemic Transport of Messenger RNA over Graft Junctions in Plants. Curr. Biol. 29 (15), 2465–2476. doi:10.1016/j.cub.2019.06.042
Yang, X., Yang, Y., Sun, B.-F., Chen, Y.-S., Xu, J.-W., Lai, W.-Y., et al. (2017). 5-methylcytosine Promotes mRNA export - NSUN2 as the Methyltransferase and ALYREF as an m5C Reader. Cell Res 27 (5), 606–625. doi:10.1038/cr.2017.55
Yu, G., Wang, L.-G., Han, Y., and He, Q.-Y. (2012). clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters. OMICS: A J. Integr. Biol. 16 (5), 284–287. doi:10.1089/omi.2011.0118
Zhang, B., Wu, Q., Li, B., Wang, D., Wang, L., and Zhou, Y. L. (2020). m6A Regulator-Mediated Methylation Modification Patterns and Tumor Microenvironment Infiltration Characterization in Gastric cancerA Regulator-Mediated Methylation Modification Patterns and Tumor Microenvironment Infiltration Characterization in Gastric Cancer. Mol. Cancer 19 (1), 53. doi:10.1186/s12943-020-01170-0
Keywords: ovarian cancer, 5-methylcytosine, RNA modification, metabolism heterogeneity, LASSO cox regression
Citation: Xu J, Liu X, Chen Y, Wang Y, Liu T and Yi P (2022) RNA 5-Methylcytosine Regulators Contribute to Metabolism Heterogeneity and Predict Prognosis in Ovarian Cancer. Front. Cell Dev. Biol. 10:807786. doi: 10.3389/fcell.2022.807786
Received: 02 November 2021; Accepted: 08 February 2022;
Published: 18 March 2022.
Edited by:
Xiaoxing Li, The First Affiliated Hospital of Sun Yat-sen University, ChinaReviewed by:
Yubin Xie, Sun Yat-sen University, ChinaQian Yang, Cedars Sinai Medical Center, United States
Copyright © 2022 Xu, Liu, Chen, Wang, Liu and Yi. 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: Ping Yi, yiping@cqmu.edu.cn
†These authors have contributed equally to this work