Drug-Induced Epigenomic Plasticity Reprograms Circadian Rhythm Regulation to Drive Prostate Cancer toward Androgen Independence

Androgen ablation has been the mainstay treatment for patients with metastatic prostate cancer ever since the direct critical link between androgens and prostate tumor progression was first described (1). The androgen receptor (AR) is the key driver of prostate cancer development and progression, and multiple therapeutic strategies have been developed over the years to effectively block the activity of this hormone-driven transcription factor (TF). Upon androgen binding, AR is associated with the chromatin at distal cis-regulatory enhancer elements, where it regulates the expression of genes through long-range chromatin interactions in three-dimensional genomic space (2, 3). AR does not operate in isolation but rather recruits a large spectrum of coregulators and other TFs to promote the expression of genes that drive cancer cell proliferation (4). Critical AR interactors in the transcription complex are HOXB13 and FOXA1, which are both upregulated in primary prostate cancer (4–6) and demarcate enhancers that drive not only primary tumorigenesis but also metastatic disease progression (7). Mechanistically, FOXA1 acts as a pioneer factor, rendering the chromatin accessible for AR to bind (8–11). FOXA1 is frequently mutated in prostate cancer (12–16), which was shown to alter its pioneering capacities, perturb luminal epithelial differentiation programs, and promote tumor growth, further highlighting the critical role of FOXA1 in human prostate tumors (17, 18).

Most patients are diagnosed with organ-confined prostate cancer, which can potentially be cured through locoregional therapies, such as surgery (radical prostatectomy), radiotherapy, and/or brachytherapy (19). However, approximately 30% of these patients experience a biochemical recurrence (BCR)—a rise in prostate-specific antigen (PSA) serum levels—indicating prostate cancer relapse (20). At this stage of the disease, suppression of androgen production is a commonly applied therapeutic intervention that can delay further cancer progression for years (21, 22). Nevertheless, the development of resistance to androgen deprivation is inevitable, resulting in castration-resistant prostate cancer (CRPC) for which there is no cure (23). Most CRPC tumors acquire molecular features that enable active AR signaling despite low circulating androgen levels, a finding that led to the development of several highly effective AR-targeted therapies. Enzalutamide (ENZ) is one of the most frequently used AR-targeting agents, which functions through a combined mechanism of blocked AR nuclear import, diminished AR chromatin binding, and decreased transcription complex formation, effectively impairing AR-driven prostate cancer growth (24). ENZ's potent antitumor activity has been demonstrated in multiple clinical trials, which led to its FDA approval in various prostate cancer disease stages—from metastatic CRPC (mCRPC; refs. 25, 26) to metastatic hormone-sensitive prostate cancer (mHSPC; ref. 27) and even nonmetastatic CRPC (28)—illustrating how AR-targeted therapies are being progressively introduced earlier in clinical practice. A clinical benefit of ENZ monotherapy as a neoadjuvant treatment prior to prostatectomy for patients with localized disease has not been established. Although effective in the CRPC setting, resistance to AR pathway inhibition will ultimately develop, and the management of advanced prostate cancer with this acquired resistance remains a major clinical challenge, especially since the underlying mechanisms are still not fully elucidated (29). Therefore, furthering our understanding of how ENZ affects prostate cancer biology may lead to the identification of acquired cellular vulnerabilities that could be therapeutically exploited.

To study global drug-induced transcriptional and epigenetic plasticity in human prostate tumors and identify cellular adaptation mechanisms to evade drug treatment, we designed a phase II clinical trial to perform multiomics studies in pre- and posttreatment samples from patients with high-risk localized prostate cancer treated with neoadjuvant ENZ monotherapy. We identified transcriptional reprogramming after treatment, with the deactivation of AR signaling and an activation of cell plasticity with neuroendocrine (NE)-like features upon 3 months of AR suppression. Posttreatment, these tumors harbored a distinct set of 1,430 de novo occupied FOXA1-positive cis-regulatory elements, positive for—yet independent of—AR activity, which are dictated by the circadian clock core regulator ARNTL to drive tumor cell proliferation instead. Using ARNTL knockout experiments, we could restore ENZ sensitivity in treatment-resistant cell line and xenograft models, revealing an unexpected biological interplay between hormonal resistance and circadian rhythm regulation, and identifying a novel, highly promising candidate drug target in the clinical management of primary high-risk prostate cancer.

To study how early ENZ intervention affects prostate tumor biology in a noncastrate environment, we performed integrative multiomics analyses as part of a single-arm, open-label, phase II clinical trial: the DARANA study (Dynamics of Androgen Receptor Genomics and Transcriptomics After Neoadjuvant Androgen Ablation; ClinicalTrials.gov identifier, NCT03297385). In this trial, 56 men with primary high-risk (Gleason score ≥7) prostate cancer were enrolled (Fig. 1A). Patient demographics and disease characteristics are summarized in Table 1, and the clinical outcomes of this study are discussed in Supplementary Fig. S1A–S1F, Supplementary Table S1, and Supplementary Data. Prior to ENZ therapy, MRI-guided core needle tumor biopsies were taken, hereafter referred to as the pretreatment setting. Subsequently, patients received neoadjuvant ENZ treatment (160 mg/day) without additional androgen deprivation therapy (ADT) for 3 months, followed by robotic-assisted laparoscopic prostatectomy. Based on baseline MRI information and palpation, additional tumor-targeted core needle biopsies were taken ex vivo—representing the posttreatment setting. This pre- and posttreatment sampling allowed us to study the epigenetic, genomic, transcriptomic, and proteomic effects of neoadjuvant ENZ therapy in individual patients (Fig. 1A). We generated chromatin immunoprecipitation [ChIP sequencing (ChIP-seq)] profiles of the prostate cancer drivers AR and FOXA1 as well as the histone modification H3K27ac before and after ENZ treatment, and we integrated these cistromic findings with pre- and posttreatment gene expression [RNA sequencing (RNA-seq)], copy number [copy-number variation sequencing (CNV-seq)], and immunohistochemistry (IHC) data from the same tumors. Stringent quality control (QC) analyses were performed on all data streams (Supplementary Fig. S2), and the following number of samples passed all QC measures (Fig. 1B): AR ChIP-seq (pre: n = 10; post: n = 12), FOXA1 ChIP-seq (pre: n = 17; post: n = 17), H3K27ac ChIP-seq (pre: n = 24; post: n = 23), CNV-seq (pre: n = 24; post: n = 24), RNA-seq (pre: n = 42; post: n = 52), and IHC (post: n = 51).

Collectively, we performed integrative multiomics analyses as part of a clinical trial that enabled us to examine ENZ-induced oncogenomic changes to identify early epigenetic steps in treatment response, but also therapy-induced resistance.

To assess how neoadjuvant ENZ treatment affects the cis-regulatory landscape in primary prostate cancer, we generated human tumor ChIP-seq profiles for the TFs AR and FOXA1 along with the active enhancer/promoter histone mark H3K27ac before and after neoadjuvant intervention. ChIP-seq quality metrics are summarized in Supplementary Fig. S3A–S3E and Supplementary Table S2. Visual inspection of known AR target genes showed high-quality data for all ChIP factors in both clinical settings (Fig. 2A). On a genome-wide scale, the H3K27ac ChIP-seq profiles were highly distinct from the TFs and divided the samples into two main clusters irrespective of their treatment status (Fig. 2B and C). Notably, the AR and FOXA1 ChIP-seq data sets were intermingled in the clustering analysis, suggesting largely comparable binding profiles, which is in line with FOXA1's role as a canonical AR pioneer factor (Supplementary Fig. S4; refs. 5, 30). As described previously (31), the highest Pearson correlation was found between H3K27ac samples, indicating comparable histone acetylation profiles among primary prostate cancer samples (Fig. 2B; Supplementary Fig. S4). Much greater heterogeneity in chromatin binding was observed for the TFs AR and FOXA1, which was further supported by the steep decrease in the number of overlapping AR and FOXA1 peaks with an increasing number of samples compared with H3K27ac (Fig. 2D; Supplementary Fig. S4). Heterogeneity was comparable when separately analyzing pre- versus posttreatment specimens and in the same order of magnitude as compared with previous reports describing TF cistromics and epigenomics in clinical samples (refs. 31, 32; Supplementary Fig. S5A and S5B) with a comparable overlap of peaks for AR and FOXA1 (Supplementary Fig. S5C and S5D). In order to maintain the high-confidence peaks that have been reproducibly identified in multiple patients without losing too much binding site heterogeneity between samples, we decided to generate consensus peak sets. To this end, we only considered binding sites that were present in at least 3 of 22 AR samples, 7 of 34 FOXA1 samples, and 13 of 47 H3K27ac samples, which corresponds to ∼25% of all binding sites identified for each factor (Fig. 2D). Genomic distribution analyses of these consensus sites revealed distinct enrichments for annotated genomic regions. Although AR and FOXA1 were almost exclusively found at intronic and distal intergenic regions, H3K27ac peaks were also enriched at promoters (Fig. 2E), which is in line with previously published genomic distributions of AR (5, 31), FOXA1 (5, 9), and H3K27ac (31, 33). In addition, motif enrichment analyses at AR and FOXA1 consensus peaks identified, as expected, androgen and Forkhead response elements among the top-ranked motifs, respectively (Fig. 2F). Analyses on correlations between factors (Fig. 2B–D), genomic distributions (Fig. 2E), and motif enrichment (Fig. 2F) were repeated for the pretreatment samples exclusively, supporting the same conclusions (Supplementary Fig. S6A–S6E).

Taken together, we generated multiple high-quality tissue ChIP-seq data streams that then allowed us to study ENZ-induced changes in patients with primary prostate cancer.

