Dissecting the effects of METTL3 on alternative splicing in prostate cancer

Although the role of METTL3 has been extensively studied in many cancers, its role in isoform switching in prostate cancer (PCa) has been poorly explored. To investigate its role, we applied standard RNA-sequencing and long-read direct RNA-sequencing from Oxford Nanopore to examine how METTL3 affects alternative splicing (AS) in two PCa cell lines. By dissecting genome-wide METTL3-regulated AS events, we noted that two PCa cell lines (representing two different PCa subtypes, androgen-sensitive or resistant) behave differently in exon skipping and intron retention events following METTL3 depletion, suggesting AS heterogeneity in PCa. Moreover, we revealed that METTL3-regulated AS is dependent on N⁶-methyladenosine (m⁶A) and distinct splicing factors. Analysis of the AS landscape also revealed cell type specific AS signatures for some genes (e.g., MKNK2) involved in key functions in PCa tumorigenesis. Finally, we also validated the clinical relevance of MKNK2 AS events in PCa patients and pointed to the possible regulatory mechanism related to m⁶A in the exon14a/b region and SRSF1. Overall, we characterize the role of METTL3 in regulating PCa-associated AS programs, expand the role of METTL3 in tumorigenesis, and suggest that MKNK2 AS events may serve as a new potential prognostic biomarker.

1 Introduction

PCa is the second most prevalent malignancy in males and the fifth leading cause of cancer mortality in men worldwide (1). Clinically, relatively indolent PCa (i.e., Gleason grade (GS)≤7) is treated with radical prostatectomy with a favorable prognosis, whereas locally progressed and metastatic PCa is typically treated with androgen deprivation therapy (ADT) (2). Notably, despite current therapies that improve prognosis, most patients will eventually fail ADT and progress to a deadly condition (androgen insensitive, AR^-) known as castration-resistant PCa (CRPC) (3). The molecular mechanisms of PCa progression, particularly CRPC, remain largely undefined.

m⁶A modification is ubiquitous in most eukaryotes (4) and plays crucial roles in a variety of bioprocesses, such as DNA damage response (5). Analogous to DNA methylation, the dynamic m⁶A is regulated by the methyltransferase complex (writers), demethylases (erasers), and reader proteins (6). Notably, METTL3 catalyzes the m⁶A modification process, which can be dramatically reduced by METTL3 deletion (7). Recently, a growing number of studies have revealed the link between imbalanced m⁶A levels and cancer progression (6). Experimentally, METTL3 acts synergistically with YTHDF2 to inhibit the expression of SOCS2, encouraging liver cancer (HCC) development and progression (8). In PCa, METTL3 enhances MYC expression thus promoting cell proliferation (9), promotes the invasion through regulating the SHH-GLI1 signaling pathway (10), and modulates the HuR-ITGB1 signaling pathway thus contributing to bone metastasis (11). Despite a few findings in PCa concentrating on specific genes on a small scale, the potential biological consequence of global m⁶A abnormalities on a genome-wide scale remains deserved.

Additionally, m⁶A alternations not only influence gene expression by changing RNA stability and mRNA export but also regulate AS and 3′-end processing (12). AS is a dynamic process coupled with transcription and involved in many physiological functions as well as disease pathogenesis (13). Evidence has demonstrated that METTL3 deletion can attenuate the inflammatory response induced by lipopolysaccharide by modifying the expression patterns of MyD88 isoforms (14). Moreover, the m⁶A reader YTHDC1 can interact with SRSF3 to regulate exon inclusion (15), while FTO regulates the process of SRSF2 binding to the m⁶A-containing RUNX1T1 transcript and results in exon skipping (16). However, a comprehensive study of the AS landscape governed by METTL3 has not been conducted in PCa.

Recently, the development of long-read sequencing (e.g., direct RNA sequencing of Oxford Nanopore Technologies) provides the possibility to study splicing variants and alternative polyadenylation events (17), which directly obtain the structural information of the entire transcript without the assembling step like standard RNA sequencing (RNA-Seq). Another significant advantage of direct RNA sequencing (DRS) over standard RNA-Seq is the use of long reads to correctly dissect transcript isoforms and simultaneously detect RNA-modified sites (e.g., m⁶A) (18).

In this study, we combined standard RNA-Seq and long-read DRS of Oxford Nanopore to comprehensively dissect the impact of METTL3 on RNA splicing in the setting of two PCa cell lines with distinct AR sensitivities. We revealed that METTL3 potentially had shared and unique mechanisms affecting alternative splicing in AR⁺ and AR^- PCa cells. The different effects of METTL3 on exon skipping and intron retention might be related to distinct splicing partners.

2 Materials and methods

2.1 Cell culture

Normal epithelial prostate cell line (RWPE-1) and PCa cell lines (LNCaP and DU145) were obtained from Prof. Andrew Chan’s laboratory in the School of Biomedical Sciences, The Chinese University of Hong Kong. RWPE-1 cells were cultured using the KSFM-1 medium kit (Thermo Fisher, No.17005042) with supplementary bovine pituitary extract (BPE) and human recombinant epidermal growth factor (rEGF). LNCaP cells were maintained in RPMI-1640 medium (Invitrogen) supplemented with 10% FBS and 1% PS, and DU145 cells were cultured in DMEM (Invitrogen) supplemented with 10% FBS and 1% PS. All cells were incubated at 37°C and 5% CO2.

2.2 siRNA transfection, RNA extraction, and polyA RNA isolation

