Molecular Characterization and Clinical Relevance of m6A Regulators in Metastatic Prostate Cancer

Qiwei Liu Third A liated Hospital of Sun Yat-Sen University Zhen Li Sun Yat-sen University Cancer Center Lizhao He Third A liated Hospital of Sun Yat-Sen University Ke Li Third A liated Hospital of Sun Yat-Sen University Jialiang Chen Third A liated Hospital of Sun Yat-Sen University Cheng Hu Third A liated Hospital of Sun Yat-Sen University Jun Wang Sun Yat-sen University Cancer Center Fangjian Zhou Sun Yat-sen University Cancer Center Yonghong Li Sun Yat-sen University Cancer Center Hengjun Xiao (  hjxiao555@126.com ) Third A liated Hospital of Sun Yat-Sen University

Accumulated studies have highlighted tight connections between m6A methylation and tumour initiation and progression [6]. In glioblastoma, downregulation of FTO or upregulation of METTL3 were involved in poor prognosis of glioblastoma by promoting the proliferation and self-renewal of glioblastoma stem cells [7]. High expression of METTL3 or METTL4 was also essential for the maintenance and self-renewal of leukaemia stem cells, thus aggravating acute myeloid leukaemia [8]. Upregulation of METTL3 and downregulation of METTL14 can both lead to progression of hepatocellular carcinoma by facilitating cell proliferation and invasion [9,10]. YTHDF2 not only enhances cell proliferation by the AKT/GSK3β/cyclin D1 signalling axis but also inhibits migration and invasion by destabilizing the m6A sites of YAP [11].
However, the mode of action of m6A methylation in prostate cancer remains largely unknown. Herein, we used published sequencing data to investigate the exact role of m6A methylation with respect to prostate cancer.

