Prostate cancer

Transcriptional profiling of matched patient biopsies clarifies molecular determinants of enzalutamide-induced lineage plasticity

Heterogenous effects of enzalutamide treatment on the tumor transcriptome across matched biopsy samples

By examining the Stand Up to Cancer Foundation/Prostate Cancer Foundation West Coast Dream Team (WCDT) prospective cohort, we identified 21 patients with CRPC who underwent a metastatic tumor biopsy prior to enza and a repeat biopsy at the time of progression and whose tumor cells underwent RNA-sequencing after laser capture microdissection. All progression biopsies were performed prior to discontinuing enza, enabling us to identify resistance mechanisms induced by continued enza treatment.

The study design is shown in Fig. 1a. Patient demographic information and prior treatments are shown in Supplementary Table 1. Bone was the most common site for both pre-treatment and progression biopsies. Eighteen of 21 patients had the same tissue type biopsied at progression. In eight patients, the exact same lesion was biopsied at baseline and progression (Fig. 1b, Supplementary Table 2). The median time on enza treatment was 226 days, shorter than previous trials conducted in this same disease state6,25. PSA response at 12 weeks and the time between biopsies for each patient are shown in Fig. 1c.

Fig. 1: Study biopsy and clinical information.
figure 1

a Study schematic. b Sankey diagram showing site of biopsy at baseline (left) and at progression (right). Values indicate number of biopsies performed on each type of tissue at baseline (left) or progression (right). c Left panel shows PSA change at 12 weeks for each patient. Right panel shows time between biopsies for each patient. Response indicates whether a patient experienced a ≥ 50% reduction in PSA level at 12 weeks vs. baseline. For subject 022, PSA information was not available. Source data are provided as a Source Data file.

To understand sample-to-sample differences, we performed unsupervised hierarchical clustering and found the nearest neighbor of 13/21 (62%) baseline samples was their matched progression sample pair (Fig. 2a). Samples did not cluster together based solely on the site of biopsy, indicating laser capture microdissection removed much of the microenvironment from these samples. Furthermore, whether the same lesion was biopsied did not impact how samples clustered.

Fig. 2: The effect of enzalutamide on tumor transcriptome is heterogenous across patients.
figure 2

a Similarity heatmap for all samples clustered by variance-stabilization transformation (vst). Hashes through biopsy site indicate that the same lesion was biopsied at baseline and progression. Bracket on right axis indicates that baseline and progression samples from the same patient are nearest neighbors. b Clinical and gene expression data for each matched pair ordered on x-axis by time between biopsies. TOS is time on study in months. AR expression is by log2(TPM + 1). Clusters for each sample were assigned based on classifications from Aggarwal et al.24 and Labrecque et al.23. AR VIPER Score is the predicted AR activity score based on the AR regulon in the VIPER package26. Source data are provided as a Source Data file.

We next examined measurements of interest in all the matched samples (Fig. 2b). To estimate AR transcriptional activity, we used Virtual Inference of Protein-activity by Enriched Regulon (VIPER) master regulator analysis26. Nine (43%) patients did not have a marked difference in inferred AR activity. Nine (43%) patients had decreased AR activity, and three (14%) patients had increased AR activity at progression (Supplementary Fig. 1a). We used a second method to measure AR activity—the ARG10 signature27. ARG10 strongly correlated with the VIPER results (Supplementary Fig. 1b). Though AR-V7 expression increased in several samples at progression, the difference in expression using the entire 21-patient cohort was not statistically significant (Supplementary Fig. 1c).

Previously, Aggarwal et al. identified five clusters of CRPC tumors by RNA-sequencing analysis24. Cluster 2 was enriched for tumors with loss of AR activity and increased E2F1 activity and contained a preponderance of tumors that had lost AR expression24, consistent with lineage plasticity. A subset of cluster 2 tumor samples was labeled NEPC based upon their morphologic appearance resembling small cell prostate cancer, though many of these tumor samples did not express canonical NEPC markers such as chromogranin A (CHGA) or synaptophysin (SYP)24.

