Introduction

Prostate cancer (PC) is a hormonally-driven disease in which androgen receptor (AR) signaling plays a fundamental role in defining a luminal epithelial cancer lineage and promoting cell survival and proliferation. Thus, androgen deprivation therapy (ADT) and AR signaling inhibitors (ARSIs) represent standard and effective frontline treatments for advanced PC. While the majority of PCs initially respond to combined treatment with ADT and an ARSI, therapeutic resistance ultimately develops over the course of months to years and leads to a disease state known as castration-resistant PC (CRPC) ^1,2. The mechanisms by which CRPC subvert ADT and ARSIs are numerous, but most are centered on aberrant reactivation of AR signaling through AR gene mutations, amplification, alternative splicing, and enhancer duplication ^2–5. Beyond reactivation of the targeted AR pathway, clinical evidence implicates lineage plasticity or a shift in cell identity as a more extreme adaptive mechanism of resistance in PC and other epithelial malignancies. A prime example is treatment-related transdifferentiation of AR-active PC (ARPC) to neuroendocrine PC (NEPC) in which the AR program is extinguished and supplanted with a NE program (NEtD) ^6,7, thereby leading to a cancer lineage state that no longer depends on AR signaling and is resistant to ADT and ARSIs.

Multiple genetic and molecular events have been associated with the process of NEtD, including amplification of MYCN, loss of RE1-silencing transcription factor (REST) activity, and loss of PTEN, TP53, and RB1 ^8–13. While our prior studies using human prostate epithelial transformation models have demonstrated that overexpression of MYCN and constitutively active myristoylated AKT1 ⁸ or the PARCB factors (dominant-negative TP53, myristoylated AKT1, short hairpin targeting RB1, MYC, and BCL2) can initiate NEPC ⁹, these experimental systems start with AR-null basal prostate epithelial cells. Thus, they do not address how AR-enforced luminal epithelial lineage commitment may be repressed and bypassed which is a critical step in the NEtD of PC. The leading biological model currently available to investigate NEtD is the LTL-331 ARPC patient-derived xenograft (PDX) model which reproducibly undergoes transdifferentiation after passage in castrated mice to the LTL-331R NEPC tumor model ¹⁴. Repeated biopsies and molecular analyses of tumors from this system have provided insights into molecular programs that are associated with NEtD ¹⁵. Yet, this model remains experimentally unwieldy for functional studies, and other more tractable approaches have been lacking.

Another challenge has been the inconsistent definition of NEPC applied in the field. Clinically, NEPC consists of a heterogeneous group of neuroendocrine tumors primarily defined by histologic morphology that includes a range of subtypes from well-differentiated carcinoid to aggressive large-cell and small-cell carcinomas ¹⁶. The expression of NE markers like chromogranin A (CHGA), synaptophysin (SYP), or CD56 (neural cell adhesion molecule 1 or NCAM1) are often assayed by immunohistochemistry (IHC) to support or confirm a diagnosis of NEPC. In general, individual NE markers expression can be heterogeneous, and combined marker scores are more robust ¹⁷. In experimental studies, phenotypic features such as NE marker expression or the emergence of dendrite-like projections from PC cells are broadly used as indicators of NE transdifferentiation of PC. The use of these sparse phenotypic features to define NEPC may be problematic as they may not reflect a bona fide cancer lineage state divergent from ARPC. A prime example is amphicrine prostate cancer (AMPC) which co-expresses AR and NE markers ¹⁸. Despite the expression of NE markers, AMPC models such as VCaP remain sensitive to the ARSI enzalutamide, indicating a functional dependence on the AR-driven lineage program.

Lineage-defining transcription factors and transcription factor networks specify cell states by regulating chromatin accessibility and gene expression programs. The basic helix-loop-helix (bHLH) proteins Achaete-scute homolog 1 (ASCL1) and neurogenic differentiation 1 (NeuroD1) are pioneer transcriptional factors that are necessary for lineage specification during neuronal development ^19,20. In cancer, ASCL1 and NeuroD1 define distinct subtypes of small cell lung cancer (SCLC) ²¹ and NEPC ²² with differential transcriptional, epigenetic, cell surface proteomic, and drug sensitivity profiles ^23–25. The role of ASCL1 as a driver of neuroendocrine transformation has been supported by the genetically engineered TCKO mouse model of high-grade neuroendocrine lung cancer where Ascl1 knockout prevents tumor formation ²¹. Further, knockdown of ASCL1 or NEUROD1 in established human SCLC or NEPC cell line models significantly impairs cell viability ^26,27, indicating that these factors are needed for sustained survival and growth. Insights into the functional role of these pioneer transcription factors in mediating the acute lineage conversion from ARPC to NEPC are limited. Recent work has shown that ASCL1 expression may be induced by ARSI therapy in PC and ASCL1 activates a neuronal stem cell-like lineage program through chromatin remodeling mediated by the polycomb repressive complex 2 (PRC2) ²⁸.

In this study, we developed a genetically defined strategy adapted from the field of cellular reprogramming to directly convert ARPC to NEPC. We identified key functional roles for ASCL1 and NeuroD1 in NE transdifferentiation, specifically new findings uncovering their ability to silence the AR program through chromatin remodeling of the AR enhancer which is somatically amplified in PC and global AR binding sites with enhancer activity. In contrast, we found that loss of REST activity induces the expression of NE markers but is insufficient to drive a NE cancer lineage state. Repeated sampling over the reprogramming time course revealed a wide range of temporal changes in the transcriptomic and epigenomic landscapes associated with the lineage conversion from ARPC to NEPC. Lastly, we found that ASCL1 and NeuroD1 expression, but not the loss of REST activity, induces the downregulation of class I major histocompatibility complex (MHC) expression in PC, providing a basis for the lineage-specific potential for immune evasion associated with NEPC.

Results

Candidate factors directly reprogram ARPC to NEPC and bypass the dependence on AR signaling

To functionally define factors required for the conversion of ARPC to NEPC, we established a NEtD assay wherein candidate factors were introduced into the human ARPC cell lines LNCaP and C4-2B after which they were propagated in media conditions permissive to the growth of NEPC and analyzed for a set of phenotypic markers of ARPC and NEPC after two weeks (Figure 1A). Initially, we introduced a pool of lentiviruses (LVs) expressing eight factors associated with NEPC including dominant-negative TP53 R175H, a short hairpin targeting RB1 (shRB1), MYCN, ASCL1, SRRM4, NR0B2, BCL2, and KRAS G12V. Loss of TP53 and RB1 as well as amplification of MYCN have been identified as common genetic features associated with NEPC ^8,9,29,30. SRRM4 is a splicing regulator that has been shown to inactivate the RE-1 silencing transcription factor (REST) which functions to suppress a major neural gene program. Loss of REST has been implicated as a common feature in NEPC that has been reported to functionally drive NEtD of PC ^12,31. In addition, the pioneer neural transcription factor ASCL1 is expressed in most NEPC and appears to define a transcriptional and epigenetic subtype of the disease ²². NR0B2, otherwise known as short heterodimer partner (SHP), is an orphan nuclear receptor that is highly expressed in NEPC that interacts with AR and inhibits AR activity ³². Lastly, BCL2 and KRAS G12V were included to potentially reduce apoptosis and enhance reprogramming efficiencies ^33,34.

