Synchronous profiling of mRNA N6-methyladenosine modifications and mRNA expression in high-grade serous ovarian cancer: a pilot study

This study aimed to synchronously determine epitranscriptome-wide RNA N6-methyladenosine (m6A) modifications and mRNA expression profile in high grade serous ovarian cancer (HGSOC). The methylated RNA immunoprecipitation sequencing (MeRIP-seq) was used to comprehensively examine the m6A modification profile and the RNA-sequencing (RNA-seq) was performed to analyze the mRNA expression profile in HGSOC and normal fallopian tube (FT) tissues. Go and KEGG analyses were carried out in the enrichment of those differentially methylated and expressed genes. MeRIP-seq data showed 53,794 m6A methylated peaks related to 19,938 genes in the HGSOC group and 51,818 m6A peaks representing 19,681 genes in the FT group. RNA-seq results revealed 2321 upregulated and 2486 downregulated genes in HGSOC. Conjoint analysis of MeRIP-seq and RNA-seq data identified differentially expressed genes in which 659 were hypermethylated (330 up- and 329 down-regulated) and 897 were hypomethylated (475 up- and 422 down-regulated). Functional enrichment analysis indicated that these differentially modulated genes are involved in pathways related to cancer development. Among methylation regulators, the m6A eraser (FTO) expression was significantly lower, but the m6A readers (IGF2BP2 and IGF2BP3) were higher in HGSOC, which was validated by the subsequent real-time PCR assay. Exploration through public databases further corroborated their possible clinical application of certain methylation regulators and differentially expressed genes. For the first time, our study screens the epitranscriptome-wide m6A modification and expression profiles of their modulated genes and signaling pathways in HGSOC. Our findings provide an alternative direction in exploring the molecular mechanisms of ovarian pathogenesis and potential biomarkers in the diagnosis and predicting the prognosis of the disease.


Clinical specimens
Both HGSOC and normal FT tissues were collected at the time of surgery in Zhejiang cancer hospital in 2021.Normal FT tissues were obtained from patients with uterine leiomyomas.Samples were quickly snap-frozen in liquid nitrogen after surgical incision and then preserved at − 80 ℃ before use.None of these patients received chemotherapy or radiotherapy before the surgery.Two independent gynecological pathologists confirmed their diagnosis.A total of eleven HGSOC and eleven FT tissues were randomly selected for this study (Supplemental Table S1).

Ethical approvement
The study was approved by the Institutional Review Board of Zhejiang cancer hospital.Each patient signed a written informed consent in the study before surgery.The study was performed in accordance with the Declaration of Helsinki.

RNA extraction from HGSOC and normal FT tissues
Total RNA was isolated and purified using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's procedure.The RNA amount and purity of each sample were quantified using NanoDrop ND-1000 (NanoDrop, Wilmington, DE, USA), and the RNA integrity was assessed by Bioanalyzer 2100 (Agilent, CA, USA).RNAs from three HGSOC and three normal FT tissues for the sequencing analyses.

MeRIP-seq and RNA-seq in HGSOC and normal FT tissues
Poly (A) RNA was purified from 50 μg total RNA using Dynabeads Oligo (dT) 25-61005 (Thermo Fisher, CA, USA) by two rounds of purification.Then the poly(A) RNA was fragmented into small pieces using Magnesium RNA Fragmentation Module (NEB, cat.e6150,USA) under 86 ℃ 7 min.The cleaved RNA fragments were incubated for 2 h at 4 ℃ with m6A-specific antibody (No. 202003, Synaptic Systems, Germany) in IP buffer (50 mM Tris-HCl, 750 mM NaCl and 0.5% Igepal CA-630).The IP RNA was reversely transcribed to cDNA by SuperScript™ II Reverse Transcriptase (Invitrogen, cat.1896649, USA), which was next used to synthesize U-labeled second-stranded DNAs with E. coli DNA polymerase I (NEB, cat.m0209,USA), RNase H (NEB, cat.m0297, USA) and dUTP Solution (Thermo Fisher, cat.R0133, Waltham, MA, USA).An A-base was added to the blunt ends of each strand, preparing for ligation to the indexed adapters.Dual-index adapters were ligated to the fragments, and size selection was performed with AMPureXP beads.After the heat-labile UDG enzyme (NEB, cat.m0280, USA) treatment of the U-labeled second-stranded DNAs, the ligated products were amplified with PCR by the following conditions: initial denaturation at 95 ℃ for 3 min; 8 cycles of denaturation at 98 ℃ for 15 s, annealing at 60 ℃ for 15 s, and extension at 72 ℃ for 30 s; and then final extension at 72 ℃ for 5 min.At last, we performed paired-end sequencing (PE150) on an Illumina Novaseq™ 6000 platform (LC-Bio Technology CO., Ltd., Hangzhou, China) following the vendor's recommended protocol.

