ATAC-seq reveals regional differences in enhancer accessibility during the establishment of spatial coordinates in the Drosophila blastoderm

Abstract

Establishment of spatial coordinates during Drosophila embryogenesis relies on differential regulatory activity of axis patterning enhancers. Concentration gradients of activator and repressor transcription factors (TFs) provide positional information to each enhancer, which in turn promotes transcription of a target gene in a specific spatial pattern. However, the interplay between an enhancer regulatory activity and its accessibility as determined by local chromatin organization is not well understood. We profiled chromatin accessibility with ATAC-seq in narrow, genetically tagged domains along the antero-posterior axis in the Drosophila blastoderm. We demonstrate that one-quarter of the accessible genome displays significant regional variation in its ATAC-seq signal immediately after zygotic genome activation. Axis patterning enhancers are enriched among the most variable intervals, and their accessibility changes correlate with their regulatory activity. In an embryonic domain where an enhancer receives a net activating TF input and promotes transcription, it displays elevated accessibility in comparison to a domain where it receives a net repressive input. We propose that differential accessibility is a signature of patterning cis-regulatory elements in the Drosophila blastoderm and discuss potential mechanisms by which accessibility of enhancers may be modulated by activator and repressor TFs.

Figure 1.

Domain-restricted ATAC-seq profiling along antero-posterior axis in the blastoderm embryo. (A) Selected AP domains are targeted by expressing a nuclear tag, UNC84-3×FLAG, under control of well-characterized enhancers of gap and pair-rule genes. All reporter constructs are integrated at the same genomic site (attP2) (Pfeiffer et al. 2008) to standardize genetic background. After homogenization of staged embryos (cellularizing blastoderm, stage 5, 2:50–3:10 h after egg laying), tagged nuclei are affinity-purified with anti-FLAG antibodies, followed by Tn5 transposase fragmentation and ATAC-seq library preparation. An ATAC-seq library representing an entire pool of nuclei from homogenized embryos (whole-embryo) serves as a control. (B) Overview of the tagged domains (D1–D7). Selected embryos immunolabeled with an anti-FLAG antibody show spatially restricted expression domains of the nuclear tag. Each domain is additionally schematized (green bars) to indicate its position along the AP axis. Embryos are positioned with anterior to the left and dorsal side up.

Figure 2.

Regional differences in chromatin accessibility. Accessibility profiles of individual tagged domains and a whole-embryo control at the locus of giant (gt), a gene of the AP patterning network of the gap class. Tracks show normalized coverage of 1- to 100-bp ATAC-seq fragments, smoothed over a sliding window of 15 bp. AP positions of the profiled domains are indicated schematically on the left (green shading). Blue bars and underlying shaded regions indicate coordinates of known giant enhancers (REDfly names and references provided in Supplemental Table S8). Spatial activity of each enhancer in blastoderm embryos is illustrated above (RNA in situ hybridization of a reporter gene) (reprinted from Schroeder et al. 2004). Note that gt_(-2)_broad was used as a driver of UNC84-3×FLAG in D7, and the additional copy in attP2 partially contributes to its elevated ATAC-seq signal (as discussed in Supplemental Methods). Genomic coordinates and gene models: FlyBase Release 5.57 (Gramates et al. 2017).

Figure 3.

Genome-wide differences in chromatin accessibility profiles along the AP axis. (A) Principal component analysis (PCA) of genome-wide accessibility variation across individual tagged domains (solid circles) and whole-embryo controls (crossed squares). Duplicates are represented as separate data points and color-coded by genotype (D1: red; D2: orange; D3: purple; D4: dark blue; D5: light blue; D6: dark green; D7: light green). PCA is based on accessibility signal (total count of Tn5 transposase cuts) over 17,345 high-confidence ATAC-seq peaks. One replicate of the D6 domain shows high similarity to whole-embryo controls, indicative of potential sample contamination with untagged nuclei. (B) Pie chart shows proportion of the accessible genome (combined size of high-confidence ATAC-seq peaks) represented by constitutive peaks that show no significant variation in their accessibility signal along the AP axis (gray), and differential peaks that are supported by a single pair-wise comparison (light blue) and multiple pair-wise comparisons (dark blue) in the DESeq2 analysis. (C) Example scatter plots show fold-change of ATAC-seq signal between selected domains against the mean normalized signal intensity. Gray: distribution of constitutive peaks, blue points: individual differential peaks (false discovery rate, FDR < 1%). The number of ATAC-seq peaks showing significant increase (up arrow) and decrease (down arrow) of their signal is indicated in the upper right corner of each plot.

Figure 4.

