Discovery of active enhancers through bidirectional expression of short transcripts

1Cellular identity and function are defined in large part by regulatory networks that determine gene expression profiles. Control of gene expression is complex, multi-faceted, and coordinated. Over the past decade, with the advent of high-throughput genomic technologies, many systems-level biological approaches have been developed to help resolve these complexities, although substantive questions remain. Recent large-scale human genetic studies have revealed that most complex disease-associated variants map to within non-coding genomic regions, providing additional impetus to expand current catalogs of gene regulatory elements and better understand cellular control of gene expression.

2The first step in gene expression is the recruitment to gene promoters of a multi-protein transcription initiation complex, which includes RNA polymerase (RNAP). Once RNAP is stably bound to the template DNA, it becomes transcriptionally engaged, and commences elongation. It was noted over two decades ago that RNAP could pause/stall at promoters, waiting for a specific signal to continue productive transcription. However, this type of regulation of transcriptional elongation was thought to be an atypical phenomenon. Three recent genome-scale approaches, employing high-throughput sequencing technologies, have revealed that promoter-proximal RNAP pausing is widespread, and likely a common mode of gene regulation.

3One of these methods, global nuclear run-on followed by high-throughput sequencing (GRO-seq), provides a density map of transcriptionally engaged RNAP across the genome by purifying, sequencing, and mapping nascent RNAs. When applied to human lung fibroblasts (IMR90), GRO-seq revealed that promoter-proximal pausing is almost always accompanied by short, divergent (anti-sense) transcription; hereafter this signature is referred to as bidirectional expression of short transcripts (BEST). Two independent methods confirmed this signature in both murine embryonic stem cells and HeLa cells, suggesting that BEST is a general feature of RNAP pausing in mammalian tissues.

4The functional consequences of promoter-proximal RNAP pausing are likely diverse. One recent study found that paused RNAP facilitates the induction and maintenance of an open chromatin conformation near a gene promoter. We reasoned that RNAP pausing, and thus BEST, may occur at specific non-promoter regions where open chromatin is present. Indeed, studies in yeast have shown that pausing also occurs at certain non-promoter sites within gene bodies.

5In this study, we first sought to assess whether BEST is a common feature in human cells. We analyzed published IMR90 GRO-seq data to define actively transcribed genes and used the GRO-seq data at the promoters of these genes for supervised training of a Naïve Bayes classifier (NBC). Using the NBC, we predicted nearly ten thousand high-confidence, non-promoter-associated BEST events genome-wide. Intriguingly, BEST significantly co-occurred with open chromatin loci (DNase hypersensitivity sites (DHSs)), and was even more strongly associated with DHSs that overlap regions enriched for histone H3 lysine 27 acetylation (H3K27ac), the most reliable chromatin marker to date of active enhancers. BEST was modest at regions bound by the CCCTC-binding factor (CTCF), which serve as either direct transcriptional modulators or insulators depending on chromatin context. Further analysis of epigenomic data revealed that several active chromatin marks, including histone H3 lysine 18 acetylation (H3K18ac) and histone H4 lysine 5 acetylation (H4K5ac), but not H3 lysine 4 tri-methylation (H3K4me3), were significantly enriched at non-promoter-associated BEST loci relative to background expectation. Overall, our findings indicate that BEST can demarcate active non-promoter regulatory elements, likely due in part to the specific association of RNAP with distal regulatory elements.

6To confirm BEST at active promoters, we first sought to define from IMR90 GRO-seq data a set of transcriptionally active genes (Materials and methods). We calculated the average reads per kilobase normalized for mapability for every known human gene longer than 3 kb, and performed a receiver operating characteristic (ROC) analysis using 1,522 expressed genes and 2,046 non-expressed genes from an IMR90 microarray expression dataset. The most accurate cutoff for transcriptional activity was determined to be 5 reads/kb/mapability (Additional file 1), yielding 14,145 active RefSeq transcripts, of which 5,213 uniquely mapped to a gene symbol. We computed the average GRO-seq signal along the length of all 5,213 transcripts, and confirmed the previously reported enrichment of sense and antisense reads near the transcription start site (TSS) (Figure 1), which is characteristic of mammalian transcriptional pausing. We also observed a similar, but substantially dampened, pausing signal at the gene end (Figure 1). As shown by Core et al., although the promoter-proximal pausing index was inversely correlated with gene transcription (Additional file 2), pausing was still detectable at the promoters of very highly expressed genes.