Figure 1.

The introduction of this candidate factor LV pool at a multiplicity-of-infection (MOI) of four for each factor in this assay led to a significant downregulation of AR and the AR target NK3 homeobox 1 (NKX3-1) as well as the upregulation of NE markers including forkhead box A2 (FOXA2), SYP, NCAM1, insulinoma-associated protein 1 (INSM1), and POU class 3 homeobox 2 (POU3F2 or BRN2) by immunoblot analysis (Figure 1B). Furthermore, immunocytochemical analysis of the C4-2B and LNCaP cell lines transduced with the candidate factor LV pool demonstrated a prominent loss of nuclear AR expression and increased SYP expression (Figure 1C). The AR^low/-NE⁺ cancer phenotype induced by the candidate factor LV pool was durable based on IHC studies in tumors obtained after subcutaneous xenografting in non-castrate male NOD.Cg-Prkdc^scid II2rg^tm1Wjl SzJ (NSG) mice (Figure S1A).

We next sought to determine whether the PC cells reprogrammed with the candidate factor LV pool were functionally independent of AR signaling. We implemented a potent AR-dependent negative selection strategy (Figure 1D) wherein C4-2B cells were transduced with LV encoding AR response elements (ARE) upstream of a FK506 binding protein-caspase 8 (FKBP-Casp8) fusion. A clonal C4-2B ARE-FKBP-Casp8 line was established and transduced with either a control LV expressing green fluorescent protein (GFP) or the candidate factor LV pool and subjected to selection with 1 nM R1881 (to induce AR signaling and expression of FKBP-Casp8 in AR⁺ cells) and AP20187 which is a synthetic FK506 analog (to activate FKBP-Casp8 and induce apoptosis). We identified no viable, adherent cells in the control LV condition after five days of selection. In contrast, we appreciated viable, adherent cells in the candidate factor LV pool condition and ∼70% of cells demonstrated functional AR bypass as they could survive stringent selection against AR signaling (Figures 1E and S1B). These finding indicated that a discrete set of factors is capable of reprogramming human ARPC, in the context of negative selection for non-reprogrammed cells (modeling highly stringent AR inhibitor therapy), to a cellular state that recapitulates critical phenotypic and functional features of NEPC.

The pioneer neural transcription factor ASCL1 suppresses the AR program during NEtD but requires additional genetic alterations to maintain a proliferative phenotype

To specify the functional contribution of each candidate factor in the LV pool to NE reprogramming, we employed a leave-one-out approach in which each factor was iteratively left out from the LV pool in the C4-2B (Figure 2A) and LNCaP (Figure S2) cell lines. Based on immunoblot analyses, we observed several established regulatory interactions that have been previously identified in cancer. These included 1) ectopic expression of N-Myc leading to elevated levels of NCAM1 in neuroblastoma ³⁵, 2) repression of POU3F2 expression by Rb1 in retinoblastoma ³⁶, 3) and enhanced SOX2 expression with Rb1 inactivation in PC ¹³. Most strikingly, we observed that the exclusion of ASCL1 diminished the decrease in AR and NKX3-1 expression and likewise abrogated the otherwise induced FOXA2, NCAM1, and INSM1 expression, implicating ASCL1 as a critical player in NE transdifferentiation. In contrast, the omission of SRRM4 restored REST expression and suppressed SYP but had little effect on AR or other NE markers. These data are consistent with prior reports in which ectopic introduction of SRRM4 in PC cell lines induced SYP and chromogranin B expression but did not directly affect AR expression ^31,37. Removal of NR0B2, BCL2, and KRAS G12V had no significant effects on AR and NE marker expression (Figures 2A and S2A).

Figure 2.

Next, we employed an orthogonal factor reconstitution strategy to identify the minimal set of factors needed to induce NEtD in the ARPC lines. Immunoblot analyses showed that ASCL1 expression alone downregulated AR markers and induced NE marker expression (Figure 2B). On the other hand, the PRNB factors (dominant-negative TP53 H175R, shRB1, MYCN, and BCL2) were insufficient to enforce NEtD (Figure 2B). SRRM4 alone imparted a partial effect on phenotype as it led to a loss of REST expression and upregulated SYP without reducing AR markers (Figure 2B). As TP53, RB1, MYCN, and BCL2 alterations are common in NEPC and other aggressive NE carcinomas we questioned why they appeared to be dispensable in our NEtD assay. To confirm whether all factors were introduced into the cells, we encoded each lentiviral construct with a unique 10-nucleotide barcode and performed single-cell DNA amplicon sequencing ³⁸ 72 hours after transduction. Of 3,870 cells analyzed, 3,785 cells (97.8%) demonstrated reads from all barcodes, while 84 cells (2.2%) were missing one barcode and 1 cell was missing two barcodes (Figure S2B). We measured cell proliferation for each of the factor reconstitution conditions as an indicator of proliferative fitness in our experimental system. PRNB significantly enhanced proliferation while ASCL1 was detrimental to cell viability over time in the ARPC lines (Figure 2C and Figure S2C). This effect of ASCL1 is consistent with prior reports that ASCL1 promotes cell cycle exit and terminal differentiation when overexpressed in neural progenitor cells and glioblastoma stem cells ³⁹. Notably, SRRM4 promoted proliferation but alone was insufficient to rescue the effects of ASCL1 without PRNB. These data support the concept that the appropriate genetic and proliferative context may be an important determinant of ASCL1’s impact on the emergence of NEPC and highlight the idea that overall fitness may be a combination of proliferation- and lineage-linked benefit. In other words, this mechanism of therapeutic escape may require a change in lineage without a substantial decrease in proliferative fitness.

NeuroD1 is competent to induce neuroendocrine lineage reprogramming

We have previously observed that most NEPCs can be assigned to ASCL1^high, NeuroD1^high, or mixed ASCL1/NeuroD1 groups based on transcriptional profiling ²⁶. These findings have also been described by other groups wherein ASCL1 and NeuroD1 appear to define molecular subtypes of NEPC ²² and the neuroendocrine subtypes of small cell lung cancer ²¹. To determine whether NeuroD1 can also drive lineage conversion from ARPC to NEPC, we performed the reprogramming assay using the C4-2B and LNCaP ARPC lines by replacing ASCL1 with NeuroD1 or combining NeuroD1 with ASCL1. The addition of NeuroD1 to PRNB and SRRM4 substantially diminished AR markers and increased NE marker expression (Figure 3A). We discovered that NeuroD1 induced ASCL1 expression and heightened overall NE marker panel expression. Bulk RNA-seq gene expression analysis of the reprogrammed cell lines with respect to established 22-gene AR signature and NE signature scores revealed a significant de-enrichment of the AR gene signature and enrichment of the NE gene signature in the PRNBSA (PRNB, SRRM4, and ASCL1) and PRNBSN (PRNB, SRRM4, and NeuroD1) conditions relative to the GFP control (Figure 3B). As C4-2B is derived from the LNCaP cell line, we also confirmed the reproducibility of these findings in the independent ARPC cell line MDA PCa 2b using this cellular reprogramming assay (Figure S3A, C).