In examining the RNA-sequencing results from the baseline tumors, four of the five Aggarwal clusters were represented (clusters 1, 3, 4, and 5) in at least one sample, while no baseline sample harbored a cluster 2 program. We also applied the Labrecque transcription-based classifier that was developed on rapid autopsy CRPC samples and identified five subsets of prostate cancer: AR-driven prostate cancer (ARPC), amphicrine prostate cancer with neuroendocrine gene expression concomitant with AR signaling, AR-activity low prostate cancer, DNPC, and NEPC23. The Labrecque classifier designated all the baseline samples in our cohort as ARPC.

To determine if any of the progression tumors in our cohort underwent lineage plasticity after enza, we determined the Aggarwal cluster and Labrecque classifier designation. Twelve of 21 matched pairs did not change their Aggarwal cluster designation. However, three of the 21 progression tumors (hereafter referred to as converters) had gene expression profiles consistent with cluster 2, suggesting enza-induced conversion to an alternate lineage. We also examined the Labrecque classifier on the progression samples. The three converter samples designated as Aggarwal cluster 2 at progression were most consistent with DNPC by the Labrecque classifier, corroborating lineage plasticity in these tumors (Fig. 2b).

We next examined additional gene signatures linked previously to lineage plasticity in progression vs. baseline biopsies. Comparing samples from the three converter patients, signature scores for genes upregulated in NEPC tumors described by Beltran et al.21 were increased (Supplementary Fig. 1d). A previously described basal stemness signature28 was also activated in these three progression samples (Supplementary Fig. 1e). We previously identified a 76 gene AR-repressed gene signature that was activated in a CRPC cell line that undergoes enza-induced lineage plasticity29. This 76 gene signature was also increased in the progression samples from the three converters (Supplementary Fig. 1f). Finally, predicted AR activity was significantly decreased in the progression samples from the converters by both VIPER and ARG10 signatures (Supplementary Fig. 1a, g). In examining pre- and post-treatment samples using the entire 21-patient cohort, none of these signatures was significantly changed, demonstrating that activation of these lineage plasticity signatures was not a generalized effect of enza treatment. Altogether, these results suggest that enza-induced lineage plasticity and conversion to an AR-independent program occurs in a subset of tumors (3/21 or 14%), similar to the frequency of cluster 2 tumors (10%) described by Aggarwal previously24.

Notably, the baseline tumors from the three converter patients did not fall into the same Aggarwal cluster (cluster 4 for sample 80 and cluster 5 for samples 135, 210). Therefore, it was not surprising that the baseline tumors from these three patients did not cluster together using unsupervised clustering (Supplementary Fig. 1h, i). These data suggest that there may be different starting points to lineage plasticity with enza treatment.

Clarification of a baseline transcriptional program linked to lineage plasticity risk

To identify genes linked with risk of lineage plasticity after enza, we examined the differentially expressed genes between the three baseline samples from converters vs. the 18 non-converters. Pathway analysis implicated activation of MYC targets, E2F targets, and allograft rejection in baseline tumors from converters (Fig. 3a). There were no significantly downregulated pathways in baseline tumors from converters. To identify differentially activated transcription factors, we performed master regulator analysis. E2F1 was the top transcription factor predicted to be activated in the baseline tumors from converters (Fig. 3b, full list Supplementary Data 1), corroborating pathway analysis and our prior work demonstrating that high E2F1 activity is linked to lineage plasticity risk29. Additionally, we found that there was an upward trend in a previously described RB1 loss signature30 in the progression samples from converters, further suggesting E2F1 activation contributes to the lineage switch (Supplementary Fig. 1j). Other highly activated transcription factors in the baseline samples from converters include MYC family members and E2F4. Conversely, TP53—whose loss has been linked to lineage plasticity27,31,32—was predicted to be the most deactivated transcription factor (Fig. 3b).

