The landscape of RNA polymerase II-associated chromatin interactions in prostate cancer

Abstract

Transcriptional dysregulation is a hallmark of prostate cancer (PCa). We mapped the RNA polymerase II-associated (RNA Pol II-associated) chromatin interactions in normal prostate cells and PCa cells. We discovered thousands of enhancer-promoter, enhancer-enhancer, as well as promoter-promoter chromatin interactions. These transcriptional hubs operate within the framework set by structural proteins - CTCF and cohesins - and are regulated by the cooperative action of master transcription factors, such as the androgen receptor (AR) and FOXA1. By combining analyses from metastatic castration-resistant PCa (mCRPC) specimens, we show that AR locus amplification contributes to the transcriptional upregulation of the AR gene by increasing the total number of chromatin interaction modules comprising the AR gene and its distal enhancer. We deconvoluted the transcription control modules of several PCa genes, notably the biomarker KLK3, lineage-restricted genes (KRT8, KRT18, HOXB13, FOXA1, ZBTB16), the drug target EZH2, and the oncogene MYC. By integrating clinical PCa data, we defined a germline-somatic interplay between the PCa risk allele rs684232 and the somatically acquired TMPRSS2-ERG gene fusion in the transcriptional regulation of multiple target genes - VPS53, FAM57A, and GEMIN4. Our studies implicate changes in genome organization as a critical determinant of aberrant transcriptional regulation in PCa.

Keywords: Epigenetics; Genetics; Oncology; Prostate cancer; Transcription.

Figure 1. Analysis of RNA Pol II–associated chromatin interactions.

(A) Pipeline for ChIA-PET data processing and identification of chromatin interactions. (B) Enhancers per gene and genes per enhancer for each cell line. For enhancers per gene, all expressed genes (fragments per kb per million mapped reads [FPKM] >1) were included; for genes per enhancer, all enhancers located in RNA Pol II peak anchor regions were included. Box plot represents the median and the 25% and 75% quantiles, with lines at 1.5 times the IQR. The significance for each pair comparison was tested using the Kolmogorov-Smirnov test (^†P < 2.2 × 10^–16).The analysis was normalized by adjusting for sequencing depth. (C) Correlation between RNA expression and chromatin interaction across the 4 cell lines. The expression level was measured by FPKM transformed by log₂, and chromatin interactions were measured by the number of promoter PETs. The correlation efficient was calculated by Spearman’s correlation, and P values are shown. (D) Significant gain and loss of RNA Pol II–associated chromatin interactions in the 4 cell lines. The 3 solid lines show the distributions of the 3 PCa cell lines compared with benign cells (RWPE-1), whereas the 3 broken lines show the distributions among the 3 PCa cell lines. The significant gain or loss interactions were obtained using the DNB model. The total number of interactions that were significantly altered between cell lines are listed accordingly. (E) Scatter plot shows the correlation between changes in gene expression and changes in promoter-associated chromatin interactions. Graphs show Spearman’s correlation coefficients for the top 10%, 20%, 40%, 60%, and 80% differentially expressed genes (right panel). FC, fold change. (F) Total number of PETs at promoter regions in the 4 cell lines. Genes are ranked by increasing number of PETs. Eight representative genes are labeled 1 through 8 in the plots.

Figure 2. Integrative analysis of chromatin interactions.

(A) Integrated genome view representing DNase-Seq; ChIP-Seq for CTCF, H3K4me3, and H3K27ac; RNA-Seq; ChIA-PET for RNA Pol II and RAD21; and Hi-C data from the LNCaP cell line for a representative region with both active and repressed genes. Chr1, chromosome 1. (B) Analysis of overlap between H3K27ac ChIP-Seq and RNA Pol II ChIA-PET data for 4 cell lines. Left: overlap in the number of total peaks; right: overlap in the number of peaks with intrachromosomal interactions. (C) Analysis of overlap between CTCF ChIP-Seq and RAD21 ChIA-PET data for LNCaP and DU145 cells. Left: overlap in the number of total peaks; right: overlap in the number of anchor peaks with intrachromosomal interactions. (D) Comparison of PET lengths between RNA Pol II ChIA-PET and RAD21 ChIA-PET data from LNCaP and DU145 cells. Box plot represents the median and the 25% and 75% quantiles, with lines at 1.5 times the IQR. ^†P < 2.2 × 10^–16, by Kolmogorov-Smirnov test. (E) Integrated genome view of the KRT8-KRT18 gene cluster and its neighborhood in LNCaP cells. The data tracks represent RNA-Seq; ChIP-Seq for CTCF, FOXA1, AR, H3K27ac, and RNA Pol II; ChIA-PET for RNA Pol II and RAD21; and Hi-C data for the LNCaP cell line. (F) ChIA-PET contact heatmap representing RNA Pol II– and RAD21-associated chromatin interactions for the KRT8-KRT18 gene cluster and neighborhood regions. The KRT8-KRT18 gene cluster is shown in light green.

