Gene Expression Signature Predictive of Neuroendocrine Transformation in Prostate Adenocarcinoma

Abstract

Neuroendocrine prostate cancer (NEPC) can arise de novo, but much more commonly occurs as a consequence of a selective pressure from androgen deprivation therapy or androgen receptor antagonists used for prostate cancer (PCa) treatment. The process is known as neuroendocrine transdifferentiation. There is little molecular characterization of NEPCs and consequently there is no standard treatment for this kind of tumors, characterized by highly metastases rates and poor survival. For this purpose, we profiled 54 PCa samples with more than 10-years follow-up for gene and miRNA expression. We divided samples into two groups (NE-like vs. AdenoPCa), according to their clinical and molecular features. NE-like tumors were characterized by a neuroendocrine fingerprint made of known neuroendocrine markers and novel molecules, including long non-coding RNAs and components of the estrogen receptor signaling. A gene expression signature able to predict NEPC was built and tested on independently published datasets. This study identified molecular features (protein-coding, long non-coding, and microRNAs), at the time of surgery, that may anticipate the NE transformation process of prostate adenocarcinoma. Our results may contribute to improving the diagnosis and treatment of this subgroup of tumors for which traditional therapy regimens do not show beneficial effects.

Keywords: estrogen signaling pathway; long non coding RNAs; neuroendocrine prostate cancer; neuroendocrine transdifferentiation; predictive gene signature; prostate cancer.

Figure 1

A Venn diagram showing the intersection between different clinical features of patients with prostate adenocarcinomas (left panel). Tumors from patients who satisfied at least one out four clinical criteria (indicated as 1/4 in the right panel) were intersected with tumors with high levels of Chromogranin A expression (right panel).

Figure 2

Volcano plots visualizing differential expressed genes (A) and miRs (B) obtained by class comparison analysis. Blue and red dots represent down- or up-regulated features, respectively. Symbols for the top 10 up- and down-regulated genes and miRNAs are displayed. (C) A circular plot of the human genome divided into three concentric circles, each representing differential expression at different molecular levels (pcRNA, ncRNA, miRNA). Red and green bars represent the position of each up- or down-regulated molecular feature along chromosomes, respectively. (D) Expression profile of the genes involved in NE prostate cancer induction and maintenance and concomitantly up-regulated in NE-like vs. AdenoPCa samples (upper part of the heatmap) and of the genes whose loss is implicated in prostate cancer invasion, metastasis formation and worse prognosis and concomitantly down-regulated in NE-like vs. AdenoPCa samples (lower part of the heatmap). Z-score values are reported (logIntensities were median-centered and divided by standard deviation).

Figure 3

(A) A selection of pathways derived from functional enrichment analysis of DEGs using Metacore (orange and light-blue bars) and David analyses (red and blue bars). The negative logarithm of the p-value (base 10) is reported on the X-axis. FDR p-values range from 3.31 × 10⁻⁵ to 0.9. The number of genes involved in each process is indicated in brackets. The complete list of pathways and processes is available in Supplementary Table S1. (B) The network of direct interactions starting from or pointing to AR and ESR1 proteins. More in detail, 70 connections started from AR, and 12 pointed to AR, while 75 links started from ESR1, and 26 pointed to ESR1. The objects that are over- or underexpressed in the data are marked with red/blue circles, respectively. Different edge colors were used to represent the activation (green), inhibition (red), and unspecified (grey) effects. (C) Transcription factor interactome analysis pointed out several enriched transcription factors that may putatively control a significant number of up- or down-regulated genes in the NE-live vs. AdenoPCa comparison.

Figure 4

(A–D) PDX model of NEPC progression evaluated from prostate adenocarcinoma before castration (LTL331) to relapsed neuroendocrine tumor (LTL331R) [16]. ESR1 (A) and CHGA (B) increased their expression throughout the transdifferentiation process, with concomitant AR down-regulation (C). (D) Correlation matrix showing a high positive correlation between ESR1 and CHGA (Pearson’s correlation coefficient = 0.86, p-value = 7.27 × 10⁻⁵) and a negative correlation between ESR1 and AR (Pearson’s correlation coefficient = −0.89, p-value = 1.54 × 10⁻⁵) and between AR and CHGA (Pearson’s correlation coefficient = –0.91, p-value = 6.11 × 10⁻⁶). (E) The overlap between up-regulated genes in NE-like vs. AdenoPCa samples and genes up-regulated in OSC11 upon treatment with estradiol (E2) or genes down-regulated in OSC11 upon treatment with DHT, with extrapolation of some genes of interest.

Figure 5

(A) Histogram depicting Pearson’s correlation coefficient for a subset of genes whose expression was positively correlated with that of HOTAIR (Pearson score > 0.3 and p-value < 0.05). (B) Histogram depicting Pearson’s correlation coefficient for a subset of prostate cancer-related genes whose expression was negatively correlated with that of HOTAIR (Pearson score < −0.3 and p-value < 0.05). (C) Venn diagram of up-regulated genes in NE-like vs. AdenoPCa samples and in hormone-deprived LNCaP cells treated with HOTAIR alone or HOTAIR plus R1881, compared to correspondent controls, with extrapolation of some genes of interest. (D) Direct interaction network between AR, ESR1, HOTAIR and miR-31-5p. Red and green edges represent inhibition and activation mechanisms, respectively. Grey edges depict undefined/unknown regulations.

Figure 6

(A) Unsupervised hierarchical clustering on 15 CRPC-NE and 34 CRPC-Adeno samples from Beltran et al. [4], using the list of up- and down-regulated genes identified in our cohort of NE-like vs. AdenoPCa samples. Pearson correlation as distance metric and average linkage method were used. (B,C) Receiver operating characteristic (ROC) curve of the score tested on Beltran dataset (B) and on GSE66187 dataset (C).