Fig. 3: Pathway and master regulator analysis implicate E2F1 in lineage plasticity risk, and a signature of lineage plasticity risk identifies tumors with poor outcomes after androgen receptor signaling inhibitor treatment.
figure 3

a Hallmark pathway analysis of activated pathways in baseline samples for the three patients whose tumors converted (underwent lineage plasticity) vs. those patients whose tumors did not upon progression. b Master regulator analysis identifies top activated and deactivated transcription factors between converters and non-converters using the baseline tumor samples. Activity scores (right) and p-values (left, calculated using a gene shuffling test of the enrichment scores) were generated in the VIPER R package26. c Dot plot showing lineage plasticity signature score for patients in this cohort, the International Dream Team dataset described in Abida et al.9 and unique patients not included in this matched biopsy cohort from Alumkal et al.17. d, e Kaplan-Meier survival curves for patients in the Alumkal et al. cohort (d) and Abida et al. cohort (e) stratified by high or low lineage plasticity risk score. p-values shown were determined using the log-rank test. f Dot plot showing lineage plasticity signature score for all castration naïve adenocarcinoma PDX models described by Lin et al.22 Source data are provided as a Source Data file.

Next, we focused on identifying genes significantly upregulated in the baseline tumors from converters vs. non-converters. We identified a 14-gene signature highly activated in the three baseline tumors from converters (Supplementary Table 3). Genes in this signature include those linked to: the Wnt pathway [RNF4333 and TRABD2A34], the spliceosome [SNRPF35], and the electron transport chain [NDUFA1236 and ATP5B37]. This signature trended downwards in the progression vs. baseline biopsies from the three converters (Supplementary Fig. 2a). These results suggest that this signature is not simply identifying tumor cells that have already undergone lineage plasticity prior to enza treatment. Rather, these genes may be markers of a transition state in cells susceptible to transcriptional conversion and lineage plasticity.

Dividing the baseline samples between converters and non-converters, we defined a cut off for this 14-gene lineage plasticity risk signature that separated the groups (Fig. 3c). Additional cohorts with matched biopsies before and after enza with lineage plasticity information are lacking. However, we hypothesized that patients whose baseline tumors had high scores for this lineage plasticity risk signature would have worse outcomes. Survival data from the time of ARSI treatment were available for several CRPC cohorts whose tumors had undergone RNA-sequencing—the International Dream Team dataset9 and a prior prospective enza clinical trial led by our group17. Because a subset of the patients in that latter enza clinical trial overlapped with the patients in this current report, we focused only on patients from that clinical trial not represented in this matched biopsy cohort. Using our pre-defined 14-gene signature score cut-off from the matched biopsy cohort, we determined that high scores were associated with worse overall survival from the time of ARSI treatment in both independent datasets (p = 0.076, p = 0.005; Fig. 3d, e). Thus, high expression of the 14-gene lineage plasticity risk signature is linked to poor patient outcomes after ARSI treatment in CRPC. To determine if the lineage plasticity risk signature was activated in primary tumors, we examined the TCGA dataset38. Importantly, only two of 495 patients had high risk scores (Supplementary Fig. 2b). The lower frequency in primary tumors vs. CRPC cohorts suggests that activation of this lineage plasticity risk program may be induced by castration.