Figure 3. Transcriptional regulation of AR and FOXA1 loci.

(A and B) Integrated genome view of the AR gene and its adjacent regions from –1400 kb to +400 kb in LNCaP and VCaP cells. RNA Pol II ChIA-PET, RNA-Seq and CTCF, FOXA1, AR, H3K27ac, and RNA Pol II ChIP-Seq data are shown. In addition, AR ChIA-PET, phospho–RNA Pol II, and ERG ChIP-Seq are shown for VCaP cells. The AR gene and upstream regions are highlighted in light blue. (C) Summary of the copy number aberrations associated with the AR and its enhancers. Heatmap shows shows the aCGH high-density probes for 27 patients with mCRPC. Gains are depicted in pink and losses in light blue, whereas amplifications are shown in red and deep deletions in dark blue. Each column is a probe in the aCGH platform, and each row represents a sample. Probes that cover the EDA2R, AR, and regions of the enhancer peaks are shown. (D) Schematic representation of a deletion between the AR gene and its enhancers. (E) Comparison of RNA Pol II–associated chromatin interactions at the FOXA1 locus and its adjacent regions in the 4 cell lines.

Figure 4. Chromatin interaction–associated transcriptional targets of the AR and FOXA1.

(A) The expression of AR/FOXA1 target genes discovered by integrating RNA Pol II ChIA-PET with AR/FOXA1 ChIP-Seq was compared with the expression of genes that were nearest to AR/FOXA1 binding peaks according to ChIP-Seq data and randomly selected control genes in VCaP and LNCaP cells. The y axis represents expression levels, measured as FPKM transformed by log₂. The box plots represent the median and 25% and 75% quantiles, with lines at 1.5 times the IQR. ^†P < 2.2 × 10^–16, by Kolmogorov-Smirnov test). (B) Gene promoters that interact with AR, FOXA1, and AR-FOXA1 co-occupied regions in the RNA Pol II ChIA-PET data sets are shown in yellow. Gene promoters that interact with AR, FOXA1, and AR-FOXA1 co-occupied enhancers in the RNA Pol II ChIA-PET data sets are shown in orange. (C) Pathway analysis for gene promoters that interact with AR-FOXA1 co-occupied regions in LNCaP and VCaP cells in the RNA Pol II ChIA-PET data sets.

Figure 5. Transcriptional regulation the KLK gene cluster and ZBTB16.

(A) Comparison of RNA Pol II–associated chromatin interactions and H3K27ac ChIP-Seq signals in the KLK3 gene and its neighborhood regions from –350 Kb to +350 Kb in 4 cell lines. (B and C) Transcriptional regulation of the ZBTB16 gene in (B) LNCaP and (C) DU145 cells.

Figure 6. Transcriptional regulation of the MYC and EZH2 genes.

(A) Comparison of chromatin interactions, H3K27ac ChIP-Seq, and RNA-Seq signals of the MYC gene and its adjacent regions from –1000 kb to +600 kb in 4 cell lines. The MYC gene is highlighted in light blue. PCa risk SNP loci located in the MYC neighborhood are shown. (B) Left: Venn diagram describes the number of PCa risk SNP loci located within the coordinates of RNA Pol II peaks for 3 cell lines. Right: Venn diagram shows the number of PCa risk SNP loci located within the coordinates of anchor regions of RNA Pol II–associated chromatin interactions for 3 cell lines. (C) Transcriptional regulation of the EZH2 gene and its neighboring genes. Comparison of RNA Pol II ChIA-PET, H3K27ac ChIP-Seq, and RNA-Seq signals at the EZH2 locus and its adjacent regions from –150 kb to +800 kb in 4 cell lines. (D) Immunoblot representing endogenous EZH2 expression in RWPE-1, LNCaP, VCaP, and DU145 cells. The EZH2 and actin blots were obtained from separate gels that were run contemporaneously.

Figure 7. Evaluation of RNA Pol II–associated peaks and interaction with 122 prostate-specific germline SNP locations.

