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.

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.

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.

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).

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).

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).

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.

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 .

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.

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.

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.

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.

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.

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.

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.

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.

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.

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).

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.

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.

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.

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.

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.