Construction of a risk prediction model using m6A RNA methylation regulators in prostate cancer: comprehensive bioinformatic analysis and histological validation

Abstract

Background

Epigenetic reprogramming reportedly has a crucial role in prostate cancer (PCa) progression. RNA modification is a hot topic in epigenetics, and N6-methyladenosine (m6A) accounts for approximately 60% of RNA chemical modifications. The aim of this study was to evaluate the m6A modification patterns in PCa patients and construct a risk prediction model using m6A RNA regulators.

Materials and methods

Analyses were based on the levels of 25 m6A regulators in The Cancer Genome Atlas (TCGA). Differentially expressed gene (DEG) and survival analyses were performed according to TCGA-PRAD clinicopathologic and follow-up information. To detect the influences of m6A regulators and their DEGs, consensus clustering analysis was performed, and tumor mutational burden (TMB) estimation and tumor microenvironment (TME) cell infiltration were assessed. mRNA levels of representative genes were verified using clinical PCa data.

Results

Diverse expression patterns of m6A regulators between tumor and normal (TN) tissues were detected regarding Gleason score (GS), pathological T stage (pT), TP53 mutation, and survival comparisons, with HNRNPA2B1 and IGFBP3 being intersecting genes. HNRNPA2B1 was upregulated in advanced stages (GS > 7, pT3, HR > 1, and TP53 mutation), as verified using clinical PCa tissue. Three distinct m6A modification patterns were identified through consensus clustering analysis, but no significant difference was found among these groups in recurrence-free survival (RFS) analysis. Six DEGs of m6A clusters (m6Aclusters) were screened through univariate Cox regression analysis. MMAB and PAIAP2 were intersecting genes for the five clinical factors. MMAB, which was upregulated in PCa compared with TN, was verified using clinical PCa samples. Three distinct subgroups were established according to the 6 DEGs. Cluster A involved the most advanced stages and had the poorest RFS. The m6A score (m6Ascore) was calculated based on the 6 genes, and the low m6Ascore group showed poor RFS with a negative association with infiltration for 16 of 23 immune-related cells.

Conclusion

We screened DEGs of m6Aclusters and identified 6 genes (BAIAP2, TEX264, MMAB, JAGN1, TIMM8AP1, and IMP3), with which we constructed a highly predictive model with prognostic value by dividing TCGA-PRAD into three distinct subgroups and performing m6Ascore analysis. This study helps to elucidate the integral effects of m6A modification patterns on PCa progression.

Introduction

Prostate cancer (PCa) is a leading malignant tumor among men [1]. PCa has primarily been treated with surgical prostatectomy or androgen deprivation therapy (ADT). However, it can become castration-resistant PCa (CRPC), and biochemical recurrence or metastasis may occur during traditional therapy, which is the main cause of cancer-specific death. Therefore, elucidating the molecular mechanisms related to PCa progression is crucial in the discovery of diagnostic biomarkers and therapeutic targets.

Epigenetic reprogramming is reported to serve a crucial role in the progression of PCa [2]. Recently, RNA modification has been regarded as a hot topic in epigenetic research, and nearly 172 different RNA modifications are present in MODOMICS [3]. Among them, N6-methyladenosine (m6A) is widespread throughout the transcriptome; indeed, m6A comprises approximately 60% of RNA chemical modifications and is present on 0.1% to 0.4% of total adenosine residues, including > 300 noncoding RNAs and 7600 mRNAs, in eukaryotes [4,5,6]. RNA m6A methylation regulates mRNA alternative splicing, stability, and intracellular localization, constituting the major posttranscriptional modification [7]. The formation of m6A is regulated by three categories of proteins: readers (which recognize m6A-modified sites, such as YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, HNRNPC, FMR1, LRPPRC, HNRNPA2B1, IGFBP1, IGFBP2, IGFBP3, RBMX, and ELAVL1), writers (methyltransferases, such as METTL3, METTL14, METTL16, WTAP, VIRMA (also called KIAA1429), ZC3H13, RBM15, RBM15B, and CBLL1), and erasers (demethylases, such as FTO and ALKBH5).

m6A regulatory genes are reported to participate in various carcinogenic and tumor progression processes [8, 9]. METTL3 is reported to advance PCa progression and is associated with poor prognosis by stabilizing the mRNAs of MYC, LEF1, and integrin β1 (ITGB1) by m6A methylation [10,11,12]. However, one study suggested that low expression of METTL3 is associated with resistance to therapy with androgen receptor antagonists via upregulation of NR5A2/LRH-1 [13]. This finding indicates the controversy and opposing functions of METTL3 in PCa. YTHDF2-induced AKT phosphorylation and MDB3B m6A modification may also promote PCa proliferation, migration, and invasion [14, 15]. FTO, an m6A demethylase, inhibits the invasion and migration of PCa cells by regulating total m6A levels [16]. Nevertheless, there are insufficient data on m6A regulators in PCa, and the role of m6A regulators remains controversial; in general, comprehensive transcriptome and genomic analysis is needed. This study fully analyzed m6A-related genes in PCa progression and prognosis.

Materials and methods

Data acquisition