As stated previously, validation datasets with matched biopsies before and after ARSI treatment that include information on lineage at time of progression are lacking. However, previously we determined the impact of surgical castration on adenocarcinoma patient-derived xenografts (PDX)22. Nine PDXs do not undergo castration-induced lineage plasticity, while one PDX—LTL331—does and converts to a resistant version called LTL331R22. Importantly, the patient from whom the LTL331 PDX is derived had evidence of lineage plasticity in his tumor when it became castration-resistant, demonstrating this model’s fidelity22,39. Our lineage plasticity risk signature was highly activated in LTL331 vs. the other hormone-naïve PDXs that do not undergo castration-induced lineage plasticity (Fig. 3f, Supplementary Fig. 2c, d). Indeed, LTL331 was the only PDX whose lineage plasticity risk score was greater than the cut-off defined in our matched biopsy cohort (Fig. 3f). Prior work demonstrates that the exome of LTL331 is strikingly similar to its castration-induced lineage plasticity derivative, strongly suggesting that transdifferentiation—rather than clonal selection—may explain conversion in this tumor22. Finally, the lineage plasticity risk score decreased in LTL331R vs. LTL331 (Supplementary Fig. 2c)—similar to the pattern we observed in the progression vs. baseline samples from converters in our matched biopsy cohort (Supplementary Fig. 2a).

Identification of transcriptional changes in tumors undergoing lineage plasticity

Next, we sought to understand changes induced by enza between the baseline and progression samples from the three converters more deeply. The top differentially expressed genes are shown in Fig. 4a. The AR, AR target genes (KLK2, KLK3, and TMPRSS2), and the AR coactivator HOXB13 had markedly decreased expression (Fig. 4a, Supplementary Data 2). In keeping with this, progression biopsies from converters had significantly reduced expression of AR target genes from the ARG10 gene signature27 (Fig. 4b). Though we found that genes from the Beltran NEPC Upregulated signature were increased in progression samples from converters (Supplementary Fig. 1d), it is worth noting that this signature contains both canonical NEPC genes and genes not explicitly associated with acquisition of neuroendocrine features that are AR-repressed. Specifically, examining canonical NEPC markers such as SYP, CHGA, and NCAM1, we found that these genes were not highly upregulated at progression (Supplementary Data 3). Importantly, other genes linked to NEPC (SYT11, CIITA, and ETV5)21 or those normally repressed by the AR (CDCA7L, FRMD3, IKZF3, and TNFAIP2)29 were more highly expressed in the progression biopsies, suggesting that these three converter tumors may be farther along the lineage plasticity spectrum than the previously described non-neuroendocrine DNPC subtype but not as far along as de novo NEPC or NEPC found at rapid autopsy by Labrecque et al.23 that harbor a more complete neuroendocrine program.

Fig. 4: Gene expression profiling identifies gene expression changes in tumors undergoing enzalutamide-induced lineage plasticity.
figure 4

a Volcano plot showing top up and down regulated genes in progression samples vs. baseline samples for the three patients whose tumors converted (n = 3 pairs, no replicates). Adjusted p-values were calculated using the Wald test in the DESeq2 R package55. b ARG10 gene signature heatmap for three converters at baseline and progression. The left shows the expression levels of individual genes in the ARG10 signature, and the right shows the ARG10 signature score. p-value shown is for a two-tailed paired t-test between baseline and progression ARG10 scores (n = 3 pairs, no replicates). c Hallmark pathway analysis shows the top up or down regulated pathways in progression vs. baseline samples for the three patients whose tumors converted. Source data are provided as a Source Data file.

Pathway analysis between baseline and progression samples from the three converters demonstrated enrichment in several pathways, including: allograft rejection, interferon gamma response, interferon alpha response, and IL6/JAK/STAT signaling (Fig. 4c). Conversely, androgen and estrogen response—both linked to luminal differentiation—were the most downregulated, confirming loss of AR-dependence. We examined differences in gene expression between baseline and progression samples from the 18 patients whose tumors did not undergo lineage plasticity. Several of the pathways activated in the converter tumors were also activated in the non-converters—namely, interferon alpha response, interferon gamma response, and TNF-α signaling (Supplementary Fig. 3). Uniquely upregulated pathways in the converters include: allograft rejection, IL6-JAK-STAT3 signaling, inflammatory response, and complement. Uniquely downregulated pathways in the progression samples from non-converters included: E2F targets, G2M checkpoint, and hedgehog signaling. The only uniquely upregulated pathway in non-converters was protein secretion, while uniquely downregulated pathways included hedgehog signaling, G2M checkpoint, and E2F targets.