To identify ENZ-induced TF reprogramming and epigenetic changes, we performed differential binding analyses comparing the pre- and posttreatment tissue ChIP-seq samples. Therefore, we first ran occupancy-based unsupervised principal component analyses (PCA) to detect whether ENZ treatment led to differences in TF chromatin binding. Although the sample size of the AR ChIP-seq data stream was not sufficient to observe significant differences in peak occupancy before versus after treatment (Supplementary Fig. S7A), the FOXA1 data did show such differences, with a clear separation of pre- and posttreatment FOXA1 samples in the second principal component (Fig. 3A). Subsequent supervised analysis (pre vs. post) revealed a total of 1,905 genomic regions [475 pretreatment-enriched (pre-enriched), 1,430 posttreatment-enriched (post-enriched); Supplementary Table S3] that showed significant differential FOXA1 binding between both clinical settings [false discovery rate (FDR) <0.05; Fig. 3B and 3C; Supplementary Fig. S7B–S7D]. Further characterization of these differential FOXA1 regions showed that both sets of binding sites were still preferentially located in intronic and distal intergenic regions (with a slight enrichment for promoters at the post-enriched sites; Supplementary Fig. S7E). In addition, Forkhead domain family motifs were the top enriched motifs at both pre- and post-enriched sites, illustrating that treatment does not alter FOXA1 motif preference and still occupies canonical FOXA1 binding sites (Supplementary Fig. S7F).

To examine whether structural variations were underlying these differential FOXA1 binding events, we performed CNV-seq on the same tumor specimens and then projected onto the differential FOXA1 cistromics the structural copy-number data. These analyses revealed a comparable level of CNV at pre- and posttreatment enriched FOXA1 sites before and after ENZ treatment, with an overall trend toward less CNV upon treatment (Supplementary Fig. S8A–S8C). However, in none of the matched sample pairs (pre– and post–CNV-seq and FOXA1 ChIP-seq; n = 15) was a strong correlation between copy-number difference and ChIP-seq signal difference observed (R = 0.11; Supplementary Fig. S8D). In total, at only 44 of 1,905 differential FOXA1 binding sites (<2.5%), we observed copy-number differences between post- and pretreatment samples that could potentially explain binding site occupancy in 3 or more patients, indicating that the vast majority of these differential binding events is based on treatment-induced TF reprogramming rather than structural variation (Supplementary Fig. S8E).

As FOXA1 dictates AR chromatin binding capacity (5), the epigenetic plasticity of FOXA1 induced by treatment may be associated with alterations in the AR cistrome. To assess this, and to explore the epigenetic landscape surrounding the differentially bound FOXA1 regions, we compared the ChIP-seq signal of all three factors (AR, FOXA1, and H3K27ac) at differential (pre-/post-enriched) and consensus (shared by ≥30 patients; n = 338) FOXA1 sites before and after ENZ therapy. Although the FOXA1 ChIP-seq signal was highest at consensus binding sites, the pre- and posttreatment enriched regions followed the expected trend and showed significantly higher signals in the corresponding settings (Fig. 3D). Notably, we also observed less binding of FOXA1 to consensus sites when treated with ENZ, although the differences were much milder compared with the effects seen at pre-enriched FOXA1 sites (P_adj = 3.62 × 10⁻²² at consensus vs. 3.76 × 10⁻¹³⁰ at pre-enriched sites, Mann–Whitney U test; Fig. 3D; Supplementary Fig. S9A). This could be explained by decreased FOXA1 gene expression levels upon ENZ treatment (Supplementary Fig. S9B). The AR ChIP-seq signal followed the same patterns as observed for FOXA1, suggesting that relocated FOXA1 upon treatment functionally drives alterations in the AR cistrome (Fig. 3D). Unexpectedly, the pre-enriched FOXA1 sites were completely devoid of any H3K27ac signal in both pre- and posttreatment samples, although the post-enriched counterparts were positive for this active enhancer/promoter mark with a significant increase post-ENZ (P_adj = 5.59 × 10⁻⁴, Mann–Whitney U test; Fig. 3D; Supplementary Fig. S9C and S9D), suggesting that pre-ENZ FOXA1 sites are inactive. To validate these observations in an independent cohort, we analyzed previously published AR (n = 87), H3K27ac (n = 92), and H3K27me3 (n = 76) ChIP-seq data from a cohort of 100 primary treatment-naive prostate cancer samples (31). Supporting our previous analyses, the vast majority of post-enriched FOXA1 sites were H3K27ac-positive and their histone acetylation status positively correlated with AR binding (R = 0.78; Fig. 3E; Supplementary Fig. S9E). The pre-enriched FOXA1 sites, however, were again H3K27ac-negative, while the repressive histone modification H3K27me3 was present, which further points toward an inactive epigenetic state of these regulatory regions (Fig. 3E).

Recently, we reported that prostate cancers can reactivate developmental programs during metastatic progression (7). These sentinel enhancers appeared to be premarked by FOXA1 from prostate gland development, and albeit inactive in normal prostate and primary tumor specimens, the sites get reactivated by AR during metastatic outgrowth. Given the inactivity of the pre-enriched FOXA1 sites, we hypothesized that FOXA1 might be decommissioned at such developmental enhancers prior to hormonal intervention. To test this, we overlapped the differential FOXA1 binding sites with the metastasis-specific AR binding sites (met-ARBS; n = 17,655), which revealed a strong enrichment for these developmental regulatory elements at pretreatment FOXA1 sites (P = 2.13 × 10⁻¹⁶, Fisher exact test; Supplementary Fig. S9F). But are the inactive preenriched FOXA1 sites solely epigenetically suppressed, or are these regions intrinsically incapable of being active in this cellular context? To address this question and to further elucidate the role of AR at these differentially bound FOXA1 sites, we integrated our tissue ChIP-seq findings with previously identified tumor-specific ARBS (n = 3,230; ref. 5) that were functionally characterized using self-transcribing active regulatory regions sequencing (STARR-seq), a massive parallel reporter assay to systematically annotate intrinsic enhancer activity (34). With this, three distinct classes of ARBS were identified (Supplementary Table S4): enhancers that were active regardless of AR stimulation (constitutively active; n = 465), ARBS with no significant enhancer activity (inactive; n = 2,479) and inducible AR enhancers that increase activity upon androgen treatment (inducible; n = 286). Interestingly, we found that posttreatment FOXA1 sites were enriched for constitutively active ARBS, which further supports the high enhancer activity and H3K27ac positivity observed at these sites, but also illustrates that this activity is constitutive and AR independent (Fig. 3F). Consistent with our postulated inactivity of the pretreatment enriched FOXA1 sites, these regions overlapped highly significantly with inactive ARBS (P = 8.60 × 10⁻⁹, Fisher exact test), which implies that these DNA elements are intrinsically inactive and incapable of acting as functional enhancers, and possibly explains why these AR-bound sites did not show active regulatory marks (Fig. 3E and F). As no enrichment of our differential FOXA1 sites was observed with inducible ARBS (pre-enriched: 4/475; post-enriched: 2/1,430), these data further support a conclusion that AR itself is not a driver at FOXA1 sites that are differentially occupied after ENZ exposure in patients.

Overall, these results suggest that prior to hormonal intervention, FOXA1 is decommissioned at inactive developmental enhancer elements, which based on their primary DNA sequence are intrinsically incapable of being active—at least in the tested hormone-sensitive disease setting. However, upon ENZ treatment, FOXA1 gets reprogrammed to highly active cis-regulatory regions, which act in an AR-independent manner.

Having assessed the cistromic and epigenomic changes in response to neoadjuvant ENZ, we next determined how transcriptional programs were affected by this hormonal intervention. PCA across both treatment states revealed that 3 months of ENZ therapy has a major effect on global gene expression profiles (Fig. 4A). Subsequently, we performed differential gene expression analysis in which we compared pre- and posttreatment RNA-seq samples. Gene set enrichment analysis (GSEA) showed that AR signaling, along with mitosis and MYC signals, was strongly decreased upon treatment (Fig. 4B and C; Supplementary Fig. S10A). As ENZ blocks the AR signaling axis, we analyzed the androgen response pathway in more detail, which revealed a strong downregulation of AR target genes in almost every patient (Fig. 4D). In contrast to this, TNFα signaling, IFNγ response, and epithelial–mesenchymal transition (EMT) signals were the most upregulated (Fig. 4B; Supplementary Fig. S10B).