7We next sought to use GRO-seq data aligning to the promoter regions of each of the active genes to train a probabilistic model (NBC) for genome-wide prediction of BEST events (Materials and methods). Specifically, for each of three categories of interest (non-transcribed, BEST, transcriptional elongation), we computed the probability distributions for six distinct features (Figure 2). Then for every non-promoter-associated 2-kb window across the entire genome, which we defined as any window at least 7 kb away from a known RefSeq transcription start site or IMR90 H3K4me3 peak (Materials and methods), we used these distributions to calculate the probability that it belongs to each of the three categories, and assigned a category based on highest probability (Materials and methods). Using a logarithm of odds (LOD) score threshold of 2.5, we predicted 9,662 high-confidence non-promoter-associated BEST loci (Figure 3a). The widespread occurrence of non-promoter-associated BEST throughout the genome, particularly within intergenic and inactive intragenic regions, indicates that RNAP specifically associates with non-promoter loci.

8Approximately 70% (n = 6,770/9,662) of the genome-wide BEST predictions overlap IMR90 open chromatin loci (DHSs; Figure 3a), which represents an approximately nine-fold enrichment over background expectation (Materials and methods). Furthermore, almost 85% (n = 5,678/6,770) of these also overlap regions significantly enriched for H3K27ac in IMR90 cells (the most reliable chromatin marker to-date of active enhancers), which represents a striking approximately 30-fold enrichment relative to background (Figure 3a). Manual inspection of several loci confirmed the tendency for non-promoter-associated BEST loci to overlap DHS regions that are enriched for H3K27ac but not for H3K4me3 (Figure 3b-d).

9To further characterize BEST at DHS regions, we compared the profiles of GRO-seq sense/anti-sense read density at non-promoter-associated DHS and non-DHS control sites within transcribed intragenic (Figure 4a), non-transcribed intragenic (Figure 4b), and intergenic loci (Figure 4c) (Materials and methods). For all three categories, we observed a significant accumulation of sense reads within DHSs and anti-sense reads immediately upstream (Figure 4a-c) - precisely the signature of BEST. The signal for BEST was even more pronounced at DHSs that overlap H3K27ac peaks (Figure 4a-c). In intragenic loci, the accumulation of GRO-seq anti-sense reads at DHSs appears more pronounced than GRO-seq sense reads (Figure 4a-b); however, this is most likely because we are normalizing read density in a particular window by the average read density in the entire region, and the average read density is always higher in the sense orientation.

10To analyze the chromatin landscape at the predicted non-promoter-associated BEST loci more comprehensively, we assessed the representation of ten different histone modifications for which IMR90 data from chromatin immunoprecipitation followed by high-throughput sequencing (ChIP-seq) and corresponding control (input) data were available for two samples from the Epigenome Atlas (Materials and methods). The most strongly enriched modification was H3K27ac (Figure 5), which is a robust discriminator between active and poised enhancers. The second-most enriched mark was H3K18ac, which is also thought to be associated with functional enhancers. The two least-enriched modifications were H3K4me3, which is associated primarily with promoters, and histone H3 lysine 4 monomethylation (H3K4me1), which is thought to be present at both active and inactive/poised enhancers. Finally, two modifications associated with repressed states, histone H3 lysine 27 trimethylation (H3K27me3) and histone H3 lysine 9 trimethylation (H3K9me3), were depleted at regions of BEST.

11We repeated the above-described analysis for non-promoter-associated DHSs that overlap H3K27ac peaks (DHS+/H3K27ac+) and for those that do not (DHS+/H3K27ac-) (Figure 5). Expectedly, DHS+/H3K27ac+ regions were most enriched for H3K27ac and other active chromatin marks, and were depleted for repressive marks (Figure 5). DHS+/H3K27ac- regions were depleted for most active and repressive marks, and were enriched only for marks often associated with poised states, H3K4me1 and H3K4me2 (Figure 5). Most importantly, regions of BEST exhibited a chromatin landscape significantly more similar to that of candidate active enhancers (DHS+/H3K27ac+) than poised enhancers (DHS+/H3K27ac-).