Figure 3.

Several neural transcription factors including FOXA2, INSM1, POU3F2, and SOX2 have also been described as documented drivers or regulators of NE transdifferentiation in PC and other cancer types ^13,40,41. We sought to determine whether these factors are independently competent to induce NE lineage reprogramming. We introduced PRNB, SRRM4, and each of the above listed neural transcription factors individually to the C4-2B line in the NEtD assay. In contrast to ASCL1 and NeuroD1, these factors were by themselves incapable of NE lineage reprogramming (Figure S3B). As ASCL1 and NeuroD1 also appear to regulate the expression of FOXA2 and INSM1, ASCL1 and NeuroD1 may sit atop a hierarchy of multiple NE-associated transcription factors required in combination to establish and maintain a NEPC transcriptional regulatory network.

To confirm the relevance of our NE transdifferentiation models in PC and to use omic-scale data to minimize the effect of individual marker heterogeneity, we performed partial least squares-discriminant analysis (PLS-DA) of RNA-seq gene expression data from previously characterized human metastatic CRPC patient specimens annotated as ARPC or NEPC ⁴² and projected the C4-2B factor-reconstituted samples onto them (Figure 3C). This analysis showed that PRNB alone did not induce a substantial shift away from the ARPC phenotype. Further, the addition of PRNB and SRRM4 (PRNBS) resulted in a minor, intermediate shift toward the NEPC phenotype. Only upon the addition of ASCL1 or NeuroD1 to PRNB (PRNBA and PRNBN) or PRNBS (PRNBSA and PRNBSN) did we observe a dramatic gene expression program-based shift from ARPC to NEPC. Similar results were observed in in MDA PCa 2b cell line (Figure S3D). These data highlight the relative lineage impacting importance of these factors in NEtD of PC.

SRRM4 promotes alternative splicing and inactivates REST through the inclusion of a neural exon that generates the REST-4 isoform ^31,43. We confirmed the presence of alternatively spliced REST mRNA in C4-2B reprogramming conditions that included SRRM4 (Figure S4A). Further, we performed PLS-DA analysis on alternative splicing data from Beltran et al., 2016 and projected the C4-2B reprogrammed samples onto them. This analysis indicated that SRRM4 drives the alternative splicing phenotype associated with NEPC with a marginal contribution from ASCL1 or NeuroD1 (Figure S4B).

Dynamic changes in cancer phenotype, transcriptional and epigenetic landscapes during NE reprogramming

Having established a defined and tractable system that recapitulates NE transdifferentiation of PC, we sought to investigate the dynamic molecular and regulatory changes underlying this shift in cancer lineage state. We applied the NE reprogramming assay to the ARPC C4-2B line and temporally sampled multiple timepoints within the 14-day reprogramming period (D2, D3, D4, D5, D8, D11, and D14) for downstream analyses (Figure 4A). Immunoblot analyses revealed a reduced expression of AR markers (AR and NKX3-1) as early as D2, with expression of the NE markers INSM1 and FOXA2 evident shortly thereafter on D4-5 (Figure 4B). These findings indicate that downregulation of the AR program precedes the establishment of the NE program in our experimental system. Transcriptome profiling by bulk RNA-seq of the C4-2B cells reprogrammed with PRNBSA or PRNBSN recapitulated the findings from the immunoblot analyses where downregulation of AR expression preceded the induction of INSM1 expression (Figures 4C and D). Notably, PCA analysis of the time course RNA-seq gene expression data from C4-2B cells reprogrammed with PRNBSA or PRNBSN, but not PRNBS, showed an arc-like trajectory during NEtD, indicating dynamic changes in gene expression programs during ARPC to NEPC reprogramming (Figures S5A-C). This arc-like trajectory is also broadly observed in normal and cancer development and cellular differentiation ^44–47.

Figure 4.

We next used PLS-DA projections from the time course RNA-seq gene expression data to further investigate the temporal nature of the phenotypic shift from ARPC to NEPC during NEtD. We distinguished three general phases of reprogramming in the C4-2B cells modified with PRNBSA and PRNBSN that could be defined as “early” on D2, “intermediate” from D3-5, and “late” from D8-14 (Figure 4E). In contrast, these phases were not apparent in C4-2B cells transduced with PRNB and SRRM4 (Figure S5D), further supporting the previous finding that ASCL1 or NEUROD1 is needed for full lineage shift via NEtD.

Lineage reprogramming of ARPC to NEPC is largely driven by epigenetic dysregulation and reprogramming as few genetic alterations apart from the enrichment of TP53 and RB1 mutations and genomic amplification of MYCN have been identified in NEPC ^30,42,48,49. To investigate global changes in the epigenetic landscape during the transition from ARPC to NEPC, we performed Cleavage Under Targets & Release Using Nuclease (CUT&RUN) analysis on PRNBSA reprogrammed C4-2B cells over the 14-day time course to examine genome-wide changes in histone modifications, including monomethylation of histone H3 at lysine 4 (H3K4me1) which is commonly associated with gene enhancers, acetylation of histone H3 on lysine 27 (H3K27ac) which indicates both active promoters and enhancers ^50,51, and trimethylation of histone H3 at lysine 4 (H3K4me3) which marks active gene promoters. Unsupervised principal component analysis (PCA) showed that these histone modifications shift significantly during the time course (Figure 4F), indicating that dynamic epigenetic changes occur during NEtD.

As ASCL1 is a pioneer transcription factor, we queried whether there may be time-dependent changes in putative ASCL1-mediated gene regulation during the acute conversion of ARPC to NEPC. We integrated the RNA-seq gene expression data and ASCL1 CUT&RUN genome occupancy data over the 14-day reprogramming period to identify genes both in proximity to ASCL1 binding and upregulated in PRNBSA compared to PRNBS. Gene set enrichment analysis (GSEA) revealed several pathways enriched early during NE reprogramming including the CDC25 pathway, TGFβ1 pathway, and estradiol response pathway (Figure 5A). In contrast, ASCL1 targets, PRC2 targets including those of the core subunits embryonic ectoderm development (EED) and SUZ12, and neuronal fate commitment pathways were significantly enriched later in the NE transdifferentiation process. The pathways enriched early may represent components of an unstable transcriptional state prior to the establishment of a stable NEPC lineage, consistent with our observation of an arc-like trajectory.

Figure 5.