Previously, we identified three distinct subtypes of primary treatment-naive prostate cancer (31), which we named clusters 1 to 3 (Cl1–3). Although Cl1 and Cl2 were mainly dominated by their ERG fusion status—with Cl1 expressing high ERG levels (ERG fusion–positive) and Cl2 expressing low ERG levels (ERG fusion–negative)—Cl3 was enriched for NE-like features, including low AR activity and a high NE gene expression score. To assess the impact of neoadjuvant ENZ therapy on these prostate cancer subtypes, we performed unsupervised hierarchical clustering in pre- and posttreatment settings using the originally identified top 100 most differentially expressed genes per cluster. Prior to hormonal intervention, we could robustly assign the samples into all three clusters (Cl1: n = 23, Cl2: n = 11, Cl3: n = 8) with highly comparable distributions, as we previously reported in another cohort of patients (ref. 31; Supplementary Fig. S11A). Our pre- and posttreatment sampling now allowed us to investigate how individual tumors were affected by neoadjuvant therapy. This revealed that three months of ENZ therapy pushed almost all of the tumors toward our NE-like Cl3 (Fig. 4E; Supplementary Fig. S11B). To ensure that the observed effects are not solely driven by the treatment-induced reduction in AR activity (Fig. 4C and D), we used a well-established NE prostate cancer (NEPC) signature (35) to calculate gene expression fold changes (FC) pre- versus post-ENZ, which confirmed an induction of NE-like signaling upon treatment (Fig. 4F). We further validated this transcriptional rewiring using gene sets that distinguish the three major lineages of prostate epithelial cells (luminal, basal, and NE; refs. 36, 37), which jointly illustrated reduced AR-driven luminal cell transcriptional activity accompanied by an enrichment of NE-like features along with a basal-type transcriptional program after treatment (Supplementary Fig. S12A). Along these lines, classic NEPC markers (38) and transcriptional disease drivers (39–41) were selectively upregulated upon treatment (CHGA, PEG10) with the acquisition of promoter-enriched H3K27ac (Supplementary Fig. S12B–S12D), while others were not affected on expression level (SYP, N-MYC) or not even expressed in primary tumors—irrespective of neoadjuvant treatment status (BRN2, encoded by the POU3F2 gene). For the classic NEPC IHC markers chromogranin A (CHGA) and synaptophysin (SYP), tissue microarrays (TMA) were stained and analyzed, showing no change (SYP) or a modest nonsignificant increase (CHGA) upon neoadjuvant ENZ treatment (Supplementary Fig. S12E).

Recently, N-MYC ChIP-seq data were reported in models of NEPC (40), which showed a limited overlap with our posttreatment FOXA1 cistrome (Supplementary Fig. S12F). Although a subset of NEPC markers were enriched at the RNA-seq level, FOXA1 reprogramming did not seem to be a crucial driver in this phenomenon based on the limited overlap of our differential FOXA1 cistromes with a recently reported NEPC FOXA1 cistrome (ref. 42; Supplementary Fig. S13A), nor was FOXA1 ChIP-seq in our study enriched for classic NEPC signature genes (Supplementary Fig. S13B). Jointly, these data suggest that altered FOXA1 cistromics after neoadjuvant ENZ treatment present a different biological state as compared with the fully developed NEPC-associated FOXA1 cistrome that presents in the advanced disease stage and may represent an early intermediate state.

FOXA1 is frequently mutated in primary prostate cancer (14) and metastatic disease, in which FOXA1 mutations were associated with loss of lineage-specific transcriptional programs and worse clinical outcomes (17, 18). Therefore, we determined the FOXA1 mutation status of our clinical samples using H3K27ac ChIP-seq and RNA-seq reads covering the FOXA1 gene (43) and tested for possible enrichment for poor clinical outcome and NE-like gene expression features specifically in the FOXA1-mutant cases. Although we observed a significant enrichment of FOXA1-mutant tumors among ENZ nonresponders (BCR ≤6 months after surgery), no such enrichment was observed at the transcriptomic level, likely affected by the almost-complete transition of all our tumor samples toward the NE-like Cl3—irrespective of FOXA1 mutation status (Supplementary Fig. S13C–S13F).

Collectively, these results demonstrate that 3 months of neoadjuvant ENZ therapy not only uniformly diminish AR signaling but also push practically all of our primary prostate cancer samples to acquire some—but not all—features of NEPC, independent of their original subtype.

Having examined the global cistromic and transcriptomic changes upon ENZ therapy, we next sought to characterize the biological consequences of the observed FOXA1 reprogramming using integrative analyses. We hypothesized that the newly acquired FOXA1 sites would be driving the expression of genes associated with tumor cell survival programs. Using H3K27ac-centric chromatin confirmation (HiChIP) data generated in LNCaP cells (44), pre- and posttreatment FOXA1 sites were coupled to their corresponding gene promoters (Supplementary Table S5). Subsequently, genome-wide CRISPR knockout screen data from Project Achilles (DepMap 20Q1 Public; VCaP) were used to identify those genes essential for prostate cancer cell proliferation (45, 46). Although genes associated with pretreatment FOXA1 sites were not enriched for essential genes (gene effect score <−1), genes under the control of posttreatment FOXA1 sites showed a significant enrichment (P = 8.66 × 10⁻⁸, Fisher exact test) for critical drivers of tumor cell proliferation (Fig. 5A), pointing toward a possible role for these sites in maintaining proliferative potential upon ENZ treatment. However, the factor regulating these genes to possibly drive proliferation remained elusive, especially since based on our STARR-seq and RNA-seq data, AR is likely not driving enhancer activity at posttreatment FOXA1 sites (Fig. 3F; Fig. 4C and D). Therefore, we sought to identify TFs involved in the activation of these regulatory regions that are selectively occupied by FOXA1 following treatment. To this end, we overlaid the genomic coordinates of the posttreatment enriched FOXA1 binding sites with those identified in publicly available ChIP-seq data sets (n = 13,976) as part of the Cistrome Data Browser TF ChIP-seq sample collection (47, 48). Besides FOXA1 and AR, which were expected to bind at these regions (Fig. 3D), we also identified the glucocorticoid receptor (GR; encoded by the NR3C1 gene), which has previously been described to be upregulated upon antiandrogen treatment and able to drive the expression of a subset of AR-responsive genes, conferring resistance to AR blockade (49–51). Unexpectedly, the second most enriched TF after FOXA1 was the circadian rhythm core component ARNTL (also known as BMAL1), which has not previously been implicated in prostate cancer biology (Fig. 5B). Interestingly, ARNTL transcript levels were upregulated upon ENZ treatment (P = 6.4 × 10⁻³, Mann–Whitney U test; Fig. 5C), which was accompanied by increased H3K27ac ChIP-seq signals at the ARNTL locus (Supplementary Fig. S14A). Consistent with this, TMA IHC analysis also revealed elevated ARNTL protein levels after treatment when comparing the prostatectomy specimens post-ENZ with those of matched untreated control patients (P = 6.89 × 10⁻¹⁹, Fisher exact test; Fig. 5D). To assess whether ARNTL levels are also associated with patient outcome, we compared the average ARNTL gene expression of patients who did not experience a BCR (responders, n = 29) with those that experienced an early BCR within ≤6 months after surgery (nonresponders, n = 8; Supplementary Table S1). Although pretreatment ARNTL levels were not significantly different between ENZ responders and nonresponders, high ARNTL levels after treatment were associated with poor clinical outcomes (P = 4.79 × 10⁻³, Mann–Whitney U test; Fig. 5E). In agreement with this observation, ARNTL levels were exclusively found upregulated in nonresponders (P = 3 × 10⁻⁴, paired Mann–Whitney U test), while overall remaining unaffected upon neoadjuvant ENZ treatment in responders (P = 0.33; Supplementary Fig. S14B).

Interestingly, while the CLOCK and NPAS2 proteins, which form a heterodimer with ARNTL to activate transcription of core clock genes, did not show differential gene expression upon ENZ treatment (Supplementary Fig. S14C), all downstream ARNTL targets were upregulated upon treatment—except for CRY1, which has recently been shown to be AR- and thus ENZ-responsive (52). In addition, the gene expression of these ARNTL dimerization partners was also not associated with clinical outcomes (Supplementary Fig. S14D), hinting toward a treatment-induced role of ARNTL that is independent of its canonical function in the circadian machinery.

Notably, in two cohorts of mCRPC (53, 54), ARNTL expression was not associated with outcome (Supplementary Fig. S15A–S15D), suggesting a context-dependent prognostic potential of this gene being associated with outcome in high-risk primary prostate cancer upon treatment with ENZ.

Taken together, these data suggest that the circadian clock regulator ARNTL may be functionally involved in ENZ resistance in high-risk primary prostate cancer by driving tumor cell proliferation processes.

To further investigate the relevance of ARNTL as a transcriptional driver at posttreatment FOXA1 sites, we performed in vitro validation experiments. To this end, we used hormone-sensitive LNCaP prostate cancer cells, which we cultured either in full medium alone (pre^LNCaP) or with ENZ for 48 hours (post^LNCaP), mimicking our clinical trial setting (Fig. 6A). Based on the ENZ-induced acquisition of NE-like gene expression profiles in our patient cohort (Fig. 4E and F), we also included the ENZ-resistant LNCaP-42D model (41) that possesses NE features (Res^LNCaP-42D; Fig. 6A), allowing us to further validate our patient-derived findings in cell lines recapitulating the transcriptional features of posttreatment clinical specimens.

We performed FOXA1 ChIP-seq experiments in all three cell line conditions (Supplementary Fig. S16A–S16D; Supplementary Table S6), which revealed highly similar FOXA1 chromatin binding dynamics as observed in our clinical samples: Although the pre-enriched FOXA1 sites identified in vivo showed less binding upon treatment, we observed that merely 48 hours of ENZ exposure was sufficient to strongly induce binding at post-enriched sites, which was further increased in the long-term exposed, treatment-resistant LNCaP-42D cell line (Supplementary Fig. S16E and S16F). Similarly, genome-wide correlation analyses indicated that short-term ENZ treatment in cell lines induced FOXA1 reprogramming to regions that are FOXA1 bound in treatment-resistant but not in treatment-naive cells (Supplementary Fig. S16G and S16H).

Having shown that differential FOXA1 chromatin binding in tumors could be recapitulated in vitro, we next sought to further assess the role of ARNTL in these preclinical models. Therefore, we first measured the intrinsic enhancer activity of our patient-derived and cell line–validated differential FOXA1 binding sites by STARR-seq for 1,209/1,905 differential regions in LNCaP cells. Notably, we identified a subset of regions (n = 968) with sustained enhancer activity upon ENZ treatment (Supplementary Fig. S17A), confirming our initial STARR-seq analysis (Fig. 3F). Although GIGGLE analyses on the inactive regions showed enrichment for FOXA1 and AR, active enhancers—irrespective of treatment—were specifically enriched for ARNTL (Supplementary Fig. S17B and S17C). These data are in full concordance with the tumor H3K27ac ChIP-seq (Fig. 3D) analyses, showing AR-independent activity at the posttreatment enriched FOXA1 sites, and uncovered once more ARNTL as a possible driver for transcriptional activity in the case of AR suppression.