Material And Methods
Prostate cancer dataset and preprocessing Public gene expression data and relative clinical information were gathered from the TCGA database (https://xenabrowser.net/datapages/). Patients without detailed survival information were removed. In addition, four eligible metastatic prostate cancer cohorts were acquired from https://www.cbioportal.org/, which include mRNA expression data, somatic mutation data and copy number variation (CNV). Clinical annotations were downloaded by the R package cgdsr, and somatic mutation data were collected using the R package TCGAbiolinks [12]. Speci c collected data are shown in Table 1, and more detailed information about the samples is presented in Supplementary Table 1. For data consistency, the original data from https://www.cbioportal.org/ were normalized by the zscore function, and the FPKM data from TCGA were transformed into the zscore normalized dataset. Finally, batch effects were corrected using the R package sva.
Gene set variation analysis (GSVA) and functional annotation To further investigate the biological signi cance of different m6A modi cation patterns, we conducted GSVA enrichment analysis with the "GSVA" R package. GSVA is a nonparametric and unsupervised technique that is commonly used to estimate changes in biological processes and signal pathways in samples [14]. The annotated gene sets of "c2.cp.kegg.v6.2.-symbols" were collected from the MSigDB database (https://www.gsea-msigdb.org/gsea/index.jsp). Adjusted P <0.5 was viewed as statistically signi cant. To carry out functional annotation for m6A-related genes, the clusterPro ler R package was used (FDR cut-off of < 0.05).
Identi cation of differentially expressed genes (DEGs) between distinct m6A phenotypes Referring to distinctly expressed m6A regulators, we classi ed four eligible metastatic prostate cancer cohorts collected from the cBioPortal website into two different m6A modi cation patterns. DEGs between the two distinct modi cation patterns were determined by the R package limma [15]. Genes with p<0.5 and 1.5<fold-change (or fold-change <0.667) were regarded as differentially expressed genes. m6Ascore calculation Redundant genes of DEGs were removed using the random forest approach[16], and the remaining genes were selected for survival analysis (p<0.05). Genes were classi ed into two clusters utilizing the Cox regression model. We then calculated m6Ascore referring to the following GGI method [17]: m6Ascore=scale(∑X-∑Y), where x or y is the gene expression value when the Cox coe cient is positive or negative, respectively. Based on the median value of m6Ascore, samples were divided into m6Ascorehigh and m6Ascore-low. Subsequently, prognostic analysis was performed between the two samples.
Correlation between the m6A gene signature and other related biological processes Mariathasan et al. constructed a series of gene sets involved in speci c biological processes, including immune checkpoints; epithelial mesenchymal transition (EMT) markers such as EMT1 and EMT2; and DNA mismatch repair [18]. We subsequently carried out correction analysis to uncover the relationships between m6Ascore and relative biological pathways.

Copy number variation (CNV) analysis
According to SNP6 CopyNumber segment data, the shared changing areas of copy number among all the samples were detected utilizing the GISTIC method. Relative parameters were set as follows: Q≤0.05, con dence level was 0.95. The above analysis was performed using the corresponding MutSigCV module of GenePattern (https://cloud.genepattern.org/gp/pages/index.jsf, an online analytical tool developed by the Broad Research Institute. Cell culture and cell transfection Human prostate cancer cell lines DU145 and PC3 were obtained from ATCC (USA). Cells were kept in RPMI1640 medium supplemented with 10% FBS at 37°C in a humidi ed incubator with 5% CO 2 .

Western Blot Analysis
Western blotting was conducted as previously described [19]. Brie y, protein concentrations were measured with a BCA Kit. Protein lysates were resolved using SDS-PAGE and transferred onto PVDF membranes (Millipore). The membrane was subsequently incubated overnight (4°C) with the following primary antibodies: anti-METTL14 (Norvus), anti-CSNK1D (Norvus), anti-SLC35E1 (Norvus) and β-actin (Invitrogen). After washing, the membranes were further subjected to the appropriate secondary antibodies (Invitrogen). Blots were visualized by a ChemiDoc XRS system, followed by quanti cation using Image Lab software (Bio-Rad).

Transwell assay
Matrigel was defrosted at 4°C overnight and diluted with serumfree medium (dilution, 1:6). Transwells were inserted in a 24-well culture plate, 40 µl of prediluted Matrigel was inoculated into each Transwell chamber, followed by 2 hours in a 37°C incubator to coagulate. Stably transfected cells were previously seeded in 6-well plates and cultured to 90% con uence. After digestion, a total of 200 µl cell suspension (8×10 4 cells/well) was dispensed to the upper chamber, and 800 µl medium containing 30% FBS was dispensed to the lower chamber. After 24 hours of incubation at 37°C, cells in the upper layer of the Transwell were removed with sterile cotton swabs, followed by PBS washing and xation with methanol for 20 min. Subsequently, cells were further stained with crystal violet dye for 5 min, washed with distilled water, imaged and counted under an inverted microscope.
Wound healing assay Transfected cells were plated into a 6well plate. Before scratching, the culture medium was replaced with serum-free medium containing 1 µg/ml mitomycin C to obtain monolayer cells. Scratches were generated using 200 µl pipette tips, followed by washing three times with PBS. Migrated cells were counted and photographed by a microscope at 0 and 24 hours after scratching.

Cck-8 assay
When the cell con uency reached 70%, drugs were added for 72 hours. DMSO was added to the control groups, and the experimental groups were administered olaparib for 72 hours. Cells were cultured to 90% con uence and then subjected to digestion, centrifugation and resuspension. Cells were further seeded in 96-well plates at a density of 4×103 cells/well. Cell proliferation was detected with a CCK-8 assay following the manufacturer instructions after culture for 24, 48 and 96 hours. The absorbance was measured at 45 nm wavelength.

Statistical Analyses
The bioinformatics differences between the two groups were analysed using the Wilcox test. Referring to the relevance between m6Ascore and patient survival, the cut-off values of different subgroups were identi ed by the survminer R package. Survival curves were generated using Kaplan-Meier analysis, and signi cant differences were determined by log-rank tests. The predictive value of m6Ascore for metastatic samples was evaluated via receiver operating characteristic (ROC) curve analysis, and the area under the curve (AUC) was calculated utilizing the pROC R package. The maftools R package was applied to plot the mutation atlas of patients with high and low m6Ascore. The R package RCircos was used to depict the location of m6A regulators on chromosomes. ns represents p > 0.05, * p ≤0.05, ** p ≤ 0.01, *** p ≤ 0.001, **** p ≤0.0001.
For the experimental data, a two-tailed t test was used with PRISM software. A P value < 0.05 was viewed as statistically signi cant.

Results
The genetic variation of m6A regulators Altogether, 21 m6A regulators (8 writers, 2 erasers and 11 readers) were identi ed. We rst analysed the mRNA expression levels of m6A regulators between metastatic and nonmetastatic samples and found that few genes were differentially expressed, such as FMR1 and FTO (Fig. 1). Subsequently, we summarized the incidence of CNV and somatic mutations of 21 m6A regulators in metastatic, nonmetastatic and NEPC samples. Except for the prevalent missing frequency of CNV in a few regulators, such as FTO, RBM15B and YTHDC2, most regulators experienced ampli cation in copy number (Fig. 1B-E, Supplementary Table 2). Among these samples, mutations of m6A regulators rarely occurred ( Fig. 1F-G). The distribution of m6A regulators on chromosomes is presented in Fig. 1H.
Unsupervised clustering for m6A regulators in metastatic prostate cancer We performed consensus clustering and univariate Cox analysis utilizing m6A gene expression pro les and survival information from the prad_su2c_2019 dataset. The m6A regulation network in Figure 2A (Supplementary Table 3) revealed the interaction and junction of m6A regulators and their impacts on the prognosis of metastatic prostate cancer. We found that not only the same functional categories of m6A regulators but also the distinct functional categories of m6A regulators displayed signi cant correlations in expression.
The above results illustrated that the interactions among distinct functional categories of m6A regulators may exert important roles in various m6A modi cation patterns. We determined the different expression patterns of 21 m6A regulators in four eligible metastatic prostate cancer cohorts downloaded from the cBioPortal website and performed unsupervised clustering analysis using the ConsensusClusterPlus R package, which led to the identi cation of two distinct subclusters (Fig. 2B, Supplementary Table 4). We termed these patterns m6A Clusters A-B.
To investigate biological behaviours among different subgroups, we conducted gene set enrichment analysis (GSEA) (Supplementary Table 5). As shown in Figure 2E, m6A Cluster A was signi cantly enriched in lysine degradation and the mTOR signalling pathway. Cluster B was mainly enriched in arachidonic acid metabolism and steroid hormone biosynthesis (Supplementary Table 6).
Furthermore, we evaluated the expression and mutation atlas of speci c genes between m6A Cluster A and m6A Cluster B (Fig. 2D-E, Supplementary Tables 7-8). Particularly, in the prad_su2c_2019 datasets, the ARV7 score and ARscore between these two clusters showed signi cant differences (Fig. 3, Supplementary Table 9). Further prognosis analysis revealed the distinct difference in the prognosis between these two clusters (Fig. 2C). Subsequently, we performed GSVA based on the gene sets constructed by Mariathasan et al. (Fig. 3C, Supplementary Table 10). The activities of matrix molecules were markedly increased in m6A Cluster B, such as the activation of epithelial mesenchymal transition, transforming growth factor-β and angiogenesis signalling pathways. In addition, the expression levels of m6A regulators in the m6A cluster.A were higher than m6A.cluster.B (Fig. 3D, Supplementary Table 11).

m6A regulators promote PCa cell metastasis and proliferation
To further investigate the function of m6A regulators in the metastasis of PCa, METTL14-overexpressing or METTL14 knockdown PC3 cell lines were established by transfecting a stable overexpressing lentivirus and shRNA, respectively. The e ciency of METTL14 knockdown and overexpression was validated by Western blotting. The results revealed that the protein levels of METTL14 were signi cantly upregulated or downregulated in PC3 cell lines ( Figure 4A). Subsequently, cell migration, invasion, wound healing, and CCK-8 assays were performed to explore the role of METTL14 in PCa cell metastasis and proliferation. Cell migration and invasion assays showed that downregulation of METTL14 decreased the migration and invasion cell numbers, while overexpressing METTL14 reversed the outcomes ( Figure 4B). Wound healing assays revealed that silencing METTL14 reduced, whereas overexpressing METTL14 increased, the wound healing of PC3 cells ( Figure 4C). Moreover, the proliferation of PC3 cells was detected by CCK-8 assay, which elucidated that METTL14 ablation inhibited the proliferation capability, while upregulating METTL14 enhanced the proliferation of PCa cells ( Figure 4D). Additionally, olaparib administration obviously reversed the cell proliferation promotion induced by METTL14 overexpression. Overall, our results implied that METTL14 played an essential role in PCa migration, invasion and proliferation.

Generation of m6A Phenotype Genes and Function
To further investigate each m6A cluster's potential biological behaviours, we determined 2330 metastatic prostate cancer-related differentially expressed genes (DEGs) using the limma package (Supplementary Table 12). The clusterPro lter package was utilized to perform KEGG analysis for DEGs, which indicated enrichment of shearing and RNA transportation (Supplementary Table 13). Then, based on the 2330 m6A phenotype-related DEGs, unsupervised clustering analyses were performed to classify patients, which similarly led to two subtypes termed the m6A gene cluster. A and B (SupplementaryTable14). We observed the poor prognosis in m6AGenecluster.A type tumour and the expression levels of most m6A regulators were higher in the m6A cluster. A than m6AGenecluster.B (Fig. S1C).

Establishment of the Prognostic Model
The above DEGs were de-redundant by the random forest algorithm to select the most category-related genes (Supplementary Table 15_sig.txt). The Cox regression model was used to uncover the relationship between these genes and patient survival. Next, we divided the above genes into two categories based on their coe cient values and scored for all the samples using the m6Ascore formula (Supplementary Table  15_nGenes.txt, Supplementary Table 15_pGenes.txt). Referring to the median m6Ascore, samples were further grouped into two categories: m6Ascore high and m6Ascore low samples (Fig. 5A, Supplementary Tables 16-17). As shown in Figure 5B, the prognosis in the m6Ascore high sample group was poorer than that in the m6Ascore low group. This means that the prognosis of samples could be characterized by the m6Ascore model. Finally, the correlation analysis of m6Ascore and feature genes selected from gene sets constructed by Mariathasan et al. revealed that m6Ascore was signi cantly associated with biological functions such as DNA repair and mismatch repair (Fig. 5C-D, Supplementary Table 18). The Wilcoxon test indicated that there was a notable difference between m6A cluster and m6A gene cluster in m6Ascore (Fig. 5E-F). m6Ascores of samples highly expressed m6A cluster.A or m6A gene cluster.A were markedly higher than those of samples with highly expressed m6Acluster.B or m6Agenecluster.B.
Similarly, we investigated the function of category-related genes in the metastasis of PCa. PC3 cell lines were stably transfected with lentiviruses expressing control shRNA, CSNK1D shRNA and SLC35E1 shRNA. The e ciency of CSNK1D or SLC35E1 knockdown was veri ed by Western blotting ( Figure 6A). Then, cell migration, invasion, wound healing, and CCK-8 assays were performed to explore the role of CSNK1D and SLC35E1 in PCa cell metastasis and proliferation. Cell migration and invasion assays showed that downregulation of CSNK1D and SLC35E1 decreased the migration and invasion cell numbers ( Figure 6B). Wound healing assays also revealed that silencing CSNK1D or SLC35E1 reduced the wound healing of PC3 cells ( Figure 6C). Moreover, proliferation of PC3 cells was detected by CCK8 assay, which showed that CSNK1D or SLC35E1 ablation inhibited the proliferation capability, and administration of olaparib further inhibited the proliferation of PCa cells in the CSNK1D or SLC35E1 ablation groups ( Figure 6D). In summary, our results revealed that CSNK1D and SLC35E1 were of great signi cance in PCa migration, invasion and proliferation.

Molecular Characteristics between high and low m6Ascore
Additional investigations of differences between m6Ascore high and low groups in prad_su2c_2019 datasets revealed that the ARscore in different groups were distinct; in the m6Ascore high group, the ARscore was also high (Fig. 7A). Then, we analysed the difference in somatic mutations between the m6Ascore high and low groups. As depicted in Figure 7B-C, the mutation numbers in the m6Ascore high groups were higher than those in the m6Ascore low group. Similarly, CNV numbers were also higher in the m6Ascore high groups than in the m6Ascore low groups (Supplementary Table 19). In the m6Ascore high groups, 18 copy number ampli cations and 31 copy number deletions were detected, while in the m6Ascore low groups, 16 copy number ampli cations and 30 copy number deletions were detected ( Fig. 7D-E).

Veri cation of m6Ascore
To further validate the prediction performance of our prognostic model, the m6Ascore of TCGA samples was calculated. The threshold of classi cation was determined by the R function surv_cutpoint. Consistently, survival analysis indicated that the prognosis in the m6Ascore-high group was poorer than that in the m6Ascore-low group (Fig. 8A, Supplementary Table 20). Furthermore, the m6Ascore also showed a signi cant difference in parts of the GLEASON_SCORE groups ( Fig. 8B-C).
In particular, we trained our prostate cancer metastasis prediction model in the prad_su2c_2019 and TCGA cohorts, which achieved an ROC AUC of 70% (Fig. 8D).

Discussion
PCa is a major malignancy affecting the male population worldwide, and effective therapeutic options for advanced-stage PCa, especially metastatic PCa, are still scarce [20]. As the most wide-ranging posttranscriptional modi cation, N6-methyladenosine (m6A) is strongly correlated with cancer cell proliferation, progression and metastasis [21]. In PCa, however, relevant studies are lacking, and there are no prediction models based on m6A regulators to evaluate the prognosis of metastatic PCa.
In this study, we found that the mRNA expression of most genes did not exhibit prominent differences between primary and metastatic samples, except for a few genes, such as FMR1 and FTO. We also performed integrative analyses on CNVs and mutation alterations and mRNA expression of m6A regulators in primary, metastatic and NEPC samples. Although few mutations were observed, their biological signi cance has been veri ed to be vital during tumour progression. A mutation in METTL14 could facilitate tumour proliferation via the AKT signalling pathway [22]. There is a paucity of studies targeting mutations in m6A regulators in PCa, but in acute myeloid leukaemia, mutations of m6A regulators were predictive of unfavourable prognosis [23]. CNVs are strongly related to mRNA expression. Speci cally, copy number gain could foster ampli cation of genes, and copy number reduction inhibits the expression of genes. Except for a few regulators, such as YTHDF2, FTO and RBM15B, most of them experienced CNV ampli cation. Ampli cation of FTO was reported to signi cantly improve the prognosis of prostate cancer [24]. However, the same m6A regulator may exert different roles in distinct tumours through diverse mechanisms. Herein, two distinct molecular subgroups of metastatic prostate cancer with obviously distinct characteristics were shown based on 21 m6A regulators related to prognosis. m6Acluster.A was signi cantly enriched in lysine degradation and the mTOR signalling pathway. While the m6Acluster.B was mainly enriched in arachidonic acid metabolism and steroid hormone biosynthesis. We know that activating mTOR signalling can enhance tumour proliferation and progression via distinct mechanisms, including the enhancement of angiogenesis, glyolytic and lipid metabolism, and inhibition of autophagy [25]. Additionally, m6A regulator expression was higher in m6A clusters. A than in m6Acluster.B. To further investigate the relationship between the expression of m6A regulators and PCa prognosis, METTL14-overexpressing or METTL14 knockdown PC3 and DU145 cell lines were constructed. Similar to previous studies, METTL14 ablation inhibited the proliferation and metastasis capability, while upregulating METTL14 enhanced the proliferation and metastasis of PCa cells [26].
Furthermore, in this study, the transcriptomic heterogeneity among distinct subgroups of metastatic prostate cancer was found to be markedly related to shearing and RNA transportation. A total of 2330 DEGs were presented as m6A phenotype-related genes. Similar to the m6A regulator clustering results, two distinct genomic subtypes were identi ed based on m6A phenotype-related genes (2330). Prognosis in m6AGenecluster.A type tumour was dismal, and most m6A regulators were expressed in the m6A cluster. A were higher than in m6AGenecluster.B. Next, we selected the most category-related genes based on the above DEGs and then constructed a prognostic model to provide a reference for treating patients with metastatic prostate cancer. We observed that m6Ascore was signi cantly correlated with biological functions such as DNA repair and mismatch repair. Similarly, m6Ascores of samples upregulated m6Acluster.A or m6Agenecluster.A were distinctively higher than those in samples overexpressing m6Acluster.B or m6Agenecluster.B. This work implied that m6A regulators play an essential role in the prognosis of metastatic PCa, and patients with high m6A scores may be more appropriate for pharmacy therapy targeted for DNA repair, such as poly (ADP-ribose) polymerase (PARP) inhibitors (PARPis).
In the high m6Ascore groups, the ARscore, mutation and CNV numbers, which were unfavourable factors for prognosis, were correspondingly elevated. In this model, CSNK1D is located on chromosome 17. Gene expression and activity changes in CSNK1D have been observed in distinct cancers [27]. In metastatic HCC, the expression level of CSNK1D was higher than that in nonmetastatic HCC [28]. SLC35E1 (solute carrier family 35, member E1) is a nucleotide sugar transporter carrier. It has been reported that for colorectal liver metastases, SLC35E1 could be predictive of the therapeutic effect of 5-uorouracil-based chemotherapy [29]. In our validation experiment, silencing CSNK1D or SLC35E1 reduced the proliferation and metastasis of DU145 and PC3 cells, which showed similar effects to the vehicle groups administered olaparib. Furthermore, KDM1A, the rst identi ed demethylase, also termed LSD1 or KIAA0601, can regulate the initiation of tumours [30]. CCCTC-binding factor (CTCF) is a well-known regulator facilitating chromatin into topologically associated domains by enhancing cohesin-mediated loop formation [31], which is strongly associated with cancer initiation [32]. RBBP4 could promote the malignant progression of colon cancer through the Wnt/β-catenin pathway [33]. CDC23 regulates the tumour cell phenotype and is upregulated in papillary thyroid cancer [34]. Cell division cycle 5-like (CDC5L) protein, a cell phase regulator of the G2/M transition, has been demonstrated to improve bladder cancer cell proliferation, migration and invasion [35]. As an RNA-binding protein, hnRNPA1 can regulate the expression and translation of several mediators involved in tumour initiation and progression [36].

Conclusions
In short, the prognostic model could be applied to guide more effective judgement of prognosis as well as treatments for metastatic prostate cancer in clinical practice. For metastatic PCa patients, a high m6A risk score indicates a dismal prognosis. Since the m6Ascore was signi cantly correlated with biological functions such as DNA repair and mismatch repair, patients with high m6Ascores may be appropriate candidates for pharmacy therapy targeted for DNA repair, such as PARPi. However, there are some pitfalls in this study. Although an independent dataset was used to validate the prognostic model and cell studies were performed to uncover the vital role of m6A-associated genes in metastatic PCa, other animal and clinical studies should be performed. Moreover, the present study is largely a bioinformatic analysis, and potential underlying mechanisms need to be further studied.

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
We all read and approved the nal manuscript.
Availability of data and material Not applicable.

Competing interests
None.  Table   Table 1 is not available with this version.  Unsupervised clustering of m6A regulator genes in metastatic prostate cancer. (A) The interaction among m6A regulator genes. The size of circle indicates the effect of each gene on survival, the larger the size, the greater the effect is; green spots inside the circle indicate risk prognostic factors, black spots inside the circle indicate factors; lines that connect genes exhibit genetic interactions, red and blue represent positive and negative associations, respectively; gene Cluster A,B and C are shown as blue, red and brown, respectively; (B) Consensus clustering m6A regulator genes in metastatic samples; (C) GSVA enrichment analysis. Heatmaps show the activation status of biological pathways, which is displayed with different m6A clusters; red denotes activation, blue denotes inhibition; (D) and (E) show the distribution of the mutation and expression of partial genes in two m6A clusters, respectively.  and low m6Ascore and the overall survival rate; (C) Pearson's correlation analysis highlighting the correlations between m6Ascore and the known genes in prad_su2c_2019 cohorts. Red, blue and X symbols represent positive, negative and nonsigni cant, respectively; the larger the circle, the more signi cant there is; (D) The distribution of enrichment scores of known genes prad_su2c_2019 cohorts between high and low m6Ascore samples; (E) and (F) show the distribution of m6Ascore among m6Aclusters and m6Ageneclusters, respectively.