Next, we evaluated the dynamic activities of active enhancers, super-enhancers, and active promoters in the NE reprogramming assay. We defined active enhancers based on overlapping H3K4me1 and H3K27ac peaks and characterized super-enhancers by applying the Rank Ordering of Super-Enhancer (ROSE) algorithm to H3K27ac peaks over the reprogramming time course ⁵². This analysis revealed a substantial temporal shift in active enhancers (Figure S6A) and super-enhancers (Figure 5B). Gene ontology (GO) analysis of downregulated active enhancer- and super-enhancer-associated genes over time revealed a highly significant enrichment of androgen response genes including the Hallmark Androgen Response gene set (Figures S6B and 5C). In contrast, genes associated with upregulated active enhancers over time were enriched for endocrine therapy resistance pathways and targets of the PRC2 complex core subunits SUZ12 and enhancer of zeste homolog 2 (EZH2) (Figure S6C). We also examined active promoters by determining the overlap of H3K4me3 and H3K27ac peaks over the reprogramming time course (Figure S7A). Genes associated with downregulated active promoters over time were enriched in androgen response and AR regulated targets (Figure S7B). On the other hand, GO analysis of upregulated active promoters over time did not highlight any gene sets of statistical significance.

ASCL1/NeuroD1 inhibits AR expression by remodeling the somatically acquired AR enhancer

To further investigate how ASCL1 and NeuroD1 may suppress AR expression and signaling during NE reprogramming, we performed single-cell (sc) multiome RNA-seq and ATAC-seq analysis on D14 after reprogramming of C4-2B cells with PRNB, PRNBS, PRNBSA, and PRNBSN. The PRNB and PRNBS conditions were included as controls as these factor combinations have no apparent effect on the AR program. Pseudotime analysis revealed a progression from PRNB to PRNBS to PRNBSA to PRNBSN which was associated with a heightened NEPC score and reduced AR expression (Figure S8A-D). In addition, we observed two distinct clusters within the PRNBSA and PRNBSN conditions indicating heterogeneous cell populations likely due to differences in expression levels of ASCL1 and NEUROD1 (Figures S8E, F).

Most strikingly, the integration analysis of RNA-seq and ATAC-seq data identified a significant reduction in chromatin accessibility of the somatically acquired AR enhancer, but not the AR gene body, that was associated with reduced AR expression in the PRNBSA and PRNBSN conditions (Figure 6A). In contrast, the reduction in chromatin accessibility of the somatically acquired AR enhancer was not seen in the control PRNB and PRNBS conditions (Figure 6A). The somatically acquired enhancer of AR and its amplification independent of the AR gene body has previously been shown to be a noncoding driver of CRPC ⁵. Attesting to the clinical relevance of this finding, we identified a hepatic metastasis from a patient with lethal mCRPC who underwent rapid autopsy as part of the University of Washington Tissue Acquisition Necropsy (UW TAN) program that showed a mixed population of AR⁺/ASCL1^- and AR^-/ASCL1⁺ tumors cells by immunofluorescence staining (Figure 6B). Similar to the findings from our NE reprogramming studies, scATAC-seq analysis of this metastatic tumor showed loss of chromatin accessibility at the enhancer of AR but not the AR gene body in AR^-/ASCL1⁺ tumor cells compared to AR⁺/ASCL1^- cells (Figure 6C).

Figure 6.

We next sought to investigate the chromatin accessibility of genome-wide AR binding sites (ARBS) with enhancer activity. ARBS with enhancer activity have previously been characterized using Self-Transcribing Active Regulatory Regions sequencing (STARR-seq) in the LNCaP cell line and classified as constitutively active enhancers, inducible enhancers, and inactive enhancers based on their responsiveness to the androgen dihydrotestosterone ⁵³. We evaluated the chromatin accessibility at these ARBS and identified a global reduction in the accessibility of ARBS with enhancer activity from least reprogrammed (PRNB) to most reprogrammed (PRNBSA or PRNBSN) cell lines (Figure 6D). In contrast, a random selection of genomic sites of similar intervals across all chromosomes, as a control, did not show reduction in chromatin accessibility. These results further demonstrate the acute epigenetic rewiring of PC during NEtD to silence an AR-driven PC lineage program.

ASCL1/NeuroD1 but not loss of REST via SRRM4 activity establish a NE transcriptional network

Lineage-defining transcription factors including neural lineage transcriptions factors have been reported to auto-activate their own expression by binding to super-enhancers to establish a positive feedback loop required to sustain a lineage program ⁵⁴. We investigated the potential for self-regulation of ASCL1 and NeuroD1 in our NE reprogramming assay based on differential RNA-seq read alignment of exogenous proviral transcripts and endogenous transcripts (Figure S9A). We specifically leveraged the presence of the synonymous single-nucleotide variant NM_004316.4(ASCL1);c.627C>G encoded in the ASCL1 lentivirus to distinguish exogenous and endogenous reads. The introduction of either exogenous ASCL1 or NeuroD1 to PRNB or PRNBS resulted in the induction of endogenous ASCL1 expression (Figure S9B). This was also confirmed by the presence of reads mapping to the 5’ and 3’ untranslated regions (UTRs) which are present in endogenous transcripts but not encoded in exogenous proviral transcripts. Using a similar strategy of mapping reads to the 5’ and 3’ UTR, we also found that exogeneous NeuroD1 could activate endogenous NEUROD1 expression (Figure S9C).

We then investigated how the different conditions in our NE reprogramming assay may impact the establishment of a transcriptional regulatory network (TRN) underlying an NEPC lineage state. We first focused on the differential expression of transcription factors due to SRRM4 (PRNBS vs. PRNB), ASCL1 (PRNBSA vs. PRNBS), and NeuroD1 (PRNBSN vs. PRNBS). We observed a significant induction of NE-associated transcription factors including INSM1, SOX2, SOX4, NK2 homeobox 2 (NKX2-2), and FOXA2 related to ASCL1 and NeuroD1 (Figure 7A). In stark contrast, not a single transcription factor demonstrated statistically significant differential expression upon the addition of SRRM4 to PRNB. Consistent with these findings, footprinting analysis of ATAC-seq data identified a significant enrichment in binding motifs associated with these NE-associated transcription factors in conditions of enforced ASCL1 and NeuroD1 expression but not SRRM4 expression (Figure 7B). Global chromatin accessibility profiles are also informative as they can differentiate lineage states in normal tissues and cancer. When we projected pseudo-bulk scATAC-seq data of reprogrammed C4-2B lines on the PCA space defined by ATAC-seq data of previously published ARPC and NEPC PDX models ²², we found that C4-2B lines reprogrammed with PRNBSA and PRNBSN clustered with the NEPC PDXs while those reprogrammed with PRNB and PRNBS were more similar to ARPC PDXs (Figure 7C). Further, we found that the expression of NEPC-associated cell surface targets (DLL3, L1CAM, CEACAM5, RET and SEZ6) were significantly increased and the ARPC-associated cell surface targets (FOLH1, PSCA, STEAP1, and TACSTD2) were significantly decreased upon the addition of ASCL1 or NeuroD1. In contrast, SRRM4 alone did not have a prominent effect on the expression of these cell surface targets (Figure 7D). Taken together, these findings indicate that, unlike ASCL1 and NeuroD1, SRRM4 activity and the resultant loss of REST repressor activity alone may be insufficient to drive a NEPC lineage due to the absence of active transcriptional and epigenetic reprogramming.