Next, we confirmed that treatment with ENZ increased ARNTL protein levels in prostate cancer models (Supplementary Fig. S18A), recapitulating the clinical observations (Fig. 5C and D). Interestingly, this treatment-induced ARNTL upregulation appeared to be FOXA1-dependent, as FOXA1 knockdown abolished the ENZ-driven increase in ARNTL levels (Supplementary Fig. S18B).

As cistromic ARNTL profiling has to date not been reported in prostate cancer models, we generated ARNTL ChIP-seq data (Supplementary Fig. S19A–S19C) to validate its binding at posttreatment FOXA1 sites. Interestingly, while we had already observed ARNTL binding to these regulatory regions in the pretreatment setting, this was strongly enhanced upon ENZ exposure (Fig. 6B; Supplementary Fig. S19D–S19F).

Functional interactions between FOXA1 and ARNTL could be further validated using ARNTL rapid immunoprecipitation mass spectrometry of endogenous proteins (RIME) experiments in ENZ-treated LNCaP-42D and LNCaP cells, confirming interactions of ARNTL not only with AR and FOXA1 but also with other classic circadian rhythm components including CLOCK/NPAS2, CRYs (CRY1, CRY2), and PERs (PER1, PER2, PER3; Fig. 6C; Supplementary Fig. S20A). As FOXA1 acts as a pioneer factor, enabling chromatin binding for other TFs including AR (9), we hypothesized that FOXA1 serves a comparable role for ARNTL. To test this hypothesis, we performed ARNTL ChIP-seq upon FOXA1 knockdown (Fig. 6D; Supplementary Fig. S20B–S20E), showing a significant decrease of ARNTL chromatin interactions exclusively for those regions co-occupied by FOXA1—highlighting FOXA1's critical role in determining ARNTL chromatin binding.

In agreement, at ARNTL consensus peaks, motifs were found to be enriched for not only CLOCK and MYC but also FOXA1 and ARNTL itself (Fig. 6E). To identify functional differences in ARNTL cistromes induced upon treatment, we overlapped the ARNTL peaks identified in all tested cell line conditions, which revealed a massive cistromic reprogramming upon ENZ treatment (Fig. 6F; Supplementary Fig. S19E and S19F). Notably, ∼70% of ENZ-gained ARNTL peaks (n = 1,752) in LNCaP cells were captured by the ARNTL cistrome in treatment-resistant cells. Interestingly, upon ENZ treatment, ARNTL binding was found to be enriched at promoter regions of key NEPC drivers, including BRN2 (POU3F2), FOXA2, EZH2, ASCL1, and SOX2 (Supplementary Fig. S21A), positioning ARNTL as a possible driver of the NE-like transcriptional program we identified. In addition, pathway overrepresentation analyses of genes coupled to Post^LNCaP–Res^LNCaP-42D shared ARNTL binding sites revealed a treatment-induced enrichment for gene sets implicated in cell-cycle progression and cell division, further supporting a possible functional involvement of ARNTL in sustaining tumor cell proliferation when AR is blocked by ENZ (Fig. 6F). To challenge this hypothesis, we assessed whether ARNTL knockdown affects the viability of hormone-sensitive and in particular of long-term ENZ-exposed cell lines. Although ARNTL targeting had minimal effect on LNCaP cell proliferation (with or without ENZ), ARNTL knockdown significantly suppressed the cell growth of ENZ-resistant LNCaP-42D cells in the absence (P = 0.031, two-way ANOVA) and even more so in the presence of ENZ (P = 7 × 10⁻⁴, two-way ANOVA), indicating that targeting ARNTL also partially restores ENZ sensitivity in this treatment-resistant cell line model (Fig. 6G). Although ARNTL was essential for sustaining cellular fitness upon ENZ treatment, exogenously introduced ARNTL did not suffice to further enhance cell proliferation when exposing LNCaP and LNCaP-42D cells to ENZ (Supplementary Fig. S21B), suggesting that ARNTL is required but not sufficient to drive the observed phenotype. Importantly, we could successfully validate the functional role of ARNTL in additional cell line models of ENZ resistance [LNCaP-ResV, originally referred to as LNCaP-ENZ^R (55), and LNCaP-ResA (56)] using ARNTL knockdown and CRISPR/Cas9-mediated ARNTL knockout (Supplementary Fig. S21C–S21E). In line with these in vitro validation experiments, ARNTL knockout also strongly inhibited the growth of LNCaP-derived ENZ-resistant xenografts (LNCaP-42D and LNCaP-ResA) in intact mice upon ENZ exposure (Fig. 6H; Supplementary Fig. S21F). Importantly, parental LNCaP cells were not affected in their proliferation potential by ARNTL knockout (Supplementary Fig. S21F), supporting the acquired dependency of ENZ-resistant cells on this circadian factor instead of a general impact on cellular fitness. Jointly, these data further highlight the treatment-induced ARNTL dependency of high-risk prostate cancer models, both in vitro and in vivo, and position ARNTL as a novel candidate therapeutic target.

GR was identified in the GIGGLE analysis as the third-most enriched factor at posttreatment FOXA1 sites, directly following ARNTL and FOXA1 itself (Fig. 5B). Given the known GR function in driving ENZ resistance in advanced CRPC (49, 57, 58), we next tested whether sustained tumor cell survival after short-term antiandrogen treatment was not only ARNTL but also GR (encoded by the NR3C1 gene) dependent. Interestingly, NR3C1 expression was upregulated upon neoadjuvant ENZ treatment in patients with primary prostate cancer (Supplementary Fig. S22A), but expression levels neither before nor after therapy were associated with clinical outcome (Supplementary Fig. S22B). Using publicly available GR ChIP-seq data from LNCaP-derived GR-positive LREX’ cells (49), we could identify GR occupancy at the majority of pretreatment FOXA1 sites and at a subset of posttreatment sites (Supplementary Fig. S22C and S22D). However, GR knockdown did not affect cellular fitness after short-term ENZ treatment in the majority of cell line models we tested (Supplementary Fig. S22E), suggesting the observed ARNTL-driven early adaptation to ENZ exposure represents a different biological entity as compared with the known GR-driven treatment resistance described in CRPC.

Overall, these data confirm the ENZ-induced FOXA1 reprogramming as observed in prostate cancer patients upon neoadjuvant antiandrogen therapy and revealed an acquired dependency on the circadian rhythm regulator ARNTL to drive tumor cell growth, positioning ARNTL as a highly promising new drug target in combination with ENZ for the treatment of high-risk prostate cancer.

In medicine, the evolutionary selection pressure imposed by drug treatment has been a well-known clinical challenge ever since the first antibiotics were discovered in the early 20th century. Also in oncology, clear escape mechanisms for both targeted therapeutics and systemic treatments are known for many years, involving ESR1 mutations in metastatic breast cancer (59), EGFR mutations in lung cancer (60), and KRAS mutations in metastatic colorectal cancer (61) but also somatic amplification of the AR locus and/or an upstream AR enhancer in CRPC (62, 63). Apart from genetic alterations, epigenetic rewiring (7, 50) or transdifferentiation are reported as mechanisms of resistance, including treatment-emergent NE prostate cancers that occur as an adaptive response under the pressure of prolonged AR-targeted therapy (64, 65).

Our unique clinical trial design with paired pre- and posttreatment biopsies of high-risk primary prostate cancer treated with ENZ monotherapy allowed us to unravel global ENZ-induced alterations in gene regulation. We report that large-scale treatment-induced dedifferentiation in prostate cancer may be a gradual process, of which the early signs are identified at the transcriptomic level within the first months of treatment onset. Although complete adenocarcinoma-to-NE transdifferentiation was not observed in any of our samples, cellular plasticity characterized by the acquisition of cistromic, transcriptomic, and proteomic features of NE disease may not only be present in primary tumors prior to treatment (31) but also become enriched upon short-term exposure to endocrine treatment, thus representing an early intermediate disease state.

In prostate cancer development (5, 32) and progression (7), AR has been reported to expose substantial plasticity in its enhancer repertoire, and we now illustrate this is also the case in primary disease upon short-term treatment. Besides AR, FOXA1 is considered a master TF and critical prostate lineage-specific regulator acting in prostate cancer, which upon overexpression during tumorigenesis gives rise to a tumor-specific AR cistrome. Also in NEPC, FOXA1 cistromes are reprogrammed (42), which indicates a direct AR-independent role of FOXA1 in prostate cancer progression. Our study confirms these observations and shows that, while co-occupied by AR, the pre- and post-ENZ enriched FOXA1 sites appeared indifferent to AR signaling.

The functional implications of the pretreatment FOXA1 sites remain unclear, as those regions were inactive, both in primary tissues and in reporter assays. A subset of these cis-regulatory elements demarcates developmental epigenomic programs that we previously reported as being occupied by FOXA1 from prostate development to tumorigenesis and metastatic progression (7), whereas others may be relevant for different physiologic processes.

