An Atlas of Accessible Chromatin in Advanced Prostate Cancer Reveals the Epigenetic Evolution during Tumor Progression

Abstract

Metastatic castration-resistant prostate cancer (mCRPC) is a lethal disease that resists therapy targeting androgen signaling, the primary driver of prostate cancer. mCRPC resists androgen receptor (AR) inhibitors by amplifying AR signaling or by evolving into therapy-resistant subtypes that do not depend on AR. Elucidation of the epigenetic underpinnings of these subtypes could provide important insights into the drivers of therapy resistance. In this study, we produced chromatin accessibility maps linked to the binding of lineage-specific transcription factors (TF) by performing assay for transposase-accessible chromatin sequencing on 70 mCRPC tissue biopsies integrated with transcriptome and whole-genome sequencing. mCRPC had a distinct global chromatin accessibility profile linked to AR function. Analysis of TF occupancy across accessible chromatin revealed 203 TFs associated with mCRPC subtypes. Notably, ZNF263 was identified as a putative prostate cancer TF with a significant impact on gene activity in the double-negative subtype (AR- neuroendocrine-), potentially activating MYC targets. Overall, this analysis of chromatin accessibility in mCRPC provides valuable insights into epigenetic changes that occur during progression to mCRPC. Significance: Integration of a large cohort of transcriptome, whole-genome, and ATAC sequencing characterizes the chromatin accessibility changes in advanced prostate cancer and identifies therapy-resistant prostate cancer subtype-specific transcription factors that modulate oncogenic programs.

Figure 1.. Chromatin accessibility changes during prostate cancer progression affect stage-specific regulatory elements.

(a) Principal component analysis (PCA) of the ATAC-seq profiles comparing different stages of prostate cancer including benign prostate, localized prostate cancer (PCa), mCRPC Adeno, and mCRPC t-SCNC. The normalized read counts of these consensus-accessible regions were used for the PCA analysis. Each dot in the plot represents an individual sample. (b) An alluvial plot demonstrating changes in accessible chromatin regions in various stages of prostate cancer. Each bar corresponds to a distinct PCa stage, with the orange and white sections indicating accessible and inaccessible chromatin regions, respectively. The shaded areas connecting the bars represent changes in the accessibility of these chromatin regions. The pink and blue shaded regions respectively represent accessible and inaccessible chromatin regions in mCRPC Adeno. (c) Heatmap representation of the chromatin between different stages. The rows are segregated by the chromatin variants in each stage. (d) The percentage of chromatin variants in mCRPC Adeno. The ATAC-seq peaks are grouped by the genomic regions (promoter, intron, or distal intergenic) to which they are mapped. (e) ATAC-seq profile plot illustrating potential regulatory regions in chromatin variants in mCRPC Adeno. The profile plot represents the overlapping region between the chromatin variants and publicly available H3K27ac ChIP-seq data from mCRPC PDX. (f) Enrichment of chromatin regions exclusively accessible in localized PCa, mCRPC Adeno, or mCRPC t-SCNC against GO Biological Processes.

Figure 2.. Chromatin accessibility in mCRPC is associated with subtypes linked to androgen signaling.

(a) Unsupervised hierarchical clustering of pairwise sample Spearman’s correlation based on the normalized read counts of consensus ATAC-seq peaks of mCRPC. (b-c) Distribution of (b) androgen receptor (AR) pathway score, (c) neuroendocrine (NE) score calculated based on RNA-seq gene expression profiles of mCRPCs classified into individual clusters in Figure 2a. Statistically significant Wilcoxon rank sum test p-values between the clusters are indicated in the plot.

Figure 3.. mCRPC transcriptional subtypes are associated with chromatin variants of prostate cancer signaling pathways.

(a) Heatmap of chromatin variants between mCRPC transcriptional subtypes. (b) ATAC-seq peaks around AR, NKX3-1, and GPR37L1 gene regions. The highlighted vertical strip illustrates the presence of ATAC-seq peaks at the enhancer region. (c) Heatmap representation of the chromatin variants between mCRPC transcriptional subtypes. The rows are segregated by the differential regions mapped to gene promoter, intron, or distal intergenic regions. (d-f) Hallmark pathways enrichment of chromatin variants between the mCRPC transcriptional subtypes mapped to (d) promoter, (e) intron, and (f) distal intergenic regions identified using GREAT enrichment analysis (see Methods section).

Figure 4.. mCRPC transcriptional subtypes are defined by DNA accessibility-guided patterns of transcription factor regulation.

(a-b) ATAC-seq TF (AR and HOXB13) footprints signal difference between TF-bound and unbound sites. (c-d) Volcano plot of differential TF footprint occupancy analysis comparing the (c) AR+NE− and AR−NE+ subtypes and (d) AR+NE− and AR−NE− subtypes. Each dot in the plot represents a TF motif. The colored dots indicate a significantly differentially bound TF motif. (Data available as Supplementary Table S4) (e) Heatmap of genome-wide active TF occupancy, determined by TF footprints, associated with different mCRPC transcriptional subtypes. Each rim of the circular heatmap represents an individual mCRPC transcriptional subtype and the sector represents TF. The darker color shade indicates the strong association of the TF with the respective mCRPC subtype. See the Methods section for details on TF occupancy phenotype score calculation.

Figure 5.. Identification of the influential transcription factors driving mCRPC transcriptional subtypes.

We hypothesized that highly active TF regulate (or influence) gene expression activity of a large fraction of target genes. The plot indicates the top influential mCRPC transcriptional subtype-associated TFs ranked by the number of target genes (based on gene expression) they influence.

Figure 6.. ZNF263 activates MYC signaling targets.

(a) The volcano plot depicts the genes that undergo activation or repression upon ZNF263 binding to their specific promoter region. Each gene is represented by a dot, and the difference in gene expression between samples with and without ZNF263 in the promoter region was measured as fold change. Additionally, the statistical significance of the difference was evaluated using the Wilcoxon rank sum test to calculate the p-value between the two groups. (b) Over-representation analysis of the predicted ZNF263 target genes against the Hallmark pathways. (c) Percentage of the predicted target genes of MYC that overlap with those of TFs associated with AR−NE− subtype. The overlap of the respective TF target genes with genes in the Hallmark MYC targets geneset is illustrated as the red line. (d) Heatmap of overlapping ZNF263 and MYC footprint sites. The red color highlights the direct overlap between ZNF263 and MYC footprints. (e) Volcano plot showing the genes that are activated when ZNF263 binds to the promoter region compared to genes that are activated with both ZNF263 and MYC are not bound to the promoter region. (f) Box plot of gene expression foldchange when different combinations of ZNF263 and MYC bind to the promoter as compared to when both ZNF263 and MYC are simultaneously absent in the promoter. Each dot represents a gene in Hallmark MYC targets geneset.