Figure 7.

Downregulation of major histocompatibility complex class I (MHC I) in NEPC is functionally attributable to ASCL1/NeuroD1-driven programs

The downregulation of MHC I antigen presentation machinery by cancer to overcome immune surveillance and avoid immune destruction is an established hallmark of cancer. In particular, this phenomenon has been associated with cancers with neuroendocrine differentiation ⁵⁵. In a recent study, examination of SCLC data sets showed that tumors with low NE scores demonstrated higher expression of MHC I genes while those with high NE scores showed decreased MHC I and reduced immune infiltration ⁵⁶. Further, cancer cells are capable of co-opting PRC2 in a lineage-specific fashion to silence the MHC I antigen processing and presentation pathway ⁵⁷. Mechanistically, mutation of the catalytic domain or pharmacologic inhibition of the PRC2 core subunit EZH2 was shown to impair MHC I repression, highlighting the critical role of PRC2 in this process.

We therefore asked whether a similar phenomenon may be operative in NEPC and could be functionally attributable to lineage-defining ASCL1 and/or NeuroD1 expression. We first examined a published RNA-seq gene expression dataset of mCRPC ⁴² including CRPC-Ad (adenocarcinoma) and CRPC-NE (NEPC) patient specimens and found that NEPC was associated with reduced MHC I antigen processing and presentation gene signature scores as well as lower B2M (beta-2-microglobulin, a critical component of MHC I) expression (Figures 8A and S10). A similar analysis using published data from tumors in the UW TAN PC rapid autopsy cohort also showed a reduced overall expression of B2M in the NEPC samples (Figure 8B). We then evaluated RNA-seq gene expression data from our reprogrammed C4-2B lines and correlated each condition and NE signature score with MHC I antigen processing and presentation gene signature scores. Consistent with prior literature, we discovered a reduction in MHC I gene expression with increasing NE differentiation (Figure 8C). Downregulation of MHC I gene expression started as early as D2 along with AR pathway repression (Figure S11). Furthermore, the addition of SRRM4 to PRNB did not have an appreciable effect on MHC I expression, indicating that loss of REST does not itself contribute to the downregulation of MHC I. In contrast, conditions where ASCL1 or NeuroD1 were added did demonstrate reduced expression of MHC I gene sets. We found that B2M expression was especially downregulated in conditions reprogrammed with NeuroD1 which were also associated with prominent NE signature scores (Figure 8C). These findings were in line with integrated scRNA-seq and scATAC-seq data from reprogrammed C4-2B cells showing an inverse correlation between NEPC score and B2M expression (Figure 8D-F) To confirm these findings at the protein level, we performed flow cytometry on reprogrammed C4-2B cells and showed that the PRNBSN condition was associated with a significant reduction in B2M positivity relative to the GFP control condition (29.2% vs. 60.9%) (Figure 8G). In addition, we observed that cells with diminished MHC I gene expression including B2M showed high EZH2 expression in both C4-2B and MDA PCa 2b cell lines (Figure 8C and Figure S12). Together, these findings functionally attribute the lineage-specific downregulation of MHC I in NEPC (relative to ARPC) to the lineage-defining transcription factors ASCL1 and NeuroD1.

Figure 8.

Discussion

AR signaling is pivotal in defining the luminal epithelial lineage of most PCs. However, potent therapies have been developed that apply stringent selective pressure against AR signaling, leading to the increased therapy-driven clinical presentation of lineage plasticity in PC. One such example is the conversion of ARPC to NEPC where the epithelial cancer state enforced by the AR program is replaced by a neural/NE cancer lineage ⁷. Alternative AR-null PC lineages have also been identified including double-negative prostate cancer (DNPC) which lacks both AR and NE marker expression ^18,58. Further, this lineage plasticity in the setting of targeted therapies has also been observed in other epithelial cancer types such as epidermal growth factor receptor (EGFR)-mutant lung adenocarcinoma in which transformation to SCLC occurs in up to 15% of patients treated with EGFR tyrosine kinase inhibitors ⁵⁹. However, the mechanisms by which AR (or EGFR) expression and signaling are effectively silenced during this transition have been poorly understood. The ability to functionally characterize these dynamic processes in a defined manner by reverse genetics has been limited by the lack of available model systems allowing for rapid manipulation and repeated sampling over time. Furthermore, due to documented heterogeneity, the application of narrow phenotypic markers to define NEPC has also been problematic and complicated our understanding of this PC lineage state.

We established an in vitro system in which we can reproducibly and durably reprogram established human ARPC cell lines to NEPC with defined genetic factors over a specified time course. This tractable system recapitulates key features of the process of NEtD including global, dynamic transcriptional and epigenetic reprogramming. Importantly, we uncovered the ability of the pioneer neural transcription factors ASCL1 and NeuroD1 to suppress AR expression and signaling by chromatin remodeling of the somatically acquired enhancer of AR and global ARBS with enhancer activity. These findings show that the activity of a new lineage-defining transcription factor can concurrently silence the prior cancer lineage program while establishing a divergent lineage. An immediate question that emerges from these findings is whether AR and ASCL1 or NeuroD1 activity is mutually exclusive in individual PC cells. Emerging single-cell sequencing data sets of human and mouse models of NEPC indicate that this may be the case ^60,61. The competence of ASCL1 and/or NeuroD1 to silence the AR program in prostate cancer may be a major constraint to NEtD, providing a potential basis for why NEPC only occurs in 15% of therapy-driven mCRPC cases. How the expression of these neural/NE transcription factors is initially induced in PC also remains to be determined and warrants further investigation.

Our studies have also distinguished the relative contributions of SRRM4 and ASCL1 or NeuroD1 to NE reprogramming of PC, and yield insight into the problematic nature of using sparse protein (IHC) or RNA analysis based on only a few markers to define NEPC lineage states. In our reprogramming studies, we find that SRRM4 activity results in the loss of REST which activates neural/NE gene expression including SYP and CHGA. However, unlike ASCL1 or NeuroD1, SRRM4 does not induce a new transcriptional regulatory network based on transcription factor expression and inferred transcription factor activity analyses, does not elicit epigenomic shifts consistent with NEtD, does not affect AR expression and signaling, and does not alter the expression of ARPC-associated cell surface targets. Thus, our findings indicate that a reliance on a sparse set of NE markers, including highly canonical markers, to define NEPC can be flawed as their expression may be independent of true NEPC lineage states. A major clinical consideration is that mCRPCs that express select NE markers but have not undergone NEtD may still be responsive to therapies intended for ARPC including ARSIs and PSMA radioligand therapy.