The treatment-induced cistromic repositioning of pioneer factor FOXA1 initiated a thus far unknown transcriptional rewiring, in which ARNTL, a classic circadian rhythm regulator and dimerization partner of CLOCK, compensates for AR inhibition and becomes essential to rescue cellular proliferation signals. Recently, it has been reported that CRY1—a transcriptional coregulator of ARNTL—is AR regulated in prostate cancer and modulates DNA repair processes in a circadian manner (52). The current data illustrate that certain components of the circadian machinery may have a potential impact on drug response, as most clock components are not only temporally regulated at the transcriptional level but also dysregulated upon exposure to hormonal therapy. Our data now show that AR blockade forces tumor cells to adapt epigenetically, upon which these cells—over time—become dependent on ARNTL as a transcriptional regulator of proliferation processes. This acquired cellular vulnerability appears to be dependent on whether or not AR activity is inhibited and cells have had time to achieve full epigenetic reprogramming, explaining the limited effect of ARNTL knockdown in hormone-sensitive prostate cancer cells as compared with the long-term ENZ-exposed treatment-resistant models.

ARNTL expression did not correlate with outcome in patients with mCRPC. Furthermore, posttreatment-induced FOXA1 profiles showed limited overlap with NEPC–FOXA1 sites, and GR action—previously reported as a driver in CRPC—did not play a decisive role in our data sets to sustain cellular fitness following short-term ENZ exposure. Jointly, these data position the clinical state as induced by short-term neoadjuvant AR-targeted therapy in primary prostate cancer as a separate biological entity, exposing already in this early clinical stage some—but not all—features of progressive therapy-resistant disease that are invoked by drug-induced epigenetic plasticity.

With the identification of ARNTL as a rescue mechanism for tumor cells to evade AR blockade, the next question is whether ARNTL could serve as a novel therapeutic target, which should be further pursued in future drug development and clinical research. With ARNTL being critically relevant for circadian rhythm regulation, it would be imperative to balance its targeting in relation to any adverse side effects. Additionally, we demonstrate that the surprisingly dynamic enhancer repertoire of FOXA1 is critical not only in prostate tumorigenesis (5) and NE differentiation (42) but also in evading AR therapy–induced growth inhibition, further supporting the rationale to intensify efforts in targeting this highly tissue-selective yet critical transcriptional regulator directly or indirectly (66).

