Comprehensive Analysis of 5-Methylcytosine Profiles of Messenger RNA in Human High-Grade Serous Ovarian Cancer by MeRIP Sequencing.

BACKGROUND: Accumulating evidence has indicated that methylation status is closely related to tumourigenesis and a few aggressive features of diverse cancers. However, as an important methylation regulation modification, the distribution of 5-methylcytosine (m5C) in high-grade serous ovarian cancer (HGSOC) remains unclear.
MATERIALS AND METHODS: We collected three pairs of human HGSOC tissues and adjacent non-tumour tissues to analyse the transcriptome-wide m5C methylation of messenger RNAs (mRNAs) by methylated RNA immunoprecipitation sequencing. Gene ontology (GO) enrichment analysis and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway analysis were performed to eExplore the potential biological functions of these genes and important cancer pathways. We used immunohistochemistry to analyse the expression of the m5C modification regulatory gene MAP2K3 in 80 HGSOC tissue samples, and their associations with clinical parameters were analyzed using the Spearman-rho test. Univariate and multivariate Cox regression analyses were performed to identify potential prognostic factors. Kaplan-Meier analysis was performed to analyze overall survival.
RESULTS: We identified 2050 dysregulated m5C peaks, 1767 of which were significantly upregulated, while 283 were significantly downregulated. GO enrichment analysis showed that genes altered by the m5C peak played a key role in system development, transporter complex, and transporter activity. KEGG pathway analysis revealed that these genes were enriched in some important pathways in cancer regulation, such as inflammatory mediator regulation of the TRP channels pathway, Wnt signalling pathway, and focal adhesion pathways. In addition, through joint analysis of MeRIP-seq and RNA-seq data, we identified 125 differentially methylated m5C peaks and synchronous differentially expressed genes. These genes play key roles in cell growth, maintenance, plasma membranes, and cell adhesion molecule activity. Immunohistochemical staining results showed that high expression of MAP2K3 was significantly correlated with CA125 level (p < 0.001), tumour size (p = 0.001), lymph node metastasis (p = 0.008), depth of myometrial invasion (p < 0.001), and FIGO stage (p < 0.001), indicating a poor prognosis.
CONCLUSION: Our results reveal the different distribution patterns of m5C in HGSOC and adjacent tissues and the possible involvement of m5C in HGSOC cell functions. Our study provides new insights into the epi-transcriptomic dysregulation of m5C in the tumourigenesis of HGSOC.

Introduction

As the second most common gynaecological malignancy, more than 80% of patients with ovarian cancer are diagnosed at an advanced stage of extensive peritoneal metastasis.1,2 High-grade serous ovarian cancer (HGSOC) is the most prevalent histotype of ovarian cancer, constituting more than 70% of ovarian cancer-related deaths.3 Despite the improvement in treatment for patients with HGSOC, the 5-year survival rate of patients diagnosed at the FIGO III or IV stage is less than 40%, making this disease the second most lethal gynaecological malignancy.4 Therefore, exploring novel and effective diagnostic and prognostic markers is urgently needed.

Increasing evidence has shown that post-transcriptional RNA modifications may serve as novel gene expression regulators, the disruption of which may lead to life-threatening diseases, including cancer.5,6 RNA epigenetic markers mainly include N6-methyladenosine (m6A), 5-methylcytosine (m5C), and pseudouridine (Ψ). RNA m5C, characterised by the addition of a methyl group at the carbon-5 position of the cytosine base, has been detected not only in rRNA and tRNA but also in mRNA by high-throughput sequencing-based transcriptome-wide mapping approaches.7 m5C is well-conserved and dynamically regulated in cellular mRNAs, and it localises to the untranslated regions (UTRs) of mRNA transcripts. It may influence a wide variety of biological functions by regulating RNA metabolism.8 Moreover, and RNA m5C modification and its regulators have been confirmed to play crucial roles in the pathogenesis of bladder cancer,9 hepatocellular carcinoma,10 glioblastoma multiforme,11 and leukaemia,12 indicating the promising prospect of m5C modification in cancer treatment. However, the quantity, distribution, and functions of m5C modification in mRNA in HGSOC remain unclear.

