Reprogramming transcription by distinct classes of enhancers functionally defined by eRNA.

Mammalian genomes are populated with thousands of transcriptional enhancers that orchestrate cell-type-specific gene expression programs, but how those enhancers are exploited to institute alternative, signal-dependent transcriptional responses remains poorly understood. Here we present evidence that cell-lineage-specific factors, such as FoxA1, can simultaneously facilitate and restrict key regulated transcription factors, exemplified by the androgen receptor (AR), to act on structurally and functionally distinct classes of enhancer. Consequently, FoxA1 downregulation, an unfavourable prognostic sign in certain advanced prostate tumours, triggers dramatic reprogramming of the hormonal response by causing a massive switch in AR binding to a distinct cohort of pre-established enhancers. These enhancers are functional, as evidenced by the production of enhancer-templated non-coding RNA (eRNA) based on global nuclear run-on sequencing (GRO-seq) analysis, with a unique class apparently requiring no nucleosome remodelling to induce specific enhancer-promoter looping and gene activation. GRO-seq data also suggest that liganded AR induces both transcription initiation and elongation. Together, these findings reveal a large repository of active enhancers that can be dynamically tuned to elicit alternative gene expression programs, which may underlie many sequential gene expression events in development, cell differentiation and disease progression.

The wide diversity of mammalian cells is determined by a large repertoire of constitutive and inducible genes, which are regulated by general and cell-type specific transcription factors and cofactors through regulatory genomic elements7,8. Recent studies reveal that gene promoters are marked by tri-methylated H3K4 (H3K4me3) and distal regulatory elements are often associated with mono-methylated H3K4 (H3K4me1)1,2. Because these H3K4me1-positive, H3K4me3-negative regions exhibit striking cell type specificity1,2, we used this signature to characterize potential enhancers in prostatic LNCaP cells in which one of key regulatory transcriptional programs is mediated by the androgen receptor (AR). We identified by ChIP-seq 14,283 H3K4me3-marked and 51,544 H3K4me1-marked loci in androgen (5α-dihydrotestosterone, DHT)-treated LNCaP cells, among which 43,565 loci are uniquely marked by H3K4me1, largely localized distal to annotated TSSs (94%), and associated with other marks linked to enhancer activities (Fig. 1a).

Figure 1

FoxA1 contributes to the enhancer code in prostate cancer cells

(a) Distribution of histone marks within ±2Kb windows around distinct genomic regions (n=43,565) marked by H3K4me1, but not H3K4me3, in androgen (DHT)-stimulated LNCaP cells. The ChIP-seq datasets for H3K4me1, H3K4me2, H3K4me3, H3K27ac, H4K5ac and p300 were each aligned with respect to the center of theH3K4me1 signal and sorted by the length of H3K4me1-marked regions. (b) Top-enriched DNA motifs with significant P-value and prospective families of DNA binding transcription factors identified by de novo motif analysis of non-promoter regions marked by H3K4me1. (c) Percentage of H3K4me1-marked regions that show FoxA1 binding events (top panel) and percentage of FoxA1 binding sites that are marked by H3K4me1 (bottom panel). Note that H3K4me1-marked regions tend to be broad, but FoxA1 binding sites are discrete, and as a result, many H3K4me1-positive regions may contain more than one FoxA1 binding site. (d) Genomic distance from FoxA1/H3K4me1-positive loci to the nearest TSS of genes in response to FOXA1 knockdown. Outliers were omitted from box plots. P-values indicate the significance in pair-wise comparisons. (e-g) Three classes of FoxA1/H3K4me1-positive loci according to the response in levels of H3K4me1 to FOXA1 knockdown: >1.5-fold decrease (e), no significant change (f), and >1.5-fold increase (g). (h) Ratio (log2) of up- and down-regulated genes in each H3K4me1 responsive category in e-g.

