Chromatin accessibility dynamics reveal novel functional enhancers in C. elegans

1Chromatin accessibility represents an essential level of genome regulation and plays a pivotal role in many biological and pathological processes, including development, tissue regeneration, aging, and cancer. However, most genome-wide chromatin accessibility studies have been in relatively simple systems to date, including cultured or purified cells as well as early embryos. Assessing chromatin accessibility directly in complex systems composed of multiple cell types could allow for high-throughput discovery of regulatory regions whose activities are restricted to rare or undefined subpopulations of cells. This is particularly relevant for enhancers, which are thought to be highly cell-type- and temporally specific.

2The primary limitation for studying chromatin accessibility in complex systems is that most assays lack the sensitivity and precision to detect regions active only in rare subpopulations of cells or require so many cells that precise temporal synchronization of samples is impractical. However, the Assay for Transposase Accessible Chromatin using sequencing (ATAC-seq) has been shown to assess native chromatin accessibility with high sensitivity and base pair resolution while requiring orders of magnitude less starting material than other assays. This approach has been used successfully in cultured or purified cells even down to single cells, though such low input relies heavily on existing knowledge. Rather than purifying specific cell types, we wondered whether ATAC-seq could be sensitive enough to detect subtle changes in chromatin accessibility in complex mixtures of tissues and, in so doing, uncover novel biological insights that would have otherwise been obscured.

3The nematode Caenorhabditis elegans is a particularly powerful model to study chromatin accessibility in a complex system and potentially identify novel regulatory regions. C. elegans has highly synchronous life stages, as well as consistent and well-characterized cellular composition throughout each stage of development. Rapid transgenesis and transparency also make C. elegans an ideal system to efficiently validate genomic regions of functional importance and visualize tissue- or cell-specificity. While there exist some preliminary reports on chromatin states in C. elegans, high-resolution, genome-wide chromatin accessibility maps throughout development have not yet been reported. In this study, we show that studying high-resolution chromatin accessibility dynamics in synchronized C. elegans populations allows us to characterize highly reproducible changes in chromatin accessibility between developmental stages and to identify functional temporal- and tissue-specific novel enhancers in vivo. Our study provides a unique resource for defining C. elegans regulatory regions as well as a guide for the interpretation of chromatin structure in complex multitissue systems in vivo.