As an emerging mechanism in gene regulation, reversible m5C modification is regulated by “writers” (methyltransferases), “readers,” and “erasers” (demethylase) and affects various biological functions, such as cell function, development, immune system, and homeostasis.13 Huang et al discovered 11 m5C regulatory genes related to TNBC neoplasia in triple negative breast cancer (TNBC), which are related to the poor prognosis of NSUN2 and NSUN6 and TNBC.14 Xue et al study in head and neck squamous cell carcinoma (HNSCC) showed that the m5C regulatory genes MAPK1, HRAS, RAC1, HIST1H4K and HIST1H3C are potential driving genes of HNSCC and are potentially related to the occurrence of HNSCC.15 In addition, the demethylation of MAP2K3 is related to inflammation and immune regulation in HNSCC patients.16 However, the gene signatures and prognostic values of m5C-related regulators in HGSOC remain unknown.

In the present study, we performed m5C-seq, RNA-seq, and in-depth bioinformatics analyses to clarify the role of m5C in the pathogenesis of HGSOC. Through comprehensive analysis of MeRIP-seq and RNA-seq data, significant differences in the number and distribution of m5C peaks were noted in human HGSOC and matched adjacent normal tissues. Bioinformatics analysis showed that genes altered by the m5C peak played a key role in cell function and were involved in important cancer pathways. Furthermore, we used immunohistochemical methods to confirm that the m5C mRNA-regulated gene MAP2K3 was significantly associated with poor prognosis in HGSOC. Our findings provide a novel viewpoint that the m5C pathway is crucial for the pathogenesis of HGSOC.

Materials and Methods

Ethical Approval

The study has been approved by the Ethics Committee of the Institutional Review Board of The Second Hospital of Hebei Medical University, and the patient’s written informed consent was obtained before the study began. And the study was conducted in accordance with the Declaration of Helsinki.

Sample Collection and RNA Preparation

Each pair of HGSOC and corresponding adjacent tissues was obtained from the same patient with HGSOC. In total, three pairs of HCC and adjacent non-tumor tissue samples were obtained. Total RNA from each sample was isolated using TRIzol reagent (Invitrogen Corporation, Carlsbad, CA). The Ribo-Zero rRNA Removal Kit (Illumina, Inc., CA, USA) was then used to reduce the content of rRNA in the total RNA. A NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) was used to measure the RNA concentration of each sample, and RNA with OD260/OD280 between 1.8 and 2.1 was marked as qualified.

RNA MeRIP-Seq Construction and Sequencing

RNA MeRIP-seq was performed based on previously reported procedures with some modifications.17 In brief, the GenSeqTM m5C RNA IP kit (China GenSeq Inc.) was used for m5C RNA immunoprecipitation according to the manufacturer’s instructions. The cDNA library was constructed using the NEBNext Ultra II RNA Library Prep Kit for Illumina (New England Biolabs, Inc.). The Agilent Bioanalyzer 2100 system (Agilent Technologies, Inc., USA) was used to evaluate the quality and purity of the cDNA library. The cDNA libraries were sequenced on a HiSeq 6000 platform (Illumina, San Diego, CA, USA) with 100 bp paired-end reads.

Sequencing Data Analysis

Briefly, quality control of the paired-end reads was performed with Q30 and reads higher than 30 were considered reliable mapping and retained. Cutadapt software (v1.9.3) was used to trim the 3′ adaptors and remove low-quality reads to harvest high-quality clean reads.18 Then, clean reads of input libraries were aligned to the reference genome (GRCh38.gencode. v32) with STAR software19 and mRNA peaks were identified using DCC software.20 Following this, clean reads of all libraries were aligned to the reference genome (MM10) using hisat2 software (v2.0.4),21 followed by identifying methylated m5C peaks on mRNAs with Model-based Analysis of ChIP-Seq (MACS) software.22 Differentially methylated m5C peaks were identified with a fold change of >2 or <0.5 (p < 0.00001) using DiffReps software.23 The MeRIP-enriched regions (peaks) identified using both software and the section of m5C that overlapped with the exon of the protein-coding genes were further annotated by scripts developed in-house.