De novo DNA motif analysis revealed several highly enriched motifs, particularly the forkhead motif (Fig. 1b). Using a specific antibody against FoxA1, a major FOX family member expressed in LNCaP cells and normal prostate gland9-11 (Supplementary Fig. 1), we identified 33,426 FoxA1-bound sites, which extensively overlap with distal H3K4me1-marked regions (Fig. 1c and Supplementary Fig. 2a; see on KLK3 enhancer12 in Supplementary Fig. 2b). RNA profiling supports the functional relevance of these FoxA1/H3K4me1 loci, as genes responsive to FOXA1 siRNA are located more proximally to FoxA1/H3K4me1-marked loci than non-responsive genes (Fig. 1d and Supplementary Fig. 3).

FoxA1 has been characterized as a “pioneer” factor to facilitate DNA binding by other sequence-specific transcription factors9,13-16 and “translate” H3K4me1/me2 into AR-mediated gene expression9. Comparing the profile of H3K4me1 and H3K27ac before and after FOXA1 knockdown, we detected three classes of FoxA1 binding sites based on the H3K4me1 signal exhibiting reduced (~22%), relatively unaffected (~74%), or even increased (~3.4%) levels over candidate enhancers (Fig. 1e-g and Supplementary Fig. 4). RNA profiling analysis agrees with the functional significance of these selective FoxA1 effects, revealing more down-regulated genes in the first class, roughly equal numbers of up- or down-regulated genes in the second, and more up-regulated genes in the third (Fig. 1h), suggesting a contribution of FoxA1 to “writing” and “reading” the “histone code” on different enhancer cohorts, in line with its critical function in prostate gland development10,11.

The rationale for our experimental strategy to use RNAi to study FoxA1-regulated enhancer network is the association of decreased FOXA1 expression with castration-resistant, poor prognostic prostate tumors (Supplementary Fig. 5). In LNCaP cells, FOXA1 RNAi enhanced cell entrance to S phase with reduced hormone (Fig. 2a). To understand the mechanistic basis for elevated hormone responsiveness, we mapped AR binding sites, identifying 3,115 high confident loci with ~65% co-incident with H3K4me1. De novo motif analysis revealed highly enriched elements for both AR and FoxA1, including a composite motif consisting of a FOX motif and AR regulatory element (ARE) half site, suggesting ternary complex formation on these sites (Fig. 2b). Indeed, 1,684 AR-bound loci (54% of total) are co-occupied by FoxA1 in DHT-treated LNCaP cells and FoxA1 appears to bind to most of these sites (~70%) before hormone treatment (Supplementary Fig. 6).

Figure 2

AR reprogramming and induced alternative hormonal response

(a)FOXA1 siRNA-induced cell progression to S phase. Relative numbers of propidium iodide (PI)-labeled cells in S-phase at different DH) concentrations were determined by FASCan. The P-value for the difference detected at each hormonal level is indicated; mean±s.e.m. is based on three independent experiments. (b) Top-enriched motifs associated with AR-occupied loci (n=3,115). (c) Comparison between genome-wide AR binding programs before and after FOXA1 knockdown in DHT-treated LNCaP cells. (d) Quantitative levels of AR binding in the “lost”, “conserved”, and “gained” programs. Outliers were omitted from box plots. (e) ChIP-qPCR validated AR binding events on randomly selected loci from the lost (n=22), conserved (n=27), and gained (n=16) programs. (f) Microarray analysis of DHT-induced genes before and after FOXA1 knockdown. (g) Quantitative analysis of gained androgen up-regulated genes based on microarray analysis in f. Outliers were omitted from box plots. (h) Genomic distance of androgen-responsive genes from TSS to the nearest AR binding site in the original and gained AR binding programs.