12To confirm this finding in another cell type, we turned to mouse embryonic stem cells (mESCs), which is the only other cell type in which both the nascent transcriptome (GRO-seq) and enhancer-related chromatin marks (ChIP-seq) have been extensively characterized. We re-trained the NBC using mESC GRO-seq data aligning to active promoters, and applied the NBC genome-wide to predict non-promoter-associated BEST loci. Using a genome-wide dataset of candidate mESC enhancers, we found that robustly active enhancers are approximately 8.5-fold enriched (P < 0.0001) for BEST relative to poised enhancers (Additional file 3). Furthermore, approximately 71% (n = 5/7) of the candidate enhancers that were validated by an in vitro reporter gene assay were predicted as BEST loci (Materials and methods). Collectively, these results indicate that BEST regions, as predicted by our classifier, are preferentially associated with active enhancer elements located within both transcribed and non-transcribed genomic regions.

13Distal regulatory elements are not limited to enhancers; another important class is target sites for CTCF, which have many known functions, including insulator activity. To assess BEST at IMR90 CTCF binding sites, we followed the same method as for DHSs to assess GRO-seq sense and anti-sense read density profiles (Materials and methods). Relative to non-CTCF control regions, we detected a robust signal for BEST at CTCF binding sites that overlap H3K27ac peaks (CTCF+/H3K27ac+), only a very modest signal at sites that overlap DHS peaks alone (CTCF+/DHS+), and no signal at sites that overlap neither (Figure 6a,b). The results are consistent with the previous finding that CTCF can sometimes recruit RNAP to CTCF binding sites. However, the pronounced signal at CTCF+/H3K27ac+ sites, together with the dampened signal at CTCF+/DHS+ sites, suggests that this recruitment may be most prevalent at CTCF binding sites that function as, or are proximal to, active enhancers.

14In this study, we analyzed the IMR90 nascent transcriptome, and developed a probabilistic model of promoter-proximal transcriptional pausing in order to identify non-promoter-associated BEST. Further computational analysis, using genome-wide IMR90 chromatin profiles (Epigenome Atlas), revealed that non-promoter BEST is significantly associated with regions enriched for chromatin marks (such as DHSs, H3K27ac and H3K18ac) that demarcate active enhancers.

15RNAP pausing is well-appreciated at promoter regions, and very recently has been discovered at specific loci within actively transcribed genes in yeast. RNAP pausing has also been observed at cohesin binding sites within a single actively transcribed human gene in human umbilical vein endothelial cells (HUVECs). However, to the best of our knowledge, our work is the first systematic, genome-scale investigation of non-promoter RNAP pausing in human cells. A recent study identified systematic biases in next-generation sequence data, such that an accumulation of GRO-seq sense reads may not necessarily reflect bona fide pausing, due to various nucleotide preferences during cDNA amplification and sequencing. Our analysis circumvents this issue by defining pausing as an accumulation of both sense and anti-sense GRO-seq reads, reflecting the widespread divergent transcription associated with promoter-proximal RNAP pausing in mammalian cells.

16Recent in silico strategies to identify RNAP pausing from nascent RNA sequencing data have used a local, deterministic approach - a minimum level of enrichment of sense read density in a particular window relative to neighboring windows. In contrast, our approach utilizes a probabilistic model trained on a reliable genome-wide dataset. The model can be trained on, and applied to, any GRO-seq dataset in order to make inferences about the most likely active enhancer elements. Applying the model to IMR90 GRO-seq data, we detected thousands of non-promoter BEST events. Some of these are located within actively transcribed regions; therefore, the BEST could be due to bound transcription factors that hinder RNAP processivity and induce pausing. However, many of the BEST events are located in intergenic regions. BEST at these loci may be due to specific recruitment of RNAP to active enhancers, as reported previously. In fact, a seminal RNA-seq-based study reported that many neuronal enhancers recruit RNAP, which then transcribes bi-directionally a novel class of transcripts termed enhancer RNAs (eRNAs). A more recent GRO-seq-based study reported significant and dynamic changes to the cellular eRNA profile upon application of an exogenous trigger. It is quite possible that predictions of BEST events using our approach coincide with regions that produce eRNAs as defined by these two studies.