Statistical Analysis of the Distribution Characteristics of m5C Peaks

The m5C peaks on the mRNA of the three samples in the HGSOC group were combined to obtain the total m5C peaks, and the adjacent tissue group was treated in the same way. The sequences of the vertices of methylated peaks, which were 50 bp on each side of the apex, were scanned using Dreme software24 to find reliable motif sequences. Subsequently, we counted the mRNA region where the m5C peak was located in each sample according to published methods25 and plotted the results as a Pie chart. MetaPlotR software26 was used to draw the distribution map of the methylation peaks on the metagene for comparison between the two groups. In addition, heatmap.2 software package was used to perform clustering based on the fold enrichment (FE) of each mRNA in six samples. The color represents the size of the logFE value. The closer the color in the figure is to red, the larger the logFE value. And volcano plots were used to visualise different m5C peaks based on the logFE value of each gene in the six samples. Red and blue represent mRNA that is up-regulated and down-regulated, respectively.

Transcriptome Sequencing Analysis and Conjoint Analysis

High-quality reads were mapped onto the genome (human GENCODE v32) with hisat2 software (v2.0.4).21 Subsequently, HTSeq software (v0.9.1) was used to obtain gene level raw counts as mRNA expression profiles.27 Edge R software (v3.16.5) was used to normalise the data and screen differentially expressed mRNAs (RNA-seq data).28 Venn diagrams were performed to screen differentially regulated methylation-regulated annotation genes and genes that overlap with RNA-seq data. FunRich software (v3.1.3) was used to perform enrichment analysis of the overlapping genes.29

Bioinformatics Analysis

Gene ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed for annotation, visualisation, and integrated discovery using the DAVID database to reveal the functions of these differentially methylated genes.30 The ontology covers three parts: cellular molecular function (MF), biological process (BP), and component (CC). KEGG pathway analyses were used to predict the functions of differentially methylated and expressed genes. A p-value <0.05 was used as the cut-off criterion.

Immunofluorescence Staining

For immunohistochemical staining, 80 paraffin-embedded HGSOC tissue samples were collected from patients at the Second Hospital of Hebei Medical University (Shijiazhuang, China). Anti-rabbit MAP2K3 (1:500, Bioss Antibodies, USA) was used as the primary antibody. IHC staining was performed using a previously described protocol.31,32 The staining intensity was evaluated as follows: 0, no staining; 1, beige; 2, darker beige; and 3, tan. The staining extent was scored as follows: 0, 0% stained; 1, 1% to 25% stained; 2, 26% to 50% stained; and 3, 51% to 100% stained. The final scores were calculated by multiplying intensity score and extent score and judged as follows: 0 to 3, low staining (-/+) and 4 to 9, high staining (++/+++).

Statistical Analysis

The association between the expression of MAP2K3 and various clinicopathological parameters was analysed using the Spearman’s Rho test. The Kaplan-Meier method was used to generate survival curves. All statistical analyses were performed using SPSS 25.0, and significance was set at p < 0.05.

Results

General Features of m5C Methylation Map in HGSOC

In total, we found 20,802 m5C mRNA peaks in normal tissues and 13,856 m5C mRNA peaks in HGSOC. We mapped up to 9969 annotated genes in normal tissues and 7089 annotated genes in HGSOC tissues. Among them, 6291 m5C mRNA peaks appeared in both normal and HGSOC tissues, corresponding to 5514 annotated genes (Figure 1A and B). Circos software was used to further investigate the distribution of m5C peaks on chromosomes and found that the quantity and distribution of m5C peaks on each chromosome were different between HGSOC and adjacent normal tissues, with the most obvious difference on the first chromosome (Figure 1C).

Figure 1