The conundrum is that, while FoxA1 is known to facilitate AR binding on several DHT-responsive genes9, FOXA1 RNAi actually markedly elevated, rather than diminished, the DHT response (Fig. 2a). We found that ~60% of the original AR binding events was “expectedly” lost in response to FOXA1 RNAi, which we refer to as the “lost” AR program (Fig. 2c,d). We refer the remaining ~40% of AR binding events to as the “conserved” AR program, which often exhibited enhanced AR binding. Strikingly, we detected a massive gain of 10,869 new AR binding loci, referred to as the “gained” AR program (Fig. 2c,d). We extensively validated each of these AR programs by conventional ChIP-qPCR (Fig. 2e). This induced AR reprogramming appears to be qualitatively and quantitatively distinct from reported AR re-targeting on androgen-resistant LNCaP-abl cells compared to parental LNCaP cells17 and is also in sharp contrast to FoxA1-dependent genomic targeting of the estrogen receptor α (ERα) in breast cancer MCF7 cells18. In concert with such massive AR reprogramming, we observed corresponding changes in gene expression in each of three AR programs (Fig. 2f,g, Supplementary Fig. 7). The newly induced AR expression program is also linked to AR binding events (Fig. 2h), suggesting a direct gain-of-function on DHT-responsive genes, as illustrated on SOX9 and other genes (Supplementary Fig. 8), which have been previously documented to play critical roles in cancer progression19,20. Because we also observed ~3-fold elevation of AR expression in FOXA1 RNAi-treated cells (Supplementary Fig. 9a), we tested the possibility that increased AR expression might trigger these effects, finding that AR overexpression alone was insufficient to induce AR reprogramming (Supplementary Fig. 9b).

To explore the mechanism for AR reprogramming, we determined FoxA1 binding on different AR programs, finding that the gained AR program is largely devoid of FoxA1, while FoxA1 is present in more than half of the lost and conserved AR programs (Supplementary Fig. 10). This raises the possibility that FoxA1 may facilitate AR binding to its original binding program, but trans-repress AR from binding to other genomic regions that lack FoxA1 binding sites in the gained program, a strategy frequently used by other transcription activators21. Indeed, as previously reported22, FoxA1 overexpression squelched AR element (ARE)-driven transcription in transfected HEK293 cells (Supplementary Fig. 11a), which is consistent with the ability of AR to directly interact with FoxA123. This mechanism appears to be exploited during tumor progression because an AR mutation identified in advanced prostate tumors lacks part of the hinge domain important for interactions with FoxA123, its ability to interact with FoxA1, and became resistant to FoxA1-mediated trans-repression (Supplementary Fig. 11b,c). Furthermore, our functional analysis indicates that the missing AR ligand-binding domain (LBD) also contributes to AR:FoxA1 interactions (Supplementary Fig. 12). Interestingly, similar AR truncations have also been reported to result from alternative splicing, gene rearrangement, and/or calpain-mediated cleavage (Supplementary Fig. 13). Based on these findings, we propose that FoxA1 regulates AR genomic targeting by simultaneously anchoring AR to cognate loci and restricting AR from other ARE-containing loci in the human genome.

To understand how reprogrammed AR binding is translated to altered hormonal response, we took advantage of the recently established global run-on strategy (GRO-seq6) to detect the functional relationship between AR binding and hormone-induced gene expression. This powerful genome-wide interrogation of on-going transcription detected a broad scope of nascent RNAs. We uncovered 28,318 transcripts with 15,656 annotated and 12,662 unannotated transcripts, among which 450 coding and 347 unannotated transcripts were induced >1.5-fold with even just 1 hr DHT treatment (Supplementary Fig. 14). The TSS of GRO-seq defined transcripts are typically marked by H3K4me3 and H3K27Ac (Supplementary Fig.15a,b). Importantly, GRO-seq also detected non-coding RNAs from a subset of H3K4me1/2 positive, H3K4me3 negative regions (Supplementary Fig. 15c). As illustrated on the enhancer of the KLK3 transcription unit (Fig. 3a), these enhancer-derived RNAs (or eRNAs)5 are largely symmetrical and bidirectional (see additional examples on other well-known hormone regulated genes, such as PMEPA1 and KLK2 in Supplementary Fig. 16). Interestingly, we often detected a large amount of nascent RNA before DHT treatment, particularly near their TSS (e.g. KLK3); DHT not only enhanced the expression of these nascent RNAs, but also permitted the extension of transcription toward the end of the gene (Fig. 3a, Supplementary Fig. 16). We estimated that ~79% of the transcription units induced by liganded AR are regulated at the level of transcriptional initiation, while ~21% appear to be primarily regulated at the level of elongation (Supplementary Fig. 17).