Protein expression analysis demonstrates switch to double negative prostate cancer in samples undergoing lineage plasticity

To understand the architecture of the tumors from the three converters, we used multiplex immunofluorescence (IF) with three luminal lineage markers (AR, NKX3.1, and HOXB13)—all downregulated at the mRNA level by RNA-sequencing (Fig. 4a)—and the NEPC marker INSM140. LuCaP PDX samples were used as positive and negative controls (Supplementary Fig. 4a). Matched tissue samples for multiplex IF were available for subjects 135 and 210 but not for subject 80 (Fig. 5). We identified one additional WCDT subject (103) with matched biopsies whose tumor underwent rapid clinical progression after enza treatment in the setting of a falling serum PSA—a clinical marker of AR-independence. Matched RNA-sequencing was not available for this subject, but his tumor exhibited evidence of lineage plasticity (Fig. 5). There was a spectrum of AR, NKX3.1, and HOXB13 expression in baseline samples with some cells expressing low levels of each marker, while other cells expressed higher levels. However, at progression, there was a convergence towards population-wide loss of AR, NKX3.1, and HOXB13 in each sample. We did not identify INSM1 upregulation in any of the baseline or progression tumors (Supplementary Fig. 4b). These results match our RNA-sequencing that failed to demonstrate upregulation of other canonical NEPC markers (Supplementary Data 3) and that characterized the three converter samples as DNPC by the Labrecque classifier, rather than NEPC (Fig. 2b).

Fig. 5: Multiplex immunofluorescence demonstrates switch to double negative prostate cancer in samples undergoing lineage plasticity.
figure 5

Patients 135, 210, and an additional West Coast Dream Team patient whose tumor converted, patient 103, were stained for AR, NKX3.1, and HOXB13 expression (n = 6 biologically independent samples with no replicates). The scale bar represents 50 µm. Signal intensity values for each marker are shown with the median value indicated. Signal intensity values were compared between each matched pair using the Mann-Whitney two-tailed test with p < 1 × 10−15 for each comparison except DTB_103 HOXB13, which is p = 1.8 × 10−14. Source data are provided as a Source Data file.

Conservation of DNA mutations in tumors undergoing lineage plasticity

Finally, to determine if the progression samples from converters represented distinct clones with unique genetic alterations vs. baseline, we performed DNA mutation and copy number analysis. For subjects 80 and 103, the same tumor lesion was biopsied at baseline and progression. DNA-sequencing of these biopsies showed identical DNA mutations. For subjects 135 and 210, matched metastatic biopsy DNA-sequencing was unavailable. However, cell-free DNA was available. DNA-sequencing of cell-free DNA samples showed that mutations and copy number alterations were conserved between baseline and progression samples (Table 1).

Table 1 DNA sequencing of matched samples from converters demonstrates conserved alterations

Loss of the tumor suppressor genes TP53, RB1, and PTEN has been linked to lineage plasticity risk in pre-clinical models31,32. However, we do not know if the presence of these genomic abnormalities in patient tumors is associated with the risk of lineage plasticity to DNPC. One of the three converter patients (subject 80) was found to have an inactivating PTEN mutation and a second patient (subject 103) had RB1 loss, but none were found to have compound TP53/RB1/PTEN loss. When available, we also examined TP53/RB1/PTEN status for tumors from the Abida et al.9 and Alumkal et al.17 cohorts that had high lineage plasticity risk scores. Of the seven high lineage plasticity risk score tumors examined from these two validation cohorts, only two tumors had loss of two or more of the genes TP53, RB1, and PTEN (Supplementary Table 4). DNA-sequencing of matched metastatic biopsies for the cohort as a whole is shown in Supplementary Table 5.