Primary prostate cancer tissues before and after ENZ treatment were acquired as part of the phase II, prospective, single-arm DARANA study (ClinicalTrials.gov #NCT03297385) at the Netherlands Cancer Institute Antoni van Leeuwenhoek hospital. The primary clinical outcome measure of the trial was the positive margin rate after neoadjuvant ENZ treatment. To allow sample size calculation, we performed a survey into the surgical margins of 1,492 in-house prostatectomy specimens (Gleason ≥7) not treated with antihormonal therapy prior to surgery, which revealed 34% not-radical resections. Earlier randomized studies on neoadjuvant androgen ablation showed reductions in the positive surgical margin rate of at least 50% (67–69). To detect a reduction of positive surgical margins from 34% to 17% with a power of 80% and an alpha level set at 0.05, 55 patients needed to be included. Inclusion criteria were over 18 years of age, Gleason ≥7 prostate cancer, and planned for prostatectomy. Prior to treatment, a multiparametric MRI scan was made to identify tumors in the prostate (cT stage) and pelvic lymph node metastasis (cN stage). Patients were treated with ENZ, once daily 160 mg orally without ADT, for 3 months prior to robotic-assisted laparoscopic prostatectomy and a pelvic lymph node dissection. The resection specimen was assessed for tumor margins, prostate tumor stage (ypT stage), and pelvic lymph node involvement (ypN stage). Secondary endpoints included the assessment of downstaging by comparison of preoperative clinical cT and cN stages with posttreatment and postoperative ypT and ypN stages, and differences in pre- and posttreatment prostate cancer cleaved caspase-3 and Ki-67 staining as markers of apoptosis and tumor cell proliferation, respectively. Moreover, various clinical time-to-event outcomes were included: time to BCR, defined as the time from trial inclusion to two consecutive rises of serum PSA with a minimal level of ≥0.2 ng/mL; ADT-free survival, defined as time from trial inclusion to the onset of ADT therapy; time to radiologic recurrence, defined as time from trial inclusion until detection of local or distant metastases by PSMA PET scanning; and time to distant metastases, defined as time from trial inclusion until the detection of distant metastases by PSMA PET scanning. The trial was approved by the institutional review board of the Netherlands Cancer Institute, written informed consent was signed by all participants enrolled in the study, and all research was carried out in accordance with relevant guidelines and regulations.

Prior to ENZ intervention, four preoperative MRI-guided 18-gauge core needle tumor biopsies were taken per patient. Directly after prostatectomy, eight additional tumor-targeted core needle biopsies (4 × 14-gauge, 4 × 5-mm) were taken from prostatectomy specimens ex vivo using previous MRI information and palpation. Biopsy and prostatectomy specimens were fresh-frozen (FF) or formalin-fixed, paraffin-embedded (FFPE) for ChIP-seq and CNV-seq or RNA-seq and IHC analyses, respectively. Prior to ChIP-seq experiments, FF material was cut in 30-μm sections, while FFPE material was cut in 10-μm sections prior to RNA extraction. Tissue sections were examined pathologically for tumor cell content, and only samples with a tumor cell percentage of ≥50% were used for further downstream analyses.

Chromatin immunoprecipitations on prostate cancer tissue specimens and cell line models were performed as previously described (70). In brief, cryosectioned tissue samples were double cross-linked in solution A (50 mmol/L HEPES-KOH, 100 mmol/L NaCl, 1 mmol/L EDTA, 0.5 mmol/L EGTA), first supplemented with 2 mmol/L disuccinimidyl glutarate (DSG; CovaChem) for 25 minutes at room temperature. Then, 1% formaldehyde (Merck) was added for 20 minutes and subsequently quenched with a surplus of 2.5 mol/L glycine. Cell lines were cross-linked using single-agent fixation. Therefore, 1% formaldehyde was added to the cell culture medium and incubated at room temperature for 10 minutes, followed by glycine quenching as described above. Tissue and cell line samples were lysed as described (71) and sonicated for at least 10 cycles (30 seconds on; 30 seconds off) using a PicoBioruptor (Diagenode). For each ChIP, 5 μg of antibody were conjugated to 50 μL magnetic protein A or G beads (10008D or 10009D, Thermo Fisher Scientific). The following antibodies were used: AR (06-680, Merck Millipore), FOXA1 (ab5089, Abcam), H3K27ac (39133, Active Motif), and ARNTL (ab93806, Abcam).

Immunoprecipitated DNA was processed for library preparation using a KAPA library preparation kit (KK8234, Roche), and generated libraries were sequenced on the Illumina HiSeq2500 platform using the single-end protocol with a read length of 65 bp and aligned to the human reference genome hg19 using Burrows-Wheeler Aligner (v0.5.10; ref. 72). Reads were filtered based on mapping quality (MAPQ ≥20), and duplicate reads were removed.

Peak calling over input controls (per tissue sample or cell line) was performed using MACS2 (v2.1.1) and Dfilter (v1.6) for tissues, and MACS2 (v2.1.2) for cell lines (73, 74). For tissue samples, only the peaks shared by both peak callers were used for downstream analyses. DeepTools (v2.5.3) was used to calculate read counts in peaks (FRiP; ref. 75). Read counts and the number of aligned reads, as well as normalized strand coefficient and relative strand correlation, which were calculated using phantompeaktools (v1.10.1; ref. 76), are shown in Supplementary Table S2 for tissue ChIP-seq data and Supplementary Table S6 for cell line ChIP-seq data. Tissue ChIP-seq samples that passed the following quality control (QC) measures were included in the final analyses; tumor cell percentage ≥50%, ChIP-qPCR enrichment, and more than 100 peaks called (Supplementary Fig. S2).

For the visualization of cell line ChIP-seq data, an average enrichment signal was generated by merging mapped reads of replicate samples using SAMtools (v1.10–3; ref. 77).

Genome browser snapshots, tornado plots, and average density plots were generated using EaSeq (v1.101; ref. 78). For snapshot overviews across multiple samples, bigWig files were generated from aligned bam files with the bamCoverage function from deepTools (v2.0), and snapshots were produced using pyGenomeTracks (v3.6; ref. 79) with the added NCBI RefSeq genome track (80, 81). Genomic distribution and motif enrichment analyses were performed using the CEAS and the SeqPos motif tools on Galaxy Cistrome (82), respectively. The CistromeDB Toolkit was used to probe that TFs and chromatin regulators had a significant binding overlap with the differential FOXA1 peak sets (48). For this, genomic coordinates of high-confidence binding sites (FC ≥1.2) were converted between assemblies (from hg19 to hg38), using the UCSC Genome Browser liftOver tool (83). The DiffBind R package (v2.10) was used to generate correlation heat maps and prostate cancer plots based on occupancy, to perform differential binding analyses using an FDR <0.05, and to generate consensus peak lists (84).

ChIP-seq signal of various data sets (FOXA1, AR, and H3K27ac from this study; AR, H3K27ac, and H3K27me3 from a previously reported study; ref. 31) at differential and consensus FOXA1 sites was investigated by counting mapped reads in FOXA1 peak regions using bedtools multicov (v2.27.1; ref. 85). Read counts were subsequently z-transformed and visualized using the aheatmap function from the R package NMF (v0.21.0; ref. 86) with a color scheme from RColorBrewer (v1.1-2; https://CRAN.R-project.org/package=RColorBrewer). To determine the significance in binding site occupancy differences between pre- and posttreatment FOXA1 sites, median z-transformed read counts were calculated per sample and compared using a Mann–Whitney U test. These median read counts per sample were also used to assess the correlation between ChIP-seq signals of AR, FOXA1, and H3K27ac at pre-enriched, post-enriched, and consensus FOXA1 binding sites.

Bedtools intersect (v2.27.1; ref. 85) was used to determine the overlap of differential FOXA1 binding sites and inactive, constitutively active, and inducible ARBS.

To assign FOXA1 and ARNTL binding regions to potential target genes, we overlapped differential FOXA1 binding sites with H3K27ac HiChIP data (44) using bedtools intersect. To assess whether or not genes coupled to FOXA1 binding sites were considered to be essential for the VCaP prostate cancer cell line, we used the DepMap (Broad 2020) 20Q1 Public gene effect data set (45) with a stringent gene effect score cutoff ≤−1. Gene set overlaps between genes linked to ChIP-seq binding sites and the Molecular Signatures Database (v7.4) were computed using GSEA (87) with an FDR q-value cutoff ≤0.05.

Prior to RNA isolation, FFPE material was pathologically assessed. The expert pathologist scored tumor cell percentage and indicated most tumor-dense regions for isolation on a hematoxylin and eosin (H&E) slide. RNA and DNA from FFPE material were simultaneously isolated from 3 to 10 sections (depending on tumor size) of 10 μm using the AllPrep DNA/RNA FFPE isolation kit (80234, Qiagen) and the QIAcube according to the manufacturer’s instructions. cDNA was synthesized from 250 ng RNA using SuperScript III Reverse Transcriptase (Invitrogen) with random hexamer primers.

Strand-specific libraries were generated with the TruSeq RNA Exome kit (Illumina) and sequenced on the Illumina HiSeq2500 platform using the single-end protocol with a read length of 65 bp. Sequencing data were aligned to the human reference genome hg38 using HISAT2 (v2.1.0; ref. 88), and the number of reads per gene was measured with HTSeq count (v0.5.3; ref. 89).

For QC purposes, total read counts per sample were determined and hierarchical clustering based on the Euclidean distance was applied. Samples with a read count ≥2 standard deviations below the mean of all sample read counts were removed, as well as samples that clustered in a separate branch.

Global gene expression differences between pre- and posttreatment samples passing QC were determined using DESeq2 (v1.22.2; ref. 90). The significance of expression level differences between pre- and posttreatment samples was determined using a paired t test.

Gene set enrichment was performed using preranked GSEA (87) based on the Wald statistic provided by DESeq2. For visualization purposes, the data were z-transformed per gene. Heat maps of gene expression values were created using the aheatmap function from the R package NMF (v0.21.0; ref. 86) with a color scheme from RColorBrewer (v1.1-2; https://CRAN.R-project.org/package=RColorBrewer).

To assign samples to previously described prostate cancer subtypes (31), the z-transformed expression levels of the top ∼100 most differentially expressed genes (n = 285) in each of the three clusters were investigated. Using these values, samples were clustered based on their Pearson correlation. The resulting tree was divided into three clusters corresponding to the previously published prostate cancer subtypes. The potential transitioning of samples from one cluster to another after treatment was visualized using a riverplot (v0.6; https://CRAN.R-project.org/package=riverplot).

To calculate FCs of NE scores upon treatment, expression of 70 NE signature genes was obtained from castration-resistant NE and prostate adenocarcinoma samples as published previously (35). The expression of 5 of the 70 NE signature genes was not included in the analysis (KIAA0408, SOGA3, LRRC16B, ST8SIA3, and SVOP) because the genes are not expressed in these samples. Expression FCs between paired pre- and posttreatment samples were calculated (n = 39), and concordance in gene expression differences (FC sign) was measured using Pearson correlation.

Low-coverage whole-genome sequencing of ChIP-seq input samples was performed on a HiSeq 2500 system (single end, 65-bp), and samples were aligned to hg19 with Burrows-Wheeler Aligner (BWA) backtrack algorithm (v0.5.10; ref. 72). Per sample, the mappability of all reads with a phred quality score of 37 and higher per 20 kb on the genome was rated against a similarly obtained mappability for all known and tiled 65-bp subsections of hg19. Sample counts were corrected per bin for local guanine–cytosine (GC) content effects using a nonlinear loess fit of mappabilities over 0.8 on autosomes. Reference values were scaled according to the slope of a linear fit, forced to intercept at the origin, of reference mappabilities after GC correction. Ratios of corrected sample counts and reference values left out bins with mappability below 0.2 or overlapping ENCODE blacklisted regions (91).

Copy-number log ratios were smoothed and segmented using the R package DNACopy (v1.50.1; https://bioconductor.org/packages/release/bioc/html/DNAcopy.html) with the parameters set to alpha = 0.00000000001, undo.SD = 2, and undo.splits = “sdundo.” Bedtools intersect (v2.27.1; ref. 85) was used to determine overlap between copy-number segments and differential FOXA1 binding sites. These data were subsequently visualized using the aheatmap function from the R package NMF (v0.21.0; ref. 86) with a color scheme from RColorBrewer (v1.1-2; https://CRAN.R-project.org/package=RColorBrewer).

To correlate FOXA1 ChIP-seq signal with copy-number status at differential FOXA1 sites, we used the z-transformed FOXA1 ChIP-seq read counts as described in the ChIP-seq section. The difference in transformed ChIP-seq read counts and the difference in normalized segmented copy-number data between matched posttreatment and pretreatment samples were calculated for every patient. Subsequently, the Pearson correlation between these two sets of differences was calculated.

For IHC analysis, we matched our ENZ-treated patient cohort (n = 51) in a 1:2 ratio to untreated control patients (not receiving ENZ prior to prostatectomy; n = 110) based on clinicopathologic parameters [initial PSA, Gleason score, tumor–node–metastasis (TNM) stage, age] using the R package MatchIt (v.4.1.0; ref. 92).

TMAs were prepared containing three cores per FFPE tumor sample. Tumor-dense areas in FFPE megablocks were marked by an expert pathologist on an H&E slide. Cores were drilled in a receptor block using the TMA grandmaster (3D Histech/Sysmex). Next, cores were taken from the donor block and placed in the receptor block using the manual tissue arrayer (4508-DM, Beecher instruments). The filled receptor block was placed in a 70°C stove for 9 minutes and cooled overnight at room temperature.

IHC was applied to TMA slides using a BenchMark Ultra autostainer (Ventana Medical Systems). In brief, paraffin sections were cut at 3 μm, heated at 75°C for 28 minutes, and deparaffinized in the instrument with EZ prep solution (Ventana Medical Systems). Heat-induced antigen retrieval was carried out using Cell Conditioning 1 (CC1, Ventana Medical Systems) for 24 minutes (cleaved caspase-3), 32 minutes (chromogranin, synaptophysin), or 64 minutes (ARNTL, Ki-67) at 95°C. The following antibodies and staining conditions were used: anti-ARNTL (ab230822, Abcam; 1:1,000 dilution; 60 minutes at room temperature), anti-Ki-67 (M7240, Agilent; 1;100 dilution; 60 minutes at 37°C), anti–cleaved caspase-3 (9661, Cell Signaling Technology; 1:100 dilution; 32 minutes at room temperature), anti-chromogranin (760–2519, Ventana Medical Systems; undiluted; 32 minutes at 37°C), and anti-synaptophysin (SYNAP-299-L-CE, Leica; 1:100; 32 minutes at 37°C). For synaptophysin, signal amplification was applied using the OptiView Amplification Kit (Ventana Medical Systems; 4 minutes). Bound antibody was detected using the OptiView DAB Detection Kit (Ventana Medical Systems). Slides were counterstained with hematoxylin and bluing reagent (Ventana Medical Systems).

Percentage of positive tumor cells or IHC staining intensity (weak, moderate, and strong) in tumor cells was scored by an expert pathologist and used for statistical analysis.

FOXA1 mutation status was assessed from H3K27ac ChIP-seq and RNA-seq reads covering the gene. We focused our search on genomic coordinates with mutations previously reported in cBioPortal (https://www.cbioportal.org). cBioPortal was queried for all somatic mutations in FOXA1 (n = 567 mutations) across all prostate cancer samples (n = 6,875 patients). Nonreference alleles were called from bam files with H3K27ac ChIP-seq or RNA-seq reads using the mpileup and call commands from bcftools (v1.9). The –prior variable for call was set to 0.05 to enhance sensitivity in the setting of low read coverage. The genomic coordinates of variants were listed in bed files and tested for overlap with FOXA1 mutations from cBioPortal.

RNA-seq data from mCRPC were processed as previously described (53, 54, 93) and converted to transcripts per million or fragments per kilobase per million reads mapped (FPKM). Patients were grouped by ARNTL expression levels as low (< median) or high (≥ median). Survival analysis was performed using the Kaplan–Meier method with endpoint overall survival from diagnosis of mCRPC, and the Wald test was used to test for statistical significance.

The LNCaP human prostate cancer cell line and HEK293T cells were purchased from the American Type Culture Collection. ENZ -resistant LNCaP-42D (41) and LNCaP-ResA (56) cells were described previously. LNCaP clones were maintained in RPMI 1640 medium (Gibco, Thermo Fisher Scientific) supplemented with 10% FBS (Sigma-Aldrich), with ENZ-resistant cell lines further supplemented with 10 μmol/L ENZ (MedChemExpress). HEK293T cells were cultured in DMEM (Gibco, Thermo Fisher Scientific) supplemented with 10% FBS. Cell lines were subjected to regular Mycoplasma testing, and all cell lines underwent authentication by short tandem repeat profiling (Eurofins Genomics). For hormone stimulation with synthetic androgen, cells were treated with 10 nmol/L R1881 (PerkinElmer) for 48 hours. For in vitro AR blockade, cells were treated with 10 μmol/L ENZ and harvested at the indicated time points.

Pooled human male genomic DNA (Promega) was randomly sheared, end-repaired, and ligated with Illumina compatible xGen CS stubby adapters (IDT) containing 3-bp unique molecular identifiers. The adapter-ligated gDNA fragments (500–800 bp) were hybridized to a custom biotinylated oligonucleotide probe (Agilent) and captured by Dynabeads M-270 Streptavidin beads (NEB). The library was designed to capture regions from clinical ChIP-seq. Any overlapping reads were collapsed using the BedTools “merge” (v2.30.0) command to eliminate possible overrepresentation. Target regions were PCR-amplified with STARR_in-fusion_F (TAGAGCATGCACCGGACACTCTTTCCCTACACGACGCTCTTCCGATCT) and STARR_in-fusion_R (GGCCGAATTCGTCGAGTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT) primers, and cloned into AgeI-HF (NEB) and SalI-HF (NEB) digested hSTARR-ORI plasmid (#99296, Addgene) by Gibson Assembly. The STARR-seq capture library was transformed into MegaX DH10B T1R electrocompetent cells (Invitrogen). Plasmid DNA was extracted using the Qiagen Plasmid Maxi Kit.

LNCaP cells (>1.6 × 10⁸ cells/replicate; three biological replicates for each cell line) were electroporated with the STARR-seq capture library (1 × 10⁶ cells: 2 μg DNA; ∼320 μg plasmid DNA/replicate) using the Neon Transfection System (Invitrogen). Electroporated LNCaP cells were immediately recovered in RPMI 1640 medium supplemented with 10% FBS, and the culture medium was refreshed 24 hours after electroporation. LNCaP cells (∼0.5 × 10⁸ cells) were treated with dimethylsulfoxide (DMSO) or 10 μmol/L ENZ for 48 hours and then either EtOH or 10 nmol/L DHT for 4 hours. All electroporated LNCaP cells were harvested 72 hours after electroporation. Cell samples were lysed with the Precellys CKMix Tissue Homegenizing Kit and Precellys 24 Tissue/Cell Ruptor (Berin Technologies). Total RNA was extracted using the Qiagen RNeasy Maxi Kit (Qiagen), and poly-A mRNA was isolated using the Oligo (dT)₂₅ Dynabeads (Thermo Fisher). FOXA1-focused STARR-seq cDNA was synthesized with the gene-specific primer (CTCATCAATGTATCTTATCATGTCTG) and amplified by junction PCR (15 cycles) with the RNA_jPCR_f (TCGTGAGGCACTGGGCAG*G*T*G*T*C) and jPCR-r (CTTATCATGTCTGCTCGA*A*G*C) primers. FOXA1-focused STARR-seq capture library plasmid DNA was extracted from 0.1 × 10⁸ transfected but untreated cells using the QIAprep Spin Miniprep Kit (Qiagen). The extracted plasmid DNA and the input plasmid DNA were PCR-amplified with the DNA-specific junction PCR primer (DNA_jPCR_f, CCTTTCTCTCCACAGGT*G*T*C) and the jPCR-r primer. After purification with Ampure XP beads, Illumina compatible libraries were generated by PCR amplification with NEBNext universal and single indexing primers (NEB) and were sequenced on an Illumina NovaSeq6000 (150-bp, paired-end).

STARR-seq data were analyzed using a custom Snakemake pipeline (https://github.com/birkiy/starr-pipe). Briefly, paired-end STARR-seq samples were aligned to the hg19 genome using BWA (v0.7.17). Raw alignment files were converted into BEDPE format using the BedTools (v2.30.0) “bamtobed -bedpe” command. The start of the first paired read and the end of its mate defined the fragments from the BEDPE file. Any fragments overlapping with hg19 blacklisted regions (https://github.com/Boyle-Lab/Blacklist) or MAPQ scores <30 were filtered from downstream analysis. Fragments containing unique genomic positions were counted using the “uniq -c” UNIX command. A count table of the unique fragment collection count was generated using a custom Julia script (v1.5.2) fragments.jl, which uses library input samples to first generate the reference fragment population and then quantifies the frequencies of each fragment.

STARR-seq aligned files were downsampled using SAMtools (v1.10) to make files with equivalent read counts across conditions. Next, count tables were generated from the downsampled files for all tested FOXA1 regions (n = 1,209) using the deepTools (v2.0) multiBamSummary function. The most correlated replicates were chosen for further analysis using the cor function in R (v3.4.4). Regions with zero counts across samples were removed, leaving 968 regions. These count tables were used as input for a differential expression analysis using DESeq2 (v1.22.2) in R. Regions with nonsignificant changes (FDR ≤ 0.05, logFC ≥ |2|) in read counts upon ENZ treatment were identified, and k-means clustering from the plotHeatmap function of deepTools was performed. To determine possible functional associations within these clusters, the sets of regions were queried using the CistromeDB Toolkit to identify factors with significant overlap.

Following treatment of LNCaP and LNCaP-42D cells with ENZ (10 μmol/L) for 48 hours, cells were fixed, lysed, and sonicated as previously described (94). The nuclear lysates were incubated with 50 μL magnetic protein A beads (10008D, Thermo Fisher Scientific) conjugated to 7.5 μg of ARNTL antibody (ab93806, Abcam) or rabbit IgG control (12-370, Merck Millipore).

Peptide mixtures were prepared and measured as previously described (4) with the following noted exceptions. Peptide mixtures (10% of total digest) were loaded directly onto the analytical column (ReproSil-Pur 120 C18-AQ, 2.4 μm, 75 μm × 500 mm, packed in-house) and analyzed by nano LC-MS/MS on an Orbitrap Fusion Tribrid mass spectrometer equipped with a Proxeon nLC1200 system (Thermo Scientific). Solvent A was 0.1% formic acid/water, and solvent B was 0.1% formic acid/80% acetonitrile. Peptides were eluted from the analytical column at a constant flow of 250 nL/minute in a 120-minute gradient containing a 105-minute step-wise increase from 7% to 34% solvent B, followed by a 15-minute wash at 80% solvent B.

Raw data were analyzed by MaxQuant (v2.0.1.0; ref. 95) using standard settings for label-free quantitation (LFQ). MS/MS data were searched against the Swissprot human database (20,397 entries, release 2021_01) complemented with a list of common contaminants and concatenated with the reversed version of all sequences. The maximum allowed mass tolerance was 4.5 ppm in the main search and 0.5 Da for fragment ion masses. FDRs for peptide and protein identification were set to 1%. Trypsin/P was chosen as cleavage specificity allowing two missed cleavages. Carbamidomethylation was set as a fixed modification, while oxidation and deamidation were used as variable modifications. LFQ intensities were log₂-transformed in Perseus (v1.6.15.0; ref. 96), after which proteins were filtered for at least three of four valid values in at least one sample group. Missing values were replaced by imputation based on a normal distribution (width: 0.3; downshift: 1.8). Differentially enriched proteins were determined using a Student t test (threshold: P ≤ 0.05 and [x − y] ≥ 1.8 | [x − y] ≤ −1.8).

Transient transfections of cell lines were performed according to the manufacturer’s instructions using Lipofectamine 2000 (Invitrogen) or Lipofectamine RNAiMAX (Invitrogen) for overexpression or siRNA knockdown experiments, respectively. ARNTL containing expression plasmid was obtained from the CCSB-Broad Lentiviral Expression Library. siRNA oligos targeting ARNTL (M-010261-00-0005), FOXA1 (M-010319-01-0020), and NR3C1 (M-003424-03-0005) and the nontargeting control (D-001206-14, D-001210-05-20) were purchased from Dharmacon. For ARNTL ChIP-seq upon FOXA1 knockdown, LNCaP and LNCaP-42D cells were reverse transfected with 50 nmol/L siFOXA1 using Lipofectamine RNAiMAX. ENZ (10 μmol/L) was added after 24 hours, and cells were fixed and harvested for ChIP-seq analysis 72 hours after transfection.

Guide RNAs targeting human ARNTL (CTGGACATTGCGTTGCATGT) and a nontargeting control guide (AACTACAAGTAAAAGTATCG) were individually cloned into the lentiCRISPR v2 plasmid (97). CRISPR vectors were coexpressed with third-generation viral vectors in HEK293T cells using polyethyleneimine (PEI; Polysciences). After lentivirus production, the medium was harvested and transferred to the designated cell lines. Two days after infection, cells were put on puromycin (Sigma-Aldrich) selection for 3 weeks, and knockout efficiency was tested using western blot analysis.

Total proteins were extracted from cells using Laemmli lysis buffer supplemented with a complete protease inhibitor cocktail (Roche). Per sample, 40 μg of protein was resolved by SDS-PAGE (10%) and transferred on nitrocellulose membranes (Santa Cruz Biotechnology). The following antibodies were used for western blot stainings: ARNTL (ab93806, Abcam), PSA (5365, Cell Signaling Technology), FOXA1 (ab5089, Abcam), GR (12041, Cell Signaling Technology), and ACTIN (MAB1501R, Merck Millipore). Blots were incubated overnight at 4°C with designated primary antibodies at 1:1,000 (ARNTL, PSA FOXA1, GR) or 1:5,000 (ACTIN) dilution and visualized using the Odyssey system (Li-Cor Biosciences).

Total RNA from cell lines was isolated using TRIzol Reagent (Thermo Fisher Scientific), and cDNA was synthesized from 2 μg RNA using the SuperScript III Reverse Transcriptase system (Thermo Fisher Scientific) with random hexamer primers according to the manufacturer’s instructions. Quantitative PCR (qPCR) was performed using the SensiMix SYBR Kit (Bioline) according to the instructions provided by the manufacturer on a Quant­Studio 6 Flex System (Thermo Fisher Scientific). Primer sequences for mRNA expression analyses were FOXA1 (forward: CGACTGGAACAGCTACTACG; reverse: TGGTGTTCATGGTCATGTAGGT) and ARNTL (forward: CTGGAGCACGACGTTCTTTCTT; reverse: GGATTGTGCAGAAGCTTTTTCG). mRNA levels are shown relative to the expression of housekeeping gene TBP (forward: GTTCTGGGAAAATGGTGTGC; reverse: GCTGGAAAACCCAACTTCTG).

For cell viability assays, LNCaP, LNCaP-42D, or LNCaP-ResA cells were seeded at 2 × 10³ cells per well in 96-well plates (Greiner) ± 10 μmol/L ENZ, and reverse transfected with 50 to 100 nmol/L siRNA (Dharmacon) using Lipofectamine RNAiMAX (Invitrogen). Cell viability was assessed 7 days after transfection using the CellTiter-Glo Luminescent Cell Viability Assay kit (Promega) per the manufacturer’s instructions. Bar charts were plotted using GraphPad Prism 9 software.

Proliferation curves for stable ARNTL knockout clones were generated using a Lionheart FX automated microscope (BioTek). Cells (LNCaP, LNCaP-42D, and LNCaP-ResA) were seeded at 2 × 10³ cells per well in 96-well plates ± 10 μmol/L ENZ. SiR-DNA (Spirochrome) live-cell nuclear stain was added 2 hours before imaging. Cell growth was recorded with a time resolution of 4 hours for a total time span of 144 hours. The microscope was maintained at 37°C in 5% CO₂, and live-cell imaging was performed using a 4× lens and a Sony CCD, 1,25-megapixel camera with two times binning (BioTek). Gen5 software (BioTek) was used to quantify cell numbers, and growth curves were plotted using GraphPad Prism 9 software.

For in vivo tumor growth studies, 7.5 × 10⁶ sgNT or sgARNTL (LNCaP, LNCaP-42D, and LNCaP-ResA) cells in PBS with 50% BME (3536-005-02, Bio-Techne) were injected subcutaneously into one of the flanks of ∼7-week-old male NOD-SCID (NSG) mice. Once tumor size reached 150 mm³, mice were randomized and treated with either 10 mg/kg ENZ (MedChemExpress) or vehicle control (1% carboxymethylcellulose sodium salt, 0.1% Tween-80, 5% DMSO; Sigma-Aldrich) through oral gavage on a daily basis. Tumor volume was monitored by caliper measurements 3 times a week. Mice were kept under standard temperature and humidity conditions in individually ventilated cages, with food and water provided ad libitum. All animal experiments were approved by the Animal Welfare Committee of the Netherlands Cancer Institute and performed in accordance with institutional, national, and European guidelines for animal research.

For differential binding and differential gene expression analyses (pre- vs. post-ENZ), an FDR cutoff <0.05 (P < 0.01) and P_adj < 0.01 was used, respectively. A Mann–Whitney U test was used to determine differences in region read counts (adjusted for multiple testing using FDR) and differences in gene expression levels before and after ENZ treatment. For peak set and gene set overlaps as well as to determine differences in IHC staining intensities, Fisher exact tests were applied. Differences in cell viability or cell/tumor growth were tested using a two-way or one-way ANOVA followed by Tukey multiple comparisons test, respectively (GraphPad Prism 9). Corresponding bar chart or growth curves show the mean with error bars representing the SD. All box plots indicate the median (center line), upper (75) and lower (25) quartile range (box limits), and 1.5 × interquartile range (whiskers). Significance is indicated as follows: ns, P > 0.05; *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001. Further details of statistical tests are provided in the figure legends.

All tissue ChIP-seq and RNA-seq raw data generated in this study have been deposited in the European Genome-phenome Archive (EGA) under the accession numbers EGAS00001006017 and EGAS00001006016, respectively. The cell line ChIP-seq, as well as all processed tissue ChIP-seq and RNA-seq data, have been deposited in the Gene Expression Omnibus (GEO) database under accession number GSE197781. The mass spectrometry proteomics (RIME) data have been deposited to the ProteomeXchange Consortium via the PRIDE partner repository with the data set identifier PXD032041. Public ChIP-seq data sets used in this study are available from GEO or EGA under the following accession numbers: GSE120738 (AR, H3K27ac, H3K27me3 ChIP-seq), GSE51497 (GR ChIP-seq), GSE117306 (N-MYC ChIP-seq), and EGAS00001003928 (FOXA1 ChIP-seq).

S.C. Baca reports other support from Precede Biosciences outside the submitted work. F.Y. Feng reports personal fees from Janssen Oncology, Bayer, Exact Sciences, Blue Earth Diagnostics, Myovant Sciences, Roivant Sciences, Astellas Pharma, Foundation Medicine, Varian, Bristol Meyers Squibb, Bluestar Genomics, Novartis, and Tempus and other support from PFS Genomics, SerImmune, and Artera AI outside the submitted work. M.L. Freedman reports other support from Precede Bio outside the submitted work. L.F.A. Wessels reports grants from Genmab outside the submitted work. N.A. Lack reports grants, personal fees, and other support from Nido Biosciences and grants from AstraZeneca outside the submitted work. H. van der Poel reports grants from Astellas during the conduct of the study. A.M. Bergman reports grants from Movember, the Dutch Cancer Society, and Astellas Pharma during the conduct of the study, as well as personal fees from Astellas Pharma outside the submitted work. W. Zwart reports grants from Astellas Pharma during the conduct of the study, as well as grants and personal fees from Astellas Pharma, and grants from AstraZeneca outside the submitted work. No disclosures were reported by the other authors.

S. Linder: Conceptualization, data curation, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. M. Hoogstraat: Conceptualization, data curation, software, formal analysis, investigation, visualization, methodology, writing–original draft, writing–review and editing. S. Stelloo: Conceptualization, funding acquisition, investigation, methodology, writing–review and editing. N. Eickhoff: Resources, formal analysis, validation, investigation, methodology, writing–review and editing. K. Schuurman: Formal analysis, investigation, methodology, writing–review and editing. H. de Barros: Resources, data curation, formal analysis, investigation, writing–review and editing. M. Alkemade: Data curation, investigation, methodology, writing–review and editing. E.M. Bekers: Formal analysis, investigation, methodology, writing–review and editing. T.M. Severson: Data curation, software, formal analysis, investigation, visualization, methodology, writing–review and editing. J. Sanders: Formal analysis, investigation, methodology, writing–review and editing. C.-C.F. Huang: Investigation, methodology, writing–review and editing. T. Morova: Software, formal analysis, methodology, writing–review and editing. U.B. Altintas: Software, formal analysis, methodology, writing–review and editing. L. Hoekman: Data curation, formal analysis, investigation, methodology, writing–review and editing. Y. Kim: Resources, software, formal analysis, investigation, visualization, methodology, writing–review and editing. S.C. Baca: Resources, data curation, software, formal analysis, investigation, methodology, writing–review and editing. M. Sjostrom: Resources, data curation, software, formal analysis, investigation, visualization, methodology, writing–review and editing. A. Zaalberg: Investigation, methodology, writing–review and editing. D.C. Hintzen: Data curation, software, investigation, methodology, writing–review and editing. J. de Jong: Formal analysis, investigation, methodology, writing–review and editing. R.J.C. Kluin: Software, investigation, methodology, writing–review and editing. I. de Rink: Data curation, software, formal analysis, methodology, writing–review and editing. C. Giambartolomei: Resources. J.-H. Seo: Resources. B. Pasaniuc: Resources. M. Altelaar: Software, supervision, methodology, writing–review and editing. R.H. Medema: Software, supervision, methodology, writing–review and editing. F.Y. Feng: Resources, supervision, methodology, writing–review and editing. A. Zoubeidi: Resources, writing–review and editing. M.L. Freedman: Resources, formal analysis, investigation, writing–review and editing. L.F.A. Wessels: Supervision, writing–review and editing. L.M. Butler: Writing–review and editing. N.A. Lack: Resources, supervision, investigation, methodology, writing–review and editing. H. van der Poel: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, methodology, project administration, writing–review and editing. A.M. Bergman: Conceptualization, resources, supervision, funding acquisition, investigation, writing–original draft, project administration, writing–review and editing. W. Zwart: Conceptualization, resources, supervision, funding acquisition, investigation, methodology, writing–original draft, project administration, writing–review and editing.

This work was supported by Movember (NKI01 to A.M. Bergman and W. Zwart), KWF Dutch Cancer Society (10084 ALPE to A.M. Bergman and W. Zwart), KWF Dutch Cancer Society/Alpe d’HuZes Bas Mulder Award (NKI 2014-6711 to W. Zwart), the Netherlands Organization for Scientific Research (NWO-VIDI-016.156.401 to W. Zwart), Department of Defense (W81XWH-19-1-0565 to W. Zwart), and Astellas Pharma (to W. Zwart, A.M. Bergman, and H. van der Poel). The authors thank the NKI-AVL Core Facility Molecular Pathology and Biobanking, the NKI Proteomics/Mass Spectrometry facility (supported by the Dutch NWO X-omics Initiative), and the NKI Preclinical Intervention Unit (MCCA) for experimental assistance; the NKI Genomics Core facility for next-generation sequencing and bioinformatics support; and the NKI Research High Performance Computing facility for computational infrastructure. We also thank all Zwart/Bergman lab members for fruitful discussions and technical advice. Finally, we thank all patients and clinical staff who were involved in the DARANA trial.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.