Figure 3

Transcriptional response on individual enhancer programs to FOXA1 down-regulation

(a) Display of nascent RNA detected by GRO-seq on the KLK3 locus. The DHT-induced AR binding is shown at bottom as a reference. (b,c) Induction of eRNA by DHT (b) or FOXA1 knockdown in DHT-treated LNCaP cells (c). The eRNA levels under different conditions (indicated at bottom) are separately displayed on three AR binding programs. (d,e) Effects of FoxA1 on binding of p300 (d) and Med12 (e) in each AR program in DHT-treated LNCaP cells. (f, g) Long-distance interaction between gene promoter and AR bound site (at ~50Kb) was determined by the 3C assay on two representative gene loci selected from the conserved and gained AR programs. Negative controls at shorter distances (~30 and ~40Kb) and a positive control with the corresponding BAC in the region are included in each case.

The ability to detect regulated eRNA expression permitted us to analyze different AR programs during transcriptional reprogramming. In the presence of FoxA1, DHT enhanced eRNA expression from AR-bound enhancers in both the lost and conserved AR programs. In contrast, a basal level of eRNAs was detectable on the gained program, but is independent of the hormone treatment, indicating that these are pre-established enhancers (Fig. 3b). In response to FOXA1 RNAi, the expression of eRNAs was diminished from the lost program, but modestly or dramatically enhanced from the conserved and gained programs, respectively (Fig. 3c). The DHT-induced nascent transcripts (detected by GRO-seq) and steady-state RNAs (detected by microarrays) best predicts direct target genes by liganded AR, as they show the shortest distance (<50kb) to nearby AR binding sites compared to genes identified by either criterion alone (Supplementary Fig. 18), indicating that AR-activated enhancers marked by increased eRNA are responsible for activation of nearby coding transcription units.

In concert with differential eRNA expression, we also observed corresponding changes in levels of another mark in the final step of enhancer activation4, specifically p300, on both conserved and gained AR programs (Fig. 3d). Interestingly, enhancers in the lost AR program continued to exhibit significant p300 binding even after AR binding and eRNA expression were diminished in FOXA1 knockdown cells (Fig. 3c,d). The transcription mediator Med12 has recently been suggested to mediate enhancer-promoter looping24. We tested Med12 binding on individual AR programs, finding that it exhibited an identical binding pattern to p300 (Fig. 3e). Enhanced Med12 binding on the conserved and gained programs after FOXA1 knockdown suggests elevated or newly activated enhancer:promoter interactions. This was demonstrated by the 3C assay on two representative genes where FOXA1 knockdown either enhanced (on the FASN locus from the conserved AR program) or create new (on the NDRG1 locus in the gained AR program) long-range interactions between AR-bound enhancers and specific gene promoters in DHT-treated cells (Fig. 3f,g and Supplementary Fig.19). These data strongly suggest that the induction of eRNAs, rather than binding of either p300 or Med12, is the most precise mark of the final, functional looping between an activated enhancer and its regulated gene promoter.

Addressing the structural basis for different functional classes of AR enhancers, we note that the distinct profiles of H3K4me1 and H3K27ac on the lost, conserved, and gained AR programs and FOXA1 RNAi had little effect on these profiles (Fig. 4a,b and Supplementary Fig. 20). The histone marks H3K4me1 and H3K27ac around the lost and conserved AR programs exhibit a bimodal distribution, which is particularly pronounced on the lost program (Fig. 4a, bottom panel). The DNA binding sites in the lost AR program are actually significantly less enriched in canonical AREs, which may render AR binding on these sites particularly dependent on FoxA1, whereas both the conserved and gained AR programs are associated with nearly perfect palindromic, canonical AREs (Supplementary Fig. 21), explaining why AR is able to target those sites in a FoxA1-independent manner. Strikingly, the gained AR binding sites are coincident with sharp H3K4me1 and H3K27ac peaks (Fig. 4a,b, middle panels), suggesting a distinct nucleosome architecture underlying the gained AR program.