Characteristics of m5C peaks in normal and HGSOC tissues. Venn diagram of (A) m5C sites identified in mRNA and (B) m5C targeted genes from normal and HGSOC tissues. (C) Visualisation of m5C sites at the chromosome level in normal and HGSOC tissues. (D) Accumulation of the region of average m5C peaks along all transcripts in normal and HGSOC tissues. (E and F) Pie chart of the source of methylated mRNA in normal and HGSOC tissues.

In addition, the analysis of the source of methylation peaks showed that m5C was distributed in all regions of mRNA and was mainly enriched in the coding sequence (Figure 1D–F). However, m5C peaks in HGSOC tissues showed a distinct pattern from peaks in normal tissues with a relative decrease in the number of m5C peaks in the coding sequence (CDS) (HGSOC: 54.74%, adjacent tissues: 55.82%) and 3′ untranslated region (3′UTR) (HGSOC: 6.39%, adjacent tissues: 8.04%), and a relative increase in the 5′ untranslated region (5′UTR) (HGSOC: 8.44%, adjacent tissues: 8.25%), start codon (startC; HGSOC: 12.13%, adjacent tissues: 11.6%), and stop codon (stopC) (HGSOC: 18.31%, adjacent tissues: 16.28%).

Motif Analysis and Differential m5C Methylation Analysis

We used the Dreme software to identify meaningful motif sequences. The results showed that CCAGSCUR (G = FFB300) was the most reliable methylation site motif (p = 7.1e-130) in normal tissues, whereas CAGSCUGG was the most reliable methylation site motif (p = 4.6e-067) in HGSOC (Figure 2A). In addition, we performed a methylation heatmap analysis and cluster analysis of the total data. The results of the cluster analysis showed the apparent expression difference between the groups and relatively consistent within group. And these differences may be closely related to the occurrence and development of HGSOC (Figure 2B). Furthermore, we used DiffReps software to identify 1474 upregulated m5C methylation peaks and 336 downregulated m5C methylation peaks (Figure 2C). We have listed the 20 mRNAs with the largest FC in Tables 1 and 2. These dysregulated m5C methylation peaks were distributed on all chromosomes and were mainly enriched on chromosome one (Figure 2D). Moreover, to determine the number of m5C peaks on each mRNA, we performed a statistical analysis of the methylation peaks and the corresponding mRNA. The results showed that most of the mRNAs with methylation sites in HGSOC tissues had only one methylation peak (65.65%), and this ratio was lower in adjacent tissues (68.98%). And the difference between the two groups has statistically significant (p < 0.05). At the same time, the number of mRNAs with three or more methylation peaks on one mRNA was higher in HGSOC tissues than that in adjacent normal tissues (p < 0.05) (Figure 2E).

Figure 2

Characteristics of altered m5C peaks in normal and HGSOC tissues. (A) The sequence motif of m5C sites in normal and HGSOC tissues. (B) Cluster analysis of methylation in normal and HGSOC tissues. (C) Volcano plots showing the significantly differential m5C peaks. (D) The distributions of altered m5C peaks in human chromosomes. (E) The number of m5C peaks in normal and HGSOC tissues on each mRNA. Most mRNAs have only one methylation peak.

Table 1

Top Ten Upregulated m5C Methylation Peaks

Chrom	PeakID	GeneName	Foldchange	P-value
chr10	diffreps_peak_289928	SLC18A2	43.5	6.10922E-06
chr19	diffreps_peak_1056117	SPTBN4	7.202531646	9.83191E-06
chr1	diffreps_peak_190032	RYR2	3.86307311	1.22851E-06
chr3	diffreps_peak_1559889	SLITRK3	3.60978836	9.87429E-06
chr2	diffreps_peak_1244125	UNC80	8.174528302	2.64657E-07
chr19	diffreps_peak_998952	ARID3A	3.949352179	7.32351E-09
chr6	diffreps_peak_1949237	TBXT	5.078212291	1.46178E-06
chr19	diffreps_peak_1022279	ZNF333	4.58560794	1.91991E-07
chr10	diffreps_peak_207329	PRKCQ	7.333333333	6.07173E-06
chr19	diffreps_peak_1010488	VAV1	5.416455696	3.09119E-07