The siRNA against METTL3 (Invitrogen, Catalog: #4392422) or control was transfected into LNCaP and DU145 cells using RNAiMAX reagent (Invitrogen). Total RNA was extracted with a Pure-link mini kit (Invitrogen) 48 h after transfection. RNA concentration was measured via a Nanodrop 2000, and RNA integrity (RIN) was detected using an RNA 6000 nanochip (Agilent Technologies) (Requirement: RIN>8.8). PolyA RNA was isolated from 75 μg of total RNA using a Dynabeads™ kit.

2.3 Real-time quantitative PCR and western blot

cDNA transcription was performed using a SuperScript III synthesis kit (Invitrogen) and then used for qPCR based on SYBR (Applied Biosystems). Two siRNAs (siM3-1 and siM3-2) targeting METTL3 were designed, and we observed that the knockdown efficiency of siM3-1 was only 50%, which we considered insufficient for further analyses. Therefore, the more effective target-2 (siM3-2) was selected for the following analysis (Figure S1). The sequences of primers and siRNAs used in this study are shown in Supplementary Table S1.

The siRNA-transfected LNCaP and DU145 cells were lysed in RIPA lysis buffer and Western blotting was conducted as previously (9). Then, the intensity of METTL3 (Proteintech, 15073-1-AP) and GAPDH (Santa Cruz, sc-25778) bands was quantified using ImageJ software and GAPDH was used as the housekeeping marker.

2.4 RNA-Seq data processing

Paired-end sequencing (2 × 150 bp) of siControl- and siMETTL3-transfected cells was performed on a HiSeq 4000 sequencing platform (Illumina). Clean reads were aligned to the human hg38 genome (UCSC) using HISAT2 (v2.2.1) software (19). PCA of three replicates was executed, with one outlier in the DU145 siControl group excluded from downstream analysis. Differential gene expression (DEG) analysis was conducted using DESeq2 (20). Gene Ontology (GO) and KEGG enrichment analyses were performed using the clusterProfiler R package (21).

2.5 TCGA data analysis

A table containing FPKM values from 551 adjacent normal, and PCa tissue samples was downloaded from the TCGA database (https://portal.gdc.cancer.gov/projects/TCGA-PRAD), wherein we also retrieved the clinical data of 499 primary PCa patients. This FPKM table was used for expression profile analysis of m⁶A regulators. The association between METTL3 expression and overall survival and disease-free survival was determined using the GEPIA2 web server (22). The functional enrichment of m⁶A regulators was analyzed using Metascape (23) and displayed using Cytoscape software (24).

2.6 Nanopore DRS data processing

As previously described, each 1 μg purified polyA RNA was subjected to library preparation following the SQK_RNA002 kit procedure and consequently sequenced on a Nanopore GridION platform for a 48-h runtime (25). Base-calling was performed with Guppy (v4.2.2) to generate passed reads with quality scores ≥7. The alignment of passed reads against human hg38 (UCSC) and transcriptome (GENCODE v32) was implemented via minimap2 (v2.17-r941) (26). Statistics of the passed reads were assessed using NanoStat (v1.4.0) (27). The sequence length distribution of each sample was generated using the wub package (https://github.com/Nanoporetech/wub).

2.7 AS analysis of DRS and RNA-Seq data

Splicing analysis of DRS data was performed with the FLAIR (v1.5) pipeline (17). Specifically, BED12 files were converted from primary BAM files and then corrected by annotated GTF files and junction files. Junction files from RNA-Seq data were obtained from the same cell lines, and splice sites (< 3 reads) were filtered out. Then, the corrected files were assembled into high-confidence isoforms with the stringent parameters of at least three supporting reads. Splicing events were determined using the FLAIR diffSplice module based on Fisher’s exact test.

Differential AS events of RNA-Seq data were identified using rMATS (v4.1.0) with a strict parameter of cstat=0.05 (28, 29) and filtered with a p-value <0.05. As an indicator of splicing events, the percent spliced in (PSI) index of differential intron retention and exon splicing events refers to the ratio of included reads to total reads (sum of included and excluded reads) (30). The PSI values of differential AS events were assessed using the Wilcoxon rank-sum test, which is used to provide a global overview of inclusion/exclusion levels between the control and knockdown groups. The global difference in PSI distributions was also evaluated by using FDR<0.05, which was supplemented in the Supplementary Figure S8.

The PSI index of MKNK2a/b isoforms in PCa was downloaded from the TCGA SpliceSeq database. The MKNK2_AA_46570 event represents the PSI index of sequencing reads covering between exon 13 (e13) and exon 14a (e14a), indicating the expression of the MKNK2a isoform. While the MKNK2_AT_46567 represents the PSI index of the sequencing reads covering between exon 13 and exon 14b (e14b), producing the MKNK2b isoform. According to the previous reference (31), we set the threshold at the 80^th percentile, dividing the PCa patients into high PSI and low PSI groups. The survival difference between high and low PSI groups was calculated by using “survival” (3.5) and “survminer” (0.4.9) R software.

2.8 Detection of differential m⁶A sites

Differential m⁶A sites between the siControl and siMETTL3 groups were detected using ELIGOS2 and DRUMMER (32, 33). Both approaches were designed to identify RNA modifications through the comparative analysis of base-calling errors in the Nanopore DRS dataset (34). In our study, ELIGOS2 was utilized to assess the global difference in m⁶A modifications, employing strict threshold values of p<0.05 and odd ratio >1.2 (https://gitlab.com/piroonj/eligos2) (32). The criterion of “odd ratio >1.2” indicated at least a 1.2-fold difference between the control and METTL3 knockdown conditions. Therefore, the differential m⁶A sites utilized in the further analysis primarily represented the decreased m⁶A sites upon METTL3 knockdown. DRUMMER can identify RNA modifications at each site on distinct transcript isoforms by using isoform mode (https://github.com/DepledgeLab/DRUMMER) (33). Therefore, DRUMMER was utilized to detect isoform-specific differences in m⁶A levels. All differential m⁶A data from ELIGOS2 and DRUMMER were then selected with a p-value <0.05 and the reference “A” site. The 5-kmer motif of differential m⁶A sites was generated via WebLogo (https://weblogo.berkeley.edu/logo.cgi). The genomic features of m⁶A sites were annotated with the ChIPseeker R package (35). The detailed differential m⁶A results were included in the Supplementary Tables S3–6.

2.9 Motif analysis of spliced genes

The regions of differential intron retention and exon skipping events were extracted for motif analysis using HOMER (http://homer.ucsd.edu/homer/motif/) software. The resulting motifs were visualized using the ggseqlogo R package.

2.10 Statistical analysis

Differences in PSI values of differential intron retention and exon skipping events between groups were assessed using the Wilcoxon rank-sum test. The significance of overlaps between m⁶A sites and DEGs, spliced genes, or junction regions was determined using a hypergeometric test by R software. The significance of candidate genes in qPCR experiments was determined by T-test.

3 Results

3.1 Increased expression of METTL3 was associated with an aggressive PCa status

Functional enrichment indicated that m⁶A regulators were more associated with mRNA stability and mRNA processing regulation (Figure 1A). To investigate how m⁶A regulates mRNA processing in PCa, we first conducted an analysis of the expression profiles of m⁶A regulators by retrieving the FPKM expression data from the TCGA database, comprising adjacent normal (n=52) and prostate adenocarcinoma (n=499). Most m⁶A “writer” and “reader” proteins were overexpressed in PCa tissues compared to normal tissues (Figure 1B), indicating high m⁶A levels in PCa.

FIGURE 1

Figure 1 High METTL3 expression is associated with an aggressive PCa status. (A) Genes and functional ontology related to m⁶A. (B) Heatmap showing the expression profile of 17 m⁶A regulators in adjacent normal (n=52) and PRAD tissues (n=499). The expression level of METTL3 was analyzed in (C) TCGA PRAD dataset, (D) GSE21034, (E) GSE74367, and (F) PCa cell lines (GSE35401). The error bars indicate the variation between the biological replicates in each dataset. (G) The difference in METTL3 expression between the low Gleason score (≤7) and high Gleason score (>7) groups. Overall survival (H), disease-free survival (I), progression-free interval (J), and biochemical recurrence (K) curves for the high and low METTL3 groups based on the median values. *p<0.05; **p<0.01; ***p <0.001; ****p <0.0001.

In particular, we focused on METTL3, the catalytic submit for m⁶A modifications (7), and its expression patterns in various stages of PCa. We observed that the increase in METTL3 expression between normal and primary groups was statically significant but relatively modest in both TCGA and GSE21034 datasets (Figures 1C, D). However, a more substantial change in METTL3 was evident when comparing the primary and metastatic groups (GSE74367) (Figure 1E), supposing a potential role of METTL3 in PCa metastasis. The overexpression of METTL3 was also observed in PCa cell lines LNCaP (androgen-sensitive, AR⁺) and DU145 (androgen insensitive, AR^-) compared to normal prostate cells (RWPE-1) (Figure 1F, GSE35401) (36). To explore the clinical significance of METTL3 in PCa, we found that the higher METTL3 level was positively associated with a high Gleason score (>7) (Figure 1G). Furthermore, overrepresented METTL3 was linked to aggressive survival stages, including overall survival (OS), disease-free survival (DFS), progression-free interval (PFI), and biochemical recurrence (Figures 1H–K). Collectively, these findings suggested a pro-oncogenic role of METTL3 in PCa progression.

3.2 Elucidating functional pathways influenced by METTL3 depletion

To gain further insight into the molecular mechanism by METTL3 in PCa, we performed RNA-Seq in DU145 (AR^-) and LNCaP (AR⁺) cells with METTL3 knockdown. First, the knockdown efficiency of siRNA against METTL3 was assessed in the two PCa cell lines by qPCR and western blotting and shown in Figures 2A, B. In parallel with qPCR and western blotting assessment, the successful knockdown of METTL3 was also confirmed by RNA-Seq data (Figure 2C). Subsequently, we evaluated genes potentially regulated by METTL3 using DESeq2. It was found that 2078 DEGs in DU145 cells between siMETTL3 and siControl groups, of which 1171 genes were upregulated and 907 genes were downregulated (padj <0.05, Figure 2D). In LNCaP cells, we found 1,263 increased genes and 958 decreased genes along with METTL3 knockdown (padj <0.05, Figure 2D). The data revealed that METTL3 depletion has a global impact on gene expression.

FIGURE 2

Figure 2 Functional annotation of DEGs of the two PCa cell lines. (A, B) The relative mRNA and protein level of METTL3 in siControl and siMETTL3 groups of LNCaP and DU145 cells. (C) The normalized counts of METTL3 in siControl and siMETTL3 groups of RNA-Seq data. (D) Volcano plot of DEGs in the comparison of the siMETTL3 against siControl groups of DU145 and LNCaP cells. The X-axis shows log2fold changes in expression and the Y-axis represents –log10 (adjusted p-value). DEGs with a padj < 0.05 & abs|log2FC| >0.5 were considered significantly differentially expressed. (E, F) Venn plot indicated the overlapped genes of upregulated and downregulated DEGs of the two PCa cells. (G, H) KEGG functional annotation of overlapped upregulated or downregulated genes. (I) The relative mRNA levels of candidate genes in the two cell lines upon METTL3 depletion. siCT, siControl; siM3, siMETTL3; *p<0.05; **p<0.01; ***p<0.001.

We next explored the consistency between LNCaP and DU145 cells in response to METTL3 knockdown. There were 458 genes commonly upregulated, and 224 genes downregulated, respectively (Figures 2E, F). KEGG enrichment of these DEGs indicated that METTL3 was negatively associated with protein processing in endoplasmic reticulum, while favorably connected to various cancer-related pathways, such as HIF-1 signaling, endocrine resistance, FoxO signaling, and proteoglycans (Figures 2G, H). We selected eight known cancer driver genes and PCa-survival genes (37, 38) and then validated the expression change using qPCR. As a result, the expression changes of these eight METTL3-regulated genes showed a consistent trend with RNA-Seq data, demonstrating the accuracy of our RNA-Seq data (Figure 2I). Moreover, these eight genes were identified as having at least one annotated m⁶A site in public miCLIP data (39, 40) (Figure S2), supposing the potential regulation by m⁶A. However, more research is required to determine the potential mechanisms involving these genes (e.g., MSI1, NPEPL1).

More interestingly, GSVA analysis also showed that many cancer hallmarks were dysregulated between siControl and siMETTL3 groups, with a noticeable unique pattern in LNCaP (AR⁺) and DU145 (AR^-) cells (Figures S3A, C). For instance, immune (e.g., IFN gamma/alpha response) and metabolism-related (e.g., cholesterol and adipogenesis) signatures were less active in the siMETTL3 group of DU145 cells (Figure S3B). In contrast, metabolism pathways such as cholesterol homeostasis and adipogenesis were more activated in the siMETTL3 LNCaP group (Figure S3D). This may be relevant to cell-specific genomic properties and AR dependence. Previous studies have demonstrated that AR controls critical genes involved in glucose metabolism and lipid metabolism (41). Taken together, the data imply that METTL3 may be engaged in multiple common cancer-related pathways, thus performing pro-oncogenic activities in PCa cells. On the other hand, METTL3 may have partially dissimilar functions in AR⁺ and AR^- PCa cells, especially in terms of metabolism.

3.3 Depletion of METTL3 impacts AS landscape

To gain an insight into the role of METTL3 in AS regulation, we applied direct RNA-sequencing (DRS) performed using the Oxford Nanopore long-read sequencing platform to generate METTL3-dependent AS profiles in LNCaP (AR⁺) and DU145 (AR^-) cells. Overall, 0.6-1.5 million reads were generated for each sample, with an average read length of 1.0 kb and average read quality of 11 (Figure S4; Table S2). Approximately 95% and 99% of the passed reads were aligned to the human genome hg38 and transcriptome (GENCODE v32), respectively (Table S2).

Then, we used FLAIR (17) software to identify the differential AS events upon METTL3 knockdown. Four types of AS events were obtained, including exon skipping (ES), alternative 5’/3’-splice site (A5SS/A3SS), and intron retention (IR). A total of 426 and 172 differential AS events were identified in DU145 and LNCaP cells, respectively (Figure 3A). The predominant AS type in the two cells was the ES type.

FIGURE 3

Figure 3 The global overview of differential AS events in DU145 and LNCaP cells. (A) Pie chart of differential AS events between siControl and siMETTL3 groups of DU145 and LNCaP cells. (B) The top 20 functional annotations and (C) ENCODE TF terms of AS-associated genes in DU145 (left) and LNCaP (right) cells. (D) Venn plot showing the overlapped genes between DEGs and AS-related genes in DU145 and LNCaP cells. (E) The intersection between differentially expressed isoforms (DEI) and AS harboring genes. Significance was assessed by a hypergeometric test. (F, H) Schematic diagram of isoforms of two known spliced genes (MKNK2 and CDKN2A) in DU145 or LNCaP cells. Stacked bar plots showing the relative proportion and the corresponding read count of each isoform and isoforms were labeled with different colors. (G, I) Bar plots displayed the relative mRNA expression of two main isoforms of MKNK2 and CDKN2A genes. siCT, siControl; siM3, siMETTL3. *p<0.05; **p<0.01.

The functional enrichment revealed that differential AS events were mainly associated with eukaryotic translation initiation, cellular response to stimuli (e.g., heat, chemical), metabolic pathways, and proliferation hallmarks (MYC targets v1, E2F targets) (Figure 3B). Additionally, transcription factor (TF) analysis showed that differential AS events were highly associated with MYC (Figure 3C), supposing the potential role of MYC in AS regulation.

3.4 METTL3 regulates isoforms switching

As we know that AS is usually coupled with transcriptional processes and thus influences gene expression (42). We found that the functional annotation of METTL3-regulated spliced genes was in line with those discovered in DEGs. However, when intersected AS genes with previously identified DEGs, we noted that only 9-15% of AS genes overlapped (Figure 3D), suggesting that AS dysregulation may have a minor influence on global gene expression but rather functionally tuned transcriptomes (43).

Then, using FLAIR (17), we inspected the isoform-level profiles to explore the functional implications of AS dysregulation on PCa transcriptome. We discovered that, as expected, the majority of differentially expressed isoforms (DEIs) were generated due to aberrant splicing events (Figure 3E). For example, two known spliced genes (e.g., MKNK2 and CDKN2A) were selected to show distinctive isoforms switching with METTL3 knockdown (Figures 3F–I) (Figure S5-6). Three isoforms of the MKNK2 gene were observed, with a proportionate rise in the shorter transcript ENST00000309340 (MKNK2b) and a comparable reduction in the longer transcript (MKNK2a) upon METTL3 depletion (Figure 3F). We validated the isoform switching of MKNK2 using qPCR and found a similar reduction of MKNK2a with METTL3 knockdown (Figure 3G). Another alternatively spliced gene, CDKN2A, has been found to express two major isoforms, of which the CDKN2A-4 isoform (ENST00000579755) was relatively decreased under the DU145 siMETTL3 condition (Figure 3H) (Figure S5). The qPCR experiment also verified a similar decrease in CDKN2A-4 and a slight rise in the CDKN2A-1 isoform upon METTL3 depletion (Figure 3I). However, more research has to be done to determine the precise mechanisms of these isoforms and how they relate to PCa development. Nevertheless, our findings suggest that METTL3-induced splicing abnormalities may influence PCa biology in part via isoform switching of key genes.

3.5 METTL3 regulates intron retention and exon-skipping events

We previously found that LNCaP (AR⁺) and DU145 (AR^-) cells have different functional enrichments upon METTL3 depletion (Figure S3). Therefore, we further investigated the overlap of AS-related genes in the two PCa cells to see if they also with divergent phenotypes. As a result, there was only limited overlap, indicating cell specificity (Figure S7). This is to be expected, given the transcriptome changes drastically as the disease progresses.

Further, we analyzed the different effects of METTL3 on AS events in DU145 and LNCaP cells and we mainly focused on ES and IR types. In the comparison of siControl versus siMETTL3 of DU145 cells, we observed a preferential increase in upregulated ES (155 up/89 down) and downregulated IR events (11 up/62 down) (Figure 4A). On the other hand, METTL3 knockdown led to the opposite effect in LNCaP cells. That is fewer upregulated ES events (34 up/44 down) and downregulated IR events (15 up/10 down) (Figure 4B). To corroborate this finding, we measured the inclusion levels (PSI values) of differential IR and ES events between siControl and siMETTL3 groups using RNA-Seq data. Consistently, the opposing changes in PSI values of IR and ES events were found in DU145 and LNCaP cells (Figures 4C, D; Figure S8). This indicated an increase in retained introns and skipped exons in the DU145 siMETTL3 group, whereas the LNCaP siMETTL3 group had the opposite tendency. For instance, the transcript (ENST00000392350) of the ORMDL1 gene with exon 2 skipping showed an increased level in the DU145 siMETTL3 group but decreased in the LNCaP siMETTL3 group (Figure 4E; Figure S9). ORMDL1 has been reported to be associated with ceramide biosynthesis and as a regulator of sphingolipid levels implicated in cell proliferation as well as migration (44).

FIGURE 4

Figure 4 Depletion of METTL3 influences intron retention and exon skipping events in DU145 and LNCaP cells. (A, B) The bar plot of up and downregulated AS events shows a significant difference between siControl and siMETTL3 groups of DU145 and LNCaP cells. (C, D) Box plots showing PSI values of differential IR and ES events in DU145 and LNCaP cells, respectively. (E) Schematic diagram of isoforms from ORDML1 gene in DU145 or LNCaP cells. Stacked bar plots showing the relative proportion and the corresponding read count of each isoform and isoforms were labeled with different colors. (F) The intersection between METTL3 pulled-down proteins and known SFs. Significance was assessed by a hypergeometric test. Several typical intersecting genes were marked in red. (G) Predicted sequence motifs of differential IR and ES events in DU145 and LNCaP cells. siCT, siControl; siM3, siMETTL3; SFs, splicing factors.

RNA splicing is a complicated process that is coregulated by multiple factors, including splice sites, cis-elements, and splicing factors (SFs) (42). Next, we intended to explore whether the distinct influence on IR and ES in AR⁺ and AR^- cells is associated with SFs. By analyzing the proteomic data of METTL3 reported by Yue et al. (45), we discovered a substantial overlap between METTL3 pulled-down proteins and known SFs (e.g., SRSF2, U2AF2, KHDRBS3, KHSRP) (P=1.2e-34, hypergeometric test), suggesting a possible interaction between METTL3 and SFs (Figure 4F). Furthermore, motif analysis revealed that SRSF1 (P=1e-13) and KHRBSH1 (P=1e-14) were enriched in differential ES and IR events of DU145 cells, whereas the SRSF2 (P=1e-8) motif was detected in differential ES and IR events for LNCaP cells (Figure 4G), indicating that different mechanisms may be involved in the two PCa cell lines. Taken together, we found that METTL3 depletion resulted in divergent phenotypes for IR and ES events in DU145 and LNCaP cells, and the regulatory effects of METTL3 may be associated with different SF proteins.

3.6 Differential m⁶A profiles

To understand whether METTL3 modulates RNA splicing through m⁶A, we first analyzed differential m⁶A sites from DRS data using ELIGOS2 software with strict criteria, focusing primarily on decreased m⁶A sites (32)(seeing Methods). In DU145 and LNCaP cells, we found 15806 and 16940 potentially decreased m⁶A sites, respectively, with 2506 m⁶A-modified genes shared by the two cell lines (Figure 5A; Tables S3, 4). Moreover, the core residues of the 5-kmer motifs around these m⁶A sites were AC (Figure 5B), which was also consistent with prior results (46). Aside from the previously reported major enrichment in the 3’UTR, we observed that identified m⁶A sites in PCa cells partially fell into exonic (16.0% or 18.1%) and intronic regions (10.3% or 11%), indicating the potential involvement of m⁶A in IR and ES events (Figure 5C). Further, cancer hallmark enrichment of m⁶A-modified genes revealed a high association with “Myc target V1/V2”, “mTORC1 signaling”, and several metabolic pathways (Figure 5D). Remarkably, this finding also coincided with previous enrichment by altered genes or AS events.

FIGURE 5

Figure 5 Decreased m⁶A profiles in DU145 and LNCaP cells with METTL3 depletion. (A) Venn diagram of the overlapping m⁶A-modified genes between DU145 and LNCaP cells. (B) The 5-kmer motif of differential m⁶A sites in DU145 and LNCaP cells after METTL3 knockdown. (C) The genomic feature distributions of decreased m⁶A sites in DU145 and LNCaP cells. (D) Cancer hallmark enrichments of m⁶A modified genes in DU145 and LNCaP cells, respectively. Only significant items with a p-value <0.05 are shown in the plot. (E–G) Examples of predicted m⁶A levels in MYC, TK1, and MTHFD2 genes. The top panel represented the proportion of predicted m⁶A levels in each gene under control (siCT) and METTL3 knockdown (siM3) condition in DU145 and LNCaP cells. Y-axis represents the m⁶A modification percentage, which is the percent error of specific bases (%ESB) identified by ELIGOS2 software. The bottom panel indicated the longest isoform of each gene.

We next investigated the relationship between decreased m⁶A sites and previously identified DEGs or splicing events. As a result, there were 443 (22.5%) and 280 (13.8%) overlapping genes between DEGs and m⁶A-modified genes in DU145 and LNCaP cells, respectively (P=1.87e-107 and P=2.76e-48, hypergeometric test) (Figure S10A). Additionally, we observed that 182 of 235 (77.4%) spliced genes in DU145 cells and 58 out of 79 (73.4%) genes in LNCaP cells intersected with m⁶A-modified genes (P=1.08e-158 and P=3.65e-56, hypergeometric test) (Figure S10B). The predicted m⁶A levels of three selected examples, namely MYC, TK1, and MTHFD2 genes were presented in Figures 5E–G. Notably, MYC has been previously shown to be regulated by METTL3 in PCa in an m⁶A-dependent manner (9). In addition, our analysis using the DRUMMER software also revealed intriguing isoform-specific differences in m⁶A levels (Tables S5, 6; Figure S11). These data suggest the function of METTL3 in the regulation of gene expression and RNA splicing in PCa cells might be mediated by m⁶A.

3.7 METTL3 modulates MKNK2 isoform expression by SRSF1 and m⁶A

Based on previous functional enrichment analysis (Figure 3B), we noted that METTL3 was primarily associated with alternative splicing of genes linked to translation initiation. Among these genes, MNK2 (encoded by MKNK2) is the known kinase responsible for phosphorylating eukaryotic translation initiation factor 4E (eIF4E) at serine 209, an essential molecule during the mRNA translation phase (47). Therefore, we further interrogate whether the spliced event of MKNK2 that was validated in PCa cell lines was clinical significance. For MKNK2, two alternative isoforms were analyzed (namely MKNK2a and MKNK2b). The whole mRNA and protein levels of MKNK2 were evaluated in PCa, suggesting its onco-prognosis in PCa (Figures 6A, B). Moreover, we explored the relationship between MKNK2a/b isoforms and the prognosis of PCa patients. The data showed that a higher MKNK2a or a lower MKNK2b had a worse prognosis (PFI survival) (Figures 6C, D) and more biochemical recurrences (Figure S12), which indicated MKNK2 isoforms might play roles in PCa progression.

FIGURE 6

Figure 6 The potential mechanism of METTL3 in modulating MKNK2 splicing in PCa. (A, B) MKNK2 mRNA and protein levels in normal and PCa tissues. The immunohistochemical results were obtained from the human protein atlas. (C, D) Survival analysis showed a significant positive/negative correlation between MKNK2a/b isoforms and progression-free survival (PFI) in PCa patients. Patients were divided into high and low PSI groups based on the 80^th percentile. The below represented the splicing pattern indicated by the red line. (E, F) The correlation between the PSI index of MKNK2a/b isoforms and METTL3 expression in PCa, respectively. (G) The violin plot of PSI index of MKNK2a/b isoforms in METTL3 high and low groups, indicating increased expression of MKNK2a isoforms and decreased expression of MKNK2b isoforms in METTL3 high group. (H) The PSI index of the MKNK2 AS event (generating the MKNK2a isoform) in siMETTL3 or STM2457 (METTL3 inhibitor) C42 cells was reduced compared to the control group. (I) The top two panels represented the proportion of predicted m⁶A levels of MKNK2 in control (siCT) and METTTL3 knockdown (siM3) conditions. Y-axis represents the m⁶A modification percentage, which is the percent error of specific bases (%ESB) identified by ELIGOS2 software. The third and fourth panels represented the read coverage of m⁶A modification levels on the MKNK2 gene in normal prostate epithelial cells (RWPE-1) and PCa cells (LNCaP) from public sequencing data. Among them, the regions highlighted by four colors are regions with differential m⁶A modification levels in the exon 14a/b region of MKNK2. The bottom shows the structural features of MKNK2a/b isoforms. (J) The correlation between the PSI index of MKNK2a/b isoforms and 1-12 SRSFs expression in PCa, respectively. The red dot represented the significant correlation. (K) The correlation between the PSI index of MKNK2a/b isoforms and SRSF1 expression in PCa, respectively. siCT, siControl; siM3, siMETTL3; *p<0.05; **p<0.01; ***p <0.001; ****p<0.0001.

We further validated the relationship between METTL3 and MKNK2a/b isoforms using TCGA Spliceseq datasets. Correlation analysis showed that there was a significant positive/negative correlation between the expression of METTL3 and MKNK2a/b AS events (Figures 6E, F). Consistently, we observed that the AS events for MKNK2a displayed a significantly higher PSI value in METTL3 high patients than in the lower group, while the PSI values were significantly lower for MKNK2b (Figure 6G). These results suggest that METTL3 participated in the splicing process of the MKNK2 gene in PCa patients. Moreover, we also verified that METTL3 depletion (METTL3 inhibitor, STM2457) led to aberrant splicing of MKNK2 using two external datasets (Figure 6H) (48), which was also consistent with our previous findings. Collectively, these data indicated the role of METTL3 in modulating MKNK2 AS events in PCa.

The previous findings indicated that METTL3 regulates the splicing process associated with splicing factors and m⁶A sites. Our nanopore DRS data showed that after METTL3 knockdown, the m⁶A modification level in the MKNK2 exon 14a/b region was downregulated (Figure 6I). At the same time, we analyzed the public MeRIP-seq data and found that compared with normal prostate cells (RWPE-1), the exon 14a/b region of the MKNK2 gene in PCa cells has a higher methylation level (Figure 6I). Therefore, we speculate that the m⁶A modification in the exon 14a/b region of the MKNK2 gene may affect the binding of splicing factors to MKNK2. We next aimed to identify the potential splicing factors targeting MKNK2. Through the correlation test, we found that several SFs were positively correlated with the level of MKNK2a and negatively correlated with the MKNK2b level (Figure 6J). It’s known that the splicing process of MKNK2 pre-mRNA is regulated by SRSF1 (49). Here, we found a significant correlation between the expression of SRSF1 and MKNK2a/b isoforms in PCa (Figure 6K). However, RNA-Seq results showed that METTL3 knockdown did not affect the expression level of most SRSFs including SRSF1 (Figure S13). Previous studies indicated that m⁶A modification near splicing sites may affect the binding of splicing factors to transcripts. Therefore, we speculate that the m⁶A modification in the exon 14a/b region of the MKNK2 gene may affect the binding of splicing factors such as SRSF1 to MKNK2. The specific regulatory mechanism remains to be studied.

4 Discussion

Although mechanisms of METTL3 in many cancers have been well documented, the role of METTL3 in PCa-related AS switch has been poorly explored. Given the inherent differences in each cancer type, a tailored analysis focusing on PCa is still worth exploring and important for drug development. Our study illustrates some findings that are consistent with those in other cancers (e.g., the association between METTL3 and MYC signaling) (50–53), but we did find many PCa-specific findings that differ from other cancer types.

Prior research has found that high METTL3 expression can promote PCa cell proliferation, survival, and invasion in vitro and in vivo (9, 10, 54–56), however, our understanding of its role in AS regulation is still limited. Thus, we performed RNA-Seq to evaluate the function of METTL3 depletion in PCa. The results showed that METTL3 knockdown resulted in alterations in many cancers and cell proliferation-related genes, such as the FOXO signaling pathway, G2M Checkpoint, E2F Targets, and DNA repair hallmark. This was consistent with results in other cancers. In liver cancer, METTL3 can stabilize FOXO3 mRNA, reducing autophagy and sorafenib resistance (57). Moreover, the METTL3-METTL14 complex participates in DNA repair processing by aggregating at UV-exposed DNA damage sites (5, 58). Moreover, we also verified the expression changes of some genes by qPCR, which showed the same trend with RNA-Seq data and was in line with existing evidence. MSI1 is characterized as an RNA-banding protein (RBP) by repressing the mRNA translation and is associated with cancer stem cell properties of various cancers (59). A recent report has implicated that MSI1 regulated the expression of m⁶A reader YTHDF1, thus involving glioblastoma cell proliferation and migration (59). Interestingly, a preprinted study by Wei et al. indicated that IGF2BP3, another m⁶A reader, destabilized the m⁶A-harboring NF1 mRNA to accelerate triple-negative breast cancer progression, which was also consistent with our data (60). A similar reduction in NPEPL1 protein has been noted in another METTL3 knockdown proteomics dataset in LNCaP cells (39), the detailed mechanism between METTL3 and NPEPL1 in PCa remains unclear and needs to be explored.

We next focused on genome-wide METTL3-regulated AS events in PCa cell lines. Our results revealed that METTL3 primarily regulated the AS switch associated with cell cycle and translation initiation in PCa, and that AS regulation may be associated with MYC signaling. C. Achour et al. reported that METTL3 may indirectly regulate AS events through MYC in breast cancer (61). In addition, studies by Yuan et al. (9) and C. Achour et al. (61) also showed that METTL3 knockdown can suppress MYC protein levels, respectively, which is also common in several other cancers, such as AML, gastric cancer, colorectal cancer, and bladder cancer (50–53). These results suggest that the regulation of carcinogenesis by METTL3 through MYC is a common phenomenon and that METTL3 regulates AS events through a similar mechanism in prostate and breast cancer. In addition, we explored the association of AS alterations with the regulation of gene expression. However, our findings showed that only 9-15% of AS-affected genes overlapped with DEGs, with the majority being related to isoform-switched genes. This suggests that in PCa cells, AS dysregulation may have little effect on global gene expression but functionally tailored transcriptomes instead (43). Interestingly, we found that functional pathways enriched by AS genes (e.g., metabolic pathways, MYC, and E2F target markers) were coordinated with those in METTL3-affected DEGs. These results indicated that METTL3 jointly regulates PCa biological processes through multi-layered gene regulation, and has a broad and complex mechanism of action, which requires further in-depth study of its coordinated mechanism.

Notably, we found that after METTL3 depletion, two PCa cells had different regulatory responses on ES and IR events. Knockdown of METTL3 in AR^- cells (DU145) resulted in more introns being preserved and more exons being skipped, whereas the opposite phenomenon was observed in AR⁺ cells (LNCaP). It has been shown that m⁶A promotes exon inclusion or exon exclusion. A breast cancer study showed that depletion of METTL3 causes exons of specific genes to be retained (61). Conversely, another study showed that m⁶A promotes exon inclusion, and depletion of METTL3 leads to exon skipping (62). The differential effects of m⁶A on AS appear to be cell- or transcript-specific and depend on specific splicing factors. Studies have shown that m⁶A reader YTHDC1 binds to m⁶A sites and recruits the splicing factor SRSF3 to promote exon inclusion (15). A parallel study showed that the accumulation of m⁶A markers following FTO depletion in mouse preadipocytes promoted the binding of another splicing factor, SRSF2, leading to increased inclusion of target exons (16). Therefore, we also hypothesized that there is heterogeneity in the regulation of AS by METTL3 in AR^- or AR⁺ cells, depending on the cellular context and specific splicing factors. Exclusion of certain exons has been shown to be controlled by antagonism of the SRSF1/SRSF2 proteins, with SRSF1 inhibiting exon skipping and SRSF2 activating exon skipping (47). Our results showed that SRSF1 or SRSF2 were enriched in DU145 and LNCaP, respectively. These studies suggest that METTL3 may play different roles in RNA splicing in DU145 and LNCaP cells by interacting with different partners. However, the regulatory network of RNA splicing is complex and dynamic, and the mechanism of METTL3 in the two PCa cell lines requires more studies.

In this study, we took the MKNK2 gene as an example to explore the clinical significance of AS dysregulation in PCa. Our previous results showed that METTL3 mainly regulates alternative splicing of genes related to translation initiation. MNK1/2 are the only known kinases responsible for phosphorylating eukaryotic translation initiation factor 4E (eIF4E) at serine 209, which is a key rate-limiting molecule in the mRNA translation phase (63). Phosphorylation of eIF4E by MNKs proteins enhances its binding to the 5’cap structure of mRNA, thereby enhancing the cap-dependent translation process (64). MNK2/eIF4E signaling pathway is associated with the translation of tumorigenesis-associated mRNAs, contributing to cell proliferation, invasion, and drug resistance (49, 63, 65–69). Our results found that MKNK2 expression was higher in PCa patients than in normal samples. Moreover, splicing variants of MKNK2a were significantly associated with poor PFI in PCa patients and significantly positively correlated with METTL3 expression. These results suggest that METTL3 was also associated with AS of MKNK2 in PCa clinical samples, and the MKNK2a subtype may play a role in promoting cancer progression. At the same time, we also verified that METTL3 knockdown led to abnormal splicing of MKNK2 and reduced expression of MKNK2a isoforms using full-length transcriptome and qPCR. Therefore, METTL3 may affect the translation process of oncogenes by regulating the AS of the translation initiation-related gene MKNK2.

Finally, we explored the possible mechanism of METTL3 in regulating MKNK2 splicing. The variable splicing process of genes is controlled by a variety of regulatory elements, mainly including splicing factors and splicing sites (42). It is known that MKNK2 pre-mRNA is one of the target molecules of the splicing factor SRSF1 (49), which is proven to modulate the MKNK2a-MKNK2b isoforms switch. In our data, a significant correlation was found between the expression of SRSF1 and MKNK2a/b isoforms in PCa, which is consistent with the literature. Furthermore, the binding motif of SRSF1 was significantly enriched in METTL3-regulated splicing events. These results suggest that SRSF1 may be involved in the regulation of MKNK2 alternative splicing by METTL3. However, RNA-Seq results showed that METTL3 knockdown did not affect the expression level of SRSF1. Existing studies have shown that m⁶A modification near the cleavage site can affect the binding of cleavage factors to transcripts, thereby changing the way the transcripts are cleaved. Our results also showed that m⁶A sites overlap significantly with splice junctions, suggesting that m⁶A is involved in the regulation of splice site recognition. Further, the data revealed that the methylation level of the exon14a/b region of the MKNK2 gene in PCa cells was significantly higher than that in normal cells. And after METTL3 knockdown, the predicted methylation level of exon exon14a/b region was significantly reduced. Therefore, the declining m⁶A level in MKNK2 caused by METTL3 knockdown potentially reduces the binding of splicing factors. In the future, further experiments are needed to verify the mechanism of METTL3 regulating MKNK2 AS events.

This study has some limitations. First, while offering advantages such as long read lengths and RNA modification analysis, the DRS approach has limitations of higher error rates and more complex bioinformatics analysis. Our results are mainly based on bioinformatic analysis, and further research in experimental cohorts is warranted. It is worth mentioning that our conclusions are supported by other independent datasets and qPCR experiments, demonstrating the validity. Second, we uncovered an intriguing mechanism of METTL3 in AR^- and AR⁺ PCa cells. However, which factors are functionally related to these responses remains unanswered, which may enhance the understanding of the carcinogenesis process of different PCa stages.

In summary, our study used genome-wide analyses to investigate METTL3-regulated AS events in AR⁺ and AR^- PCa cell lines via a combination of standard RNA-Seq and long-reads direct RNA-Seq of Nanopore. We revealed both common and phenotypic heterogeneity of different PCa cells regulated by METTL3, which were potentially associated with m⁶A modification and SFs. The study identified the clinical relevance of MKNK2 AS events in PCa and explored the potential mechanism of METTL3 in modulating MKNK2 AS events. These findings provide new insights into RNA modifications in PCa and could serve as a favorable molecular basis for novel treatment strategies in androgen-sensitive and androgen-resistant prostate cancer.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement

Ethical approval was not required for the studies on humans in accordance with the local legislation and institutional requirements because only commercially available established cell lines were used.

Author contributions

LW: Study design, Data processing, and Writing-Original draft preparation. LYW, YW, and CF: discussed the data and assisted in the sequencing library preparation. LS, YL, JN and CY: validated the differentially expressed or spliced genes. JH, PC, and CN: Writing-Reviewing and Editing. ST: Supervision, Writing-Reviewing and Editing. All authors contributed to the article and approved the submitted version.

Funding

This work was made possible by the General Research Fund from Research Grants Council of Hong Kong (Reference numbers: 14175617, 14119219, 14119420, R4012-18F) and Health and Medical Research Fund from Food and Health Bureau of Hong Kong (Reference numbers: 06171016, 07181266).

Acknowledgments

We thank Professor Andrew Chan’s laboratory from the School of Biomedical Sciences, The Chinese University of Hong Kong, for kindly providing PCa cell lines.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2023.1227016/full#supplementary-material

References

1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: Cancer J Clin (2021) 71(3):209–49. doi: 10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Sharifi N, Gulley JL, Dahut WL. Androgen deprivation therapy for prostate cancer. Jama (2005) 294(2):238–44. doi: 10.1001/jama.294.2.238

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Feldman BJ, Feldman D. The development of androgen-independent prostate cancer. Nat Rev Cancer (2001) 1(1):34–45. doi: 10.1038/35094009

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Roundtree IA, Evans ME, Pan T, He C. Dynamic RNA modifications in gene expression regulation. Cell (2017) 169(7):1187–200. doi: 10.1016/j.cell.2017.05.045

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Xiang Y, Laurent B, Hsu C-H, Nachtergaele S, Lu Z, Sheng W, et al. RNA m6A methylation regulates the ultraviolet-induced DNA damage response. Nature (2017) 543(7646):573–6. doi: 10.1038/nature21671

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Chen X-Y, Zhang J, Zhu J-S. The role of m6A RNA methylation in human cancer. Mol Cancer (2019) 18(1):1–9. doi: 10.1186/s12943-019-1033-z

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Geula S, Moshitch-Moshkovitz S, Dominissini D, Mansour AA, Kol N, Salmon-Divon M, et al. m6A mRNA methylation facilitates resolution of naïve pluripotency toward differentiation. Science (2015) 347(6225):1002–6. doi: 10.1126/science.1261417

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Chen M, Wei L, Law CT, Tsang FHC, Shen J, Cheng CLH, et al. RNA N6-methyladenosine methyltransferase-like 3 promotes liver cancer progression through YTHDF2-dependent posttranscriptional silencing of SOCS2. Hepatology (2018) 67(6):2254–70. doi: 10.1002/hep.29683

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Yuan Y, Du Y, Wang L, Liu X. The M6A methyltransferase METTL3 promotes the development and progression of prostate carcinoma via mediating MYC methylation. J Cancer (2020) 11(12):3588. doi: 10.7150/jca.42338

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Cai J, Yang F, Zhan H, Situ J, Li W, Mao Y, et al. RNA m6A methyltransferase METTL3 promotes the growth of prostate cancer by regulating hedgehog pathway. OncoTargets Ther (2019) 12:9143. doi: 10.2147/OTT.S226796

CrossRef Full Text | Google Scholar

11. Li E, Wei B, Wang X, Kang R. METTL3 enhances cell adhesion through stabilizing integrin β1 mRNA via an m6A-HuR-dependent mechanism in prostatic carcinoma. Am J Cancer Res (2020) 10(3):1012.

PubMed Abstract | Google Scholar

12. Adhikari S, Xiao W, Zhao Y-L, Yang Y-G. m6A: signaling for mRNA splicing. RNA Biol (2016) 13(9):756–9. doi: 10.1080/15476286.2016.1201628

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Wang Y, Liu J, Huang B, Xu YM, Li J, Huang LF, et al. Mechanism of alternative splicing and its regulation. Biomed Rep (2015) 3(2):152–8. doi: 10.3892/br.2014.407

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Feng Z, Li Q, Meng R, Yi B, Xu Q. METTL 3 regulates alternative splicing of MyD88 upon the lipopolysaccharide-induced inflammatory response in human dental pulp cells. J Cell Mol Med (2018) 22(5):2558–68. doi: 10.1111/jcmm.13491

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Xiao W, Adhikari S, Dahal U, Chen Y-S, Hao Y-J, Sun B-F, et al. Nuclear m6A reader YTHDC1 regulates mRNA splicing. Mol Cell (2016) 61(4):507–19. doi: 10.1016/j.molcel.2016.01.012

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Zhao X, Yang Y, Sun B-F, Shi Y, Yang X, Xiao W, et al. FTO-dependent demethylation of N6-methyladenosine regulates mRNA splicing and is required for adipogenesis. Cell Res (2014) 24(12):1403–19. doi: 10.1038/cr.2014.151

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Tang AD, Soulette CM, van Baren MJ, Hart K, Hrabeta-Robinson E, Wu CJ, et al. Full-length transcript characterization of SF3B1 mutation in chronic lymphocytic leukemia reveals downregulation of retained introns. Nat Commun (2020) 11(1):1–12. doi: 10.1038/s41467-020-15171-6

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Workman RE, Tang AD, Tang PS, Jain M, Tyson JR, Razaghi R, et al. Nanopore native RNA sequencing of a human poly (A) transcriptome. Nat Methods (2019) 16(12):1297–305. doi: 10.1038/s41592-019-0617-2

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods (2015) 12(4):357–60. doi: 10.1038/nmeth.3317

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol (2014) 15(12):1–21. doi: 10.1186/s13059-014-0550-8

CrossRef Full Text | Google Scholar

21. Yu G, Wang L-G, Han Y, He Q-Y. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics: J Integr Biol (2012) 16(5):284–7. doi: 10.1089/omi.2011.0118

CrossRef Full Text | Google Scholar

22. Tang Z, Kang B, Li C, Chen T, Zhang Z. GEPIA2: an enhanced web server for large-scale expression profiling and interactive analysis. Nucleic Acids Res (2019) 47(W1):W556–W60. doi: 10.1093/nar/gkz430

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun (2019) 10(1):1523. doi: 10.1038/s41467-019-09234-6

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res (2003) 13(11):2498–504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Garalde DR, Snell EA, Jachimowicz D, Sipos B, Lloyd JH, Bruce M, et al. Highly parallel direct RNA sequencing on an array of nanopores. Nat Methods (2018) 15(3):201–6. doi: 10.1038/nmeth.4577

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Li H. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics (2018) 34(18):3094–100. doi: 10.1093/bioinformatics/bty191

PubMed Abstract | CrossRef Full Text | Google Scholar

27. De Coster W, D’hert S, Schultz DT, Cruts M, Van Broeckhoven C. NanoPack: visualizing and processing long-read sequencing data. Bioinformatics (2018) 34(15):2666–9. doi: 10.1093/bioinformatics/bty149

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Shen S, Park JW, Z-x Lu, Lin L, Henry MD, Wu YN, et al. rMATS: robust and flexible detection of differential alternative splicing from replicate RNA-Seq data. Proc Natl Acad Sci (2014) 111(51):E5593–E601. doi: 10.1073/pnas.1419161111

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Jiang M, Zhang S, Yin H, Zhuo Z, Meng G. A comprehensive benchmarking of differential splicing tools for RNA-se q analysis at the event level. Briefings Bioinf (2023) 24(3):bbad121. doi: 10.1093/bib/bbad121

CrossRef Full Text | Google Scholar

30. Katz Y, Wang ET, Airoldi EM, Burge CB. Analysis and design of RNA sequencing experiments for identifying isoform regulation. Nat Methods (2010) 7(12):1009–15. doi: 10.1038/nmeth.1528

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Zhang B, Yang L, Wang X, Fu D. Identification of survival-related alternative splicing signatures in acute myeloid leukemia. Bioscience Rep (2021) 41(7):BSR20204037. doi: 10.1042/BSR20204037

CrossRef Full Text | Google Scholar

32. Jenjaroenpun P, Wongsurawat T, Wadley TD, Wassenaar TM, Liu J, Dai Q, et al. Decoding the epitranscriptional landscape from native RNA sequences. Nucleic Acids Res (2021) 49(2):e7–e. doi: 10.1093/nar/gkaa620

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Abebe JS, Price AM, Hayer KE, Mohr I, Weitzman MD, Wilson AC, et al. DRUMMER—rapid detection of RNA modifications through comparative nanopore sequencing. Bioinformatics (2022) 38(11):3113–5. doi: 10.1093/bioinformatics/btac274

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Hong A, Kim D, Kim VN, Chang H. Analyzing viral epitranscriptomes using nanopore direct RNA sequencing. J Microbiol (2022) 60(9):867–76. doi: 10.1007/s12275-022-2324-4

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Yu G, Wang L-G, He Q-Y. ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization. Bioinformatics (2015) 31(14):2382–3. doi: 10.1093/bioinformatics/btv145

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Lee BH, Taylor MG, Robinet P, Smith JD, Schweitzer J, Sehayek E, et al. Dysregulation of cholesterol homeostasis in human prostate cancer through loss of ABCA1Loss of ABCA1 promotes prostate cancer. Cancer Res (2013) 73(3):1211–8. doi: 10.1158/0008-5472.CAN-12-3128

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Zhang B, Zhang Z, Hu H, Dong J-T. Novel gene signatures predictive of patient recurrence-free survival and castration resistance in prostate cancer. Cancers (2021) 13(4):917. doi: 10.3390/cancers13040917

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Yarmishyn AA, Yang Y-P, Lu K-H, Chen Y-C, Chien Y, Chou S-J, et al. Musashi-1 promotes cancer stem cell properties of glioblastoma cells via upregulation of YTHDF1. Cancer Cell Int (2020) 20(1):1–15. doi: 10.1186/s12935-020-01696-9

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Cotter KA, Gallon J, Uebersax N, Rubin P, Meyer KD, Piscuoglio S, et al. Mapping of m6A and its regulatory targets in prostate cancer reveals a METTL3-low induction of therapy resistance. Mol Cancer Res (2021) 19(8):1398–411. doi: 10.1158/1541-7786.MCR-21-0014

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Zheng Y, Nie P, Peng D, He Z, Liu M, Xie Y, et al. m6AVar: a database of functional variants involved in m6A modification. Nucleic Acids Res (2018) 46(D1):D139–D45. doi: 10.1093/nar/gkx895

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Massie CE, Lynch A, Ramos-Montoya A, Boren J, Stark R, Fazli L, et al. The androgen receptor fuels prostate cancer by regulating central metabolism and biosynthesis. EMBO J (2011) 30(13):2719–33. doi: 10.1038/emboj.2011.158

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Naftelberg S, Schor IE, Ast G, Kornblihtt AR. Regulation of alternative splicing through coupling with transcription and chromatin structure. Annu Rev Biochem (2015) 84:165–98. doi: 10.1146/annurev-biochem-060614-034242

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Zhang D, Hu Q, Liu X, Ji Y, Chao H-P, Liu Y, et al. Intron retention is a hallmark and spliceosome represents a therapeutic vulnerability in aggressive prostate cancer. Nat Commun (2020) 11(1):1–19. doi: 10.1038/s41467-020-15815-7

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Zhu T, Chen Y, Min S, Li F, Tian Y. Expression patterns and prognostic values of ORMDL1 in different cancers. BioMed Res Int (2020) 2020. doi: 10.1155/2020/5178397

CrossRef Full Text | Google Scholar

45. Yue Y, Liu J, Cui X, Cao J, Luo G, Zhang Z, et al. VIRMA mediates preferential m6A mRNA methylation in 3′ UTR and near stop codon and associates with alternative polyadenylation. Cell Discovery (2018) 4(1):1–17. doi: 10.1038/s41421-018-0019-0

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Dominissini D, Moshitch-Moshkovitz S, Schwartz S, Salmon-Divon M, Ungar L, Osenberg S, et al. Topology of the human and mouse m6A RNA methylomes revealed by m6A-seq. Nature (2012) 485(7397):201–6. doi: 10.1038/nature11112

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Pandit S, Zhou Y, Shiue L, Coutinho-Mansfield G, Li H, Qiu J, et al. Genome-wide analysis reveals SR protein cooperation and competition in regulated splicing. Mol Cell (2013) 50(2):223–35. doi: 10.1016/j.molcel.2013.03.001

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Haigh DB, Woodcock CL, Lothion-Roy J, Harris AE, Metzler VM, Persson JL, et al. The METTL3 RNA methyltransferase regulates transcriptional networks in prostate cancer. Cancers (2022) 14(20):5148. doi: 10.3390/cancers14205148

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Liu H, Gong Z, Li K, Zhang Q, Xu Z, Xu Y. SRPK1/2 and PP1α exert opposite functions by modulating SRSF1-guided MKNK2 alternative splicing in colon adenocarcinoma. J Exp Clin Cancer Res (2021) 40(1):75. doi: 10.1186/s13046-021-01877-y

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Su R, Dong L, Li C, Nachtergaele S, Wunderlich M, Qing Y, et al. R-2HG exhibits anti-tumor activity by targeting FTO/m6A/MYC/CEBPA signaling. Cell (2018) 172(1-2):90–105.e23. doi: 10.1016/j.cell.2017.11.031

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Yang D-D, Chen Z-H, Yu K, Lu J-H, Wu Q-N, Wang Y, et al. METTL3 promotes the progression of gastric cancer via targeting the MYC pathway. Front Oncol (2020) 10:115. doi: 10.3389/fonc.2020.00115

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Xiang S, Liang X, Yin S, Liu J, Xiang Z. N6-methyladenosine methyltransferase METTL3 promotes colorectal cancer cell proliferation through enhancing MYC expression. Am J Trans Res (2020) 12(5):1789.

Google Scholar

53. Cheng M, Sheng L, Gao Q, Xiong Q, Zhang H, Wu M, et al. The m6A methyltransferase METTL3 promotes bladder cancer progression via AFF4/NF-κB/MYC signaling network. Oncogene (2019) 38(19):3667–80. doi: 10.1038/s41388-019-0683-z

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Chen Y, Pan C, Wang X, Xu D, Ma Y, Hu J, et al. Silencing of METTL3 effectively hinders invasion and metastasis of prostate cancer cells. Theranostics (2021) 11(16):7640. doi: 10.7150/thno.61178

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Li J, Xie H, Ying Y, Chen H, Yan H, He L, et al. YTHDF2 mediates the mRNA degradation of the tumor suppressors to induce AKT phosphorylation in N6-methyladenosine-dependent way in prostate cancer. Mol Cancer (2020) 19(1):1–18. doi: 10.1186/s12943-020-01267-6

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Lothion-Roy J, Haigh DB, Harris AE, Metzler VM, Alsaleem M, Toss MS, et al. Clinical and molecular significance of the RNA m6A methyltransferase complex in prostate cancer. Front Genet (2023) 13. doi: 10.3389/fgene.2022.1096071

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Lin Z, Niu Y, Wan A, Chen D, Liang H, Chen X, et al. RNA m6A methylation regulates sorafenib resistance in liver cancer through FOXO3-mediated autophagy. EMBO J (2020) 39(12):e103181. doi: 10.15252/embj.2019103181

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Yu D, Horton JR, Yang J, Hajian T, Vedadi M, Sagum CA, et al. Human MettL3-MettL14 RNA adenine methyltransferase complex is active on double-stranded DNA containing lesions. Nucleic Acids Res (2021) 49(20):11629–42. doi: 10.1093/nar/gkab460

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Yarmishyn AA, Yang Y-P, Lu K-H, Chen Y-C, Chien Y, Chou S-J, et al. Musashi-1 promotes cancer stem cell properties of glioblastoma cells via upregulation of YTHDF1. Cancer Cell Int (2020) 20(1):597. doi: 10.1186/s12935-020-01696-9

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Ding Q, Zhang X, Shi L, Sun H, Wang Z-W, Xu F, et al. IGF2BP3 Mediates the mRNA Degradation of NF1 to Promote Triple-Negative Breast Cancer Progress via an m6A-Dependent Manner. Res Square (2022). doi: 10.21203/rs.3.rs-1752002/v1

CrossRef Full Text | Google Scholar

61. Achour C, Bhattarai DP, Groza P, Román Á-C, Aguilo F. METTL3 regulates breast cancer-associated alternative splicing switches. Oncogene (2023) 42(12):911–25. doi: 10.1038/s41388-023-02602-z

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Xu K, Yang Y, Feng G-H, Sun B-F, Chen J-Q, Li Y-F, et al. Mettl3-mediated m6A regulates spermaTogonial differentiation and meiosis initiation. Cell Res (2017) 27(9):1100–14. doi: 10.1038/cr.2017.100

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Furic L, Rong L, Larsson O, Koumakpayi IH, Yoshida K, Brueschke A, et al. eIF4E phosphorylation promotes tumorigenesis and is associated with prostate cancer progression. Proc Natl Acad Sci (2010) 107(32):14134–9. doi: 10.1073/pnas.1005320107

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Yang X, Zhong W, Cao R. Phosphorylation of the mRNA cap-binding protein eIF4E and cancer. Cell Signalling (2020) 73:109689. doi: 10.1016/j.cellsig.2020.109689

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Ueda T, Sasaki M, Elia AJ, Chio II, Hamada K, Fukunaga R, et al. Combined deficiency for MAP kinase-interacting kinase 1 and 2 (Mnk1 and Mnk2) delays tumor development. Proc Natl Acad Sci USA (2010) 107(32):13984–90. doi: 10.1073/pnas.1008136107

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Prabhu SA, Moussa O, Gonçalves C, LaPierre JH, Chou H, Huang F, et al. Inhibition of the MNK1/2-eIF4E axis augments palbociclib-mediated antitumor activity in melanoma and breast cancer. Mol Cancer Ther (2023) 22(2):192–204. doi: 10.1158/1535-7163.MCT-22-0092

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Prabhu SA, Moussa O, Miller WH Jr., Del Rincón SV. The MNK1/2-eIF4E axis as a potential therapeutic target in melanoma. Int J Mol Sci (2020) 21(11):4055. doi: 10.3390/ijms21114055

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Wang J, Da C, Su Y, Song R, Bai Z. MKNK2 enhances chemoresistance of ovarian cancer by suppressing autophagy via miR-125b. Biochem Biophys Res Commun (2021) 556:31–8. doi: 10.1016/j.bbrc.2021.02.084

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Xie J, Shen K, Jones AT, Yang J, Tee AR, Shen MH, et al. Reciprocal signaling between mTORC1 and MNK2 controls cell growth and oncogenesis. Cell Mol Life Sci (2021) 78(1):249–70. doi: 10.1007/s00018-020-03491-1

PubMed Abstract | CrossRef Full Text | Google Scholar