Differential ATAC-seq peaks display strong features of axis patterning enhancers. (A) Proportional distribution of genomic annotations among different classes of accessible regions: all high-confidence ATAC-seq peaks (all peaks), constitutive peaks, differential peaks, top quarter of differential peaks (highest values of maximum log₂ fold-change reported in DESeq2) (Supplemental Fig. S10). (UTR) 5′ and 3′ untranslated regions, (CDS) coding sequence. (B) Bar plot shows the proportion of different classes of ATAC-seq peaks that map to intronic and intergenic regions (numbers of intervals in the legend) and colocalize with ChIP signal of different classes of proteins. (ORI complex) ChIP-seq peaks of ORC2 (origin recognition complex subunit 2) (Eaton et al. 2011), (Insulator proteins) ChIP-chip peaks of BEAF-32, CP190, CTCF, and Su(Hw) (Celniker et al. 2009), (DV TFs) ChIP-chip peaks of four maternal and zygotic DV TFs: Dorsal, Mothers against dpp, Snail, and Twist (MacArthur et al. 2009), (AP TFs [broad set]) ChIP-chip peaks of 14 maternal, gap, and pair-rule AP TFs: bicoid, caudal, giant, hunchback, knirps, Kruppel, huckebein, tailless, Dichaete, fushi tarazu, hairy, paired, runt, and sloppy paired 1 (Li et al. 2008; MacArthur et al. 2009), (AP TFs [narrow set]) ChIP-seq peaks of six maternal and gap AP TFs: Bicoid, Caudal, Giant, Hunchback, Knirps, and Kruppel (Bradley et al. 2010). Asterisks indicate significant differences between constitutive and differential peaks; (***) P < 0.0001, Fisher's exact test. (C) Bar plot shows the proportion of annotated CREs that overlap differential peaks. Total number of CREs from each category is indicated above the bars: Vienna Tiles (Kvon et al. 2014) active at stage 4–6, REDfly CREs (Gallo et al. 2011) active in blastoderm embryos, and AP enhancers driving patterned expression specifically along the AP axis (Supplemental Table S8). (D) Box plots show distributions of maximum log₂ fold-change of accessibility signal reported for all differential peaks (blue) and differential peaks overlapping the three categories of annotated CREs (green).

Figure 5.

Elevated accessibility of AP enhancers in tagged domains coincides with their activity. (A) Coverage of 1- to 100-bp ATAC-seq fragments (mean over two replicates) of four selected AP enhancers. Comparison between a whole-embryo control (control: gray), a tagged domain that encompasses the enhancer's activity pattern (IN: blue), and a tagged domain from which the enhancer's activity is excluded (OUT: green). Activity pattern along the AP axis is indicated schematically in dark blue, with color-coded outlines representing positions of respective domains. RNA in situ hybridization images of a reporter gene: gt_(−10) and nub_(+2) (reprinted from Schroeder et al. 2004), Dfd_(−13) (reprinted from Fisher et al. 2012), and Antp_(−16) (reprinted from Kazemian et al. 2010). (B) Violin plots show distribution of log₂ fold-changes of ATAC-seq signal (total count of Tn5 transposase cuts, mean over two replicates) between a given tagged domain and its corresponding whole-embryo control, over two classes of AP enhancers: IN (blue) and OUT (green). Asterisks indicate significant differences between IN and OUT enhancers (Student's t-test); (*) P-value < 0.05, (***) P-value < 0.0001. Numbers of enhancers in each class are indicated below individual plots. Domains D3 and D6 are not presented; due to their limited size, none of the AP enhancers had its activity pattern fully included in the domains.

Figure 6.

Accessibility of AP enhancers is highly correlated with their transcriptional output. (A) Coverage of 1- to 100-bp ATAC-seq fragments over an enhancer of runt, run_(+17), measured in D1–D7 domains and a whole-embryo control (mean over replicates). Domains are color-coded as in the schematic that shows their positions along the AP axis and with respect to the activity pattern of the enhancer (dark blue). (B) ATAC-seq signal (mean number of Tn5 transposase cuts per bp) plotted against the proportion of an embryonic domain in which the enhancer is active (active nuclei). Each point represents an individual replicate of D1–D7 samples and whole-embryo controls (pooled replicates from multiple strains), color-coded as in panel A. D6 replicate 1 is excluded due to its close similarity to whole-embryo controls (Fig. 3A). (C) Box plot represents distribution of correlation coefficients and linear regression slopes across 88 AP enhancers. (D) Scatter plots of example enhancers, each representing a different quarter of correlation coefficients (all enhancers in Supplemental Fig. S14). (E) Box plots show ATAC-seq signal distribution (mean number of Tn5 transposase cuts per bp) of active AP enhancers (signal from tagged domains with 100% active nuclei; active state), inactive AP enhancers (signal from tagged domain with 0% active nuclei; inactive state), and 9309 background regions of the genome (mean signal across all tagged domains; background regions).

Figure 7.

Local accessibility modulation within axis patterning enhancers. (A) Normalized frequency of Tn5 transposase cleavages, (B) normalized coverage of 1- to 100-bp nucleosome-free ATAC-seq fragments, and (C) predicted nucleosome occupancy along Ubx_(−10) enhancer in D7 domain (blue; 100% active nuclei) and D1 domain (green, 0% active nuclei). Inset in a gray frame represents magnification of the shaded region. Nucleosome occupancy was predicted with NucleoATAC (Schep et al. 2015). r = Pearson correlation coefficient of the compared profiles. (D) Box plots show distribution of correlation coefficients across all AP enhancers. Distribution of transposase cuts, coverage of 1- to 100-bp ATAC-seq fragments and predicted nucleosome occupancy was compared between tagged domains with 100% active and 0% active nuclei. Coverage of 1- to 100-bp fragments was compared with the predicted nucleosome occupancy in tagged domains representing active enhancer states.