Identification of Molecular Markers Associated with Prostate Cancer Subtypes: An Integrative Bioinformatics Approach

1. Introduction

In a recent demographic study, it has been observed that the increasing population ageing, due to demographic and social transitions, contributes to a rapid increase of new cancer cases (24 million in 2035), among which prostate cancer (PCa) still remains the leading cause of oncologic death in men worldwide [1]. As this pathology affects predominantly older men, cancer management can be complicated by comorbidity and age-related variations with significant social and economic implications. Considering the high impact that the burden of PCa has on families and health services, there is an urgent need for improving cancer surveillance and treatment, thus ensuring adequate disease management.

In the early stages of PCa development, tumour cell growth is dependent on circulating testosterone, providing the rationale for using androgen deprivation therapy (ADT) for localised disease. Under genetic and epigenetic alterations, PCa clones with a marked malignant phenotype evolve into the metastatic state (mPCa), which precedes the insurgence of castration resistance. The PCa clonal heterogeneity and gene instability, amplified by antineoplastic therapy effects, can be associated with the ineffectiveness of conventional ADT [2]. The early stages of "castration-resistant" prostate cancer (CRPC) are characterised by heterogeneous PCa cell populations generated under anti-androgens pressure, where the stress conditions allow adaptive cellular reprogramming to low circulating and tissue levels of testosterone. This condition leads to tumour progression and proliferation of malignant clones through either androgen receptor (AR) pathway reactivation or AR-independent pathways activation [3,4]. As a consequence, in a single patient, it is possible to highlight extensive phenotypic heterogeneity characterised by the co-presence of cell populations with different evolutionary histories and drug susceptibility. Interestingly, it has been proposed that this diversity may originate from non-mutational mechanisms which result in an expansion of isogenic populations differing by their sensitivity to ADT [5]. Based on recent observations, it is hypothesised that ADT resistance results from the interplay between redundant genetic and epigenetic mechanisms engaged in complex crosstalk with cellular plasticity that facilitates adaptation to prolonged drug exposure [6,7,8].

Although multimodal approaches and new therapeutic resources improved metastatic CRPC (mCRPC) patients' care, reliable prognostic and predictive criteria for selecting adequate therapy are missing, especially for high-risk patients. The lack of biomarkers able to guide clinicians to the most appropriate therapeutic choice led to weak prognosis improvement after treatment of metastatic PCa (mPCa) and mCRPC [9]. The precision medicine inadequacy shown by several neoplasm treatments would depend on the absence of guidelines on target molecule selection based on clinical evidence, as defined by the ESMO Scale for Clinical Action of Molecular Targets-ESCAT [10]. In addition, although many studies have highlighted genetic alterations associated with PCa evolution, a significant role of epigenetic regulation has emerged in controlling the cancer cell plasticity involved in androgen-resistance acquisition [7,11,12]. The improvement of clinical practices, then, requires further studies to identify candidate target key genes and pathways associated with different PCa and CRPC cell populations.

The possibility of performing genome-wide gene deletion experiments, such as CRISPR-Cas9 and RNAi, allowed investigation of the gene essentiality in hundreds of cancer cell lines [13]. These studies highlighted the conditional nature of gene essentiality, a dynamic property that can change at the modification of genetic and/or environmental context. In particular, diverse sets of essential genes have been observed in different cancer tissues and even in cell lines deriving from the same tissue. These findings have great relevance in the view of precision cancer therapy. For instance, cancer-specific essential genes represent ideal candidates to target cancer cells sparing healthy ones, in which those same genes are dispensable.

The increasing availability of public datasets and independent data derived from multiple studies allows the validation and generalisation of methods and findings, increasing their reliability. Large-scale data integration from multiple experiments, although permitting to increase the robustness of statistical tests by increasing the sample size, still presents challenging tasks related to the exact genetic and phenotypic correspondence of samples. Batch effect correction methods, which are often applied, risk overcorrecting the data removing true but unknown biological differences. To overcome these issues, an alternative approach can be the integration of the results obtained from the analysis of different datasets, searching for a consensus among the independent studies.

Furthermore, while generally, the investigation of the essentiality starts from the definition of the context-specific essentialome to identify candidate biomarkers, here we present an inverse approach. We first identified gene sets relevant to PCa subtypes from a functional and clinical point of view, and then, we characterised the selected genes in terms of drug targeting and essentiality.