4To sensitively measure high-resolution chromatin accessibility at different life stages in C. elegans, we used the Assay for Transposase Accessible Chromatin using sequencing. We optimized the ATAC-seq protocol for C. elegans by including a step of native nuclei isolation by mechanical homogenization before the transposition step (see Methods; Supplemental Extended Protocol). The low input requirements of ATAC-seq (several orders of magnitude less than standard ChIP-seq) allowed us to grow C. elegans in standard conditions (i.e., plates, whereas most high-throughput assays in C. elegans require growth in liquid) and to generate three independent biological replicates that are tightly synchronized at three key life stages—early embryo, larval stage 3 (L3), and young adults—thereby limiting variation within stages (see Methods; Fig. 1A). We generated and sequenced ATAC-seq libraries, as well as an input control, to a median depth of over 17 million unique, high-quality mapping reads per sample (Supplemental Table S1). The insert size distribution of each C. elegans ATAC-seq library displays a stereotypical 147-bp periodicity that is consistent with the expected nucleosome occupancy of chromatin (Supplemental Fig. S1A), indicative of ATAC-seq library quality. We designed a computational framework to integrate the input control (Supplemental Fig. S1B) and emphasize single base pair resolution (Fig. 1B), resulting in the identification of 13,000–27,000 high-confidence, accessible peaks per developmental stage and more than 30,000 consensus ATAC-seq peaks found in at least one of the three stages (Supplemental Tables S2, S3; see Methods). The high correlation of ATAC-seq signal between each of the three biological replicates (Spearman's ρ > 0.837) demonstrates the high reproducibility of this approach (Fig. 1C). The ability to cluster samples by their developmental stage also shows that chromatin accessibility is strikingly different between these three life stages (Fig. 1C). These differences between developmental stages are likely due to both changes in accessibility within cells, as well as the organisms’ changing cellular composition throughout development. Finally, ATAC-seq signal is enriched at transcription start sites (TSSs) (Fig. 2B), another indication of the quality of ATAC-seq data. Together, these results indicate that reproducible high-resolution chromatin accessibility can be obtained from low amounts (at least an order of magnitude less than standard histone ChIP-seq or DNase-seq) of complex, multitissue samples.

5To investigate the changes in chromatin accessibility between life stages, we identified and characterized the ATAC-seq peaks that significantly changed accessibility between early embryo and L3 (12,193 peaks) (Supplemental Fig. S1C) and between L3 and young adult (783 peaks; FDR < 0.05) (Supplemental Fig. S1D; see Methods). The larger number of differentially accessible peaks (both decreased and increased) observed during the transition from early embryo to L3 versus L3 to young adult could not simply be explained by differences in sequencing depth (see Methods) and is likely due to the massive changes in cell number and tissue composition that occur during this transition.

6An example of a decrease in chromatin accessibility from embryo to L3 can be seen in the promoter region of the cav-1 gene, which is expressed during embryogenesis but not larval development (Fig. 1D; ). Conversely, several ATAC-seq peaks drastically increase from embryo to L3 in the promoter and regions upstream of the daf-12 gene, which is a key regulator of stage-specific developmental programs, particularly at L3 (Fig. 1E;, ). Confirming these specific examples, the most enriched gene ontology (GO) terms for genes with decreased chromatin accessibility from embryo to L3 include early-development terms such as embryonic morphogenesis and cell fate specification (Fig. 1F; Supplemental Table S4), while the most enriched terms for genes with increased chromatin accessibility from embryo to L3 include larval development and locomotion (Fig. 1G; Supplemental Table S5). Similarly, we observe strong enrichments of GO terms reflecting the major phenotypic changes occurring between L3 and adult, including terms like larval development and reproduction (Supplemental Fig. S1E,F; Supplemental Tables S6, S7).

7Together, these results indicate that ATAC-seq in whole organisms can identify changes in DNA accessibility that represent key biological differences between stages, regardless of whether these changes are due to activation/repression of specific regions within a cell type or to changes in cell-type composition.

8Accessible chromatin encompasses several key features of the epigenome, including active and poised regulatory regions. To verify that our ATAC-seq data correctly identify regulatory regions throughout the epigenome, we used multiple histone modification ChIP-seq data sets from modENCODE and ChromHMM to build predictive models of the epigenome. ChromHMM is a hidden Markov model that classifies regions of the genome into chromatin states (e.g., heterochromatin) using the co-occurrence of multiple histone modifications from each life stage—in our case, ChIP-seq data sets characterizing eight distinct histone modifications across all three stages (Supplemental Fig. S3A–E; Supplemental Tables S8–S11; see Methods). ATAC-seq peaks from all three stages were significantly enriched in active and poised regulatory chromatin states (e.g., promoter), as defined by this ChromHMM model, and significantly depleted in heterochromatic states (Fig. 2A; Supplemental Fig. S3B,C). ATAC-seq peaks were also enriched in H3K27me3-repressed regions, which is likely due to H3K27me3 marking inactive, yet accessible poised sites. ATAC-seq signal was correlated with individual active histone modifications at TSSs (Fig. 2B) and genome-wide (Supplemental Fig. S3F). Thus, ATAC-seq correctly identifies poised and active regulatory regions at both specific loci and genome-wide.

9To further assess the relevance of whole organism chromatin accessibility, we compared our ATAC-seq data to publicly available gene expression data (GRO-seq, and RNA-seq). Both GRO-seq (Supplemental Fig. S2A,B) and RNA-seq (Supplemental Fig. S2C,D) were positively correlated with ATAC-seq signal near the TSS. These observations support the relationship between chromatin accessibility at the TSS and gene expression, despite the numerous other layers of gene regulation (e.g., RNA stability).

10Beyond simply identifying regulatory regions important for individual life stages, chromatin accessibility dynamics should highlight regulatory regions critical for transitions from embryo to larval stages, and from larval stages to adulthood. We examined whether genomic regions that showed accessibility changes from one life stage to another were enriched for specific chromatin state transitions. Regions that lost chromatin accessibility from embryo to L3 or from L3 to young adult were enriched for transitions from active regulatory chromatin states (especially predicted enhancers) to repressed or heterochromatic states (Fig. 2D; Supplemental Fig. S3G). Conversely, the regions that gained accessibility during development were enriched for transitions from inactive chromatin states to active regulatory states (again, especially predicted enhancers) (Fig. 2C,E; Supplemental Fig. S3H). Collectively these results show that ATAC-seq performed on an entire organism, even a complex multitissue adult, is sensitive enough to detect important global changes in chromatin structure that are consistent with predicted chromatin state transitions. Thus, ATAC-seq as a single assay constitutes an attractive alternative to performing multiple histone modification ChIP-seq experiments for identifying key regulatory regions, especially considering that ATAC-seq requires orders of magnitude less input than a single ChIP-seq.

11Enhancers are key regulators of temporal- and tissue-specific gene expression that play important and conserved functions during development. However, the identification and characterization of novel enhancers remains a challenge because they can be specifically active in rare cell populations, some of which may not have even been characterized. The extent and functional importance of distal regulatory regions in the C. elegans genome has been particularly underexplored. While several methods to identify potential enhancers genome-wide have recently been developed, these methods suffer from notable drawbacks. For example, the co-occurrence of multiple histone modifications (e.g., H3K27ac, H3K4me1) or RNA polymerase II (PolII) through ChIP-seq experiments lacks the sensitivity to detect enhancers active only in rare subpopulations of cells and the resolution to precisely identify the active enhancer region. We therefore investigated whether whole-organism ATAC-seq could overcome these challenges and facilitate the identification of tissue- and stage-specific enhancers. To this end, we examined distal noncoding ATAC-seq peaks (defined as ATAC-seq peaks at least 1 kb upstream of or 0.5 kb downstream from a TSS and not in an exon). Distal noncoding ATAC-seq peaks in all three stages were highly and significantly enriched for active and repressed enhancers (Supplemental Fig S4A), as defined by the ChromHMM model and distinguished by the presence of H3K27ac (active enhancers) or H3K27me3 (repressed enhancers) (Supplemental Fig. S3A). Furthermore, these distal noncoding ATAC-seq peaks were significantly more conserved than expected by chance (Supplemental Fig. S4B,C), a defining feature of enhancers. Finally, distal noncoding ATAC-seq peaks also exhibited significant overlap with regions previously predicted to be enhancers by overlapping short capped RNA, transcription factor binding, and a H3K4me3 dearth (P < 1 × 10−323, one-sided Fisher's exact test). These analyses support the notion that distal noncoding ATAC-seq peaks include enhancers.

12In C. elegans, a small number of functional enhancers have been experimentally mapped, including four enhancers in the upstream regulatory region of hlh-1, the C. elegans MyoD ortholog that regulates muscle development. ATAC-seq peaks overlap three of these four regions, and despite its lack of statistical significance, the fourth region still exhibits noticeable ATAC-seq signal in embryos (Supplemental Fig. S4D). Given that hlh-1 is exclusively expressed in muscle, a tissue comprising <10% of C. elegans’ cellular composition, these observations indicate that ATAC-seq performed on a whole organism is sensitive enough to identify functional tissue-specific enhancers.

13We next sought to determine whether ATAC-seq dynamics could be leveraged to identify novel functional enhancers. Previous work has demonstrated that enhancers are precisely activated/inactivated at very specific times in development. We hypothesized that the temporal resolution of our ATAC-seq data (due to the precise synchronization of C. elegans populations) would enable us to capture distal regulatory dynamics of development. To experimentally test the functional enhancers predicted by our ATAC-seq data, we selected 13 distal ATAC-seq peaks (i.e., at least 1 kb away from the nearest TSS, defined on the UCSC Genome Browser) that exhibited the largest fold changes in accessibility between any two stages. We generated multiple transgenic C. elegans strains with these 13 putative enhancer regions upstream of a minimal promoter (pes-10) driving expression of green fluorescent protein (GFP) containing a nucleolar localization signal (NoLS) (Fig. 3A; Table 1; Supplemental Tables S12, S13). To control for specificity, we also tested regions flanking 10 of the 13 putative enhancer regions (Table 1; Supplemental Tables S12, S13). We created numerous (8–11) independent genetic strains for each validated regulatory site to ensure that no artifact of transgenesis (e.g., hybridization of the rol-6 marker in the extrachromosomal array) could be driving spurious GFP expression (Table 1). By fluorescence microscopy, we examined the spatiotemporal GFP pattern in these transgenic strains to assess the temporal- and tissue-specificity of the putative enhancers and control regions. Using stringent criteria for defining enhancer activity (see Methods), we found that six of the 13 putative enhancer regions (identified by dynamic distal noncoding ATAC-seq peaks) led to specific and consistent spatiotemporal GFP pattern and, for all but one of these, this pattern was observed regardless of the genomic orientation of the region—an important characteristic of enhancers (Table 1; Fig. 3B–G; Supplemental Fig. S5A–E; ). In contrast, zero of the 10 flanking regions led to such a GFP pattern, indicating enrichment for functional enhancers in distal noncoding ATAC-seq peaks that change between stages (P = 0.038, one-sided Fisher's exact test) (Table 1). Thus, dynamic chromatin accessibility on its own can mark functional enhancers and can be used to successfully identify novel distal regulatory regions.

14Importantly, the enhancers we experimentally identified display diverse spatiotemporal activity patterns. Three of the enhancer regions (putatively associated with gei-13, mlt-8, and nhr-25) are active in the head or tail hypodermis during development, while others are active specifically in the pharynx (C54G6.3) or along the flank of the worm (swip-10). In addition, the enhancers are located at a wide diversity of genomic positions relative to their putatively associated gene: upstream of the TSS (at distances varying from 1 to 9 kb), within introns, as well as downstream from the coding sequence. An interesting example is the nhr-25-associated enhancer, which is downstream from the 3′ UTR, more than 5 kb downstream from the TSS (as defined on the UCSC Genome Browser; this enhancer is also downstream from npr-10 3′ UTR on the other side, about 10 kb from the TSS) (Fig. 3D). nhr-25 is a conserved nuclear receptor primarily expressed in the hypodermis (and somatic gonad) during larval development. The 3′ nhr-25 enhancer specifically drives GFP expression in approximately 20 hypodermal cells in the head and tail of the worm during larval development (Fig. 3E). The limited expression pattern of the 3′ nhr-25 enhancer (as well as the other enhancers we identified) indicates that ATAC-seq performed in whole worms is sensitive enough to identify regulatory regions active only in specific cell types, though we cannot rule out the possibility that these regions are also accessible but their activity repressed (e.g., by H3K27me3) in other cells.

15Collectively, these findings demonstrate the power of using ATAC-seq dynamics as an unbiased approach to identify functional and conserved enhancers active in a small subset of cells. When applied to whole organisms or complex samples composed of diverse cellular populations, this approach could capture regulatory regions that may have been missed in studies of isolated cell populations.

16To explore the regulatory underpinnings of chromatin accessibility dynamics, including that of enhancers, we examined the occurrence of experimentally defined C. elegans transcription factor (TF) binding motifs in the dynamic chromatin accessibility regions identified by ATAC-seq. In regions that change chromatin accessibility between early embryo and L3 or between L3 and young adult, we observe significant enrichment of motifs associated with TF homologs and orthologs that have previously been connected to chromatin accessibility dynamics (Fig. 4A; Supplemental Fig. S6A). For example, the motif for BLMP-1, the C. elegans ortholog of human BLIMP-1/PRDM1, a TF that recruits chromatin-remodeling complexes in human B cells, is enriched in peaks that are more accessible in L3 when compared to embryos or adults. Furthermore, the motif for ELT-3, a GATA TF, is enriched in peaks more accessible in L3 than embryos (the GATA family was recently shown to be an important regulator of chromatin accessibility during human hematopoiesis).

17The DNA binding motif for EOR-1, which resembles a dimeric version of the canonical GAGA motif, was significantly enriched in distal noncoding ATAC-seq peaks (Supplemental Fig. S6B) and in ATAC-seq peaks that gained in accessibility in L3 versus embryo (Fig. 4A). This GAGA/EOR-1 motif was also present in two of the five functional enhancers that we experimentally validated, including the nhr-25-associated enhancer described above. TFs binding GAGA motifs have been identified as modulators of chromatin structure dynamics from plants to humans, and the GAGA motif itself is required for full functionality in two previously defined C. elegans enhancers. Finally, GAGA/EOR-1 factors might play an important role in the regulation of chromatin accessibility as EOR-1 genetically interacts with components of two nucleosome remodeling complexes, SWI/SNF and RSC, and a dimeric GAGA motif nearly identical to the EOR-1 motif is highly enriched in ChIP-seq peaks for two separate SWI/SNF components in C. elegans.

18To independently verify the importance of the GAGA/EOR-1 motif, we developed a machine learning model to identify the TF binding motifs that are predictive of regions that gained or lost accessibility between early embryo and larval stage 3. This machine learning method is unbiased and allowed for the integration of the largest available set of previously defined C. elegans TF binding motifs (107 high-confidence motifs from 98 TFs out of more than 750 TFs in the C. elegans genome) as well as 59 motifs discovered de novo in regions that changed in accessibility between early embryo and L3 (Fig. 4B; Supplemental Fig. S6C; see Methods). Using this machine learning framework, we identified the motifs that are the most predictive of ATAC-seq dynamics; the top three most informative motifs in our model are the GAGA motif that we identified above and motifs for known chromatin regulators (BLMP-1 and ELT-3) (Fig. 4C). Thus, machine learning independently supports the GAGA/EOR-1 motif as a potential important regulator of chromatin accessibility.

19We next investigated whether EOR-1 may play a role in chromatin accessibility changes. An EOR-1 ChIP-seq data set at L3 revealed enrichment for a dimeric GAGA motif (Supplemental Fig. S7A) and was significantly enriched for regions that overlapped ATAC-seq peaks that gain accessibility in L3 compared to early embryo (Supplemental Fig. S7B), indicating that EOR-1 is indeed bound to the GAGA motif and that EOR-1 may be involved in regulating chromatin accessibility. To examine if EOR-1 is found at closed chromatin, we quantified not only ATAC-seq signal but also ATAC-seq fragment size, as larger ATAC-seq fragments correlate with higher nucleosome occupancy and less accessible chromatin regions. EOR-1 ChIP-seq peak summits at L3 have significantly larger fragment sizes (indicative of less accessible chromatin regions) than all 24 canonical TFs and the histone deacetylase HDA-1 ChIP-seq peak summits (FDR < 0.05) (Fig. 4D; Supplemental Fig. S7C,D; Supplemental Table S14). In fact, the only factor with larger ATAC-seq insert sizes was the heterochromatin-associated protein HPL-2. The larger ATAC-seq insert sizes at EOR-1 binding sites are unlikely to be due to the binding of other TFs hindering the ability of transposons to interact with the genomic DNA because EOR-1 ChIP-seq peaks did not show unusual overlap with the other TFs assessed (Supplemental Fig. S7E). These observations are consistent with the possibility that EOR-1 is uniquely present at regions with less accessible chromatin compared to other TFs.

20To more closely investigate the chromatin accessibility landscape of EOR-1 peaks at the L3 stage, we generated aggregated histograms of the median ATAC-seq fragment sizes for the ChIP-seq peak summits of EOR-1, all canonical TFs, and the chromatin-associated factors HPL-1 and HDA-1 (Fig. 4E; Supplemental Fig. S7F). As expected, canonical TFs are predominantly found at regions with short ATAC-seq insert sizes (<147 bp, the length of DNA wrapped around a nucleosome), whereas chromatin factors, especially HPL-2, are also present at regions with larger ATAC-seq fragment sizes (with many >147 bp). Strikingly, EOR-1 was found both at regions with short ATAC-seq insert sizes as well as at regions with larger ATAC-seq fragment sizes (>147 bp). Together, these results suggest that EOR-1 may, in some cases, bind to or at least be immediately adjacent to nucleosome-occupied DNA.

21We next quantified nucleosome occupancy dynamics in the summits of ChIP-seq peaks using publicly available histone H3 ChIP-seq data. EOR-1 peak summits exhibited a large decrease in nucleosome occupancy from embryo to L3, and this was distinct from the canonical TFs assessed (Fig. 4F). In contrast, chromatin factors such as HPL-2 and HDA-1 exhibited an increase in nucleosome occupancy. These results suggest that EOR-1 binding sites exhibit changes in chromatin accessibility during the transition from early embryos to L3 that are unique from the other factors analyzed.

22Taking these observations collectively, we propose four potential models to explain the unique binding profile of EOR-1: (1) EOR-1 is capable of binding to both open and closed chromatin (depending on the loci or the tissue); for example, EOR-1 could bind open sites in one tissue but closed sites in another tissue (Fig. 5A); (2) EOR-1 binds immediately adjacent to nucleosomes (resulting in larger ATAC-seq fragment sizes), including those nucleosomes that shift location between early embryo and L3 (explaining the change in accessibility and nucleosome signal) (Fig. 5B); (3) EOR-1 binds open chromatin as part of a larger complex, including factors not assayed here (e.g., EOR-2), thereby hindering the ability of transposons to interact with the genomic DNA (Fig. 5C); and (4) EOR-1 binds closed chromatin as a pioneer factor and contributes to its opening (Fig. 5D). Future experiments, such as nucleosome binding assays and chromatin profiling of eor-1 mutants, will be needed to distinguish between these models and to fully elucidate the mechanism underlying the unusual binding profile of EOR-1. Collectively, these analyses suggest that ATAC-seq on a complex sample can identify factors with unusual chromatin binding patterns that could regulate chromatin dynamics during development.

23Here, we show for the first time that sensitively measuring chromatin accessibility in a whole metazoan can detect dynamic changes throughout development and even identify novel functional enhancers active in only a small subset of the whole organism. Among the three developmental stages surveyed, we identified over 30,000 accessible sites which could serve as a catalog to facilitate the discovery of previously unknown distal regulatory loci (Supplemental Table S3) such as insulators and enhancers. This approach, developed here for C. elegans, should be readily applicable to other complex samples in vivo such as whole mammalian organs or tumor samples.

24In C. elegans, enhancer identification has historically been limited, with most previous studies employing a single gene approach and focusing on promoter-proximal regions. Other groups have identified enhancers genome-wide but have not functionally validated predicted enhancer activity. In this study, we identify and functionally characterize distal ATAC-seq peaks as novel, active enhancers. These enhancers have a range of spatiotemporal activity patterns that are orientation-independent and are found at a diversity of genomic locations, suggesting that enhancers may be more prevalent in C. elegans gene regulation than previously appreciated. We note that some of these regions could also be promoters, especially when they are closer to the TSS, since promoters could also function in a bidirectional manner. Indeed, two of the five regions we identified as enhancers (gei-13 and mlt-8) have previously been suggested to be promoters (despite their ChromHMM enhancer annotation). Thus, further experimental work will be required to verify the exact nature of the regulatory regions defined by ATAC-seq.

25An interesting example is the enhancer downstream from the conserved transcription factor nhr-25. This putative nhr-25 enhancer is located 5 kb downstream from the nhr-25 TSS and acts in an orientation-independent manner. Thus, this enhancer might be looping to the nhr-25 TSS to enhance transcription, although the limited spatial resolution of current Hi-C data in C. elegans does not provide information on the three-dimensional chromatin structure of the nhr-25 locus. This enhancer model contrasts with the promoter-proximal model most often thought to regulate gene expression in C. elegans but agrees with the recent Hi-C study which identified insulator-like loci throughout the C. elegans genome. Insulators are important regulators of three-dimensional chromatin architecture. Interestingly, nhr-25 orthologs in fly (Ftz-F1), mouse (Nr5a1), and human (NR5A1) exhibit a consistent signature of three-dimensional chromatin architecture downstream from the 3′ UTR (insulator class I in flies and CTCF-binding sites in mice and humans). While C. elegans does not possess a CTCF ortholog, our results raise the intriguing possibility that the nhr-25 enhancer we have identified is evolutionarily conserved.

26Furthermore, we have uncovered a potential role for a likely GAGA factor in C. elegans, EOR-1. The EOR-1 motif we initially detected closely resembles a dimeric version of the canonical GAGA motif bound by Trl/GAGA-Associated Factor (GAF) in Drosophila. GAF is a multifaceted transcription factor that can associate with heterochromatin, remodel chromatin in concert with nucleosome remodelers, and act as a transcriptional activator, in part due to its ability to increase chromatin accessibility. Like Drosophila GAF, C. elegans EOR-1 is a transcriptional activator in the Ras/ERK signaling pathway that has been suggested to also repress gene expression. GAF and EOR-1 are also similar in that both proteins have a BTB/POZ domain on their N terminal as well as C2H2 zinc-fingers and polyQ domains on their C terminals. This dual action (gene activation and repression) is particularly interesting considering the bimodal binding pattern we found when examining ATAC-seq fragment size in EOR-1 binding sites.

27In C. elegans, EOR-1 genetically interacts with at least two chromatin remodeling complexes (SWI/SNF and RSC). This is noteworthy given our findings that EOR-1 is found in less accessible and potentially even nucleosome-occupied regions of the genome and that the EOR-1 motif is predictive of increased accessibility in development. GAGA-factors are conserved regulators of gene expression, and a GAGA motif has been found to be necessary for enhancer functionality in C. elegans. Indeed, two out of five of the novel validated enhancers identified in this study (nhr-25 and swip-10) contain an EOR-1/GAGA motif. Collectively, these data suggest that EOR-1 might regulate chromatin accessibility at enhancers in C. elegans and potentially other species.

28Using C. elegans as a paradigm, we have shown that ATAC-seq performed on complex, heterogeneous samples can reveal novel, spatiotemporally specific genetic regulators and that measuring chromatin accessibility across a developmental time course can identify important dynamic regions. We have highlighted important applications of this approach: discovering functional distal regulatory regions active in only a small subset of the total sample and identifying candidate regulators of genome-wide chromatin dynamics. This data set, which represents an initial atlas of genome-wide chromatin accessibility and candidate distal regulatory sites in C. elegans, should be a valuable resource for the community. The fact that a genome-wide chromatin accessibility assay performed in whole organisms can sensitively identify previously undiscovered functional enhancers in vivo raises the exciting possibility that distal regulation plays a more important role than previously believed in the nematode.

29Brief methods can be found below, but in all cases more details can be found in the Supplemental Methods.

30Three sets of completely independent biological replicates were prepared by harvesting and then flash-freezing tightly temporally synchronized samples at three different stages: early embryo (in utero), larval stage 3 (36 h post-egg lay), and young adult (57 h post-egg lay). Native nuclei were purified from frozen samples using mechanical homogenization as previously described. The purified nuclei were immediately used for the ATAC-seq protocol. An input control was also generated by using 10 ng of genomic DNA. Sequencing was performed using 101-bp paired-end sequencing on an Illumina HiSeq 2000.

31The ATAC-seq libraries were sequenced to a median depth of over 17 million unique, high-quality mapping reads per sample (Supplemental Table S1). Prior to mapping, standard next-generation sequencing quality control steps, as well as ATAC-seq-specific quality control steps, were performed. Individual replicate ATAC-seq peaks were called using a custom pipeline which used MACS (v2.1) to call the peaks. To identify the most high-confidence set of peaks, we employed a strategy that emphasized peaks which were consistently observed across replicates. All peaks were masked for regions with significant signal in the input control and for previously identified regions known to give spurious results in next generation sequencing assays. The result was a set of 30,832 consensus ATAC-seq peaks.

32For the enhancer screen, each putative regulatory region was cloned in the pL4051 plasmid (a gift from Andrew Fire) upstream of a minimal promoter (pes-10) driving expression of a C. elegans intron- and photo-stability-optimized GFP containing an N-terminal nucleolar localization signal. Putative regulatory regions were chosen by selecting ATAC-seq peaks that exhibited the largest differential accessibility between two stages and that were at least 1 kb from a transcription start site (as defined on the UCSC Genome Browser). Flanking negative control regions were chosen by selecting regions within 2 kb of the putative regulatory regions that were not in peaks of accessibility. Primers were designed to amplify each region as well as 50–500 bp flanking either side (Supplemental Table S13).

33Stable extrachromosomal transgenic lines for putative enhancer regions (in both orientations), negative control regions, and the no-insert controls were generated. Multiple independent lines were generated per construct (Table 1). For each line, mixed-staged worms were screened for GFP signal distinct from the no-insert control background signal, which is 1–2 nuclei near the pharynx in all larval and adult stages (Supplemental Fig. S5E). To quantify the consistency of GFP expression pattern, all lines (including the negative control lines) (Table 1) were scored for GFP expression pattern in a blinded manner.

34To predict changes in accessibility between early embryo and larval stage 3, the number of each mapped C. elegans motif from cisBP (v1.02) and each de novo discovered motif found in the dynamics ATAC-seq peaks between early embryo and L3 (59 in total) was counted to create a matrix of 166 motif-counts and 30,832 ATAC-seq peaks. The ATAC-seq peaks were split into a training set (70%) and a testing set and several classification models evaluated. We found that a generalized boosting model (GBM) (https://CRAN.R-project.org/package=gbm) had the best accuracy while still allowing for interpretation of which motifs were the most informative. Given the unbalanced classification problem, we used balanced accuracy as our primary metric of classification success (Supplemental Fig. S6D).

35All analysis source code is freely available at https://github.com/brunetlab/CelegansATACseq, as well as in the Supplemental Material.

36The raw data as well as stage-specific peaks for all three stages and the input control have been submitted to the NCBI Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) under accession number GSE89608.