Data analysis
Cutadapt and perl scripts in-house were used to remove reads that contained adaptor contamination, low quality bases, and undetermined base.The sequence quality was verified by Fastp.The high-quality clean IP and input reads were aligned to the genome sequences of homo sapiens (Version: v96) with default parameters using HISAT2 software.The aligned reads were used for peak calling of the MeRIP regions by R package exomePeak, while IGV software was used to visualize identified m6A peaks with bed or bam formats (http:// www.igv.org/).MEME and HOMER were used for de novo and known motif finding, followed by localization of the motif involving peak summit by perl scripts in-house.ChIPseeker was used to conduct called peaks annotation and intersection with gene architecture.
The levels of methylated peaks were normalized by RPM (Reads per million reads).The expression levels for all genes from input libraries were calculated as FPKM (Fragments per kilobase of exon per million fragments mapped = total_exon_fragments/mapped_reads (millions) X exon_length (kB)) using StringTie.Based on the log2 (fold change) > 1 or log2 (fold change) < − 1, the differentially expressed methylation peaks were identified by the Poisson Test (P < 0.05), and their related genes by the Exact test (P < 0.05).The levels of mRNAs in RNA-seq data were presented as FPKM.The differentially expressed mRNAs were identified by the same fold change criteria and the Exact test (P < 0.05) using R package edgeR.Go and KEGG 23 analyses were performed in differentially methylated and expressed genes. www.nature.com/scientificreports/

Real-time (RT)-PCR validation of m6A regulator genes
The analyses were performed in 11 HGSOC and 11 FT tissues (Supplementary Table S1).The primer of the target genes and β-actin in RT-qPCR were designed and synthesized by Sangon Biotech Co., Ltd (Shanghai, China) (Supplementary Table S2).The SuperScript IV cDNA synthesis Kit (Thermo Fisher, Waltham, MA, USA) was used for cDNA synthesis.SYBR Green master mix from Bimake (Munich, German) and the ABI 7500 system (Thermo Fisher, Waltham, MA, USA) was used for RT-PCR.All PCR experiments were conducted in triplicate and the average value was calculated.The gene expression values were normalized to β-actin.
The intensity and distribution patterns of IHC staining were evaluated using a semi-quantitative Immunoreactive Score (IRS).The IRS is calculated based on two parameters: the staining intensity (0, no staining; 1, weak staining; 2, moderate staining; and 3, strong staining) and the percentage of positive cells (0, no positive cells;1, less than 10% of cells positive; 2, 10-50% of cells positive; 3, 51-80% of cells positive; 4, more than 80% of cells positive).The IRS was calculated by multiplying these two scores.

Correlation of gene expression with clinical significance using public databases
We then explored the expression of genes between ovarian cancer and normal ovary tissues using the TCGA and GTEx public databases (https:// xenab rowser.net) 24 .The Kaplan-Meier Plotter database was applied to examine the association between the expression of genes and survival in serous ovarian cancer patients (https:// kmplot.com/ analy sis/) 25 .

Statistical analysis
The Mann-Whitney U Test was performed to compare gene expression from the RT-PCR, and the different IRS scores of IHC staining between HGSOC and FT groups.Statistical analyses were performed using SPSS (version 22.0, SPSS Inc., Chicago, IL, USA).A P value less than 0.05 was considered statistically significant.

Global m6A modification patterns in HGSOC and FT
Reads data and quality testing for MeRIP-seq were shown in Supplementary Table S3.Venn diagrams were constructed to show the m6A peaks (Fig. 1A) and their modified genes (Fig. 1B) in the two groups.A total of 28,491 m6A peaks (Fig. 1A) and their modified 4629 genes (Fig. 1B) were detected in only the HGSOC group.As many as 26,515 m6A peaks (Fig. 1A) and their modified 4372 genes (Fig. 1B) were identified only in the FT group.In addition, 25,303 m6A peaks and their modified 15,309 genes were shared by both groups (Fig. 1A,B).
The abundance of m6A methylated genes between the two groups was compared.A total of 2444 genes were hyper-methylated and 3265 genes were hypo-methylated in the HGSOC group compared with the FT group (Fig. 1C).The top 30 methylated genes with altered m6A methylated peaks are listed in Table 1.There was no significant difference in the number of m6A peaks per gene between the two groups (Fig. 1D).The majority of genes (75%) had 1-3 m6A peaks and about 36% of genes had one modified m6A peak (Supplementary Fig. 1A).Among those differentially m6A modified genes, 77% of them harbored one m6A methylated peak and approximately 98% of genes had 1-3 m6A modified sites (Fig. 1E and Supplementary Fig. 1B).

M6A topological patterns in HGSOC
The distribution of m6A peaks in the entire epitranscriptome was analyzed in both HGSOC and FT groups.The m6A peaks were identified in the start and codons, 3′ untranslated regions (3′UTR), 5′UTR, and coding sequence (CDS) on transcripts.In general, these m6A peaks were more enriched in CDS and 3′UTR (Fig. 2A).The distributions of differential and common m6A peaks between the two groups were displayed using pie graphs (Fig. 2B).The m6A hypermethylated peaks in the HGSOC group were found more in CDS (64.8% vs 57.6%) but less in 3′UTR (23.4% vs 30.8%) compared to those in the FT group.The common m6A peaks showed a relatively higher proportion in 3′UTR than in CDS (45.4% vs 31.6%).The classical consensus GGAC motif structures were identified and thus confirmed the authenticity of m6A peaks in both HGSOC and FT groups (Fig. 2C).

Differentially methylated m6A genes were enriched in important signaling pathways
Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) analyses were performed in differentially methylated genes to investigate the biological significance and potential functions of m6A modification.In GO analysis, genes were classified into three functional domains: biological process, cellular component, and molecular functions.The top 25 biological process terms, 15 cellular component terms and 10 molecular function terms are listed in Fig. 3A.The differentially methylated genes were mainly enriched in the regulation of transcription, signal transduction, and regulation of transcription by RNA polymerase II, et al.In terms of cellular components, genes were mainly involved in cell cytoplasm, membrane, and nucleus, et al.In molecular function terms, genes played a part in protein binding, metal ion binding, and DNA binding, et al.GO enrichment analysis revealed that genes were significantly enriched in DNA binding, metal ion binding, regulation of transcription, G protein-coupled receptor signaling pathway, et al. (Fig. 3B).KEGG pathway analysis revealed that transcriptional misregulation, FoxO signaling pathway, Rap1 signaling pathway, IL-17 signaling pathway, cAMP signaling pathway, Hedgehog signaling pathway, MAPK signaling pathway and Wnt signaling pathways were significantly altered in HGSOC (Fig. 3C).

Functional analysis of differentially expressed genes between HGSOC and FT tissues
The RNA-seq was used to examine differentially expressed genes between HGSOC and FT tissues.The heatmap was used to present the top differentially expressed genes between the two groups (Fig. 4A).There were 4807 www.nature.com/scientificreports/differentially expressed genes between two groups, including 2321 upregulated and 2486 downregulated in HGSOC.The expression of 29,884 genes was not significantly different between groups (Fig. 4B,C).

Conjoint analysis of m6A methylation and gene expression profiles
By conjoint analysis of MeRIP-seq and RNA-seq data, genes were classified into four groups: 659 hyper-methylated m6A peaks with 330 genes upregulated (hyper-up) and 329 genes downregulated (hyper-down), and 897 hypo-methylated m6A peaks with 475 genes upregulated (hypo-up), and 422 genes downregulated (hypo-down) (Fig. 4D).Results of GO and KEGG analyses in these four groups of genes are presented in Supplementary Fig. S2A-H.The top differentially modulated genes in each group are presented in Tables 2 and 3. We performed GO and KEGG analyses in those differentially m6A methylated and expressed genes.The top 10 terms of biological process, cellular component and molecular functions are listed in Fig. 5A.KEGG pathway analysis revealed that these genes were enriched in cGMP-PKG pathway and calcium signaling pathway, et al. (Fig. 5B).

Expression of genes related to m6A regulators between HGSOC and FT tissues
RNA-seq data revealed the expression levels of genes related to m6A regulators in HGSOC and FT tissues (Fig. 6A-J).Our data showed that IGF2BP2 and IGF2BP3 were significantly increased, whereas FTO was significantly decreased in HGSOC.In contrast, there was no significant change in other m6A regulators.Expression of these genes was further examined by RT-PCT in 11 HGSOC and 11 FT tissues (Fig. 7A-J).The protein levels of these m6A regulators were further examined using IHC staining (Supplementary Fig. S3A,B).Both PCR and IHC staining results matched the expression patterns observed from profiling analyses.

Results of m6A methylation regulators and differentially expressed genes in public databases
We also explore differentially expressed genes in the TCGA and GTEx public databases.Consistent with our results, IGF2BP2 and IGF2BP3 were significantly increased in ovarian cancer, whereas FTO was significantly decreased in ovarian cancer.PDE10A and PTGS2 were significantly decreased, but ITGB7 was significantly increased in ovarian cancer (Fig. 8A-K).Data from the Kaplan-Meier Plotter database indicated the association between the expression of genes and survival in ovarian cancer patients (Fig. 9A-G).

Discussion
Abnormal m6A modifications have been associated with proliferation, drug resistance, and prognosis of ovarian and other cancers 22,26,27 .The global m6A profile has previously been screened in several cancers, such as endometrioid ovarian cancer 28 , invasive malignant pleomorphic adenoma 29 , and colorectal cancer 30 .For the first time, this study examined the epitranscriptome-wide m6A methylation and synchronously determined the gene expression profile in HGSOC and FT tissues.Besides the feature of m6A modification sites, our study revealed the differential global m6A methylation and differentially expressed gene expression profile between HGSOC and FT tissues.The conjoint analysis further identified those differentially m6A modified and expressed genes.Enrichment analysis suggests the involvement of these genes in multiple complex functions and signaling pathways correlated with ovarian development and progression.Among them, m6A eraser FTO and readers IGF2BP2 and IGF2BP3 were differentially modulated between HGSOC and FT tissues.Parts of our findings were validated with subsequent RT-PCR and IHC staining.The expression levels and clinical roles of significant genes were further examined in ovarian cancer patients using public databases.Previous studies have demonstrated the critical role of M6A regulators FTO, IGF2BP2 and IGF2BP3 in ovarian and other cancers.IGF2BP2 and IGF2BP3 were positively associated with the expression of C3AR1, a member of the G protein coupled receptor 1 family, in ovarian cancer 31 .C3AR1 overexpression significantly increased the proliferation of ovarian cancer SKOV3 cells and was associated with immunosuppression and poor prognosis of ovarian cancer 31 .Overexpression of C3AR1 was associated with a poor prognosis in patients with www.nature.com/scientificreports/gastric adenocarcinoma 32 .IGF2BP2 contributed to colorectal cancer pathogenesis and progression by SOX2 m6A methylation and preventing its degradation 33 .IGF2BP3 was shown to modulate cell cycle, cell proliferation and angiogenesis through m6A modification to increase both expression and stability of cyclin D1 (CCND1)   34 .Downregulation of m 6 A erasers FTO and ALKBH5 increased FZD10 m 6 A modification and its mRNA expression, and reduced PARPi sensitivity in BRCA-mutated epithelial ovarian cancers 35 .Overexpression of FTO led to hypomethylation and reduced mRNA stability of two phosphodiesterase genes (PDE1C and PDE4B) but augmented cAMP signaling and dampen stemness features of ovarian cancer cells 19 .FTO decreased the m6A level and stability of SNAI1 mRNA, causing IGF2BP2 dependent downregulation of SNAI1 and inhibiting epithelial-mesenchymal transition 36 .Overexpression of FTO inhibited tumor growth in nude mice by facilitating the oxidative stress response and apoptosis of ovarian cancer cells via activation of the P53 signaling pathway 37 .Our study found that m6A eraser FTO was significantly lower, whereas IGF2BP2 and IGF2BP3 were significantly higher in HGSOC than in normal FT.The profiling data was confirmed by subsequent RT-PCR analyses in more samples.This change pattern of these genes between normal and ovarian cancer tissues is observed in data from the TCGA database.Data from public database Kaplan-Meier Plotter revealed that the expression of FTO and IGF2BF3 was significantly correlated with the prognosis of ovarian cancer patients.Other m6A regulators, such as m6A writers METTL3, METTL14, WTAP, m6A eraser ALKBH5, and m6A readers YTHDF1, YTHDF2, YTHDF3 and IGF2BP1, have been shown their important roles in the ovarian and other cancers 13,22,38,39 .In an m6A-dependent way, IGF2BP1 inhibited the miRNA-mediated degradation and thus increased the expression of SRF, FOXK1, and PDLIM7 40,41 .Increased expression of IGF2BP1 is associated with a poor prognosis of ovarian and other cancers [40][41][42] .METTL3 was shown to promote ovarian cancer cell proliferation, invasion, and migration through targeting pri-miRNA 126-5p, lncRNA RHPN1 antisense RNA 1 (head to head) (RHPN1-AS1) and AXL receptor tyrosine kinase (AXL) mRNA 18,43,44 .ALKBH5 erased the m6A modification of Janus kinase 2 (JAK2) and stabilized JAK2 mRNA, leading to cisplatin resistance in ovarian cancer 45 .The m6A reader YTHDF1 promoted ovarian cancer progression via augmenting EIF3C translation 46 .In addition, data from the TCGA database showed the differential expression of tIGF2BP1, YTHDF1 and YTHDF2 genes between normal and ovarian cancer tissues.This study did not observe the difference in expression of these genes between HGSOC and normal FT tissues.Though more studies are needed to corroborate the result, our study implies that FTO, IGF2BP2 and IGF2BP3 may play more important roles than other regulators in the m6A regulation in HGSOC in Chinese women.It is worthwhile to further study mechanisms in modulating these m6A regulators in ovarian cancer.
GO analysis revealed that the G protein-coupled receptor (GPCR) signaling pathway was one of the functions related to m6A modification.GPCRs are the largest family of cell membrane receptors involving signal transduction related to many key physiological processes.Recent studies have shown that the aberrant GPCRs contribute to many aspects of tumorigenesis, including proliferation, invasion, angiogenesis, metastasis, and survival 47,48 .Systematic analysis of large-scale genomes revealed that GPCRs were mutated in about 20% of all human tumors, including ovarian cancer 49 .The alterations of gene expression and promoter methylation of GPCRs have also been reported previously 50 .An estrogen receptor GPER, a member of the GPCR family, was reported to promote the migration and invasion of ovarian cancer cells 51 .In contrast, other studies suggested that GPER was correlated with the prolonged overall survival of epithelial ovarian cancer 52,53 , and the agonist of GPER can inhibit ovarian cancer cell proliferation 54 .
KEGG analysis showed that m6A hyper-methylated genes were enriched in the FoxO signal pathway.The forkhead box O (FoxO) family of transcriptional factors consists of four members: FoxO1, FoxO3, FoxO4 and FoxO6.They participate in the regulation of various steps of cancer initiation and metastasis, though their regulation patterns have not been fully illustrated.Sisci et al. 55 reported that in breast cancer, FoxO3a overexpression decreased the motility, invasiveness, and anchorage-independent growth in estrogen receptor α-positive (ERα+)  This study also identified many other differentially m6A methylated and expressed genes involving signaling pathways related to ovarian pathogenesis and progression 37 .Our data showed that the Wnt singling pathway is one of the m6A modulated pathways in HGSOC.The previous study indicated that upregulation of the Wnt signaling pathway by elevated m6A contributed to Poly (ADP-ribose) polymerase inhibitors (PARPi) resistance.The Wnt signaling pathway is a potential therapeutic target in epithelial ovarian cancer 35 .All m6A modulated genes and their related complex signaling pathways broaden our understanding of the molecular mechanism for pathogenesis, drug resistance and progression of ovarian cancer.More studies are warranted to utilize these genes and pathways as biomarkers for early diagnosis and targets for treatment.
Limitations of this study include that both MeRIP-seq and RNA-seq were conducted in a small number of samples.As a retrospective cross-sectional study, there is a potential bias in sample selection.This study is limited to examining profiling data at one time point.The changes in gene methylation and expression during ovarian cancer development and progression are not determined.The MeRIP-seq has its limits in the detection of m6A methylation 59 .The functions of m6A methylation regulators are fulfilled through their proteins.Their protein levels were determined with IHC staining but with a limited number of samples.The target genes are important in mediating the functions of these regulators.Further knock-out or knock-in studies are needed to reveal the functions of these regulators and their target genes in ovarian cancer development.More studies were needed to develop these identified genes as biomarkers in diagnosing HGSOC.The strength of this study is the www.nature.com/scientificreports/synchronous profiling of both m6A methylation and gene expression in HGSOC and FT tissues and the validation using RT-PCR, IHC staining, and public databases.

Conclusions
These differentially m6A methylated and expressed genes and their related signaling pathways may help to better understand molecular mechanisms for the development and progression of ovarian cancer.These genes may have the potential to diagnose and predict the prognosis of ovarian cancer.More in-depth studies are warranted to further explore their application in novel diagnosis and treatment of the disease.

Figure 1 .
Figure 1.Overview of m6A distributions in genes between HGSOC and FT tissues (n = 3 each).Venn diagrams showing (A) m6A methylation peaks and (B) their related genes only in HGSOC, or only in FT tissues, or in both.(C) Volcano plots of the significantly differentially m6A methylated genes between HGSOC and FT tissues.(D) Number of m6A peaks per gene in HGSOC and FT tissues.(E) Number of m6A peaks per differentially methylated gene.FT, fallopian tube; HGSOC, high grade serous ovarian cancer.

Figure 2 .
Figure 2. Topological patterns of m6A distributions in genes between HGSOC and FT tissues (n = 3 each).(A) Distribution site of the m6A peaks on all the transcripts in HGSOC and FT.Each transcript was divided into 5ʹ-UTR, CDS, and 3ʹ-UTR regions.(B) Pie charts showed the distribution proportion of unique and common m6A peaks between HGSOC and FT groups.(C) The representative m6A motifs displayed GGAC conserved sequence from altered m6A peaks in HGSOC and FT groups.CDS, coding sequence; FT, fallopian tube; HGSOC, high grade serous ovarian cancer; UTR, untranslated region.

Figure 3 .
Figure 3. Result of GO function and KEGG pathway analyses in differentially methylated genes.(A) Bar plot of the top 25 biological process terms, 15 cellular component terms and 10 molecular function terms of differentially methylated genes in HGSOC in GO analysis.(B) Scatter plotter of the top 20 significant GO terms.(C) A scatter plot of the top 20 significant KEGG enrichment pathways.

Figure 4 .
Figure 4. Conjoint analysis result of MeRIP-seq and RNA-seq data.(A) Heatmap plots of the differentially expressed genes between HGSOC and FT tissues.(B) The number of differentially expressed genes.(C) Scatter plot of the differentially expressed genes.(D) The four-quadrant plot of differentially m6A methylated and expressed genes.FT, fallopian tube; HGSOC, high grade serous ovarian cancer.

Figure 5 .
Figure 5. Result of GO terms and KEGG pathways of differentially methylated and expressed genes between HGSOC and FT tissues (n = 3 each).(A) Bar plot of the top 10 biological processes, cellular component and molecular function terms in GO analysis.(B) Scatter plot of the top 10 KEGG pathways.FT, fallopian tube; HGSOC, high grade serous ovarian cancer.

Figure 9 .
Figure 9. Survival curves of ovarian cancer patients with high or low expression of genes using the Kaplan Meier plotter database.HR, Hazard ratio.

Table 2 .
The top hypermethylated and differentially expressed genes in HGSOC vs FT.FT, fallopian tube; HGSOC, high grade serous ovarian cancer.