(A) Peak analysis. The red dashed lines indicate the observed number of peaks containing SNPs for each cell line. The histograms illustrate the results from 10,000 simulations that assessed the expected number of peaks containing SNPs. The mean of the simulations is shown with a dashed black line. ***P < 0.001, for significant differences between the expected and observed values. RWPE-1, LNCaP, VCaP, and DU145 have 17, 13, 14, and 4 peaks that overlap SNPs, respectively. The black dashed lines indicate the expected number of overlapping peaks (the mean of all the simulations). The expected values for RWPE-1, LNCaP, VCaP, and DU145 cells are 4.88, 4.09, 4.03, and 3.44, respectively. (B) Interaction analysis. The same procedure was repeated except using only the peaks involved in interactions. The red dashed lines indicate the observed number of SNPs, and the black dashed lines show the expected values. *P < 0.05 and ***P < 0.001, for significant differences between the expected and observed values. RWPE-1, LNCaP, VCaP, and DU145 cells had an observed value of 10, 9, 5, and 1, respectively, as indicated by the red dashed lines. The expected values for RWPE-1, LNCaP, VCaP, and DU145 cells are indicated by the black dashed lines and were 1.81, 2.17, 2.60, and 1.30, respectively. (C) Enrichment analysis. Fisher’s enrichment analysis was performed to compare the number of SNP-positive peaks with the rest of the genome as well as to compare the number of SNP-positive and interaction-positive peaks with the rest of the genome. ***P < 0.001, for significant enrichment. Fisher’s exact test was used to determine the P values for C. Empirical method was used to determine the P values for A and B.

Figure 8. Transcriptional regulation by the PCa risk SNP rs684232.

(A) Integrated genome view of RNA Pol II–associated chromatin interactions in the genomic region harboring the PCa risk SNP rs684232. (B) rs684232 was significantly associated with mRNA abundance of FAM57A, VPS53, and GEMIN4 in the CPC-GENE, TCGA, and Porto cohorts. Box plots represent the median and the 0.25 and 0.75 quantiles, with whiskers at 1.5 times the IQR. mRNA abundance was measured in FPKM. The numbers below the genotypes indicate the number of samples in each group. P values and effect size are from a linear model. (C) Epigenetic features of the PCa risk SNP rs684232 locus in VCaP cells are described using ChIP-Seq analysis. (D) rs684232 falls in an active enhancer region, and the alternative allele was found to be significantly associated with decreased H3K27ac binding. Heatmap shows H3K27ac ChIP-Seq signal within chr17:614900-622900 (x axis) for 92 patients (y axis). Color indicates ChIP-Seq signal intensity, and the black bar in the covariate along the top indicates the location of rs684232. Box plot shows H3K27ac signal intensity stratified by genotype in the Porto cohort (Mann-Whitney U test for the recessive model). The y axis represents the number of H3K27ac ChIP-Seq read counts mapped to the SNP rs684232 region, which were normalized by the trimmed mean of M values (TMM). (E) Box plots show H3K27ac signal intensity in the promoter regions of FAM57A and GEMIN4 stratified by genotype in the Porto cohort (Mann-Whitney U test for the recessive model). P values and effect size are from a linear model. (F) Sequence analysis to confirm the cloning of the WT and risk (rs684232) alleles in the pGL2 luciferase reporter plasmid. (G) Luciferase reporter assays in LNCaP and VCaP cells. Cells were cotransfected with pSV-Renilla and the luciferase reporter encoding the WT or risk (rs684232) allele and processed 48 hours after transfection. Firefly Luc/Renilla Luc activity was determined (mean ± SD, n = 6; ****P < 0.0001, by 2-tailed Student’s t test).

Figure 9. Transcriptional regulation and clinical correlates of the chromatin interaction targets of the PCa risk SNP rs684232.

(A) qRT-PCR validation of AR knockdown and the expression of VPS53, GEMIN4, and FAM57A genes upon treatment of LNCaP cells with AR siRNA. **P < 0.01, and ***P < 0.001, by 2-tailed Student’s t test. Error bars indicate the SD of 3 technical replicates. (B) Box plot shows AR ChIP-Seq signal intensity stratified by genotype in the Porto cohort (Mann-Whitney U test for the recessive model). The y axis shows the number of AR ChIP-Seq read counts mapped to the SNP rs684232 region, which were normalized by the TMM method. Box plot represents the median and the 0.25 and 0.75 quantiles, with whiskers at 1.5 times the IQR. (C) The regulatory impact of rs684232 was enhanced in the presence of the TMPRSS2-ERG fusion. Box plots show the mRNA abundance (FPKM) of each gene stratified by genotype and further split by ERG status in the CPC-GENE cohort. β_positive and β_negative, and the associated P values quantify the eQTL within ERG-positive and -negative patients, respectively (linear model). (D) Box plots show the mRNA abundance of FAM57A, GEMIN4, and VPS53 genes in PCa specimens from TCGA cohort. Tumors were classified into various ISUP grade groups. Relationship between mRNA abundance and ISUP Grade Group was quantified using Spearman’s correlation and represented as Spearman’s rho and corresponding P values. The P values for Spearman’s correlation were computed using the algorithm AS 89 (47). (E) BCR-free survival curves for PCa patient groups defined by transcript abundance for FAM57A, GEMIN4, and VPS53 genes in the CPC-GENE cohort. P values in E were determined by log-rank test.