The aim of this study was to identify and investigate key genes involved in PCa progression and CRPC evolution toward AR-dependent (AR+) or AR-indifferent (AR-) subtypes in order to derive candidate markers with functional driver roles useful in PCa patients' clinical management. To this extent, we adopted an integrative bioinformatics approach for the analysis of multi-sources data, identifying two gene sets associated with differential expression patterns between primary and metastatic PCa (PCa-gene set) or between CRPC AR+ and AR- (CRPC-gene set). In addition to chemotherapeutic agents already in clinical use, computational drug prediction of both gene sets identified novel compounds for the treatment of advanced variants of PCa and CRPC.

2. Materials and Methods

2.1. Datasets

The present work was based on the integration of several published datasets in order to get a consensus list and a source of validation of genes with a crucial role in PCa progression and CRPC pathogenesis.

In particular, the comparison of metastatic versus primary PCa and mCRPC versus primary PCa have been performed by using the following 8 datasets from Gene Expression Omnibus (GEO) portal: GSE3325 [14], GSE3933 [15], GSE68882 [16], GSE32269 [17], GSE6811 [18], GSE70770 [19], GSE6752 [20], GSE35988 [21].

The datasets published by studies in which gene lists involved in CRPC have been provided and here investigated and integrated with our findings are summarised in Table S1 and are: GSE101607 [22], GSE104786 [23] and GSE77930 [24] from GEO portal and SU2C/PCF (Dream Team, Cell 2015) study [25] from cBioPortal (https://www.cbioportal.org/).

Tissue-specific expression was independently validated in the following studies from cBioPortal: i) PCa-gene set: Prostate Adenocarcinoma (PRAD-TCGA PanCancer Atlas) and Prostate Adenocarcinoma (MSK, Cancer Cell 2010) [26]; ii) CRPC-gene set: SU2C/PCF (Dream Team, PNAS 2019) [27] and Neuroendocrine-PC (Multi-Institute, Nat Med 2016) [28].

2.2. Differential expression analysis

Differential expression analysis of the GEO datasets has been performed exploiting the GEO2R application of the GEO portal for the following comparisons:

• metastatic versus primary PCa groups for GSE3325, GSE3933, GSE68882 datasets;
• mCRPC versus PCa groups for GSE32269, GSE6811, GSE70770, GSE6752, GSE35988 datasets;
• AR-driven versus non AR-driven groups for GSE101607;
• CHGA negative versus CHGA positive/SYP positive/SR negative groups were compared for GSE77930.

Genes with ǀlog2fold-changeǀ ≥1 and adjusted P value<0.05 have been considered significant.

A separate mention is needed for the SU2C/PCF (Dream Team, Cell 2015) dataset from cBioPortal. We selected the samples belonging to NE positive and negative groups based on the RPKM values of 3 well-known markers: CHGA, SYP and AR. We considered the mean value as the threshold and the values above and below it as positive and negative, respectively. This way, we obtained 30 CHGA-positive/SYP-positive/AR-negative and 11 CHGA-negative/SYP-negative/AR-positive samples. We downloaded the ​​RNA-Seq expression median normalised values and performed a t-test to evaluate the differences between the two groups.

2.3. Over-representation analysis