In our NE reprogramming assay, we have found that ASCL1 and NeuroD1 but not other NE-associated transcription factors are competent to induce NEtD of PC. This supports a transcription regulatory network model where ASCL1 and NeuroD1 sit atop a hierarchy of NE-associated transcription factors necessary to establish self-sustaining NEPC transcriptional convergence. We have shown that ASCL1 and NeuroD1 are capable of self-regulation but the full set of regulatory interactions between these and other NE-associated transcription factors are unknown in NEPC. Future studies to define these relationships through genetic loss-of-function screening using this and other NEtD assays elucidate these functional interactions, identify potential vulnerabilities, and promote therapeutic development to inhibit NEtD under the pressures of otherwise effective targeted therapies.

Prior work identified suppression of MHC I antigen processing and presentation in neuroendocrine cancers and a functional association with lineage-specific PRC2 activity ^55–57. Because our studies directly implicated ASCL1 and NeuroD1 in lineage reprogramming of ARPC to NEPC as well as the induction of PRC2 activity, we examined the expression of MHC I pathway gene sets and discovered that MHC I pathway genes were suppressed in the effective ASCL1- or NeuroD1-driven reprogramming conditions. In contrast, the addition of SRRM4 and subsequent loss of REST activity had little to no effect on the MHC I pathway. These findings support that downregulation of the MHC I antigen processing and presentation pathway in NEPC is functionally attributable to ASCL1 and NeuroD1 activity. Immune checkpoint inhibitors have demonstrated limited utility in mCRPC except in the case of DNA mismatch repair deficiency ⁶². Thus, our results support therapeutic efforts to induce immunogenic responses in NEPC using combinations with epigenetic inhibitors such as those targeting EZH2, histone deacetylases, or demethylases to enhance MHC I expression ^63,64.

Overall, this study provides important mechanistic details of NEtD of PC by highlighting key lineage-defining activities of neural pioneer transcription factors in this process. These findings have important implications on the molecular definition of NEPC and therapeutic considerations related to the application of targeted agents and immune-based therapies.

Methods

Cell lines