Table 2

Top Ten Downregulated m5C Methylation Peaks

Chrom	PeakID	GeneName	Foldchange	P-value
chr22	diffreps_peak_1405002	LRP5L	2.443803461	7.36179E-06
chr7	diffreps_peak_2083662	RBM33	2.20000869	2.48523E-06
chr19	diffreps_peak_1066761	ELSPBP1	9.132893497	2.26821E-07
chr20	diffreps_peak_1284252	PAK5	56.65506329	6.48132E-08
chr22	diffreps_peak_1394029	MICAL3	8.106723319	7.87452E-07
chr17	diffreps_peak_913277	ITGB4	5.617481203	8.24225E-07
chr5	diffreps_peak_1806753	CAMK2A	5.044444444	3.79784E-07
chr14	diffreps_peak_643710	BDKRB2	3.272917393	5.2717E-07
chr2	diffreps_peak_1237629	ALS2CR12	8.163206181	4.23215E-07
chr10	diffreps_peak_240544	ANXA8	14.26625534	4.36573E-07

GO and KEGG Pathway Analysis

To examine the biological function of differentially methylated mRNAs in HGSOC, we performed GO and KEGG pathway analyses. For up-methylated m5C sites in HGSOC tissues, GO analysis indicated that these genes were mainly enriched in system development (GO term: BP), transporter complex (GO term: CC), and transporter activity (GO term: MF) (Figure 3A–C). For down-methylated m5C sites in HGSOC tissues, GO analysis indicated that these genes were mainly enriched in the regulation of signalling (GO term: BP), synapse (GO term: CC), and protein serine/threonine kinase activity (GO term: MF) (Figure 3D–F).

Figure 3

GO and KEGG enrichment analyses of HGSOC and adjacent tissues. (A) Biological processes, (B) cellular components, and (C) molecular functions of genes annotated by upregulated m5C peaks. (D) Biological processes, (E) cellular components, and (F) molecular functions of genes annotated by down-methylated m5C peaks. KEGG pathways of genes annotated by (G) up-methylated and (H) down-methylated m5C peaks. The top ten most significant terms are shown for each analysis.

Pathway analysis showed that the upregulation of mRNA in HGSOC tissues was mainly involved in the calcium signalling pathway, GABAergic synapse, and circadian entrainment pathway (Figure 3G). Down-methylated mRNAs were significantly enriched in aldosterone synthesis and secretion, inflammatory mediator regulation of TRP channels, and the gastric acid secretion pathway (Figure 3H).

Conjoint Analysis of MeRIP-Seq and RNA-Seq Data

We detected the transcriptome profiles of HGSOC tissues versus adjacent normal tissues using hisat2 software (v2.0.4). Differentially expressed genes were detected using edgeR software (v3.16.5). Compared with normal tissues, HGSOC tissues had 541 significantly upregulated genes and 319 significantly downregulated genes. The top 20 differentially expressed genes are listed in Tables 3 and 4. Subsequently, all differentially methylated m5C peaks with differential mRNA levels (125) were divided into four groups by conjoint analysis of the MeRIP-seq and RNA-seq data. We identified 101 hypermethylated m5C peaks in mRNAs that were significantly upregulated (65; hyper-up) or downregulated (36; hyper-down), while 27 hypomethylated m5C peaks in mRNAs were significantly upregulated (15; hypo-up) or downregulated (12; hypo-down) (Figure 4A) ( for details).

Table 3

Top Ten Upregulated Differentially Expressed mRNAs