Figure 4

Distinct classes of AR enhancers in the human genome

(a,b) Profiles of H3K4me1 (a) and H3K27ac (b) associated with the lost (bottom panels), conserved (top panels), and gained (middle panels) AR programs in DHT-treated LNCaP cells in response to FOXA1 knockdown. (c, d) Profiles of H3K4me2 around AR binding loci at the nucleosomal resolution in response to DHT stimulation in control siRNA-treated (c) or FOXA1 siRNA-treated (d) LNCaP cells. (e) Profiles of the histone variant H2A.Z on the three different AR programs. (f) Model for FoxA1-mediated AR targeting and reprogramming in LNCaP cells. In Class I (the lost AR program), FoxA1 licenses liganded AR to bind to ARE in relatively nucleosome-free regions. AR binding does not induce nucleosome remodeling in this class of enhancers. In Class II (the conserved AR program), AR binds independently of FoxA1 to ARE, inducing nucleosome remodeling. In Class III (the gained AR program), FoxA1 restricts AR binding, despite the presence of strong AREs. Although pre-established, these gained loci exhibit a strong central nucleosome and are associated with H2AZ, which is not affected by AR binding. FOXA1 knockdown converted these sites to androgen-responsive sites. In all these three classes, eRNAs were generated or increased after AR binding.

A recent study suggests that AR binding leads to dynamic dismissal of a central, H2A.Z-containing nucleosome, replacing by two flanking H3K4me2-marked nucleosomes25. We found that the lost AR program was largely devoid of a “central” nucleosome even prior to AR binding (Fig. 4c, bottom panel). The conserved AR program exhibited DHT-induced switch from the central H3K4me2-marked nucleosome to two flanking H3K4me2-marked nucleosomes, which is largely independent of FoxA1 (Fig. 4c,d top panels). The gained program showed a strong H3K4me2-marked central nucleosome both before and after AR binding (Fig. 4c,d, middle panel). Thus, this gained AR program represents a new type of enhancer topography that requires no nucleosome remodeling for enhancer recognition and subsequent enhancer-promoter interactions. H2A.Z is prevalently associated with the gained AR program, modestly with the conserved AR program, and absent in the lost AR program (Fig. 4e). Together, these findings establish distinct chromatin structures underlying functionally distinct classes of AR enhancers.

In summary, our findings infer a general principle for establishing cell type-specific transcription programs. Cell lineage-specific factors (such as FoxA1) coupled with other general transcriptional factors “create” a cell type-specific enhancer network, allowing other regulated factors (such as AR) to “activate” these pre-established enhancers (Fig. 4f). The enhancer activation process is tightly linked to eRNA production, which appear to serve as a more robust indicator of enhancer activities than any enhancer-bound transcription activators or chromatin marks. On the current biology model, AR reprogramming dramatically altered the androgen responsive pathway, which, according to GO analysis (Supplementary Fig. 22 and Fig.23), may contribute to enhanced cell growth and the establishment of appropriate microenvironment in advanced prostate cancer26-28. Together, these findings provide a conceptual framework to understand complex gene expression switching events, as occurred during disease progression and development.

METHODS SUMMARY

Experiments were performed on LNCaP cells, LNCaP-AR cells (generous gift of Dr. C. Sawyers), and HEK293 cells. Chromatin immunoprecipitations (ChIPs) were done as previously described29 and global run-on (GRO) was performed as described6,30. Control siRNA was purchased from Qiagen (Cat# 1027280). FOXA1 siRNA #1 (M-010319) and #2 (sense 5′-GAGAGAAAAAAUCAACAGC-3′; antisense 5′-GCUGUUGAUUUUUUCUCUC-3′) (ref9) were purchased or synthesized from Dharmacon. Full Methods and associated references are available in the online version of the paper at www.nature.com/nature.