Gene Set Enrichment Analysis (GSEA) of the 32 genes (Table S2) identified in our previous papers [6,29] was performed with the hallmark (H), chemical and genetic perturbations (CGP) and canonical pathways (CP) gene set collections from MSigDB (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp), considering statistically significant the overlapping results showing FDR q-value <0.05 (Table S3). Exploiting the GEPIA database (http://gepia.cancer-pku.cn), we performed ANOVA test to evaluate differences between tumour and paired normal human prostate samples (ǀlog2fold-changeǀ ≥1, adjusted P value < 0.05) of the TCGA-Prostate Cancer dataset (PRAD) and to extract the genes associated with significant disease-free survival (log-rank < 0.05) (Table S4).

Enrichment analysis of the gene-sets was performed by using the enrichR R package (v.3.2). Only terms having an adjusted p-value < 0.05 were considered significant.

The correlation between drug sensitivity and the gene expression was performed through the Gene Set Cancer Analysis platform (GSCA) (http://bioinfo.life.hust.edu.cn/GSCA) by using data from the Genomics of Drug Sensitivity in Cancer (GDSC) database and Cancer Therapeutics Response Portal (CTRP).

2.4. Statistical analysis

Principal component analysis was performed by PCAtools R package v.2.6.0 (https://github.com/kevinblighe/PCAtools) [30] using the following arguments: center = TRUE, scale = TRUE, removeVar = 0.1.

Pearson correlation of samples was obtained by using corrr R package v.0.4.3 (https://CRAN.R-project.org/package=corrr) [31].

Receiver Operating Characteristics (ROC) curves were constructed using OriginPro 2021b software for statistical evaluation of the PCa and CRPC gene sets as potential molecular classifiers. The area under the ROC curve (ROC-AUC) was evaluated, for each gene, according to its deregulated expression (Table S7 and S8) using test directions Positive versus High (up-regulated) and Positive versus Low (down-regulated).

Kaplan-Meier analysis and log rank test were performed to estimate disease free survival and overall survival of PCa patients with high or low expression level of ADAMTS1, SPON2 and EDN3. The cutoff value for each gene was defined according to the Youden index [32], which defines the optimal value as the one maximising the difference between true positive rate and false positive rate over all possible cut-point values.

3. Results

3.1. Identification and characterization of clinically relevant gene sets in PCa progression to metastatic and castration-resistant phenotypes.

In our recent papers, by using in vitro PCa models, we have identified 32 genes (Table S2) associated with resistant phenotype acquisition after androgen deprivation treatments [6,29]. Here, additional computational analyses were conducted to assess the functional and clinical relevance of the 32-gene set and investigate the mechanisms involved in the PCa progression and acquisition of castration resistance (Figure 1).

To this extent, we performed a GSEA of the 32 genes overlapping with the hallmark (H), chemical and genetic perturbations (CGP) and canonical pathways (CP) gene set collections from MSigDB. We found 43 genes significantly (FDR q-value <0.05) associated with processes or pathways enriched by the 32 genes (Table S3). Performing an investigation on TCGA-Prostate Cancer dataset (PRAD) using GEPIA database, we selected, among the above 43 genes, a subset of 18 genes (Table S4) associated with significant disease-free survival and/or with differential expression levels in tumour vs normal human prostate samples, thus extending our set from 32 to 50 genes (Table S5).

Interestingly, by analysing data collected from the literature [29,33,34,35,36,37,38,39,40], among the 50 genes, 17 (34%) are associated with AR as co-regulators of its transcriptional activity and/or defined as AR-controlled genes, and 24 (48%) are associated with prostatic neoplasms, castration resistance, invasiveness, metastasis or carcinogenesis according to Gene Set to Diseases web platform (http://cbdm-01.zdv.uni-mainz.de/~jfontain/cms/?page_id=592) (Table S5).

3.2. Evaluation of 50-gene set alterations in PCa tissue samples at different evolutive stages

To identify a consensus among different experiments and to generalise our findings as much as possible, we performed differential expression analyses of the 50-gene set by using mRNA abundances of patients showing the three different PCa stages from eight public datasets (Table S6). Comparing mPCA to primary PCa, we found 27 (58%) differentially expressed genes (DEGs), including 15 up- and 12 down-regulated genes (molecular profile1, MP1). Two genes (KLK3, QKI) showed a discordant expression pattern and consequently were discarded. Similarly, 25 DEGs (50%) were obtained comparing mCRPC versus PCa, with 12 up- and 13 down-regulated genes (molecular profile 2, MP2). One gene (IL4) was excluded due to a discordant expression profile.

Among the 34 genes included in MP1 and/or MP2 groups (hereinafter called PCa-gene set, Table S7), 15 genes showed the same expression pattern (Figure 2A), and 16 were differentially expressed only in one of the two groups. Of note, an opposite expression profile for ADAMTS1, ETV1 and SMAD2 genes was observed between MP1 and MP2 (Figure 2A).

To investigate the potential application in clinical settings of the PCa-gene set, we evaluated their performance as molecular classifiers in PCa patients by performing a ROC curve analysis. To this extent, we exploited two multidimensional datasets containing both transcriptomics and clinical data: the PRAD-TCGA dataset with 493 samples of primary PCa and the Prostate Adenocarcinoma dataset (MSK, Cancer Cell 2010) including 181 primary PCa and 37 metastatic PCa samples. Significant AUC values (P< 0.05) associated with disease-specific survival, progression, recurrence, metastasis presence and/or high recurrence risk group 1 [41] were shown for 20 genes (Figure 2B and C, Table S8). Kaplan-Meier analysis and log-rank test confirmed that low expression of ADAMTS1, SPON2 and EDN3 was associated with high-risk PCa-specific mortality (Figure 1D).

3.3. Identification of genes associated with different CRPC phenotypes

To investigate the genes with key roles in CRPC, we exploited the findings of two other studies [22,25] that have previously identified genes associated with different CRPC phenotypes (Table S1). The union of their genes with our PCa-gene set gave rise to a new list of 88 genes (Figure 3C).

The association between this new list and AR+ or AR- phenotype was investigated on data provided by Ylitalo et al. (GSE10167) [22], in which AR+ and AR- samples have been identified and labelled as AR-driven and non-AR-driven, respectively, with a strong unbalancing in favour of AR+ (32 AR+; 8 AR-).

From PCA representation and sample correlation analysis, it is evident that the mRNA abundance values of these 88 genes separated the two groups (Figure 3D-F) in a manner comparable to all the genes of the dataset (Figure 3A-B).

The genes having higher loading scores in PC1, and thus giving a higher contribution to the variance, were STEAP2 and CD40 (Figure 3E).

We further included two additional studies [23,24] in our analysis and searched for a consensus representative of CRPC AR+ or AR- phenotypes, where the latter is often associated with neuroendocrine (NE) features and identified as NE+, giving rise to three different profiles: adeno-CRPC (AR+/NE-), NE-CRPC with neuroendocrine features (AR-/NE+) or double-negative CRPC (DNPC) (AR-/NE-).

To this extent, we extracted lists of differentially abundant genes from the four datasets listed in Table S1 as described in the Material and Methods section.

We first extended the 88 gene-set by adding the marker genes identified by [24] (referred to as Nelson list), obtaining a final list of 104 genes. These genes have been searched in the four differentially abundant gene lists to find a consensus. We considered concordant genes having the same expression pattern (up or down) and significantly different abundances between AR+ and AR- (P-value < 0.05) in at least two datasets. This way, we obtained a consensus list of 29 genes: ALDH1A3, AR, CDKN2A, CHGA, CHGB, EHF, ENO2, EZH2, FKBP5, FOXA1, HES6, HNRNPK, HOXA13, HOXB13, KLK2, KLK3, NKX3-1, PCSK1, PIK3CB, PLPP1, PMEPA1, RB1, SCG3, SCN3A, SORD, SPON2, STEAP2, STEAP4, TMPRSS2. In the following, this list will be referred to as CRPC-gene set (Table S9).

The CRPC-gene set showed its discriminative power on the AR+ and NE+ samples of the GSE77930 dataset (Figure 4), as well as on the AR+ and AR- samples of the GSE101607 dataset (Figure 5 C-D), determining a clearer separation than using the Nelson list (Figure 5A-B).

3.4. Validation of CRPC-gene set on independent datasets

The CRPC-gene set was validated by verifying the association of different CRPC phenotypes with transcriptomic alterations using two independent datasets from cBioPortal: the SU2C/PCF (208 mCRPC samples/201 patients) and the Neuroendocrine-PC (49 samples/35 patients). For this purpose, mRNA expression of the CRPC-gene set was evaluated using ROC analysis.

In the SU2C/PCF study, tissue samples can be stratified into AR+ (AR-score ≥ 0.25), AR- (AR-score <0.25), NE+ (NE-score ≥ 0.4) and NE- (NE-score < 0.4) [28]. The significant AUC values (Asymptotic Probability <0.05) obtained for the differential expression level of the 29 CRPC genes (Table S9, “Deregulated expression” column), indicated a clear discrimination of all four groups (Figure S1). Interestingly, using the NE-score cutoff of 0.4, the CRPC-genes down-regulated in AR- were able to classify with greater accuracy the NE+ samples .

Considering the pathological classification available in the SU2C/PCF and the Neuroendocrine-PC studies, significant AUC values were obtained for all CRPC-genes. More specifically, in both examined studies, the genes with positive LogFC (Table S8) efficiently identified NE-CRPC samples, while the negative ones were associated with Adeno-CRPC samples (Table S10; Figure 6A and B).

3.5. Functional enrichment of gene-sets

The functional role of the genes composing the three gene sets was assessed by performing the enrichment analysis of Gene Ontology Biological Processes (GO-BP) and KEGG pathways terms and selecting the top 50 in terms of overlapping genes.

Most of the terms were shared between MP1 and MP2 genes, while different top-enriched terms were obtained by using the CRPC set, confirming their important distinction (Figure 7). All sets enriched “Regulation of apoptotic process” from GO-BP and “Prostate cancer”, “Pathways in cancer”, and “Pancreatic cancer” from KEGG.

MP1 and MP2 genes were particularly involved in the regulation of transcription and processes related to tumoral transformation, such as angiogenesis, proliferation and differentiation. MP1 genes were more involved in biosynthetic processes, regulation of cell adhesion and response to stimuli, MP2 genes instead in signalling pathways. CRPC genes enriched hormone-related pathways, inflammation response and membrane traffic, terms characterising the advanced status of the tumoral processes.

The KEGG pathways confirmed the shared and distinct roles highlighted from GO-BP terms, as MP1 and MP2 appeared to be involved in regulation and interaction processes, MP1 alone in adhesion molecules, cell cycle and proteoglycans, while MP2 in the acquisition of stem cells properties as alteration of NOTCH signalling. CRPC genes enriched many cancer pathways.

3.6. Investigating the context-specific essentiality

The genes we identified as key players in the different phenotypes under study have been investigated in terms of context-specific essentiality.

To this extent, we exploited DepMap, the most curated portal containing data from experiments for the identification of essential genes in cancer.

We downloaded the Gene Effect scores (GEs) from CRISPR and RNAi gene deletion experiments on cell lines of Prostate Carcinoma [42]. The more the score is negative, the more the gene is essential for the survival and growth of the cell. Considering the above results, we focused on FOXA1, HOXB13, ETV1, ADMATS1, EZH2 and AR (Figure S2). All the genes showed context-specific essentiality as having a different behaviour in the different cell lines, remarking their potential as markers for targeted therapies.

According to the CRISPR scores (shown in the left panel of the figure), all the genes showed essentiality in VCaP and/or LNCaP cell lines, particularly marked values were observed in the case of FOXA1 and HOXB13. FOXA1 showed highly negative values also in the case of 22Rv1 cells according to CRISPR and NCIH660 and VCaP according to RNAi (middle panel). HOXB13 resulted strongly essential also in PC3 (RNAi). ADAMTS1 showed slight essentiality for VCaP, PC3 and 22Rv1 cells in CRISPR experiments and in DU145 and VCaP in the case of RNAi. EZH2 and ETV1 also showed no strong essentiality, but the most negative values were obtained in LNCaP and MDAPCA2B cells for both genes and SHMAC5 only for EZH2. AR showed essentiality, particularly in LNCaP (RNAi) and 22Rv1 (CRISPR) cells.

It is worth noticing that an inverse proportionality can be observed between GEs and mRNA expression (right panel), as a confirmation that essential genes are strongly overexpressed in cancer.

3.7. Computational drug identification based on the PCa- and CRPC-gene sets

In order to predict candidate compounds to treat aggressive subtypes of PCa, we explored the drug sensitivity of the PCa- and CRPC-gene sets using GSCA.

Correlating mRNA expression of the 20 genes up-regulated in MP1 and/or in MP2 and drug sensitivity data from both GDSC and CTR resources, we highlighted 60 compounds with potential inhibitory activity on PCa progression (Figure S3).

Significant positive correlations were observed for PPARG, ATP1B1, ASPH, ADAMTS1, VEGFA and NOTCH3 genes, involved in different stages of PCa progression (Figure 8A). Interestingly, 17-AAG, a drug candidate in human clinical trials (phase 3), potentially exerted an inhibitory activity on EZH2, PAX5 and KHDRBS1 genes associated with the risk of developing metastatic phenotypes (Figure 8A).

Similarly, performing the analysis for the CRPC-gene set (Figure S4), we identified 60 candidate chemicals to be used for Adeno-CRPC in combination with or alternatively to anti-androgen therapy (Figure S4A) and 25 for NE-CRPC treatment (Figure S4B). The overexpressed genes STEAP2, PMEPA1, PLPP1, ALDH1A3, FOXA1 and EHF likely are drug targets in the case of Adeno-CRPC treatment, while ENO2, HES6 and EZH2 for treating the NE-CRPC (Figure 8B).

It is worth noticing a significative relationship between the expression of genes found to be essential in the previous analysis (FOXA1, HOXB13, ETV1, ADMATS1, EZH2 and AR) and 60 chemicals, among which several were already identified in the above analyses (Figure S5).

In Table S11, we reported 20 FDA-approved drugs for the treatment of different cancer histotypes and potentially to be evaluated also for the advanced forms of PCa. Of note, 14 compounds have already been tested in preclinical and clinical studies on their potential efficacy for PCa and CRPC treatment, and 11 were potentially active on essential genes. Drugs that exerted an inhibitory effect on epigenetic regulation of transcription seem to have increased activity in androgen-dependent forms while those which interfere with MAPK-regulated signal transduction pathways could be more effective on the NE-CRPC.

4. Discussion

The significant efforts made to study the molecular mechanisms involved in PCa progression, identified molecular alterations potentially relevant for improving the clinical management of the patients. Over the past decade, the use of new drugs inhibiting the AR axis improved the treatment of mPCa and some forms of CPRC, but these chemicals have led to increased CRPC with neuroendocrine features that are still incurable [3,46]. Hence the hypothesis that PCa progression may depend not only on genetic alterations but also on epi- and/or non-genetic factors caused by drug-stress pressure. For this reason, researchers are currently focusing on identifying the key mechanisms and the master genes driving PCa cells' fate toward metastatic and resistant phenotypes to develop novel and more effective therapeutic approaches [7,8,9,11,66].

Heterogeneity is an intrinsic characteristic of PCa that accentuates during neoplastic evolution, particularly in the advanced phases of disease and in the acquisition of resistance to anti-androgen treatments. Over the years, numerous in vitro and ex-vivo models have been developed [42] in an attempt to mimic human PCa evolution and study the biological mechanisms and the key proteins involved in cell growth and proliferation leading to neoplastic progression. Studies performed using preclinical models have identified molecular alterations as potential markers for a specific pathological stage and/or as therapeutic targets [10]. Nonetheless, due to the intrinsically static and homogeneity of the in vitro models, frequently, the results obtained were not suitable to be translated into clinical practice.

Starting from 32 proteins previously identified in in vitro models mimicking the early stages of androgen-resistance [6,29] and using an integrative bioinformatics approach (schematically illustrated in Figure 1), we derived two gene panels that dynamically track the evolution of phenotypic changes : the PCa-gene set and the CRPC-gene set (Table S7 and S9). The 34 genes included in the PCa-gene set, according to their expression in PCa tissues, were divided into two groups: i) MP1, containing genes altered during primary PCa progression toward metastatic phenotype and ii) MP2, involving genes modified by anti-androgen resistance acquisition (Figure 2A). The CRPC-gene set, instead, was based on the expression profile of 29 genes discriminating between the AR-driven (Adeno-CRPC) and the AR-indifferent (NE-CRPC) resistant phenotypes.

Nine genes (AR, EZH2, FOXA1, HOXB13, HOXA13, KLK3, EHF, SORD and SPON2) were shared by both PCa and CRPC gene sets suggesting their potential role in the multiple evolutive phases of PCa progression (Figure 8). The percentage of genes participating directly or indirectly in the AR activity increased from 47% for those shared between MP1 and MP2 (7 over 15 genes) to 78% (7 over 9 genes) for genes common to PCa- and CRPC-gene lists. These results confirmed the essential role of AR in regulating PCa cells growth and proliferation during the early stages of the oncogenesis, as well as in the resistance acquisition [3]. It has been reported that the anti-androgen therapy pressure, exerted on the intrinsic cellular heterogeneity of PCa tissues, led to the selection of tumour cells showing AR gene alterations or epigenetic perturbation of its regulated pathways [8,11]. Consequently, during the tumour progression, two different phenotypes can be developed: androgen-responsive, where AR always acts as a driver gene, or androgen-insensitive, independent of AR activity [66]. In agreement with this statement, we observed that AR was highly expressed in the stages leading to mPCa and Adeno-CRPC, while was poorly expressed in NE-CRPC (Figure 8).

The ROC curve analysis of multidimensional datasets including transcriptomic and clinical data of PCa tissues from the Prostate Adenocarcinoma study (MSK, Cancer Cell 2010), revealed that the transition to the metastatic phenotype was significantly associated with over-expression of AR, EZH2, NOTCH3, ZIC2 and SOX4 genes, and down-regulation of ALDH3A2, EHF, EPHA3, KLK3, SORD, XBP1, ADAMTS1 and PRDM5 genes (Figure 2B). The correlation between NOTCH3-MMP3 axis activation and bone metastases induction [68], the pro-metastatic activity of EHF knockdown in PCa cells [69] and the involvement of SOX4 and EZH2 in enhancing the PCa cells invasiveness toward the activation of AKT and β-catenin pathways [11] have been experimentally demonstrated. According to the above, high levels of SOX4 and EZH2 were observed in the patients included in the high recurrence group 1 of the TGCA-PRAD study (Figure 2C).

The PCa progression toward the androgen-refractory state was exclusively correlated, in agreement with previous studies [70,71,72], with low expression levels of ALDH3A2, increased expression of NOTCH1, NOTCH4 and QKI (Figure 8), that significantly identified patients with poor prognosis (Figure 2B and C).

Our investigation interestingly highlighted that ADAMTS1 may act, in the context of PCa tissue, as tumour suppressor or as a pro-tumorigenic factor in agreement with previous experimental evidence [73]. Furthermore, the down-regulation of ADAMTS1 was associated with the enhancement of tumour toward metastasis, while high levels were observed in CRPC. Recently, it has been reported that, depending on tumoral contexts, ADAMTS1 can induce proteolytic extracellular matrix modification activating cell plasticity, a biological process involved in the acquisition of therapeutic resistance [74]. It is not surprising that the low level of ADAMTS1 in PCa tissues had a significant prognostic value as a disease-specific marker (Figure 2D).

Using two datasets reporting tissue-specific transcriptomic data and histo-pathologic information of CRPC patients, we verified the high discrimination power of the CRPC gene-set in identifying the expression patterns associated with Adeno- or NE-phenotype (Figure 6). Six genes regulating AR activity (AR, TMPRSS2, HOXB13, NKX3-1, FKBP5, ALDH1A3, PMEPA1) and PLPP1, a gene involved in sphingolipid metabolism and linked to androgen signalling [75], resulted over-expressed in Adeno-CRPC tissues of both datasets (Figure 6), in agreement with the described role of AR in these patients [3]. In NE-CRPC samples, the six molecular markers of NE phenotype CHGA, CHGB, ENO2, PCSK1, SCG3 and SCN3A [67,76] were up-regulated, as well as HES6 and EZH2 genes, involved in the regulation of cell fate decision. Ramos-Montoya et al. [77] reported that HES6 had a driving role in androgen independence acquisition by activating AR-independent pathways for sustaining the survival of PCa cells treated with anti-androgen therapies. Interestingly, EZH2 and CDKN2A were included in the HES6-associated gene signature strongly connected with unfavourable outcomes of PCa patients. In our study, with respect to genes involved in AR signalling, an inverse behaviour was observed for the expression of EZH2. Indeed, it gradually increased as the tumour progressed toward metastatic and castration-resistant forms (Figure 8), in agreement with experimental evidence previously reported [11,78]. It has been suggested that EZH2, acting as an epigenetic regulator, controls cell cycle and cell stemness, and promotes epithelial-mesenchymal transition (EMT), metastatic progression in PCa and neuroendocrine trans-differentiation in CRPC [12,79,80].

Our work, in agreement with recent experimental evidence, highlighted that, in addition to the above-mentioned genes, FOXA1, HOXB13 and ETV1 have a pivotal role in supporting cell growth and proliferation during the stages of PCa evolution (Figure 8). FOXA1 and HOXB13 encode for a chromatin remodeler that acts as a pioneer factor (PF) by reprogramming AR transcriptional activity [81]. In normal prostate cells, these PFs regulate AR accessibility to activate the transcriptional processes while, in tumoral cells, their altered availability triggers aberrant transcription programs leading to the different oncogenic phenotypes of the evolutive forms of PCa [11,80,81]. Consequently, in accordance with the AR behaviour, the expression of FOXA1 and HOXB13 genes was inversely correlated with the aggressiveness of PCa forms such as mPCa, mCRPC and NE-CRPC (Figure 8).

The ETS transcription factor ETV1 cooperates with AR in regulating the transcription of androgen-driven genes involved in cell growth and proliferation. Beana et al. [82] have found that ETV1 overexpression activates, in concert with AR, an oncogenic program in PCa cells leading to metastatic phenotype and patient's worse outcome. Although in our analysis, the ETV1 expression was up-regulated in localised vs metastatic PCa, lower expression level of this gene was detected in mCRPC with respect to PCa (Figure 8). The role of ETV1 in CRPC subtypes is still unclear, although it has been observed that the interplay between TGF-β/SMAD signalling and the ETV1 oncogenic activity could be regulated depending on cellular genetic context [83]. We hypothesised that SMAD2, an intracellular signal transducer of TGF-β, is likely to be involved in this mutual interaction, as we observed a concordant expression pattern for SMAD2 and ETV1, overexpressed in mPCa and down-expressed in mCRPC (Figure 8). In support of this hypothesis, it has been pointed out that, in the early stage of PCa development, the TGF-β signalling maintains cell homeostasis acting as a tumour suppressor, while, in the advanced stages of the disease, promotes cell proliferation and de-differentiation [84].

The enrichment analysis performed on MP1 and MP2 gene sets evidenced common processes and pathways associated with the dynamic switching of PCa cells between proliferative, metastatic and resistant phenotypes such as the regulation of signal transduction, transcription, differentiation and proliferation (Figure 7). This analysis also revealed other biological processes and pathways promoting cell plasticity that, under androgen deprivation therapies, can facilitate the acquisition of stem cells properties. In this regard, MP2 genes were associated with NOTCH, TGF-β and BMP pathways that controlled cell proliferation and differentiation in CRPC unlike the normal prostate tissue in which they maintain homeostasis [3,13,60,85].

As expected, among the CRPC genes, the ones up-regulated in Adeno- vs NE-CRPC contributed to androgen resistance toward aberrant AR signalling pathway, while neuroendocrine transdifferentiation was sustained by increased EZH2 and HES6 expression promoting androgen-independence and the neuroendocrine phenotype acquisition of PCa cells [2,78,80].

From the previously reported data, both gene sets were involved in sustaining the dynamic evolution of PCa and could be considered as potential targets to develop new therapeutic strategies aimed at intercepting the different stages of PCa. To this purpose, using the GSCA platform, we extracted the chemicals having a positive and significant correlation with expression of PCa- and CRPC-gene sets. Globally, our results showed that inhibitors of regulators controlling epigenetic transcription and MAPK signal transduction (Table S11) can represent promising and alternative drugs to counteract PCa progression and resistance acquisition. The findings of our study, in agreement with other experimental evidence [66,78,82], provide the rationale for testing the efficacy of epigenetic drugs in the treatment of PCa and Adeno-CRPC, as the molecular alterations observed in PCa evolution toward metastatic and resistant phenotypes mainly affect transcription factors and coregulator involved in AR activity (Figure 8, Table S7 and S9). It has also been reported that the histone deacetylase inhibitors panobinostat and vorinostat simultaneously block the AR expression and inhibit the transcription of genes under its control [45]. Conversely, treatment with molecules targeting proteins regulating BRAF and MEK signal cascade, and TGF-β pathway can be used for NE-CRPC, suggesting that these pathways had a pivotal role in sustaining neuroendocrine phenotype. In this regard, Vellano et al. [13] observed that treatment of melanoma with BRAF/MEK inhibitor improved recurrence-free survival more efficiently in female than in male patients and, using preclinical models, attributed this significant difference to AR inactivity. This intriguing result had a relevant clinical implication of considering drugs such dabrafenib, dasatinib, selumetinib and trametinib for the treatment of CRPC as monotherapy or in combination with ADT [8,60,61]. A recent study provided evidence for repositioning the antifungal ciclopirox (CPX), an iron chelating compound, for cancer therapy [65]. CPX exerted its anti-tumoral effect inducing apoptosis and inhibiting cell proliferation, migration and angiogenesis. The PCa and CRPC cells treatment with CPX impaired the (WNT)/β-catenin pathway that, in advanced tumoral forms, regulated cell plasticity inducing cancer cell stemness [66]. Therefore, it is not surprising that the “iron ion transport” was included in the enriched GO-BP terms of CRPC-gene set (Figure 7) and the CPX was listed among the drugs to be considered for PCa and Adeno-CRPC (Table S11).

Taken together, these results suggested that there is a set of genes (AR, EZH2, FOXA1, HOXB13, HOXA13, KLK3, EHF, SORD, SPON2, ADAMTS1, ETV1 and SMAD2) involved in modulating different transcriptional program that may determine phenotype-specific expression profile during PCa progression. It should be noted that among those, six genes showed context-specific essentiality, involved in sustaining PCa cell growth and differentiation (Figure 8). The therapeutic implication of these findings was clear by analysing data from the GSCA platform (Table S11), according to which 55% of drugs can inhibit these essential genes.

He et al., in an exhaustive review [46], reported the recent therapeutic advances for improving the clinical management of PCa patients, describing several chemicals, targeting cell signalling or pathways associated with specific molecular alterations, in addition to AR-signalling inhibitors representing the gold standard for the treatment of androgen-sensitive PCa forms. Interestingly, among the new drugs proposed by authors, panobinostat, vorinostat, dasatinib and temsirolimus are included in Table S9 of the present study. The other molecules we have identified can likely provide further helpful knowledge for designing or repositioning drugs, thus contributing to achieving the goal of precision medicine in PCa treatment.

5. Conclusions

This study provided a proof-of-concept for demonstrating the benefits of using an integrative bioinformatic approach, that, by joining the recent literature and data was able to highlight the master genes involved in non-mutational mechanisms essential to support cell growth and differentiation toward metastatic and androgen-independent PCa phenotypes. In our opinion, these findings may provide useful insights to develop therapeutic strategy supporting personalised medicine.