Transcriptid	p-value	LogFC	GeneName
ENSG00000149043	0.000254916	8.250880408	SYT8
ENSG00000110203	6.48026E-06	8.173523376	FOLR3
ENSG00000166736	2.1651E-08	8.103963073	HTR3A
ENSG00000125285	0.00010903	7.88815242	SOX21
ENSG00000167608	0.000360554	7.810911589	TMC4
ENSG00000138100	3.48501E-05	7.61848835	TRIM54
ENSG00000167754	0.005993373	7.577582095	KLK5
ENSG00000179593	7.67683E-05	7.539994514	ALOX15B
ENSG00000186847	0.041739418	7.522791842	KRT14
ENSG00000162592	0.000656117	7.257742134	CCDC27

Table 4

Top Ten Downregulated Differentially Expressed mRNAs

Transcriptid	p-value	LogFC	GeneName
ENSG00000254126	8.55723E-07	−10.53087034	CD8B2
ENSG00000163501	4.05136E-06	−10.1102667	IHH
ENSG00000075388	1.97585E-05	−9.567486124	FGF4
ENSG00000177984	2.24728E-05	−9.519507451	LCN15
ENSG00000165556	4.825E-05	−9.222425608	CDX2
ENSG00000180347	5.04864E-05	−9.204099823	CCDC129
ENSG00000175065	6.00178E-05	−9.133359088	DSG4
ENSG00000131686	7.25797E-05	−9.054127258	CA6
ENSG00000168875	9.01074E-05	−8.672964422	SOX14
ENSG00000123454	0.000168206	−8.655541281	DBH

Figure 4

Joint analysis of MeRIP-seq and RNA-seq data. (A) Venn diagram displaying differentially expressed mRNAs with significant changes in both the m5C modification and mRNA levels. The top 10 GO terms of (B) biological processes, (C) cellular components, and (D) molecular functions were significantly enriched for mRNA-annotated genes with significant changes in both the m5C modification and mRNA levels. (E) The top ten KEGG pathways of mRNA-annotated genes with significant changes in both the m5C modification and mRNA levels.

GO and KEGG pathway analyses were performed to investigate the biological significance of these genes (125). GO analysis indicated that these genes were mainly enriched in cell growth and maintenance (12.6%) (GO term: BP), plasma membrane (46.6%) (GO term: CC), and cell adhesion molecule activity (4.5%) (GO term: MF) (Figure 4B–D). KEGG pathway analysis revealed that these genes were primarily enriched in transmission across chemical synapses (10.3%), transmembrane transport of small molecules (17.9%), and integrin cell surface interactions (7.7%) (Figure 4E).

Gene Expression Regulated by m5C Modification Correlates with Clinical Parameters of Patients with HGSOC

Finally, we found there are 15 hypomethylated m5C peaks with upregulated in mRNAs, of which MAP2K3 is an oncogene, we observed that the m5C regulatory gene MAP2K3 was enriched in several important pathways, such as the “inflammatory mediator regulation of TRP channels”, “GnRH signalling pathway”, “rap1 signalling pathway” and “MAPK signalling pathway”. Therefore, we first verified the expression of MAP2K3 in 3 pairs of samples. As expected, MAP2K3 has low expression in the paired normal tissues, and the expression in HGSOC tissues ranges from weakly positive to strong positive expression (Figure 5A). Then, we obtained 80 clinical samples to investigate the potential relationship between the abnormal m5C regulatory gene MAP2K3 and the prognosis of patients with HGSOC. Figure 5B shows representative images of MAP2K3 staining in HGSOC tissue. MAP2K3 staining was primarily in the nucleus, but the staining scores varied from low to high among HGSOC tissues. Table 5 summarises the association between the expression of MAP2K3 and the clinicopathological characteristics of patients with HGSOC. There were 28/80 cases with low MAP2K3 expression and of 52/80 cases with high MAP2K3 expression. MAP2K3 expression was correlated with CA125 levels (p < 0.001), tumour size (p = 0.001), lymph node metastasis (p = 0.008), depth of myometrial invasion (p < 0.001), and FIGO stage (p < 0.001). Cox regression analysis revealed that high expression of MAP2K3 indicated a poor prognosis (p = 0.005, Figure 5A).

Figure 5