Transcriptome profiling and single nucleotide variation data for prostate adenocarcinoma in The Cancer Genome Atlas (TCGA) were downloaded from the GDC Data Portal (https://portal.gdc.cancer.gov/). Copy number and clinical phenotype data were downloaded from the University of California Santa Cruz Xena (https://xena.ucsc.edu/). Gene expression matrices were extracted and obtained through Practical Extraction and Report Language (Perl) (version 5.34.0) and R software (4.0.3) (R Development Core Team, Vienna, Austria). The R package “RCircos” was used to generate Circos plots.

Differentially expressed gene (DEG) analysis

mRNA levels were analyzed with TPM (transcripts per kilobase of exome per million mapped reads) data, which were transformed from the HTSeq-FPKM transcriptome profiling data of TCGA-PRAD. The R packages “limma” and “ggpurb” were then used to identify DEGs between the normal and tumor groups and for further statistical analysis. The Wilcoxon or Kruskal−Wallis test was performed to determine DEG levels, and P < 0.05 was identified as statistically significant. The R packages “clusterProfiler” and “enrichplot” were used to analyze Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment.

Survival and correlation analyses

Survival data for “biochemical_recurrence”, “days_to_first_biochemical_recurrence”, and “days_to_last_follow_up.diagnoses” were obtained from clinical phenotype data. When “days_to_first_biochemical_recurrence” was indicated, we regarded these patients as having “recurrence status”, and the time notated was used as the “recurrence follow-up time”. Other cases were regarded as having “no recurrence”, and “days_to_last_follow_up.diagnoses” was used as “recurrence_follow_up_time”. The R packages “survival” and “survminer” were used for survival analysis. Survival curves were evaluated through Kaplan–Meier and log-rank tests. Correlation analysis was performed through Pearson or Spearman correlation analysis, and a prognostic network map was drawn using the R packages “igraph”, “psych”, “reshape2”, and “RColorBrewer”. To build a PCa prognostic model using DEGs of m6A clusters, univariate Cox regression analysis was conducted, and P < 0.05 was used for later gene consensus clustering analysis.

Consensus clustering analysis and principal component analysis (PCA)

To determine whether m6A regulators are related to PCa prognosis, the cohort from TCGA was allocated into different groups based on the consensus level of m6A regulators. The process was performed using the R package “ConsensusClusterPlus” and resulted in cluster consensus and item-consensus results. The graphical output consisted of heatmaps, consensus cumulative distribution function (CDF) plots and delta area plots. The cluster number was determined through a high consistency of clusters, a low coefficient of variation, and no significant increase in the CDF curve. The chi-square test or Fisher’s exact test was used to analyze clinicopathological characteristics and clustering. The heatmap was generated through the R package “pheatmap”. Recurrence-free survival (RFS) was detected among groups using Kaplan–Meier and log-rank tests. PCA was performed to judge the fitness of the classification with the prcomp function of R software.

Gene set variation analysis (GSVA)

GSVA is used for assessing KEGG gene set enrichment, which is a nonparametric and unsupervised method [17]. GSVA could comprehensively score the gene sets of interest and translate them into the pathway level. In this study, to evaluate the potential pathway changes of different clusters, we downloaded “KEGG gene sets as Gene Symbols” from the GSEA website (http://www.gsea-msigdb.org/gsea/) and used the GSVA algorithm to comprehensively score each gene set.

Tumor mutational burden (TMB) estimation

TMB was calculated by the total number of mutated/total covered bases [18]. PCa was classified into two groups based on the optimum threshold segmentation of TMB status population. We analyzed the relationship between the m6A score (shown as m6Ascore below) and the TMB and then performed survival analysis comparing prognosis between the high TMB and low TMB groups.

Tumor microenvironment (TME) cell infiltration

TME infiltration levels were calculated through single-sample gene set enrichment analysis (ssGSEA) and quantified using enrichment scores [17]. The gene set of each TME infiltrating immune cell type was obtained as previously reported [19]. Correlation analysis between m6Ascore and immune-associated genes was performed to illustrate the relationship.

Clinical PCa samples

Forty-five pathologically diagnosed PCa patients (15 with Gleason score (GS) < 7, 15 with GS = 7, and 15 with GS > 7) were recruited from Beijing Tongren Hospital and Beijing Chaoyang Hospital in accordance with the Ethics Committee of Beijing Tongren Hospital and Beijing Chaoyang Hospital, affiliated with Capital Medical University. All patients underwent prostatectomy between 2016 and 2021; PCa and adjacent normal tissues were removed and stored in liquid nitrogen. The clinicopathological characteristics of PCa patients are shown in Table 1.

Table 1 Clinicopathological characteristics of PCa patients

Quantitative real-time PCR (qPCR) analysis

qPCR analysis was performed as previously described by our group [20, 21]. Total RNA was isolated using TRIzol™ reagent (Invitrogen), and complementary DNA (cDNA) was synthesized through One-Step gDNA Removal and cDNA Synthesis SuperMix (TransGen Biotech, Beijing, China). qPCR was performed using Top Green qPCR SuperMix (TransGen Biotech) on an SDS 7500 FAST Real-Time PCR system (Applied Biosystems, Foster City, CA, USA). GAPDH or 18S ribosomal RNA was used as an endogenous reference gene. The relevant primer sequences are shown in Additional file 1: Table S1.

Immunophenoscore (IPS) analysis

IPS analysis was performed as previously reported [22]. IPS determines immunogenicity by referring to effector cells, immunosuppressive cells, MHC molecules and immunomodulators. IPS (ranging from 0 to 10) was calculated according to the gene expression of the representative cell types without bias using machine learning methods. The IPS results of TCGA-PRAD patients were downloaded from The Cancer Immunome Atlas (TCIA) (https://tcia.at/home).

Statistical analysis

Statistical analyses were conducted using R software (4.0.3), SPSS software version 23 (IBM, Armonk, New York, USA), and GraphPad Prism 7.0 (GraphPad Software, La Jolla, CA, USA). Incomplete data were excluded. Chi-square or Fisher’s exact tests were used for categorical variables, and Wilcoxon or Kruskal–Wallis tests were used for continuous data. Correlations between two levels were assessed through Pearson or Spearman correlation analysis. Survival analysis was applied through the log-rank test of Kaplan–Meier survival analysis and hazard ratio (HR) with the 95% confidence interval (CI) of univariate Cox proportional hazard models. P < 0.05 was indicated as significant.

Results

Characteristics of m6A regulators in PCa

In the TCGA-PRAD dataset, we analyzed copy number variation (CNV) analysis alterations, DEGs, and the mutation frequency of m6A regulators of PCa through comparison with normal samples. For CNV events, approximately 76% (19/25) of m6A regulators lost DNA copy number, with ZC3H14 having the highest degree of copy number loss (28.49%) (Fig. 1A). Six m6A regulators gained copy number, among which VIRMA had the highest such percentage (3.19%) (Fig. 1A). The m6A regulator CNV alterations and locations on chromosomes are shown in the lower panel of Fig. 1A. Of the 499 tumor and 52 normal samples in the dataset from TCGA, DEGs of m6A regulators were statistically estimated using TPM data. The levels of METTL3, HNRNPA2B1, RBM15B, and IGFBP2 were higher, while those of ZC3H13, FTO, and IGFBP3 were lower in PCa tissues than in normal tissues (P < 0.001) (Fig. 1B and Table 2). In mutation frequency analyses, 16 m6A-related genes were mutated among 19 of 484 (3.93%) samples; mutations in VIRMA (KIAA1429), YTHDC2, RBM15B, YTHDF2, and IGFBP1 were detected in one sample (Fig. 1C). ZC3H14 exhibited the highest mutation rate (4/484), and all 25 m6A regulators exhibited low mutation rates (< 1%) in TCGA samples. This suggests the relatively conserved and stable expression of m6A regulators during PCa progression.

Fig. 1

Characteristics of m6A regulators in prostate cancer. A (top) The CNV variation frequency of m6A regulators in TCGA-PRAD. Blue dot: deletion frequency; red dot: amplification frequency. (bottom) The location of CNV alterations of m6A regulators on 23 chromosomes. Blue dot: copy number loss, red dot: copy number gain. B Box plot of the expression of the 25 m6A regulators in PCa and normal tissues. Median values ± interquartile ranges are shown in the graph. ns P > 0.05; * P < 0.05; ** P < 0.01; ***P < 0.001. C The mutation frequency of 25 m6A regulators in 484 TCGA-PRAD patients is shown on the waterfall plot. Columns represent individual patients, the upper bar plot shows the TMB, the right bar plot shows the proportion of each variant type, and the stacked bar plot below shows the transformed fraction of each patient

Table 2 The expression levels of 25 m6A regulators in TCGA-PRAD and normal tissues

Expression of m6A regulators in PCa prognosis and different clinicopathological characteristics

The expression of m6A regulators in different GSs and pathological T (pT) stages was estimated using TCGA-PARD TPM data. As PCa with GS ≥ 7 is associated with worse prognosis [23, 24], the PCa samples were divided into three groups: GS < 7, GS = 7, and GS > 7. Compared with the GS < 7 group, IGFBP3, HNRNPA2B, RBMX, RBM15B, YTHDF1, HNRNPC, VIRMA, and FMR1 were significantly highly expressed in the GS > 7 group (P < 0.001) (Fig. 2A and Table 3). In addition, T2 (188), T3 (295), and T4 (10) pT stages were compared. Among the 25 m6A regulators, IGFBP3, HNRNPA2B1, VIRMA, and RBMX were more highly expressed in the T3 stage than in the T2 stage (P < 0.001) (Fig. 2B and Table 4). Based on Kaplan–Meier curves in survival analysis, high expression of HNRNPA2B1, IGFBP1, and ELAVL1 was associated with poor RFS (P < 0.001) (Fig. 2C). Moreover, TP53 mutation in PCa was associated with shorter radiographic progression-free survival (rPFS) and time to CRPC [25]. Compared with the TP53 wild-type group of PCa samples, the levels of VIRMA and IGFBP3 were higher and the level of IGFBP2 was lower in the TP53 mutation group (P < 0.001) (Fig. 2D and Table 5). Risk factors and relationships of m6A regulators are summarized in the prognosis network shown in Fig. 2E.

Fig. 2

Expression of m6A regulators in relation to different PCa clinicopathological characteristics and prognoses. A–D Distribution of 3–4 representative m6A regulators (3–4 with the lowest P values) in TCGA-PRAD data stratified by GS (A), pT (B), RFS (C), and TP53 mutation (D). The box plots show the median ± interquartile range of values, and P values are presented above each pair of comparisons. E Prognosis network of interactions between m6A regulators in PCa. The P values of each regulator with respect to the prognosis are shown as circles of different sizes. Purple in the right hemisphere: risk factors for RFS; green in the right hemisphere: favorable factors for RFS. The erasers, readers, and writers of the m6A regulator are shown as blue, orange, and red colors, respectively on the left. Positive or negative correlations of m6A regulators are linked with lines of different colors (pink: positive, blue: negative), and the correlation strength between them is shown as different line thicknesses. F Intersecting DEGs of 25 m6A regulators in GS (GS > 7 vs. GS < 7), pT (T3 vs. T2), RFS (P value of univariate Cox regression analysis < 0.05), TN (tumor vs. normal PCa tissues), and TP53 (mutation vs. wild type) are shown as Venn diagrams. G PCa tissues and adjacent normal tissues were divided into four groups according to GS score (15 of normal tissue, 15 of GS < 7, 15 of GS = 7, and 15 of GS > 7). mRNA levels of intersecting genes were compared in different groups of PCa and adjacent normal tissues (PCa vs. adjacent normal; adjacent normal vs. GS < 7 vs. GS = 7 vs. GS > 7). Relative HNRNPA2B1 and IGFBP3 mRNA levels were assessed through qPCR analysis and normalized to adjacent normal tissues or those of endogenous reference genes. Means ± SEMs are shown in the graphs. ns P > 0.05 and * P < 0.05

Table 3 Expression levels of m6A regulators in TCGA-PRAD data stratified by GS

Table 4 Expression levels of m6A regulators in TCGA-PRAD data stratified by pT

Table 5 Expression levels of m6A regulators in TCGA-PRAD data stratified by TP53 mutation

We analyzed DEGs of m6A regulators in GS (GS > 7 versus (vs.) GS < 7), pT (T3 vs. T2), RFS (P value of univariate Cox regression analysis < 0.05), TN (tumor vs. normal PCa tissues), and TP53 (mutation vs. wild type) comparisons. Intersecting genes were identified via a Venn diagram, and HNRNPA2B1 and IGFBP3 were differentially expressed in all comparisons (Fig. 2F). The level of HNRNPA2B1 was higher in the PCa group than in the normal group and high in association with advanced-stage parameters: namely, GS > 7, pT3, HR > 1, and TP53 mutation. However, the level of IGFBP3 was lower in PCa tissues but higher in the presence of the above advanced-stage indicators. This suggests that different molecular mechanisms and expression patterns may exist in PCa and throughout progression. We then assessed the expression of HNRNPA2B1 and IGFBP3 in 45 PCa tissues (including 15 with GS < 7, 15 with GS = 7, and 15 with GS > 7) and 15 adjacent normal prostate tissues. Levels of HNRNPA2B1 and IGFBP3 were higher in 66.7% (10/15, P < 0.05) and 40.0% (6/15, P > 0.05) of PCa tissues than in adjacent normal tissues, respectively, and the GS > 7 group exhibited elevated expression of HNRNPA2B1 and IGFBP3 compared with the GS < 7 group (both P < 0.05) (Fig. 2G).

Consensus clustering analysis according to m6A regulators

To further explore the phenotypes of m6A regulators in different prognoses and clinicopathological characteristics of the TCGA-PRAD cohort, we performed consensus clustering analysis to identify subgroups of 495 PARD cases according to 25 m6A regulators (m6Acluster). Relatively high consistency, a low coefficient of variation, and an appreciable increase in the area under the CDF curve were determined as the criteria of cluster number. With regard to the relative change in the area under the CDF curves for the cluster number, k = 3 was determined to be the best category number of clusters (Fig. 3A). Subclasses were evaluated via PCA, and all three clusters were significantly distinguished, especially in clusters A and C (Fig. 3B), indicating correct prediction of the m6Acluster (these clusters are shown as “clusters 1–3 or A-C”). Next, we compared different clinicopathological characteristics of PRAD among the clusters, and the results showed relatively higher pathological N (pN) stage, GS 8–10, and biochemical recurrence rates in cluster A (Fig. 3C). The different KEGG pathways between clusters A and C were analyzed through GSVA (Fig. 3D). Cancer- and apoptosis-related pathways were significantly dysregulated between clusters A and C. We also analyzed changes in immune cell infiltration in different m6Aclusters and found 12 of 23 subpopulations of immune cells to be differentially expressed in the three clusters (Fig. 3E and Table 6). Kaplan–Meier survival curves of RFS based on the m6Aclusters demonstrated no significant differences among the groups (P = 0.475) (Fig. 3F).

Fig. 3

Consensus clustering analysis based on m6A regulators. A (left) Relative change in area under the CDF curve from k = 2 to 9. (right) Color-coded heatmap of the consensus matrix for k = 3. Color gradients represent values from 0–1 (white: 0, dark blue: 1). B PCA of the transcriptome profiles of three m6Acluster patterns, showing a marked difference in the transcriptome between different m6Aclusters. Different dots of red, orange, and blue in the scatter diagram represent m6Aclusters A to C. C The heatmap shows unsupervised clustering of 25 m6A regulators in TCGA-PRAD. PSA grade, pT, pN, GS, and biochemical recurrence were used for patient annotation. Red represents high expression, and blue represents low expression. D Heatmap of KEGG enrichment analysis showing the activation states of biological pathways in distinct m6Aclusters. Red represents activation, and blue represents inhibition. E Abundance of each infiltrating immune cell among the three m6Aclusters is shown in the boxplot. The median ± interquartile range of values is shown in the graph. ns P > 0.05; * P < 0.05; ***P < 0.001. F Survival analysis for RFS among three m6Aclusters based on 495 TCGA-PRAD patients. Kaplan–Meier curves and log-rank P values are shown in the graph, and the numbers at risk are shown at the bottom

Table 6 Abundance of each infiltrating immune cell among the three m6Aclusters

Consensus clustering analysis based on DEGs of m6Aclusters

DEGs (P value of Bayes test < 0.001) were analyzed in every pairwise comparison according to the three clusters based on m6A regulators, and the 74 intersecting genes were identified via a Venn diagram (Fig. 4A). We then analyzed the biological functions of these 74 genes, which were categorized into GO terms of biological process (BP), cell component (CC), and molecular function (MF). Under the stringent threshold of P-adjust < 0.05, only 1 specific CC (proton-transporting V-type ATPase, V0 domain) was enriched (Additional file 2: Fig. S1A-B). Additionally, KEGG signaling pathway analysis of 74 genes indicated significant enrichment in the oxidative phosphorylation KEGG pathway (Additional file 2: Fig. S1C-D).

Fig. 4

Consensus clustering analysis according to DEGs of m6Aclusters. A DEGs were analyzed in every pairwise comparison of m6Aclusters and intersecting genes were identified via a Venn diagram. B Univariate Cox regression analysis was performed to evaluate RFS in relation to 74 intersecting genes. The statistical significance (P < 0.05) of 6 genes (BAIAP2, TEX264, MMAB, JAGN1, TIMM8AP1, and IMP3), P value, and HR value (with 95% CI) are shown in the plot. C Consensus clustering analysis according to DEGs of m6Aclusters was performed. (left) Relative change in area under the CDF curve from k = 2 to 9. (right) Color-coded heatmap of the consensus matrix for k = 3. Color gradients represent values from 0–1 (white: 0, dark blue: 1). D The heatmap shows unsupervised clustering of the 6 genes in TCGA-PRAD. PSA grade, pT, pN, GS, and biochemical recurrence were used for patient annotation. E Survival analysis for RFS among three geneClusters based on 495 TCGA-PRAD patients. Kaplan–Meier curves and log-rank P values are shown in the graph, and the numbers of patients at risk are shown at the bottom. F Expression of m6A regulators among three geneClusters is shown in the boxplot. The median ± interquartile range of values is shown in the graph. ns P > 0.05; * P < 0.05; ** P < 0.01; ***P < 0.001. G As described in Fig. 2F, intersecting DEGs of GS, pT, RFS, TN, and TP53 are shown in a Venn diagram. H As described in Fig. 2G, PCa tissues and adjacent normal tissues were divided into four groups, and mRNA levels of intersecting gene (MMAB) in different groups were compared. * P < 0.05

Univariate Cox regression analysis was also performed on the 74 intersecting genes, and 6 (BAIAP2, TEX264, MMAB, JAGN1, TIMM8AP1, and IMP3) related to PCa recurrence were selected (P < 0.05) (Fig. 4B). The characteristics of the 6 or 5 DEGs (TIMM8AP1 was excluded because of being a processed pseudogene and a lack of copy number data) regarding CNV alterations, DEGs, and mutation frequency were analyzed in PCa compared with normal samples (Additional file 3: Fig. S2A-D and Additional file 1: Table S2). We found that the levels of 4 (BAIAP2, IMP3, JAGN1, MMAB) of the 6 genes were significantly upregulated in PCa tissues. Interestingly, JAGN1 exhibited the highest degree of copy number loss (5.18%), but its expression was higher in PCa tissues than in normal prostate tissues (P < 0.01). When assessing the expression of the 6 genes with respect to the clinicopathological characteristics GS, pT, and TP53, we found that most were downregulated in advanced stages (GS > 7, pT3, and TP53 mutation). Survival analysis also indicated that low expression of these 6 genes was significantly associated with poor RFS (Additional file 4: Fig. S3A-D and Additional file 1: Table S3–S5). The risk factors and relationships of the 6 genes are summarized in the prognosis network illustrated in Additional file 4: Fig. S3E.

Consensus clustering analysis was performed for these 6 genes (geneCluster), and k = 3 was determined to result in the best classification with respect to the delta area results (cluster 1–3 or A-C) (Fig. 4C). PCA verified the significance of the three subgroups (Additional file 5: Fig. S4A).

A heatmap of clinicopathological features in the three geneClusters indicated that pT stage, N stage, GS, and biological recurrence were significantly higher in geneCluster A than in the other clusters (Fig. 4D). Kaplan–Meier survival analysis revealed significant differences between the three m6Aclusters, among which geneCluster A had the poorest RFS (P = 0.012) (Fig. 4E). Next, we analyzed the expression of m6A regulators in different geneClusters and found that 21 of 25 m6A regulators were differentially expressed among them (Fig. 4F and Table 7). Moreover, 16 of 23 subpopulations of immune cells were differentially expressed among the three clusters (Additional file 5: Fig. S4B and Additional file 1: Table S6). In the GSVA of different KEGG pathways, a small cell lung cancer-related pathway was dysregulated between geneClusters A and C (Additional file 5: Fig. S4C).

Table 7 Expression of m6A regulators among three geneClusters

We also analyzed DEGs of the 6 genes regarding GS, pT, RFS, TN, and TP53, as illustrated in Fig. 2G. The intersecting genes MMAB and BAIAP2 were identified via a Venn diagram (Fig. 4G), and both were highly expressed in PCa tissues compared with normal tissue but were low with respect to the four parameters, indicating advanced stage. This suggests opposite expression patterns of MMAB and PAIAP2 in PCa occurrence and progression. The mRNA expression level of MMAB was analyzed through qPCR using 45 PCa tissues and 15 adjacent normal prostate tissues, as mentioned above. It was upregulated in 60.0% (9/15, P < 0.05) of PCa tissues compared with adjacent normal tissues, with the GS > 7 group exhibiting reduced expression compared with the GS < 7 group (P < 0.05) (Fig. 4H). The differential expression of MMAB in PCa occurrence and progression was histologically verified.

Low m6Ascore based on the 6 genes was associated with poor prognosis in RFS

The clustering analysis above was based on the population of the patient and could not be used to quantitatively assess the patterns of m6A regulators. Considering the complexity of m6A regulators and individual heterogeneity, the m6Ascore was calculated through PCA of the 6 gene levels in TCGA-PARD, and statistical analysis was performed to examine the relationship between m6Ascore and geneCluster/m6Acluster. The results revealed a significant difference in both clusters: geneCluster A had the lowest m6Ascore, and geneCluster C had the highest (Fig. 5A). PRAD patients were divided into two groups based on the optimum threshold segmentation of m6Ascore in RFS analysis, and according to Kaplan–Meier survival analysis, poor RFS was associated with the low m6Ascore group (Fig. 5B). The alluvial diagram shows the attribute changes in m6Acluster, geneCluster, m6Ascore, and recurrence status (fustat) (Fig. 5C). All patients in the low m6Ascore group were classified into geneCluster A, which was relevant to the worse RFS outcome (recurrence rate: 27.6%). The correlation of m6Ascore and immune-associated genes was evaluated, resulting in 16 of 23 immune-related cell infiltrates being found to be negatively related to m6Ascore (Fig. 5D). This suggests the negative regulation of the immune reaction during m6Ascore elevation. A previous study showed that low immune infiltration was potentially related to prolonged survival [26]. Consistent with our results, the high m6Ascore group had a tendency toward favorable prognosis in terms of RFS (Fig. 5B).

Fig. 5

The low m6Ascore group based on the 6 genes was associated with poor prognosis in RFS. A The m6Ascore was calculated through PCA of the 6 gene levels in TCGA-PARD, and m6Ascores among geneClusters (left) and m6Aclusters (right) were compared. The median ± interquartile range of values and P values are shown in the box plots. B PRAD patients were divided into low and high m6Ascore groups based on the optimum threshold segmentation of m6Ascores in relation to RFS. The prognosis of the two groups was evaluated through survival analysis of Kaplan–Meier curve and log-rank P value. C Alluvial diagram showing changes in m6Aclusters, geneClusters, m6Ascore, and recurrence status (fustat of 0: no biochemical recurrence, 1: biochemical recurrence). D Correlations between m6Ascore and immune cell infiltration in TCGA-PRAD. Red: positive correlation, blue: negative correlation

Characteristics of m6Ascore status for TCGA-PRAD tumor mutation and subtypes

A study suggested that the TMB score was positively associated with tumor-related oncogenic mutations, which led to tumor progression but was easier to detect by immune cells, making them more likely to be sensitive to immunotherapy [27, 28]. As shown in Fig. 5B, TCGA-RRAD patients were divided into two groups (58 in the low m6Ascore group and 437 in the high m6Ascore group), and the association of m6Ascore and TMB was assessed. The results showed no significant difference in TMB scores in the m6Ascore groups (P = 0.13) (Fig. 6A, left). Spearman correlation analysis also showed no significant correlation between the expression levels of m6Ascore and TMB (P > 0.05) (Fig. 6A, right). This is discordant with other study in gastric cancer, which showed that m6Ascore and TMB exhibited a significant negative correlation [19]. The PRAD patients were then divided into two groups based on the optimum threshold segmentation of TMB value in RFS analysis (242 in the low TMB group and 231 in the high TMB group). In Kaplan–Meier analysis, a tendency toward poor RFS was found in the high TMB group compared with the low TMB group, but the difference was not significant (P = 0.051) (Fig. 6B, left). This suggests that the higher malignancy of PCa patients potentially correlated with the high TMB score group. When analyzing m6A and TMB together, we found a great RFS advantage for the combination of low TMB with high m6Ascore (Fig. 6B, right). Consistent with our previous result, the low TMB and high m6Ascore groups had better prognoses in terms of RFS (Fig. 5B and Fig. 6B, left).

Fig. 6

Characteristics of m6Ascore status in TCGA-PRAD tumor mutation and subtypes. A (left) TMB values in low and high m6Ascore groups. The P value is shown in box plot. (right) Correlation analysis of m6Ascore and TMB value in PCa was performed through Spearman correlation analysis. Different dots of red, orange, and blue in the scatter diagram represent geneClusters A to C, respectively. B (left) PRAD patients were divided into low and high TMB groups based on the optimum threshold segmentation of TMB values in relation to RFS. Survival analysis of the two groups is shown in the graph. (right) PRAD patients were divided into four groups based on TMB and m6Ascore status (high (H)-TMB + H-m6Ascore; H-TMB + low (L)-m6Ascore; L-TMB + H-m6Ascore; and L-TMB + L-m6Ascore), and survival analysis was performed. C The waterfall plot of tumor somatic mutations in PRAD with low m6Ascore (left) and high m6Ascore (right). D (left) The proportion of biochemical recurrence status in the m6Ascore groups (fustat of 0: no recurrence, 1: recurrence). (right) m6Ascore with respect to different biochemical recurrence statuses. The P value is shown in box plot. E Survival analysis of different m6Ascore groups among TRAD patients with GS 8–10 (left) and pT3 (right). F Box plot of PD-L1 expression in the low and high m6Ascore groups. The P value is shown in box plot. (G) Violin plot of IPS scores in ips_ctla4_neg_pd1_neg (CTLA4 negative response and PD1 negative response) in respective m6Ascore groups. The P value is shown in violin plot

The distribution differences of somatic mutations were analyzed between m6Ascore groups, and the results indicated a more extensive mutation frequency of TP53 (18% vs. 8%) but a reduced mutation frequency of SPOP (2% vs. 12%) in the low m6Ascore group (Fig. 6C).

For the relationship between m6Ascore and fustat, the low m6Ascore group experienced a higher rate of recurrence (Fig. 6D left), which was consistent with the RFS analysis in Fig. 5B. Biochemical recurrence was also associated with lower m6Ascore (Fig. 6D, right). We then analyzed prognosis in the m6Ascore groups with regard to different clinicopathological features through RFS analysis (Fig. 6E and Additional file 6: Fig. S5A-B). The low m6Ascore group was associated with poor RFS in GS 8–10 and pT3 stages (P < 0.05) (Fig. 6E).

The expression of PD-L1 was detected to elucidate a potential response to immunotherapy, and the low m6Ascore group showed relatively high levels of expression (P = 0.023) (Fig. 6F). Finally, we further predicted the value of the risk score for immune checkpoint blockade (ICB). The IPS results of TCGA-PRAD were downloaded from the TCIA website. PD1 and CTLA4 were enrolled for IPS analysis and further classified into four parts: ips_ctla4_neg_pd1_neg (CTLA4 negative response and PD1 negative response), ips_ctla4_neg_pd1_pos (CTLA4 negative response and PD1 positive response), ips_ctla4_pos_pd1_neg, and ips_ctla4_pos_pd1_pos (Fig. 6G, Additional file 6: Fig. S5C). The average IPS in the comparison of the low and high m6Ascore groups showed no significant difference in the four parts of the negative or positive response of PD1 and CTLA4 (Fig. 6G, Additional file 6: Fig. S5C). These results indicate that m6Ascore may lack efficacy in the risk score model for predicting the response to PD1 and CTLA4 therapy.

Discussion

In this study, we comprehensively evaluated the m6A modification patterns in TCGA-PRAD and found 6 DEGs (BAIAP2, TEX264, MMAB, JAGN1, TIMM8AP1, and IMP3) based on m6Aclusters. A model of optimal clusters was established, and m6Ascore analysis was performed based on the above 6 genes.

In general, elucidating the molecular mechanism of PCa progression remains crucial. It has been reported that epigenetic reprogramming is key for PCa progression [29, 30] and that inhibition of the epigenetic regulator EZH2 might effectively overcome ADT resistance [2]. Our previous study found that epigenetically activating AR/NDRG1 signaling through histone methylation and DNA methylation significantly suppresses CRPC progression [20, 21].

Nevertheless, the effect of mRNA posttranscriptional modification in PCa is still unclear. m6A RNA methylation may be among the most extensive RNA modifications [31,32,33]. Yabing Chen et al. highlighted that total RNA m6A modification levels are significantly increased in PCa tissues due to upregulation of METTL3 [34]. In addition, knockdown of METTL3 significantly reduces PCa cell migration and invasion. YTHDF2, a reader of m6A modification, synergistically induces PCa progression with METTL3 by regulating AKT phosphorylation [15]. However, low levels of METTL3 were also found to lead to advanced metastatic PCa that is resistant to androgen receptor antagonists [35]. These controversial results reveal the complexity of m6A regulators and individual heterogeneity of m6A regulators in PCa, and various m6A regulators may form a complex network structure and interact with each other to affect PCa progression.

In this study, we examined three representative groups of m6A regulators: 14 “readers”, 9 “writers”, and 2 “erasers”. Bioinformatic analysis was then performed to comprehensively investigate the association between m6A modification patterns and PCa clinicopathological characteristics and prognosis.

Compared with normal prostate tissues in TCGA database, METTL3 was significantly highly expressed in PCa, and negligible CNV and mutation rate were found in this analysis. ZC3H13 exhibited the highest degree of copy number loss (28.49%) and highest mutation rate (4/484) among all 25 m6A regulators. The mRNA of ZC3H13 was also expressed at low levels in the TN comparison (P < 0.001). Therefore, a potential mechanism including ZC3H13 mutation may exist in PCa progression. Insufficient study has been conducted to elucidate this point, which warrants further verification.

In survival and Cox regression analyses, we evaluated RFS instead of overall survival (OS) because only 10 of 493 TCGA-PARD patients died, making it difficult to obtain statistically significant results by assessing OS. In our RFS analysis, “biochemical_recurrence” and “days_to_first_biochemical_recurrence” data did not coincide, and we ultimately used the “days_to_first_biochemical_recurrence” data to confirm “recurrence status” and “recurrence follow-up time”.

The landscape of m6A variation in PCa was recently reported [26, 36,37,38,39,40], although no scholars have analyzed sufficient m6A regulators and verified them using clinical PCa tissues. To overcome these shortcomings, we explored all 25 widely-acknowledged m6A regulators and assessed relevant associations of various clinicopathological characteristics, such as TN, GS, pT, TP53 mutation, and survival analysis of RFS using TCGA-PRAD data. Additionally, we analyzed DEGs of m6A regulators in relation to the five factors, and the intersecting genes (HNRNPA2B1 and IGFBP3) were verified with clinical PCa tissues. The level of HNRNPA2B1 was higher in PCa tissues than in normal prostate tissues and was highly correlated with the other four advanced factors. However, the level of IGFBP3 was lower in PCa tissues but higher in the presence of the above advanced-stage indicators. This suggests that independent molecular mechanisms and expression patterns exist in PCa occurrence and progression. Consistent with the bioinformatic analysis, the mRNA levels of HNRNPA2B1 and IGFBP3 were elevated in PCa tissue of the high GS group (GS > 7) compared with the low GS group (GS < 7). In the TN comparison, the expression of HNRNPA2B1 was higher in PCa, but no significant difference in the level of IGFBP3 expression was found.

It has been reported that upregulation of HNRNPA2B1 by PCAT6 promotes PCa progression and neuroendocrine differentiation [41]. HNRNPA2B1 may also be an independent prognostic factor and contribute to cancer progression [42]. Coexpression network analysis using clinical data from the GSE70768 dataset as well as quantitative proteomic mass spectrometry profiling and gene enrichment analysis using LNCaP and PC3 cell lines suggest that HNRNPA2B1 is associated with PCa progression and prognosis [43, 44]. The level of IGFBP3 in PCa is controversial, with one meta-analysis indicating that the CC genotype of the IGFBP3–202A/C polymorphism is associated with an increased risk of PCa [45,46,47]. Overall, tissue verification and previous studies support a certain level of accuracy of our bioinformatic analysis using TCGA-PRAD data.

To analyze the overall effects of m6A regulators in PCa, we performed consensus clustering analysis and divided the TCGA-PRAD cohort into three subgroups based on 25 m6A regulators. As no significant difference in the three clusters for RFS was observed, this cluster analysis with m6A regulators may not be suitable for inclusion in a prognostic risk prediction model.

As epigenetic regulators, posttranscriptional modification of target genes is the primary function of m6A regulators. To analyze the downstream genes of m6A regulators, we focused on the DEGs of the m6Acluster and screened six risk genes (BAIAP2, TEX264, MMAB, JAGN1, TIMM8AP1, and IMP3) through univariate Cox regression analysis. Only IMP3 has previously been reported in PCa studies, and it appears to be increased in PCa tissues and associated with higher GS, PCa metastasis, and PCa-specific survival [48,49,50]. Regardless, according to our bioinformatic analysis of TCGA-PRAD, the expression of IMP3 was decreased in three groups related to advanced PCa (HR (RFS) < 1, TP53 mutation < wild type, pT3 < pT2, and tumor > normal). This finding indicates discordance between TCGA-PRAD data and previous results, which warrants further confirmation.

Consensus clustering analysis was then performed based on the 6 genes to construct a risk prediction model, and TCGA-PARD patients were divided into three subgroups. In these clusters, the clinicopathological parameters of PSA grade 3 (grade 1: > 0 and < 1, grade 2: 1–10, grade 3: > 10), pT 4, pN 1, and GS > 7 were gathered in cluster A, which was related to significantly poor RFS, suggesting an accurate risk prediction model based on 6 genes. We then analyzed the DEGs of the 6 genes in correlation with the above five factors, and MMAB and BAIAP2 were identified as intersecting genes. Levels of both were higher in PCa tissues than in adjacent normal prostate tissues but lower in tissues with advanced disease factors. We next verified MMAB based on different GSs of PCa tissues and adjacent normal tissues and found that its level of mRNA was elevated in PCa tissues but decreased in the GS > 7 group, consistent with the bioinformatic analysis.

To quantitatively illustrate the m6A signature, we calculated m6Ascore according to the expression of the 6 genes in TCGA-PARD: m6Ascore correlated positively with geneCluster, and lower m6Ascore was associated with poor prognosis in RFS. Furthermore, m6Ascore was negatively correlated with 16 of 23 immune-associated cells. Thus, a potential mechanism by which the m6A signature protects against PCa progression is by negatively regulating immune cell infiltration. A high TMB score corresponds with tumor-related mutations and sensitivity to immunotherapy [27, 28], and Yue Zhao et al. reported that m6A modification in PCa may contribute to immunotherapy strategies [26]. In our study, the lower m6Ascore group tended to have a higher proportion of TP53 mutations (18%), and the high m6Ascore group tended to have a higher proportion of SPOP mutations (12%). However, no significant correlation with TMB was found in the two m6Ascore groups in our study, and there were also no differences in PD-L1 expression or IPS score with response prediction of negative or positive PD1 and CTLA4 therapy. In summary, the results of this study do not clearly support an influence of m6A modification patterns on PCa immune cell infiltration and immunotherapy.

Conclusions

In this study, we comprehensively analyzed m6A regulators in TCGA-PRAD through bioinformatic analysis and tissue validation. First, we screened DEGs of m6Aclusters and found 6 genes (BAIAP2, TEX264, MMAB, JAGN1, TIMM8AP1, and IMP3), through which we divided PCa patients into three subgroups and calculated m6Ascore to construct a risk model with high predictive value for recurrence. This study may contribute to determination of the effects of m6A signaling on the progression and prognosis of PCa.