LNCaP (Cat# CRL-1740, RRID: CVCL_0395), C4-2B (Cat# CRL-3315, RRID: CVCL_4784), MDA PCa 2b (Cat# CRL-2422, RRID: CVCL_4748) cell lines were purchased from ATCC. All cell lines were validated by short tandem repeat analysis after receipt and on a yearly basis. LNCaP and C4-2B were maintained in RPMI medium supplemented with 10% FBS, 100 U/mL penicillin and 100 μg/mL streptomycin, and 4 mmol/L GlutaMAX. MDA PCa 2b was maintained in F-12K medium supplemented with 20% FBS, 25 ng/ml cholera toxin, 10 ng/ml mouse epidermal growth factor, 0.005 mM phosphoethanolamine, 100 pg/ml hydrocortisone, 1X insulin-transferrin-selenium (ITS-G), 100 U/mL penicillin and 100 μg/mL streptomycin, and 4 mmol/L GlutaMAX.

Lentiviral constructs and lentivirus production

Double-barcoded lentiviral vectors used as a lentiviral backbone vector were generated from FU-CGW by sequentially inserting matched 10-nucleotide barcodes into the Pacl site distal to the HIV FLAP using the Quick Ligation Kit (New England Biolabs) and PCR amplification of the WPRE sequence and barcode with insertion into the Clal site proximal to the 3’ LTR by HiFi DNA Assembly (New England Biolabs). U6 promoter and shRNA cassettes were isolated by digesting pLKO.1 TRC shRNA clones with PspXI and EcorRI and were inserted into the digested double-barcoded plasmid as previously described ^8,9. Lentiviral vectors expressing dominant-negative TP53, BCL2, MYC, SRRM4, NR0B2, KRAS, ASCL1 and NeuroD1 were generated by PCR amplification of the open reading frames and subcloning into the EcoRI site of double-barcoded lentiviral vectors by NEBuilder HiFi DNA Assembly (New England Biolabs). Lentiviruses were prepared and titered as previously described ⁸.

Immunoblotting

Whole cell extracts were collected in 9 M urea lysis buffer and quantified using the Pierce Rapid Gold BCA Protein Assay Kit (Thermo Fisher Scientific). Protein samples were fractionated by SDS-PAGE using Bolt Bis-Tris Plus gels and transferred to a nitrocellulose membrane using Bolt Transfer Buffer according to the manufacturer’s instructions (Invitrogen). Membranes were blocked with 5% non-fat milk in PBST (DPBS + 0.5% Tween 20) for 30 minutes while shaking at room temperature, then incubated with primary antibodies at 4°C overnight. Membranes were washed three times for 5 minutes with PBST and incubated with horseradish peroxidase (HRP) conjugated anti-mouse or anti-rabbit secondary antibodies for 1 hour at room temperature. Blots were washed three times for 5 minutes each with PBST and developed with Immobilon Western Chemiluminescent HRP Substrate (EMD Millipore) for 3 minutes at room temperature. Blot images were acquired using a ChemiDoc Imaging System (Bio-Rad).

Immunohistochemistry (IHC) and immunofluorescence (IF) studies

Cells and xenografts were fixed in 10% buffered formalin for 12 hours and then embedded in HistoGel (Thermo Fisher Scientific) and paraffin, sectioned, and mounted on glass slides (Thermo Fisher Scientific). IHC was performed as previously described ²⁶. In brief, slides containing formalin-fixed paraffin-embedded sections were deparaffinized in xylene and rehydrated in 100%, 95%, and 70% ethanol. Slides were then heated in antigen retrieval buffer (0.2 mol/L citric acid and 0.2 mol/L sodium citrate) within a pressure cooker followed by PBS wash. Slides were blocked with 2.5% horse serum for 30 minutes and then incubated with primary antibody diluted in 2.5% horse serum overnight at 4°C. HRP was detected with ImmPRESS-HRP anti-mouse or anti-rabbit IgG peroxidase detection kits (Vector Laboratories) and staining was visualized with DAB peroxidase substrate (Dako). Slides were counterstained with hematoxylin and dehydrated for mounting.

Immunofluorescence studies using specific antibodies were carried out on archival formalin fixed paraffin embedded tissues. In brief, 5 µm paraffin sections were de-waxed and rehydrated following standard protocols. Antigen retrieval consisted of steaming for 40 minutes in Target Retrieval Solution (Agilent). Slides were then washed and equilibrated in TBS-Tween buffer (Sigma) for 10 minutes. Primary antibodies were applied at at 37°C for 60 minutes. Sequential dual-immunofluorescence labeling was carried out using Tyramide SuperBoost kits (Thermo Fisher).

Neuroendocrine transdifferentiation assay

Cells were seeded in 6-well tissue culture plates at a density of 3 x 10⁵ cells per mL in 3 mL of RPMI medium supplemented with 10% FBS, 100 U/mL penicillin and 100 μg/mL streptomycin, and 4 mmol/L GlutaMAX. Cells were transduced approximately 4-6 hours after seeding at a defined multiplicity of infection (MOI) of 4 for each lentivirus. 72 hours after transduction, cells were trypsinized, washed, and transferred to 100 mm tissue culture plates in 15 mL of neural stem cell media (N-SCM) consisting of Advanced DMEM/F12 medium supplemented with 1X serum-free B27, 10 ng/mL recombinant human bFGF, 10 ng/mL recombinant human EGF, 100 U/mL penicillin and 100 μg/mL streptomycin, and 4 mmol/L GlutaMAX. Media were replenished every 3-4 days. Cells were collected 14 days post-transduction for analysis.

Mouse studies

All animal care and studies were performed in accordance with an approved Fred Hutchinson Cancer Center Institutional Animal Care and Use Committee protocol and Comparative Medicine regulations. Eight-week-old male NSG (NOD-SCID-IL2Rγ-null, RRID:BCBC_4142) mice were obtained from the Jackson Laboratory. A total of 1 x 10⁶ cells from each C4-2B reprogramming condition were suspended in 100 μL of cold Matrigel (Corning) and implanted by injection subcutaneously into NSG mice.

Single-cell DNA amplicon sequencing library preparation and sequencing

A custom panel was designed for the Mission Bio Tapestri to amplify segments of ten mouse genes at two exons each, the 5’ and 3’ lentiviral barcodes, and lentiviral GFP. Libraries were generated from C4-2B cells transduced with PRNBSA LV pool at MOI of 4 for each lentivirus using the Mission Bio Tapestri Single-cell DNA Custom Kit according to the manufacturer’s recommendations. Single cells (3,000 to 4,000 cells per μl) were resuspended in Tapestri cell buffer and encapsulated using a Tapestri microfluidics cartridge, lysed, and barcoded. Barcoded samples were subjected to targted PCR amplification and PCR products were removed from individual droplets, purified with KAPA Pure Beads (Roche Molecular Systems), and used as a template for PCR to incorporate Illumina P7 indices. PCR products were purified by KAPA Pure Beads, and quantified by Qubit dsDNA High Sensitivity Assay (ThermoFisher Scientific). Sample quality was assessed by Agilent TapeStation analysis. Libraries were pooled and sequenced on an Illumina MiSeq or HiSeq 2500 with 150 bp paired-end reads in the Fred Hutchinson Cancer Center Genomics Shared Resource.

Cell proliferation assay

For cell proliferation assays, the CellTiter-Glo® 2.0 Cell Viability Assay system (Promega) was used according to the manufacturer’s instructions. 5×10⁴ cells were plated into each well in a 96-well and 25μL per well of CellTiter-Glo® 2.0 reagent was added at each endpoint for analysis. After a 30 minute incubation at room temperature, luminescence was measured using a BioTek Synergy H1 Multimode Reader (Agilent). Six replicate wells per time point were used to obtain measures of cell proliferation.

RNA-seq

Total RNA was extracted with a Purelink^TM RNA mini kit (Invitrogen). cDNA libraries were prepared from isolated RNA using the Illumina Stranded mRNA prep kit. High-throughput sequencing with 50 bp paired-end reads was performed using an Illumina NovaSeq 6000 SP. Sequencing reads were mapped to human genome reference hg38 and gene expression was quantified and normalized using the Toil RNA-seq pipeline (v4.1.2). Sequence adapters were trimmed using CutAdapt (v1.9) ⁶⁵. Trimmed sequences were then aligned to human reference genome GRCh38 using STAR (v2.4.2a) ⁶⁶ and gene expression quantification was performed using RSEM v1.2.25 ⁶⁷. RSEM expected counts were subsequently normalized by upper-quartile (UQ) normalization for downstream analysis.

The enrichment scores were calculated in R using the GSVA package using a previously published set of 8 AR-regulated genes for AR-activity scores and the 14 genes in the NEURO I and NEURO II gene sets for NEURO scores¹⁸. Differential gene expression analysis was performed using the R package DEseq2 (v1.36.0) RNA-seq data from Beltran et al., 2016 was processed through the Toil RNA-seq pipeline (v4.1.2). Gene expression data was projected onto the Beltran et al., 2016 PCA space to examine reprogramming towards the NEPC phenotype. Gene signature scores were calculated based on RNA-seq expression data using the R package GSVA ⁶⁸. The neuroendocrine score was previously described based on CRPC patient data ⁴². Briefly, 70 genes were used to construct the score based on differential expression or methylation between CRPC-NE and CRPC-Ad or relevance to CRPC-NE based on the literature. The mean expression of the 70 genes was calculated to form the NE score. The scores were visualized in heatmap format using ComplexHeatmap (v2.12.0) ⁶⁹.

Differential RNA alternative splicing analysis

Differential RNA alternative splicing (AS) between vector transduction conditions was calculated from bulk RNA-sequencing using the rMATS (version 4.1.1)⁷⁰. Gene and transcript annotations were sourced from Ensembl GRCH38 release 104⁷¹. The program was run under the statoff mode to obtain the percent spliced in (PSI) matrix for all C42B samples. AS events that were defined by both junction counts and exon counts were used for downstream analysis. The PSI values of all five AS events: skipped exon (SE), alternative 5’ splice site (A5SS), alternative 3’ splice site (A3SS), mutually exclusive exons (MXE) and retained intron (RI) were aggregated together to form a large matrix for PCA. CRPC-Adeno and CRPC-NE Beltran et al 2016 patient RNAseq data was processed using the same rMATS workflow for AS calling. C42B was projected onto Betlran et al 2016 PSI PCA space for examining the remodeling of AS by oncogene transduction towards the NE direction.

CUT&RUN

Cells were harvested by centrifugation at 1200 rpm for 3 minutes in a swinging bucket rotor and washed twice in ice cold wash buffer (20 mM HEPES pH 7.5; 150 mM NaCl; 0.5 mM Spermidine; 1 tablet of Roche complete EDTA-free protease inhibitor). Concanavalin A-coated beads (Bangs Laboratories) were activated by resuspension and washed in binding buffer (20 mM HEPES-KOH at pH 7.9; 10 mM KCl, 1 mM CaCl₂; 1 mM MnCl₂). Cells were mixed with activated beads at a ratio of 9:1 at room temperature for 10 minutes. The cell and bead mixtures were resuspended in antibody solution (2 mM EDTA; 0.05% digitonin in wash buffer). Corresponding antibodies were added at desired concentrations to 150 µl of cell and bead mixture and rotated overnight at 4°C. Libraries were prepared using the AutoCUT&RUN protocol (https://www.protocols.io/view/autocut-run-genome-wide-profiling-of-chromatin-pro-6qpvre6zblmk/v1). High-throughput sequencing with 50 bp paired-end reads was performed on CUT&RUN libraries using an Illumina NextSeq.

Single-cell multiome RNA-seq and ATAC-seq

Cells were collected 14 days post transduction and processed per manufacturer’s instructions using the 10x Chromium Next GEM Single Cell Multiome ATAC + Gene Expression platform to generate libraries. Single-cell RNA-seq libraries and single-cell ATAC-seq libraries were sequenced on an Illumina NovaSeq 6000 targeting 35,000 read pairs/nuclei for RNA-seq and 40,000 read pairs/nucleus for ATAC-seq.

Single-cell RNA-seq and ATAC-seq analysis

Single-cell multiome ATAC- and RNA-sequencing reads were processed through the CellRanger ARC analysis pipeline. RNA-sequencing filtered feature count matrices and ATAC-sequencing fragment paths from each sample were loaded individually loaded into R with Signac (version 1.7.0)⁷². The data was filtered for total RNA count between 1000 and 50000, total feature RNA count between 500 and 10000, percent mitochondria RNA less than 20%, total ATAC fragment count between 1000 and 50000, nucleosome signal score less than 2, and a TSS enrichment score above 1.5. ATACseq peaks were called using MACS2 ⁷³. Peaks on nonstandard chromosomes and in genomic blacklist regions of hg38 were removed. ATAC counts were quantified for all peaks across all samples.

Single-cell trajectory analysis

Filtered Seurat object was loaded into Monocle3 (version 1.2.9) to calculate Single-cell trajectory ^74,75. First 50 PCA dimensions were used in the preprocessing step. Cells were aligned by regressing on percent mitochondria RNA.

CUT&RUN data processing and analyses

Peak calls for H3K4me1, H3K4me3 and H3K27ac from CUT&RUN data at different timepoints were performed using SEACR (v1.3) ⁷⁶. Bedgraph files were prepared using bedtools (v2.27.1) to serve as the inputs for SEACR. IgG control conditions at each timepoint were used to define empirical thresholds for peak calling. The “norm” parameter was used in SEACR to normalize the data. Peaks were called in “stringent” mode. Active promoters were defined by checking if H3K4me3 peaks overlap with H3K27ac signals using the bedtools intersect -wa option. The resulting regions within +/- 5 kb of transcription start sites (TSS) were retained. TSS regions were calculated by applying GTFtools (v0.8.5) ⁷⁷ to the gencode v39 (GRCh38.p13) GTF file. Active enhancers were obtained by overlapping H3K4me1 peaks with H3K27ac peaks, followed by subtracting the +/- 5 kb TSS regions using the bedtools subtract –A option. Super-enhancers were called using ROSE (v0.1) ^52,78 on H3K27ac CUT&RUN data. STITCHING_DISTANCE was set to 12.5 kb. TSS_EXCLUSION_ZONE_SIZE was set to 2.5 kb. hg38 was used as the reference genome. The resulting genomic regions of active promoters, active enhancers and super-enhancers were loaded into DiffBind ⁷⁹ for differential binding analysis. Counts that mapped to the genomic regions were normalized using the Relative Log Expression (RLE) method. For super-enhancers, the counts of member regions that belong to each parental broad genomic region were aggregated together to form the overall count. The correlations between time and genomic regions were calculated using Kendall test followed by Benjamini-Hochberg P-value adjustment. Genomic regions were mapped to genes using ChIPseeker (v1.32.0) ⁸⁰. Genomic regions and genes demonstrating significant changes along with time were visualized in heatmap format using ComplexHeatmaps (v2.12.0)⁶⁹ and analyzed by Gene Ontology using hypergeometric test with MSigDB gene sets ⁸¹.

AR binding site analyses

Inducible, constitutively active, and inactive AR binding sites (ARBS) were obtained from a previously published study ⁵³. A random set of genomic regions with the same length of ARBS (700 bp) was generated by randomly selecting 50 start sites for each chromosome. Fragments of scATAC-seq data were quantified against these sites. The significance of the difference between conditions was calculated based on the log transformed counts using Welch’s t-test.

AR enhancer analyses

The genomic coordinates of the somatically acquired AR enhancer were obtained from a previously published study ⁵. scATAC-seq data were examined at the same overall genomic region to determine the location of the AR enhancer in our model system and in the prostate cancer hepatic metastasis. The correlation between scATAC chromatin accessibility and scRNA gene expression in reprogrammed cells was visualized using CoveragePlot in the R package Signac (version 1.7.0) ⁷².

Statistical analyses

Data analysis was performed on GraphPad Prism 9 (GraphPad Software, Inc.). Quantitative PCR results were analyzed in Excel. Statistical significance was determined using the unpaired two-tailed Student t test, unpaired two-tailed Welch t test where the variances are shown to be different via F-test, one-way ANOVA, or two-way ANOVA. Only two-tailed tests were used. Results are depicted as mean + SD unless stated otherwise. All P values of <0.05, <0.01, <0.001, and <0.0001 were considered significant. Pearson correlation coefficient was used to determine correlation between genes. The symbols used to represent the P values were: ns, nonsignificant for P > 0.05; *, P ≤ 0.05; **, P ≤ 0.01; ***, P ≤ 0.001; ****, P ≤ 0.0001. The test used in each statistical analysis is specified in the figure legends.

Data availability

All raw and analyzed sequencing data including RNA-seq, scRNA-seq, and scATAC-seq have been uploaded to the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) under accession SuperSeries GSE225026 and are publicly available for download.

Acknowledgements

We would like to thank the Fred Hutch Genomics Shared Resource, Comparative Medicine Shared Resource, and Experimental Histopathology Shared Resource (supported by NIH/NCI Cancer Center Support Grant P30CA015704). We also acknowledge support by the Pacific Northwest Prostate Cancer SPORE (NCI P50CA97186) and the Institute for Prostate Cancer Research. SL was supported by a Department of Defense Prostate Cancer Research Program Early Investigator Research Award (W81XWH-20-1-0083). The work was also supported by NIH DP2 CA271301 (J.K.L.); NCI R01 CA234715 (J.K.L., M.C.H., P.S.N.); NCI R01 CA266452 (J.K.L., P.S.N., M.C.H.); NCI R01 CA280056 (P.S.N.); NCI R01 CA222877 (T.G.G.); NCI P50CA092131 (T.G.G., J.K.L.); NCI P01CA163227 (P.S.N., H.W.L.); the W.M. Keck Foundation (T.G.G.); the V Foundation (M.C.H.); the UCLA Eli and Edythe Broad Center of Regenerative Medicine and Stem Cell Research (T.G.G.); Hal Gaba Director’s Fund for Cancer Stem Cell Research (T.G.G.); and Prostate Cancer Foundation Challenge Awards (J.K.L., P.S.N.).

Additional information

Author contributions

Conceptualization, S.L., K.S., H.S., T.G.G., and J.K.L.; Methodology, S.L. and J.K.L.; Investigation, S.L., K.S., H.S., Y.T., A.H., V.B., B.H., R.A.P., H.W.L., C.M., M.C.H.; Writing – Original Draft, S.L. and J.K.L.; Writing – Review & Editing, P.S.N., T.G.G., and J.K.L.; Resources, C.M., P.S.N.; Supervision, T.G.G. and J.K.L.

Additional files

Supplementary Materials