Relationship between MAP2K3 and HGSOC clinical prognosis. (A) The expression of MAP2K3 in three pairs of HGSOC tissue specimens (×200). (B) Low and high expression of MAP2K3 was detected by IHC in 80 clinical HGSOC samples (×200). (C) Association of overall survival with MAP2K3 expression.

Table 5

Clinicopathological Variables and the Expression Status of MAP2K3 in HGSOC Patients

Characteristics			MAP2K3
			Low (%)	High (%)	p
Age	<50 years	31	10 (32.3)	21 (67.7)	0.683
	>50 years	49	18 (36.7)	31 (63.3)
CA125*	<150 U/L	23	18 (78.3)	5 (21.7)	<0.001
	>150 U/L	57	10 (17.5)	47 (82.5)
Tumor size*	<3 cm	52	25 (48.1)	27 (51.9)	0.001
	>3 cm	28	3 (10.7)	25 (89.3)
Lymph node metastasis	No	30	16 (53.3)	14 (46.7)	0.008
	Yes	50	12 (4.0)	38 (4.0)
Depth of myometrial invasion*	<1/2	21	15 (71.4)	6 (28.6)	<0.001
	>1/2	59	13 (22.0)	46 (78.0)
FIGO stage*	I or II	24	15 (62.5)	9 (37.5)	0.001
	III	34	11 (32.4)	23 (32.4)
	IV	22	2 (9.1)	20 (90.9)

Notes: Spearman-rho test was used. *P < 0.05.

Discussion

Previous studies on m5C and ovarian cancer were at the level of DNA methylation, and the focus was more on the production of 5-hydroxymethylcytosine (5-hmC) after m5C was oxidized by TET family enzymes.33 In this study, we first utilised the MeRIP-seq approach to characterise m5C methylomes in human HGSOC transcriptomes. We discovered 2050 aberrant m5C peaks in HGSOC tissues, of which 1767 m5C peaks were significantly upregulated, while 283 m5C peaks were significantly downregulated. In addition, we further revealed that the differentially methylated transcripts were enriched in various important cellular biological processes and cancer-related pathways. Moreover, we discovered genes with differentially methylated m5C peaks and synchronous differential expression by combining MeRIP-seq and RNA-seq data and explored the potential functions by GO and KEGG pathway analyses. Finally, the differentially expressed m5C modified regulatory gene MAP2K3 was further proven to be associated with the prognosis of patients with HGSOC. Taken together, our data provide a comprehensive description of m5C methylation in HGSOC and provide novel insights into the biological processes and important pathways of m5C methylation in the pathogenesis of HGSOC.

Considering that m5C modification is involved in various aspects of the mRNA life cycle, including mRNA export, translation, mRNA stability, and mRNA pre-splicing,34 increasing evidence has shown that m5C methylation modification plays an important role in the pathogenesis of many diseases including tumours.35 Mei et al showed that the RNA methyltransferase NSUN2 inhibits the expression of p57Kip2 in an m5C-dependent manner to promote the proliferation of gastric cancer cells, providing a new target for the treatment of gastric cancer.36 Previous studies have revealed that m5C modification maintains the stability of its target mRNA by recruiting ELAVL1 through the recognition protein YBX1. In addition, YBX1 drove the pathogenesis of human urothelial carcinoma by targeting the m5C methylation site in the 3′-untranslated region of HDGF.9 Recent studies have shown a wide range of abnormal m5C modifications in hepatocellular carcinoma.37 Our research identified that m5C sites were distributed in all regions of mRNA and occurred on all chromosomes, especially CD regions (HGSOC: 54.74%, adjacent tissues: 55.82%) and chromosome one, which was consistent with the results of Bohnsack et al.38 In addition, we found thousands of methylation peaks in mRNA, and found significant differences in the number and distribution of methylation peaks between the two samples. The results indicate that the methylation frequency and methylation genes of HGSOC are far less than the frequency of adjacent tissues proves that there is a clear connection between m5C and HGSOC. This may be caused by demethylation of normal tissues.39,40 The cluster analysis results showed that the degree of methylation can clearly distinguish HGSOC from adjacent tissues, which further confirms the potential relationship between m5C and HGSOC. Since previous studies on m5C and 5-hmc in ovarian cancer are more based on DNA methylation level, Zhang et al in 170 tissue samples also confirmed that compared with normal ovarian tissue, the expression of m5C in epithelial ovarian cancer tissue was significantly increased and was related to pathological stage, tumor grade, lymph node metastasis and vascular thrombosis.33 As the function of mRNA m5C in HGSOC is not completely clear, we can only make some suggestions on clinical relevance or application: compared with DNA methylation, mRNA m5C methylation has a higher degree of methylation modification. In the clinic, the combination of CA125 and other tumor markers has better sensitivity and specificity in the early diagnosis of HGSOC. These hypothesis could better guide the clinical.