17Two major unanswered questions are whether RNAP-mediated BEST occurs at all active enhancers, and whether it contributes to the maintenance of open chromatin, or is just a consequence of the presence of other factors at those sites. A critical related question is whether the short bi-directional transcripts produced by BEST at these sites have functional relevance, or are simply transcriptional noise tolerated by evolution because of relatively minor metabolic cost. Either way, the ability to detect BEST by analysis of GRO-seq data contributes another important approach for the dissection of genomic regulation in higher eukaryotes.

18Long-range regulatory elements are important modulators of gene expression, but they remain poorly annotated. Recent approaches for genome-wide identification of regulatory elements have focused on analyzing the chromatin state. This study contributes an alternative, complementary strategy. We developed a probabilistic model to capture the transcriptomic signature, BEST, of promoter-proximal polymerase pausing. We used this model to predict non-promoter-associated BEST regions, which were significantly enriched for chromatin marks (such as H3K27ac) that are associated with active long-range regulatory elements.

19All human RefSeq transcripts (n = 35,983) were downloaded in hg18 coordinates from the UCSC Table Browser, build 36. Only validated mRNA transcripts ('NM_' prefix) were retained (n = 30,326). For each transcript longer than 3 kb (n = 27,863), we defined the body of the transcript as 1 kb downstream of the transcription start site to the annotated gene end. The level of transcription for each transcript body in IMR90 cells was determined by computing the average GRO-seq sense reads/kb/mapability using the previously published IMR90 GRO-seq data and the 'Duke Uniq 35' mapability data downloaded from the UCSC Table Browser, build 36. A reads/kb/mapability cutoff for transcriptional activity was chosen according to the maximal accuracy measure - the reads/kb/mapability that achieves the optimal combination of sensitivity (true positive rate) and specificity (true negative rate) using high-confidence true positive genes (expressed genes; n = 1,522) and true negative genes (non-expressed genes; n = 2,046) from a published IMR90 microarray dataset. Accuracy was measured by calculating the following: (Number of expressed genes identified + Number of non-expressed genes identified)/All genes. At the most accurate cutoff of 5 reads/kb/mapability, 14,145 transcripts were called active, of which 5,213 uniquely mapped to a gene symbol.

20GRO-seq data aligning to (i) active promoters (n = 5,213), (ii) active transcript bodies (n = 5,213), and (iii) randomly selected intergenic regions (with a similar length distribution as active RefSeq transcripts, but non-overlapping with any known RefSeq transcript; n = 5,213) were used to train a NBC to identify BEST, transcriptional elongation, and non-transcribed regions, respectively. First, 2-kb windows were centered at all active start sites of transcription, mid-points of active transcripts, and mid-points of the randomly selected intergenic regions. Then, for all 2-kb test windows in each class, the values for the following six features were calculated: (i) sense strand reads/kb/mapability (hereafter 'read density'), (ii) antisense strand read density, (iii) ratio of the sense strand read density to that in a 2-kb window immediately 3' of the test window, (iv) ratio of the antisense strand read density to that in a 2-kb window immediately 3' of the test window, (v) ratio of the sense strand read density to that in a 2-kb window immediately 5' of the test window, and (vi) ratio of the antisense strand read density to that in a 2-kb window immediately 5' of the test window. The values of these features were used to compute a probability distribution for each feature for each class. These distributions were utilized by the NBC, according to the following equation, to differentiate BEST events from productive transcriptional elongation and transcriptional noise:

21where the three classes B, E, and N represent BEST, transcriptional elongation, and non-transcribed, respectively. The prior probabilities, Pr(ci), were set to be equal for all classes. The six described features are represented by f1 to f6. The NBC was applied genome-wide, on both strands, avoiding regions associated with promoters or gene ends (within 7 kb of known transcription start sites, annotated gene ends, and IMR90 H3K4me3 peaks). On the plus strand, plus strand GRO-seq reads are interpreted as 'sense' and minus strand GRO-seq reads are interpreted as 'antisense'; on the minus strand, plus strand GRO-seq reads are interpreted as 'antisense' and minus-strand GRO-seq reads are interpreted as 'sense'. On each strand, for each non-overlapping 2-kb test window in the search space, a LOD score was computed comparing the probability of BEST with that of the other two classes:

22Test windows with LOD scores > 2.5 on one or both strands were set as high-confidence BEST loci.

23IMR90 DNase-seq read data for four biological replicates were downloaded from the Epigenome Atlas, release 3. MACS version 1.4 was run on each dataset, using the parameter values described previously, to identify genomic regions of enrichment for DNase-seq reads. Regions called as enriched in all four replicates were defined as 'DHS peaks'. IMR90 H3K27ac ChIP-seq read data for two biological replicates, and corresponding control (input) data, were downloaded from the Epigenome Atlas, release 3. MACS version 1.4 was run on each dataset, using the default parameter values, to identify genomic regions enriched for H3K27ac. Regions called as enriched in both replicates were defined as 'H3K27ac peaks'. Finally, IMR90 ChIP-chip-derived CTCF peaks were downloaded from the Ren laboratory website and converted to hg18 coordinates using the command line liftOver program with the -minMatch parameter set to 0.9.

24DHS/CTCF peaks were categorized as located within actively transcribed intragenic regions, inactive intragenic regions, or intergenic regions, with respect to the RefSeq dataset used in this study (see the 'Identifying actively transcribed genes' section of the Materials and methods). To avoid promoter-associated peaks, DHS/CTCF peaks + 5 kb flanking regions that were within 2 kb of known transcription start sites, annotated gene ends, or IMR90 H3K4me3 peaks were discarded. For each of the remaining DHS/CTCF peaks within each category, GRO-seq sense and anti-sense reads/kb/mapability were computed in 150-bp windows from the start of the DHS/CTCF peak to the end of 5-kb flanking regions on either side. Then, for each DHS/CTCF peak and flanking region, nucleotide distance was converted to proportional distance. For example, for a DHS/CTCF peak that is 300 bp in length, the first 150 bp immediately upstream of the peak corresponds to '-0.5 to 0', the first 150 bp within the peak corresponds to '0 to 0.5', the second 150 bp within the peak corresponds to '0.5 to 1', the first 150 bp immediately downstream of the peak corresponds to '1 to 1.5', and so on.

25IMR90 ChIP-seq read data for ten different histone modifications, each with at least two biological replicates, and corresponding control (input) data, were downloaded from the Epigenome Atlas, release 3. For each histone modification dataset, the read density (reads/bp) was computed at predicted, high-confidence BEST loci, and then divided by the read density at randomly generated background (control) regions (2 kb in length and drawn from the same genomic locations as BEST loci), to yield an enrichment value. The enrichment value was then divided by the enrichment value for input, to yield a normalized enrichment value.

26To perform genome-wide prediction of BEST loci in an additional cell type, the NBC was trained and applied on publicly available mESC GRO-seq data in exactly the same manner as was done using GRO-seq data from IMR90 cells. Genome-wide candidate mESC enhancers (poised, weak, and strong) were downloaded from Zentner et al. and in vitro validated mESC enhancers were downloaded from Schnetz et al.. In both cases, only those not within 7 kb of known transcription start sites, annotated gene ends, and mESC H3K4me3 peaks were retained for further analysis.

27BEST: bidirectional expression of short transcripts; ChIP-seq: chromatin immunoprecipitation followed by high-throughput sequencing; CTCF: CCCTC binding factor; DHS: DNase hypersensitive site; eRNA: enhancer RNA; GRO-seq: global nuclear run-on assay followed by high-throughput sequencing; H3K18ac: histone H3 lysine 18 acetylation; H3K27ac: histone H3 lysine 27 acetylation; H3K4me1: histone H3 lysine 4 mono-methylation; H3K4me3: histone H3 lysine 4 tri-methylation; IMR90: human lung fibroblasts; LOD: logarithm of odds; mESC: mouse embryonic stem cell; NBC: Naïve Bayes classifier; RNAP: RNA polymerase.

28The authors declare that they have no competing interests.

29PS conceived of, designed, and co-supervised the study, carried out computational analyses, and wrote the manuscript. MFM carried out computational analyses and participated in manuscript preparation. FSC co-supervised the study and edited the manuscript. All authors read and approved the final manuscript.