Novel mRNAs 3′ end-associated cis-regulatory elements with epigenomic signatures of mammalian enhancers in the Arabidopsis genome

Genome-wide identification of putative enhancer-like elements in Arabidopsis

Transcriptional enhancers activate gene expression independent of orientation and distance. Due to the scarcity of knowledge on plant transcriptional enhancers, we based our criteria on (i) information known about metazoan enhancers and (ii) our analysis on the well experimentally characterized EASE enhancer in Arabidopsis. We then used these criteria to isolate genomic regions that carry these canonical chromatin signatures from data sets produced from Arabidopsis seedling of similar ages and conditions (schematic in Fig. 1B and details in SI Text). We first examined the sites of nucleosome depletion in the Arabidopsis genome. To assess the nucleosome-free regions, 38,290 DNase I hypersensitive sites (DHSs) were extracted from DNase-seq data set (Zhang et al. 2012) for young Arabidopsis seedlings (Fig. 1B, step 1).

Open chromatin regions are also marked by specific histone modifications, including H3K4 methylation. Monomethylation on H3K4 generally marks bona fide canonical enhancers in metazoans, although different chromatin signatures exist in metazoans; additionally, poised enhancers are marked by H3K4me1 before they acquire the proper H3K27 marks (Heintzman et al. 2007; Ernst et al. 2011; Zentner et al. 2011; Bonn et al. 2012; Djebali et al. 2012; Calo and Wysocka 2013; Shlyueva et al. 2014; Deaton et al. 2016). Recent data in Arabidopsis, pea, and maize showed a positive correlation between active enhancers and H3/H4 acetylation, and a correlation between inactive enhancers and H3K27me3, in agreement with the chromatin signatures of metazoan enhancers, suggesting that some plant enhancers may exhibit chromatin features similar to those of animals (Weber et al. 2016). Additionally, using data sets from chromatin immunoprecipitation-microarray (ChIP-chip) experiments on Arabidopsis seedlings (Zhang et al. 2007; Charron et al. 2009; Chodavarapu et al. 2010), we determined that the experimentally characterized Arabidopsis EASE (Yang et al. 2005) and its neighboring chromatin are enriched for H3K4me1 compared to H3K4me3 (Fig. 5A; P < 0.05; see below for more for details). This result further indicates that H3K4me1 is informative for the identification of putative enhancers in Arabidopsis. In other species, more than 50% of potential enhancers were predicted to be outside of intergenic regions (Heintzman et al. 2007; Kim et al. 2010; ENCODE Project Consortium et al. 2012; Kowalczyk et al. 2012; Arnold et al. 2013). Since H3K4me1 is found predominantly in the gene body of less than one-third of all protein-coding genes in Arabidopsis (Supplemental Fig. S1A; Zhang et al. 2009), the use of H3K4me1 enrichment as an enhancer signature also allows to examine all DHSs genome-wide. Therefore, we next set out to isolate DHSs that colocalize with H3K4me1 enrichment, using H3K4me1 ChIP-chip data sets from the high-density whole-genome tiling microarray for young Arabidopsis seedlings (Zhang et al. 2009). We found that 6001 (P < 0.0001, Bedtools Fisher's exact test was used to test the statistical significance of data associations with a set of randomized control) DHSs colocalized with areas enriched for H3K4me1 (Fig. 1B, step 2).

Since the vast majority of TF binding sites are located in the vicinity of promoters, we had to discriminate against promoters, which are typically marked with H3K4me3 (as observed in metazoans and Arabidopsis) (Zhang et al. 2009; Spitz and Furlong 2012; Calo and Wysocka 2013; Shlyueva et al. 2014). To discriminate between enhancers and promoters, we set a threshold for the ratio of mono- and trimethylation of H3K4 of >2 (Fig. 1B,C; HMM posterior probability scores represent the signal intensities per specific histone modification) using H3K4me1 and H3K4me3 ChIP-chip data sets for young Arabidopsis seedlings (Zhang et al. 2009). This analysis yielded 1919 regions with the canonical signatures of mammalian enhancers and EASE (Fig. 1B, step 3). These 1919 regions, which varied in size, could harbor the candidates for PEs in the Arabidopsis genome. We also found that 65% of these regions (n = 1249, P < 0.0001) overlap with the nucleosome-free regions identified by ATAC-seq in root tissues by Maher et al. (2018).

To uniformly analyze the 1919 regions, further explore their genomic features, and consider the nucleosomes up- or downstream from the DHS, we first added 1.5-kb up- and downstream from the DHS midpoints converting each DHS into a 3-kb region. We termed these 3-kb regions PEs and the center of the DHS the PE midpoint (Fig. 1; Supplemental Dataset S1).

PEs are highly H3K27 acetylated, bound by Pol II, and transcribed

Metagene analysis was used to profile the distribution of chromatin modifications across PE regions. Analysis of mono-, di-, and trimethylated H3K4 (Zhang et al. 2009) in the PE regions revealed patterns very similar to those reported for mammalian enhancers (Heintzman et al. 2007; Kim et al. 2010; Calo and Wysocka 2013; Shlyueva et al. 2014): a significant enrichment of H3K4me1, very low H3K4me3, and no fluctuations in H3K4me2 throughout the PE region (Fig. 1D; Supplemental Fig. S2). Notably, examination of H3K27 (Zhang et al. 2007; Charron et al. 2009), which can distinguish between active and poised enhancers (Calo and Wysocka 2013; Shlyueva et al. 2014), revealed that PEs have high H3K27ac and little H3K27me3 (Fig. 1D; Supplemental Fig. S2), suggesting that PEs could be active enhancers.

Recent studies showed that active enhancers are transcribed in mammals (Kim et al. 2010; Li et al. 2016). In accord, using a Pol II ChIP-chip data set for young Arabidopsis seedlings (Chodavarapu et al. 2010), we found widespread Pol II occupancy throughout the entire PE region (Fig. 1D; Supplemental Fig. S2), peaking around the PE midpoints. Given that plants have two additional plant-specific RNA polymerases, Pol IV and V, we confirmed that very few PE regions overlap with Pol IV/Pol V-associated loci, suggesting that Pol IV/V are not responsible in transcribing PEs (Wierzbicki et al. 2012; Law et al. 2013) (details in SI Text). All together, the epigenomic profiles of PEs replicate the canonical signatures of mammalian enhancers and a well experimentally characterized Arabidopsis enhancer EASE, and also indicated that PEs may represent active enhancer-like elements in Arabidopsis.

To examine transcriptional activity in PE regions, we profiled RNA expression in WT seedlings. To determine whether transcripts in PE regions are exosome-sensitive, we compared the RNA expression profile in WT with the RNA expression profile in inducible RNAi lines of the RNA exosome core complex subunit RRP41 (rrp41-i). Since eRNAs can vary in length and polyadenylation status and are transcribed uni- or bidirectionally (Li et al. 2016), we sequenced polyadenylated RNAs [poly(A)⁺] and ribo-minus RNAs depleted of rRNAs to increase our detection of nonpolyadenylated RNAs (Supplemental Dataset S2). Our metagene analysis of the regions ±1.5 kb from the PE midpoint in all RNA-seq data sets found that 99.8% (P < 0.0001, χ² test) of PEs are transcribed in both directions (Fig. 1E,F; Supplemental Fig. S2). The transcriptional profile of poly(A)⁺ and ribo-minus RNA around PE midpoints appears very different from typical bidirectional transcription from TSSs of metazoan protein-coding genes (Andersson et al. 2015). Instead, the peaks of transcriptional activity appear upstream of PEs on both DNA strands and end near the PE midpoint (Fig. 1F).

In examining the effect of the RNA exosome, we found small but discernable differences in RNA abundance in the rrp41-i mutant lines compared to the WT (Fig. 1G), and this difference was particularly pronounced for the ribo-minus RNAs (Fig. 1G, right panel). Of the PEs, 207 express polyadenylated RNAs affected by depletion of the exosome subunit RRP41, with 67 PEs exclusively expressed only in rrp41-i (Fig. 1E; Supplemental Fig. S3; Supplemental Dataset S3; details in SI Text).

PEs are hotspots for TF binding and binding sites specifically cluster around the PE midpoints

Enhancers serve as cis-regulatory modules providing a platform for TF binding (Spitz and Furlong 2012; Core et al. 2014). To examine the occurrence of TF binding motifs in PEs, we extracted the binding motifs of >1700 experimentally characterized TFs and >50 TF families from the Arabidopsis regulatory information server (AGRIS) (Yilmaz et al. 2011). The presence of these known TF motifs were scanned using find individual motif occurrences (FIMO) (Grant et al. 2011) on both strands of each PE region, and only the statistically significant motif occurrences were subjected to further analysis (statistical threshold, P < 0.0001). Metagene analysis revealed significant enrichment of TF binding motifs in PE regions preferentially clustered in the vicinity of PE midpoints (Fig. 2A), as would be expected for transcriptional enhancers. Based on scanning for the presence of their motifs with a statistical threshold of P < 0.0001, MADs-box, MYB, HSF, HB, and bZIP TFs were the most prevalent families of TFs with binding sites in the PE regions (Fig. 2B), suggesting that PEs could be bound by a wide range of TFs.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 2.

PEs are hotspots for TF binding and binding sites specifically cluster around PE midpoints. Panels (A–D) show the distribution of known TF binding motifs in PE regions. The experimentally characterized TF motifs were extracted from the AGRIS database and their occurrences in PE regions were determined with MEME-FIMO (P < 0.0001). Panels (E–G) show the experimentally identified binding of TFs at PE regions. The peaks of TF binding at PE regions were extracted from DAP-seq data sets (A) Distribution of all TF motifs in identified PEs. (B) Abundance of different TF families with binding motifs at PEs, shown as a total number of binding motifs present in the PE region for each TF family. (C) Distribution of specific TF motifs bearing PAS-like (in red) and TATA-box like (in green) sequences in PE regions. (D) The top 20 most frequent TF binding motifs identified in the PE regions; their distribution is shown in Supplemental Figure S4. (E) Distribution of TF binding sites in the identified PEs. The peak enrichment of TF binding was observed close to the PE midpoints (orange dotted line), defined as PE midpoint-proximal region, ±500 nt flanking the PE midpoint (in yellow). (F) Numbers of TFs binding in PE midpoint-proximal regions (***) P < 0.0001. The five specific TFs represented by the red bar bind at a total of >600 PEs (30%–40% of total PEs); and (G) their information (P < 0.0001, Bedtools Fisher's exact test).

Surprisingly, the most enriched motifs near PE midpoints included polyadenylation sequence (PAS)-like motifs (Fig. 2C,D). Moreover, most of the TF motifs are A-rich, closely resembling PASs and the sequences located at the 3′ ends of protein-coding genes (Fig. 2D; Supplemental Fig. S4). Additionally, our analysis of TF binding motif occurrences showed that the PAS-like motifs (AAAATAAA and AAATAAAA) present in the majority of PE regions (76%; P < 0.0001, χ² test) are bound by the Arabidopsis TF BELLRINGER (BLR), which belongs to the conserved eukaryotic Homeobox TF family. BLR has been shown to function as repressor of AGAMOUS (AG) by binding to AG CREs located in a 3-kb intronic region (Bao et al. 2004). A total of 1460 PEs have at least one BLR-binding-PAS-like motif (76%; P < 0.0001, χ² test) and the majority of these PEs harbor more than one motif (n = 1083). The homeobox TF family includes a wide range of well-conserved eukaryotic TFs; specifically, one specific homeobox TF, NANOG, was shown to bind and activate primed enhancers in mice to promote the pluripotent state in mouse embryonic stem cells (Murakami et al. 2016). The Uniport sequence search was used to examine whether BLR shares sequence similarity to any proteins involved in enhancer binding. BLR shared significant, high sequence similarity with the homeobox protein Meis1 in human (e = 3 × 10⁻¹³) and mice (e = 2 × 10⁻¹³), and Pknox1 in human (e = 4 × 10⁻¹³) and mice (e = 3 × 10⁻¹³). Pknox1 and Meis I bind to the Hoxb2 enhancer and mediate Hoxb2 hindbrain enhancer activity (Jacobs et al. 1999). Additionally, Meis1 is a homolog of homothorax (hth) in Drosophila, which (in complex with other TFs) can function as a molecular switch in regulating rho CREs (enhancers) that are essential for Drosophila peripheral nervous system (PNS) development (Li-Kroeger et al. 2008).

To identify how TF binding sites clustered in PEs, we extracted peaks of TF binding in PE regions from DNA affinity purification sequencing (O'Malley et al. 2016) (DAP-seq) data sets (P < 0.0001, Fisher's exact test was used to test the statistical significance of data associations with a set of randomized control) and subjected them to metagene analysis (Supplemental Datasets S4, S5). The Cistrome project produced a genome-wide catalog of the binding peaks for 1812 TFs (composing 80 TF families) in Arabidopsis using DAP-seq (O'Malley et al. 2016), and 529 (30%) of these Arabidopsis TFs were further analyzed for their binding peaks and motifs. Similar to the results from mapping of TF motifs (Fig. 2A), we found TF binding peaked in the region close to the PE midpoints (Fig. 2E).

To determine which specific groups of TFs bind near PE midpoints, we restricted the analysis window to the PE midpoint-proximal region (yellow area in Fig. 2E), ±500 nt of the PE midpoints. We then analyzed all TF binding peaks determined by DAP-seq in the midpoint-proximal region of 1919 PEs (Supplemental Dataset S6). The TFs with binding sites within ±500 nt of PE midpoints can be subdivided into two groups, depending on the number of PE regions they bound (P < 0.0001, χ² test). The first group included a wide variety of TFs, but each TF bound to only a few PEs (Fig. 2F, green bars). For example, 391 different TFs bind to <10% of all PEs (Fig. 2F, green bar at far right). The second group comprised only five TFs (Fig. 2F, red bar; Fig. 2G), but these TFs bind a large number of >600 PEs (30%–40% of total PEs).

TF binding motifs in the midpoint-proximal regions are distinct from general consensus binding motifs

We further analyzed the five specific TFs that bind the most PEs (Fig. 2G). A study that analyzed the TF motifs in ATAC-seq also implicated one of these five TFs, MYB77, in coregulation of common gene sets in root tips across four different plant species (Maher et al. 2018). We first examined the binding profile of these five TFs and the distribution of their binding sites over the PE regions individually (Supplemental Fig. S5A–E, panel i; Supplemental Dataset S7). We found that binding sites for these TFs (with the exception of the G2-like TF) are highly enriched at the midpoint-proximal regions, with a sharp binding peak around the midpoint. We then extracted the binding sites for each of these five TFs from the entire 3-kb PE region and the 1-kb midpoint-proximal regions and subjected them to motif discovery analysis separately to identify the three most enriched binding motifs in these regions (E-value ≤0.05, MEME-ChIP [Machanick and Bailey 2011]). The motifs determined from TF binding sites across the entire 3-kb region closely resembled the consensus motifs determined by interrogating all peaks of TF binding genome-wide by O'Malley et al. (2016) (Supplemental Fig. S6). In contrast, the motifs from TF binding sites within midpoint-proximal regions are extremely AT-rich and tend to cluster preferentially around PE midpoints (Supplemental Fig. S5A–E, panel ii; Supplemental Dataset S7). The difference between the motifs from the PE midpoint-proximal regions and the general consensus sequences suggests that midpoint-proximal regions might be bound by TFs through a distinct set of binding motifs.

Spatial relationship between PEs and the TSSs and TESs of protein-coding genes

Metazoan enhancers occur in inter- or intragenic regions across the genome (Heintzman et al. 2007; Kim et al. 2010; ENCODE Project Consortium et al. 2012; Kowalczyk et al. 2012; Arnold et al. 2013). To characterize the genomic locations of the PEs, we mapped the PE regions on all Arabidopsis chromosomes and found that PEs are distributed across euchromatic regions on all chromosomes, but are notably absent from centromeric regions (Supplemental Fig. S7). Examination of PE positions relative to annotated genes showed that nearly all PEs (99%, P < 0.0001, Fisher's exact test was used to test the statistical significance of data associations with a set of randomized control) overlap with protein-coding gene units and are thus intragenic PEs (Fig. 3A; Supplemental Datasets S8, S9; P < 0.0001, χ² test). Only 21 PEs (1%) were located between protein-coding genes. Therefore, in contrast to the previously identified intergenic enhancer candidates (Zhu et al. 2015; Hetzel et al. 2016), 99% of the putative enhancers-like elements (PEs) identified here are intragenic.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 3.

Genomic architecture of the identified PEs and the spatial relationship between PEs and the TSSs and transcription end sites (TESs). (A) Nearly all PEs are intragenic, overlapping with protein-coding genes. Only 21 PEs (1%) are located in intergenic regions. (***) P < 0.0001. (B) Pie charts show the proportion of genomic features (protein-coding genes and non-protein-coding features) overlapping with intragenic PEs (n = 1898) and intergenic PEs (n = 21) (left and right, respectively). (C) Gene Ontology (GO) analysis of the protein-coding genes overlapping exclusively with intragenic PEs (n = 1348) in the molecular function category (P < 0.0001; Supplemental Dataset S10). (D) Analysis of distances between PE midpoints (n = 1348) and TSSs of all 2450 protein-coding genes overlapping with PEs, conducted in a strand-specific manner. Regions beyond the PE borders extended to ±10 kb from the PE midpoints were analyzed. These protein-coding genes were sub-divided into two groups based on where their TSSs are located relative to the PE midpoints, either upstream (i) or downstream (ii). Upstream genes: TSSs of genes located 5′ to PE midpoint; downstream genes: TSSs of genes located 3′ to PE midpoint. The bars represent the total number of TSSs in 1-kb intervals from PE midpoints. (E) Analysis of intragenic PEs based on the numbers of protein-coding genes overlapping with them. (F,G) Heatmaps illustrating protein-coding gene structures relative to PE midpoints for Group 1 and 2 PEs, respectively. A red-to-blue color scale represents the protein-coding gene structures, with red referring to the 5′ end (TSSs) and blue referring to the 3′ end (TESs) of protein-coding genes. The three panels in figure C and D differ by the region examined, ranging from ±1500 nt, ±500 nt, and ±100 nt from the PE midpoint (from left to right). (F) Group 1 PEs overlap with one gene; (G) Group 2 PEs overlap with two protein-coding genes per PE region. For both Group 1 and 2 PEs, TESs (blue) are located preferentially in close proximity to the PE midpoint in contrast to TSSs (red). The same trend is found for TESs of protein-coding genes encoded on both DNA strands.

Most of the intragenic PEs (n = 1348) identified here overlap exclusively with protein-coding genes (Fig. 3B; Supplemental Datasets S8, S9). GO analysis of the 2450 protein-coding genes overlapping these PEs found that regulation of transcription is among the most frequent GO terms in the molecular function category (Fig. 3C; Supplemental Fig. S8; Supplemental Dataset S10; P < 0.0001, Fisher's exact test). The remaining 550 intragenic PEs overlap with protein-coding genes and non-protein-coding genes, including transposons (TEs) (Fig. 3B). The intergenic PEs (n = 21) located between protein-coding genes can also overlap with non-protein-coding features, including TEs (Fig. 3B). The integrative genomics view (IGV) visualization of three intragenic PEs and their overlapping protein-coding genes is shown in Figure 4.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 4.

Integrative genomics viewer (IGV) visualization of three PE loci and their surrounding protein-coding gene annotations. The PE region (indicated by orange boxes) and the protein-coding genes (indicated by black boxes with the direction of transcription) in the same genome region are shown. The IGV genome-browser view of H3K4 methylation (mono-, di-, and tri-), acetylation and trimethylation of H3K27, and RNA Pol II occupancy are shown for PE_000087 (A), PE_001121 (B), and PE_000107 (C). Each genomic-browser view is centered around the PE midpoint (indicated by orange dotted lines).

In the compact genome of Arabidopsis, gene bodies and intergenic regions occur about every 2–3 kb (The Arabidopsis Genome Initiative 2000), indicating that CREs, like enhancers, in Arabidopsis are likely to be located close to protein-coding genes. To understand how the intragenic PEs overlap with specific protein-coding gene structures, we first used a metagene-like analysis and mapped the TSSs of all protein-coding genes that overlap with PE regions, in a strand-specific manner. To account for the TSSs of all genes overlapping with PEs, we expanded the window of analysis to ±10 kb from the PE midpoint. To consider the distance between PE midpoints and the TSSs of all protein-coding genes that overlap with PE regions, all genes overlapping with PEs were grouped based on the distance of their TSSs to PE midpoints (shown as bars on Fig. 3D). For the genes that initiate downstream from PE midpoints (top right and bottom left on Fig. 3D, panel ii), the average distance of their TSSs from the midpoint was 454–458 nt, indicating that their promoters are located away from PE midpoints. In contrast, the average distance of TSSs located upstream of PE midpoints was >3.2 kb (top left and bottom right on Fig. 3D, panel i). Notably, twice as many genes overlapping PEs have TSSs located far upstream of the PE midpoint (n = 1647) compared to downstream from the midpoints (n = 803) (P < 0.0001, χ² test; Fig. 3D; Supplemental Fig. S9). This suggests that most TSSs of protein-coding genes are not preferentially located near PE midpoints. Additionally, based on chromatin modification signatures (shown in Fig. 1D), PEs do not resemble conventional promoters, when compared to conventional promoters previously characterized by Hetzel et al. (2016).

We then analyzed the connection between transcription end sites (TESs) of protein-coding genes and the PE midpoints. To visualize where different protein-coding gene structures are enriched relative to the PE midpoints, we depicted the positions of the TSSs and TESs in different colors using a heatmap-like approach. Since some intragenic PEs overlap with more than one protein-coding gene, we subdivided the PEs into groups, based on the number of genes they overlap and analyzed them separately (Fig. 3E). Most PE regions (93%) overlap with one or two protein-coding genes (Groups 1 and 2, respectively). We mapped the positions of protein-coding genes relative to gene structures and the PE midpoint for the Group 1 and 2 PEs separately (Fig. 3F,G), taking into account the orientation of each gene structure. Three different ranges of distance from PE midpoints were analyzed, ±1500, ±500, and ±100 nt, to zoom in on the regions closest to the midpoint, which most likely represent the regulatory modules. This analysis indicates that the midpoint-proximal regions of PEs are dominated by the TESs of protein-coding genes.

To confirm that there is no bias between the chromatin signature present at PEs and the 3′ ends of all other protein-coding genes that do not harbor PEs, we further analyzed the 3′ ends of the group of all other protein-coding genes that do not house PEs. We profiled the distribution of H3K4 methylation (mono-, di-, and tri-) and acetylation and trimethylation of H3K27 in regions flanking TESs of this group of protein-coding genes (that do not harbor PEs). As expected, we found that the chromatin signature of the 3′ ends of all protein-coding genes that do not harbor PEs is very different from the chromatin signature of PEs (Supplemental Fig. S10). The additional profiling of the distribution of H3K4 methylation at the gene bodies of protein-coding genes that harbor PEs at their 3′ UTRs showed that the H3K4me1 signal dramatically declines toward the 3′ ends of genes (Supplemental Fig. S1B). These comparisons further indicate that there is no bias between the chromatin signature present at PEs and the 3′ ends of the protein-coding genes (details in SI Text). We also confirmed that there is no bias between convergent gene pattern and PEs’ identification (details in SI Text).

Together with the TF binding data, these observations suggest that most of the PE midpoint-proximal regions overlap with the 3′ ends of some genes, where PASs are also located, and cleavage/polyadenylation takes place. Therefore, it is possible that these regulatory regions may serve dual purposes in Arabidopsis, by also functioning as transcriptional enhancers.

Connection between PEs and the experimentally characterized egg apparatus-specific enhancer

Only five enhancers have been extensively experimentally characterized in Arabidopsis (Weber et al. 2016) and we wanted to find out if any of these enhancers share similarities with PEs. Using Arabidopsis seedling ChIP-chip data sets (Zhang et al. 2007; Charron et al. 2009; Chodavarapu et al. 2010), we determined that one of them, the egg apparatus-specific enhancer (EASE) and its neighboring chromatin, are enriched for H3K4me1 compared to H3K4me3 (Fig. 5A; P < 0.05, two-tailed paired t-test; IGV visualization of EASE locus: Fig. 5C), providing further support for our use of H3K4me1 as a genomic property of a certain group of Arabidopsis enhancers.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 5.

Chromatin signatures in the region flanking the EASE genomic locus and the IGV visualization of the EASE locus. The accumulation of HMM posterior probability (as a measurement of total ChIP-chip signal intensity) for each ChIP-chip data set (±200-bp the EASE genomic locus). (A) Histone H3K4 mono- and trimethylation. Statistically significant (*) P < 0.05, two-tailed paired t-test. (B) Histone H3K27 acetylation and trimethylation. Statistically nonsignificant (n.s.), P = 0.43269, two-tailed paired t-test. (C) IGV visualization of the EASE locus and the ±200 bp flanking region. The signal intensity of H3K4 methylation (mono-, di-, and tri-), acetylation and trimethylation of H3K27, and RNA Pol II occupancy are shown.

EASE is located in the intergenic regions on chromosome IV and its enhancer activity has been characterized experimentally (Yang et al. 2005). Previous work demonstrated that a 77-bp module located at the EASE locus can work independently of position and orientation and is sufficient to drive spatially restricted β-glucuronidase (GUS) expression mainly in the egg apparatus; this module is well conserved in different accessions of Arabidopsis (77-bp EASE module, GenBank accession no. AX100536). We did not observe the EASE locus among the 1919 PEs, since it was not located in the nucleosome-free region in the seedling stage, which is most likely the result of EASE being an egg apparatus-specific enhancer. Likewise, we did not find significant differences between H3K27ac and H3K27me3 enrichment at the EASE locus (Fig. 5B; n.s., P > 0.05, two-tailed paired t-test), which might be because at the seedling stage, EASE is only primed by H3K4me1, and has yet to be activated as a functional enhancer in the egg cells.

Next, we examined whether the PE regions harbor modules bearing similarities to the EASE 77-bp module that has been shown to be sufficient to control specific gene expression in the egg apparatus. However, there is generally no sequence conservation between different transcriptional enhancers, so we used the 77-bp EASE module as a model to search for significant matches based on position-specific scoring matrix (FIMO [Grant et al. 2011], P < 0.0001) on both strands of each PE region. This analysis allows the identification of modules with tolerable variations to the EASE 77-bp sequences that are also statistically significant. The identified modules were termed EASE-like. We found that 77% of identified PEs harbored at least one EASE-like module (n = 1468; statistical threshold: P < 0.0001). Moreover, more than half of these PEs carry more than one EASE-like module (n = 899; P < 0.0001, χ² test). All together, these data, combined with the ability of PEs to engage in chromatin looping and being bound by a wide range of TFs, provide additional evidence that PEs have the potential to function as enhancer-like elements in Arabidopsis.

The intragenic nature of the PEs, and their potential function in seedlings, makes classical functional tests, usually performed using minimal promoters and a reporter gene, challenging to conduct and interpret. However, the relationship between PEs and known enhancers that function in seedlings can be examined using enhancer trap lines. To this end, we examined enhancer trap lines that are expressed in germinating seedlings and found that among the 10 positionally mapped enhancer trap insertion lines, (Liu et al. 2005) two PEs are located in the same genomic location as the reporter gene insertion (Supplemental Fig. S11A,B), and one PE is located in close proximity (<6 kb) from the insertion site (Supplemental Fig. S11C; details in SI Text). This result indicates that the PEs overlapping with or located close to the enhancer trap reporter gene insertion could be active enhancers in Arabidopsis. Also, the enhancer trap lines have the advantage of providing evidence of the enhancer's activity in its native genomic context, which is a crucial point given the intragenic nature of PEs.

Intriguingly, we also found that in addition to occurring in PEs, the EASE-like modules can be found in the 3′ UTRs of a much wider group of Arabidopsis genes (53% of all 3′ ends of Arabidopsis protein-coding genes have one or more EASE-like modules; P < 0.0001; details in SI Text). The presence of the EASE-like modules in 3′ ends of many Arabidopsis protein-coding genes may suggest the possibility of a broader connection between transcriptional enhancers and the sites of mRNA 3′ end processing. It also suggests that different groups of transcriptional enhancers bearing chromatin signatures that differ from those of PEs might be discovered in the future at the 3′ ends of Arabidopsis genes.

PEs can form intrachromosomal chromatin loops and physically interact with other protein-coding genes

The hallmark of transcriptional enhancers is that they can function over large genomic distances by forming chromatin loops (independently of distance and orientation) to target promoters (Sanyal et al. 2012; Spitz and Furlong 2012; Shlyueva et al. 2014). Enhancers can also interact with regions outside of promoters and form physical interactions with a target gene body by associating with elongating Pol II (Lee et al. 2015). To examine possible intrachromosomal interactions involving PE regions, we utilized Arabidopsis Hi-C data at 2-kb resolution in young seedlings (Wang et al. 2015), and identified 5773 pairs of interacting loci (2-kb each). We found that 452 PEs (24%; P = 9.2 × 10⁻¹⁵, Bedtools Fisher's exact test was used to test the statistical significance of data associations with a set of randomized control) can form 1076 intrachromosomal loops with 879 loci, termed PE-interacting loci, eILs (Fig. 6A–C; Supplemental Fig. S12; Supplemental Dataset S11). These data also suggest that these interactive PE regions are located at genomic loci that have the high probability of forming intrachromosomal interactions with another loci.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

FIGURE 6.

PEs can physically interact with other protein-coding genes intrachromosomally. (A) Illustration of the chromatin looping between PEs and eILs. (B) 24% (P = 1.7 × 10⁻¹⁴) of all PEs can form chromatin loops with other loci (eILs); these interactive PEs can be grouped based on the number of eILs they interact with (shown as bars). (C) The identified 1076 PE–eIL interactions fall into three groups: neighboring loci (red bars), intermediate-range interactions (2–10 kb, blue bars), and long-range interactions (>10 kb, green bars). (D) Most eILs occur in protein-coding gene regions (73%; P < 0.0001). eILs also occur in other genomic features, including TEs, intergenic regions, and other non-protein-coding genes. (E) Over 50% of eILs are mapped to 5′ UTRs and 3′ UTRs of protein-coding genes; and the rest are mapped to gene body or entire gene regions (from 5′ to 3′ UTRs of genes with length <2-kb). (F,G) GO analysis of protein-coding genes that interact with PEs (Fisher's exact test [*] P < 0.05, [***] P < 0.0001; Supplemental Dataset S13). (F) Regulation of transcription is the most frequent GO term in the molecular function category. (G) GO terms for subcellular localization show that many of the genes encode for membrane-associated functions. (H) 79% (P < 0.0001, χ² test) of genes encoding for membrane-associated functions interact with PEs via their 5′ or 3′ UTRs.

In addition, interactive PE regions can form multiple chromatin loops with different eILs (Fig. 6B); about 23% of PEs can form more than three loops, with some PEs forming 10–45 loops. We also found that some PEs overlap with pairs of interreacting loci that consist of two neighboring loci (Fig. 6C, panel i). Based on the size of the chromatin loops, PE–eIL interactions also involve (ii) intermediate-range interactions (2–10 kb) and (iii) long-range interactions (>10 kb), with the largest loop being 27 Mb (Fig. 6C). We also found that the intermediate- and long-range interactions occur in similar numbers, 38% and 30% of total interactions, respectively, and that the more chromatin loops a single PE locus can form, the larger the size of the loops (Supplemental Dataset S12). For example, the average size of chromatin loops for PEs forming 3–4 loops is 7.8 kb, but for the PEs forming 10–45 loops, the size is 2.74 Mb. Moreover, the more extensive loop-forming PEs tend to interact more with TEs, rather than protein-coding genes.

The eILs occur in protein-coding gene regions (73%), TEs (23%), and other non-protein-coding regions (4%) (Fig. 6D; P < 0.0001, Fisher's exact test). This indicates that most interactive PEs form chromatin loops with protein-coding genes. More than half of the protein-coding genes interact with PEs via 5′ UTRs (23%) and, intriguingly, 3′ UTRs (32%) (Fig. 6E). The remaining genes interact via the gene body or entire gene structure (Fig. 6E), which is similar to results reported in mice (Lee et al. 2015). GO analysis of the protein-coding genes interacting with PEs showed that the majority of these genes are involved in the regulation of transcription (Fig. 6F; Supplemental Fig. S13; Supplemental Dataset S13; P < 0.0001, Fisher's exact test). Half of the genes interacting with PEs are annotated as membrane function associated (Fig. 6G,H), and the majority (79%; P < 0.0001, χ² test) of these genes interacts with PEs via either their 5′ or 3′ ends (36% and 43%, respectively) (Fig. 6H).

Intriguingly, many of the genes interact with PEs through their 3′ UTRs. In yeast and mammals, the 3′ and 5′ ends of protein-coding genes often interact forming self-gene loops, which are widely involved in the regulation of these genes’ expression (Hampsey et al. 2011; Grzechnik et al. 2014). Self-gene loops are also common in Arabidopsis (Liu et al. 2016). Using the Arabidopsis Hi-C data set at gene-level resolution (Liu et al. 2016), we found that ∼22% of the eILs are located in the 3′ UTRs of genes that also form self-gene loops between their 3′ and 5′ UTRs (P < 0.0001, Fisher's exact test). It is possible that the interaction of the PEs with the 3′ ends of their target genes may be one of the mechanisms influencing transcription at the 5′ end via targeting the 3′ end of self-looped genes.

We also used a panel of 16 epigenetic data sets from young Arabidopsis seedlings (Wang et al. 2015) to examine the epigenomic signatures of PEs as well as whether different categories of eILs associated with specific genomic features (TEs or 5′ UTRs, 3′ UTRs and gene bodies of protein-coding genes) have distinct chromatin signatures. (Supplemental Fig. S14A,B, details in SI Text). Remarkably, the pool of interactive PEs has more-pronounced 3′ UTR chromatin signatures than the entire population of PEs (Supplemental Fig. S14B, bottom panel). Additionally, this pool of PEs had a noticeable enrichment of the histone H3.3 variant (Supplemental Fig. S14B). Histone H3.3 has been shown to associate with the 3′ end of genes in Arabidopsis and animals (Ooi et al. 2010; Stroud et al. 2012; Wollmann et al. 2012), associate with less stable and easier to displace than canonical nucleosomes in animals, and also actively mark transcriptional enhancers (Deaton et al. 2016). The epigenetic signatures based on comprehensive chromatin profiling confirmed our findings that PEs are closely associated with the 3′ UTRs of protein-coding genes, where cleavage/polyadenylation also takes place (Fig. 3F,G), yet bear the signatures of canonical mammalian enhancers (Fig. 1C,D).

PE regions are conserved in 1135 Arabidopsis accessions

DNA sequence conservation is often associated with important functions (Wittkopp and Kalay 2011; Burgess and Freeling 2014). Therefore, we applied comparative genomics to the collection of 1135 recently sequenced Arabidopsis accessions (1001 Genomes Consortium 2016). Enhancers have higher sequence variability than promoters within intra-species because enhancers are responsible for the recruitment of multiple TFs that control gene expression as shown in animal studies (Wittkopp and Kalay 2011). Yet, the genomic locations of functionally homologous enhancers are often conserved between species and also exhibit some degree of intra-species sequence conservation.

To gain insight into the conservation of PE positions and sequence, we first interrogated the high-diversity panel of 10 accessions (Kawakatsu et al. 2016; 1001 Genomes Consortium 2016) selected based on their geographic and genetic diversity (Supplemental Dataset S14). The calculated average of allelic variants (AVs) across the entire Arabidopsis genome in a context-independent manner was reported to be 300 nt per 3-kb region across all accessions. Analysis of AVs at each genomic position in all PEs indicated that PE regions are highly conserved within the high-diversity panel (Supplemental Dataset S14). Furthermore, analysis of all 1135 accessions showed that the average AV at PEs is significantly lower than the expected average AV across the entire genome (Supplemental Fig. S15; the expected average AV was calculated in a genomic context-independent manner based on a random set of 3-kb genomic regions, see SI methods for detailed methods). More than 72% of total PEs had AVs <300 nt per 3-kb, with the average being 256 nt per 3-kb PE region, P < 0.0001, suggesting that the sequences of PEs are conserved in the 1135 Arabidopsis accessions.