It is well known that dynamic modification of m5C methylation is involved in many biological processes, such as cell proliferation,41 mitosis,42 senescence,43 migration,44 and differentiation,45 and these phenotypes are all signs of cancer. In our study, we found that the differentially expressed m5C mRNA-annotated genes were involved in system development, cell metabolic stress, and cell signal transduction by GO functional enrichment analysis. In animal models, m5C RNA modification promotes murine leukaemia virus replication in an ALYREF reader protein-dependent manner,46 and the mutations in the NSUN7 reader protein can cause sperm motility defects, leading to infertility or complete infertility in male mice.47 Notably, we found that m5C mRNA-regulated genes were enriched in some important cancer pathways, such as inflammatory mediator regulation of TRP channels pathway, the Wnt signalling pathway, and focal adhesion pathways. The regulation of inflammatory mediators in the TRP channel pathway makes an important contribution to cell proliferation, migration, inflammatory reactions, pain, and vascular muscle tension.48,49 The wnt signalling pathway and focal adhesion pathways, as classic carcinogenic pathways, have been extensively studied in ovarian cancer.50,51 In summary, these results indicate that differentially expressed m5C mRNA modification exerts precise functions in regulating the occurrence of HGSOC.

By combining the MeRIP sequence and RNA sequence data, we identified 125 genes with different m5C methylation peaks and synchronous expression in HGSOC. Similarly, these genes were enriched in some important functional processes, such as cell growth or maintenance, as well as in plasma membrane and cell adhesion molecule activity. These results indicate that differentially expressed mRNAs also contribute significantly to HGSOC. Messenger RNA as a coding RNA and mRNA modification may affect most of the post-transcriptional steps in gene expression.52 Therefore, research on the mechanism of m5C mRNA regulation has become a hot topic. Fukumoto et al showed that methylation modification of FZD10 mRNA contributes to PARPi resistance by up-regulating the Wnt/β-catenin pathway in BRCA-mutated ovarian cancer cells.53 In ovarian cancer cell models, overexpression of the mRNA methyltransferase METL3 significantly increased cell proliferation, focal formation, invasion, and tumour formation. The high expression of METL3 was significantly related to tumour grade, FIGO stage, and overall survival rate,54 which was consistent with our study results. We further detected the expression of the differentially expressed m5C mRNA regulatory gene MAP2K3 in 80 human HGSOC tissues and found that MAP2K3 was highly expressed in HGSOC tissues and was correlated with poor prognosis in patients with HGSOC. Research by Park et al showed that abnormal methylation of MAP2K3 is related to the occurrence of lung cancer.55 However, the specific mechanism by which the m5C modified gene MAP2K3 leads to the poor prognosis of HGSOC needs further study. Taken together, our research shows that the regulation of m5C mRNA modification is a promising therapeutic strategy for HGSOC in the future.

In conclusion, our study revealed, for the first time, the difference in the number and distribution of RNA m5C peaks between HGSOC tissue and adjacent tissues. We further conducted a bioinformatics analysis showing the potential function of m5C modification in the pathogenesis of HGSOC at the level of epitranscription. Targeting the regulation of m5C mRNA represents a new therapeutic approach to guide personalised treatment for patients with HGSOC, but the specific mechanisms involved need to be further explored.