Article

Chromatin Architecture Emerges during Zygotic Genome Activation Independent of Transcription Graphical Abstract

Authors Clemens B. Hug, Alexis G. Grimaldi, Kai Kruse, Juan M. Vaquerizas

Correspondence [email protected]

In Brief Chromosome architecture gains complexity during development in concert with the onset of zygotic genome activation but without reliance on transcription.

Highlights d

Chromatin conformation emerges during Drosophila zygotic genome activation

d

TAD boundaries are established at early transcription loci enriched in housekeeping genes

d

TADs and TAD boundary formation are transcription independent

d

Zelda depletion results in locus-specific loss of insulation

Hug et al., 2017, Cell 169, 216–228 April 6, 2017 ª 2017 Elsevier Inc. http://dx.doi.org/10.1016/j.cell.2017.03.024

Article Chromatin Architecture Emerges during Zygotic Genome Activation Independent of Transcription Clemens B. Hug,1 Alexis G. Grimaldi,1 Kai Kruse,1 and Juan M. Vaquerizas1,2,* 1Max

Planck Institute for Molecular Biomedicine, Roentgenstrasse 20, 48149 Muenster, Germany Contact *Correspondence: [email protected] http://dx.doi.org/10.1016/j.cell.2017.03.024 2Lead

SUMMARY

Chromatin architecture is fundamental in regulating gene expression. To investigate when spatial genome organization is first established during development, we examined chromatin conformation during Drosophila embryogenesis and observed the emergence of chromatin architecture within a tight time window that coincides with the onset of transcription activation in the zygote. Prior to zygotic genome activation, the genome is mostly unstructured. Early expressed genes serve as nucleation sites for topologically associating domain (TAD) boundaries. Activation of gene expression coincides with the establishment of TADs throughout the genome and co-localization of housekeeping gene clusters, which remain stable in subsequent stages of development. However, the appearance of TAD boundaries is independent of transcription and requires the transcription factor Zelda for locus-specific TAD boundary insulation. These results offer insight into when spatial organization of the genome emerges and identify a key factor that helps trigger this architecture. INTRODUCTION In the interphase nucleus, chromatin is organized in a compartmentalized hierarchical fashion consisting of chromosome territories and globular domains of active and inactive chromatin (Bickmore and van Steensel, 2013; Lieberman-Aiden et al., 2009; Sexton and Cavalli, 2015). Within this chromatin organization, topologically associating domains (TADs) have emerged as a critical organizational feature of eukaryotic genomes (Dixon et al., 2012; Hou et al., 2012; Nora et al., 2012; Sexton et al., 2012). TADs—regions of high physical self-interaction displaying significant insulation from neighboring regions—are essential for the correct three-dimensional arrangement of enhancer-promoter contacts and deployment of gene expression programs (Dixon et al., 2015). When elements that determine TAD boundaries, such as binding motifs for architectural proteins, are disrupted, this can result 216 Cell 169, 216–228, April 6, 2017 ª 2017 Elsevier Inc.

in distinct developmental disorders and cancer (Flavahan et al., 2016; Hnisz et al., 2016; Lupia´n˜ez et al., 2015; Narendra et al., 2015). Chromatin organization at the TAD level seems to be evolutionarily conserved in eukaryotic species, although a significant number of regulatory loops are cell type specific (Dixon et al., 2012, 2015; Rao et al., 2014; Vietri Rudan et al., 2015). Furthermore, the overall chromatin topology seems to be rather invariant to perturbations of the cell’s environment, and enhancers and promoters interact structurally even before gene expression starts (Ghavi-Helm et al., 2014; Li et al., 2015). However, despite its fundamental role in regulating gene expression, it is unclear when during eukaryotic development this organization is established and how it is affected by the onset of transcriptional activity. To address this question, we examined chromatin conformation during early development in Drosophila melanogaster, when most of the zygotic genome is transcriptionally inert (De Renzis et al., 2007). Early stages of development in the fruit fly are characterized by a series of 13 synchronous nuclear divisions until the zygotic embryo is cellularized and undergoes zygotic genome activation (ZGA) at nuclear cycle (nc) 14 (Foe and Alberts, 1983). Before ZGA, at around nc8, a small subset of genes is already actively transcribed (Lott et al., 2011). ZGA then occurs as a stepwise process consisting of a large-scale de novo recruitment of RNA polymerase II (RNA Pol II) at nc13 and the subsequent activation of transcriptional elongation genome-wide at nc14 (Blythe and Wieschaus, 2015). By probing chromatin conformation at stages surrounding ZGA, we aim to decouple the effects of transcriptional activity on chromatin organization in a system where transcription is naturally not occurring. This allows us to examine how chromatin is organized in a naive state, when the zygotic genome is not being transcribed, and thereby assess how the subsequent transcriptional activation affects chromatin conformation. In this study, we generate in situ Hi-C maps for embryos at times surrounding ZGA. Analysis of these datasets allows us to uncover drastic changes in chromatin conformation concomitant with ZGA. Before ZGA, the genome is mostly unstructured, except for
180 regions enriched in housekeeping genes that resemble TAD boundaries. De novo recruitment of RNA Pol II co-occurs with the establishment of TADs throughout the genome and the three-dimensional co-localization of housekeeping gene clusters, features which remain mostly invariant throughout development. Pharmacological inhibition of RNA

Figure 1. Chromatin Architecture Is Established during Zygotic Genome Activation

A

B

D

C

Pol II activity does not preclude TAD formation, but has a marked effect on the contact properties of TADs and prevents the colocalization of housekeeping gene-enriched TAD boundaries. This suggests that while the mechanism of TAD formation is independent of transcription, transcription is necessary to maintain a correct organization of chromatin. Finally, we show that the transcription factor Zelda is necessary for TAD boundary formation at a subset of loci highly enriched in Zelda binding. Overall, our study provides insight into the temporal dynamics of the establishment of chromatin organization in vivo during development and the mechanisms associated with its emergence. RESULTS Chromatin Conformation Is Established at the Onset of ZGA To determine chromatin conformation at different stages during development with nuclear cycle resolution, we performed in situ Hi-C (Rao et al., 2014) for hand-sorted non-mitotic (S-phase) embryos at nc12, nc13, and nc14 and for embryos at 3–4 hr post fertilization (hpf) (Figure 1A). To achieve nuclear cycle resolution, we used a fly line containing an eGFP-PCNA transgene (Blythe and Wieschaus, 2016) that allowed us to accurately classify developing embryos based on their cell-cycle status and their developmental stage (Figure S1). Visual inspection of normalized chromatin contact probability maps for all samples at 10 kb resolution revealed a striking difference in chromatin organization at nc12 and nc13 when compared with later embryos (Figure 1B). Overall, later-stage embryos displayed a chromatin topology similar to that previ-

(A) Schematic representation of the experimental design. Arrows represent embryo collection time points for in situ Hi-C experiments. (B) Chromatin interaction probability maps at 10 kb resolution. (C) log2 fold-change of interaction probabilities (red, increased; blue, decreased) for each stage relative to nc12. (D) Representative fluorescence microscopy images of manually staged eGFP-PCNA embryos used for in situ Hi-C. White bars represent 200 mm. See also Figures S1 and S2 and Table S1.

ously described for older embryos and Kc167 cells (Figure S2C). In contrast, early-stage embryos were characterized by a low level of chromatin organization, displaying unusually uniform chromatin contact probabilities across long genomic distances. This observation was supported by unsupervised classification of our Hi-C data for each replicate that resulted in clustering according to developmental timing (Figure S2A). TADs, TAD boundaries, and the previously reported plaid-like pattern of longrange inter-TAD interactions (Dixon et al., 2012; Hou et al., 2012; Lieberman-Aiden et al., 2009; Nora et al., 2012; Sexton et al., 2012) were increasingly apparent during ZGA. To validate our Hi-C data, we compared our chromatin contact maps from 3–4 hpf embryos with chromatin contact maps for Drosophila generated previously for embryos at 16–18 hpf and the Kc167 cell line (Li et al., 2015; Sexton et al., 2012), which resulted in a highly significant correlation (Pearson’s correlation coefficients: 0.93 for Kc167 cells; 0.87 for 16–18 hpf embryos; p value < 2.2 3 1016 in both cases; Figures S2C and S2D). Furthermore, we were able to detect enhancer-promoter loops in our Hi-C data previously reported using 4C-seq for 3–4 hpf and 6–8 hpf embryos (Figure S2E; Ghavi-Helm et al., 2014). Importantly, in situ Hi-C on mitotic embryos up to nc14 indicated that the differences in contact probabilities observed for early embryos were not related to changes in chromatin compaction due to cell cycle (Naumova et al., 2013) and suggested instead that they were specific to the developmental stage studied (Figures 1B and S2B). Together these results indicate that our in situ Hi-C data are robust and of high quality. To examine changes in chromatin conformation during development, we calculated the difference in contacts for all time points relative to nc12 (Figure 1C). This analysis revealed that, starting from nc13 and continuing during ZGA, embryos displayed an increasingly marked enrichment of chromatin contacts within TADs. At the same time, embryos displayed a significant depletion of contacts surrounding TADs, suggesting an establishment of insulation at TAD boundaries. Together, these data indicate that chromatin structure is dramatically

Cell 169, 216–228, April 6, 2017 217

A

B

D

E

C F

Figure 2. TAD Boundaries Are Associated with Active Transcription (A) Number of nuclear cycle-specific boundary calls at each developmental time point (dots represent individual biological replicates; bars represent mean values). (B and C) Overlap of boundaries between different developmental stages (B) and previously published data (C). (D) Corrected Hi-C contact map and insulation scores for nc12 embryos at 10 kb resolution. Negative insulation scores (red) correspond to areas with depleted contacts between neighboring domains and are therefore indicative of TAD boundaries. log2 fold enrichment of RNA Pol II ChIP-seq signal over input (bottom). (E) Heatmaps of insulation scores around RNA Pol II-bound loci. Rows are sorted according to the strength of the RNA Pol II peak. (F) Data from (E) averaged over all RNA Pol II peaks. ***p value < 1 3 104; permutation test; n = 10,000; mean value (dashed lines) and 95th-percentile intervals (gray area) of permutation tests are shown. See also Figure S3 and Tables S2 and S3.

remodelled from a naive unordered state into a complex, feature-rich state in the developmental window between nc12 and nc14. Developmental Maintenance of TAD Boundaries and RNA Pol II Binding To systematically assess differences in chromatin organization at the TAD level during development, we used a recently developed insulation score metric (Crane et al., 2015), a measure of the overall level of contacts that occur over a given locus (Figures S3A–S3C). Insulation scores displayed a significantly lower variance throughout the genome for nc12 and nc13 embryos (Figures S3D and S3E), supporting our previous observation that

218 Cell 169, 216–228, April 6, 2017

at early stages, prior to nc14, chromatin has a limited higherorder structure. Next, we used local minima of the insulation scores to detect TAD boundary-like regions across the genome. The number of boundary calls increased progressively from nc12 and nc13 embryos, before reaching a plateau after the activation of the zygotic genome, which is then maintained throughout later stages of development (Figure 2A). To assess the level of chromatin remodelling during development, we compared the position of these boundaries for the different samples, which resulted in a remarkable level of overlap across stages (Figures 2B and S3F). Boundaries detected in 3–4 hpf embryos were also maintained in later-stage embryos and Kc167 cells (Figure 2C), suggesting

A

C

Figure 3. TAD Boundaries Are Enriched in Co-localizing Housekeeping Genes

B

E

(A and B) Enrichment of housekeeping genes and enhancers in TAD boundaries. NS: not significant; **adjusted p value < 1 3 104; ***adjusted p value < 2.2 3 1016; Fisher’s exact test with Bonferroni correction. (C) Schematic representation of boundaryboundary contacts analysis. Numbers represent ordinal boundary-boundary distances from reference. Arrows indicate boundary regions. (D) Aggregate Hi-C plots of boundary-boundary contact regions. Maps show the average contact probability in a 100 kb window around the crossing point of two boundaries in the contact matrix. Columns represent increasing distances between boundaries. ND: not detected. (E) Distribution of contact probabilities in boundary-boundary contact regions compared with the contact density in mirrored control contact regions. NS: not significant; **adjusted p value < 1 3 104; ***adjusted p value < 2.2 3 1016; MannWhitney U test with Bonferroni correction. See also Figure S4.

D

that, once established, the majority of these boundaries are invariant. These results suggest that concomitantly with the activation of the zygotic genome, there is a marked reorganization of chromatin structure that is subsequently stably maintained at later stages of development. This is in agreement with findings in other organisms that have shown a significant degree of stability at the TAD level between cell types and regions of synteny (Dixon et al., 2012, 2015; Vietri Rudan et al., 2015). We next investigated the genomic features present at TAD boundaries at early stages of development. Previous studies have shown that some genes are transcribed before the main wave of zygotic genome activation (Lott et al., 2011). To test whether these early expressed loci are associated with establishing chromatin conformation, we determined sites of active transcription using stage-specific RNA Pol II chromatin immunoprecipitation followed by sequencing (ChIP-seq) data for nc12 embryos (Blythe and Wieschaus, 2015). Visual inspection of the data suggested a strong correlation between TAD boundary-like regions and RNA Pol II occupancy (Figure 2D). A genome-wide examination of this correlation revealed a striking

RNA Pol II dose-dependent enrichment of insulation for these sites, which was also apparent at later nuclear cycles (Figures 2E and 2F). 39%, 43%, 55%, and 62% of all RNA Pol II peaks were located within TAD boundaries at nc12, nc13, nc14, and 3–4 hpf, respectively. Consistent with this, 70% (nc12) to 91% (3–4 hpf) of boundaries at each developmental stage were located within 20 kb of a RNA Pol II peak (p value < 1 3 104; permutation test; n = 10,000). Both stalled and elongating RNA Pol II-bound loci displayed negative insulation, although this was more pronounced for elongating RNA Pol II (Figure S3J). These results suggest that RNA Pol II-bound loci are associated with the establishment of inter-TAD insulation. TAD Boundaries Are Enriched in Housekeeping Genes and Co-localize in 3D The association between RNA Pol II binding and TAD boundaries prompted us to examine the functional properties of genes located at TAD boundaries. To do so, we classified loci into housekeeping or developmental/tissue-specific genes based on their expression patterns (Graveley et al., 2011). Using a stringent threshold, a total of 2,891 genes were selected as housekeeping, which displayed a strong enrichment in TAD boundary regions and a significant depletion inside TADs throughout the developmental nuclear cycles analyzed (Figure 3A). Similar results were obtained when housekeeping genes were selected based on the presence of 5,956 housekeeping enhancers as determined by STARR-seq (Figure 3B; Zabidi et al., 2015). Therefore, we concluded that TAD boundaries are enriched in housekeeping genes.

Cell 169, 216–228, April 6, 2017 219

A

B

D

C

Figure 4. Dynamics of Transcription and TAD Boundaries during Zygotic Genome Activation (A and B) Insulation scores and RNA Pol II binding at newly established boundaries. (C) Metaregion plot summaries for the data in (A) and (B). (D) Chromatin contact probability maps (5 kb resolution), insulation index heatmaps, and RNA Pol II binding for the Bsg25A/Elba3 locus (shaded area).

Our Hi-C maps displayed a prominent enrichment of interactions between adjacent TAD boundaries (Figures 1B and 3C). Because of the enrichment of housekeeping genes at TAD boundaries, the increased level of interaction could be indicative of spatial clustering of housekeeping genes in the three-dimensional nuclear space. To test this hypothesis, we computed the average level of inter-TAD-boundary interactions within 1 Mb distance and compared it with the level of interactions at equivalently spaced mirrored regions, to account for the strong bias due to decay of contacts away from the Hi-C diagonal (Figure 3C). This analysis revealed a significant enrichment of inter-TAD-boundary contacts starting at nc13 and increasing from nc14 onward, which are maintained at later stages of development and in Kc167 cells (Figures 3D, 3E, and S4). Therefore, our results support the emergence of three-dimensional spatial clustering of housekeeping gene-rich boundaries in the nucleus

220 Cell 169, 216–228, April 6, 2017

upon RNA Pol II recruitment and transcription initiation at the main wave of ZGA. Dynamics of TAD Formation We next focused on the dynamics of the establishment of TAD boundaries in relation with binding of RNA Pol II. To investigate this, we examined the presence of TAD boundaries and RNA Pol II binding during the transition between nc12 and 3–4 hpf embryos, which revealed that newly established TAD boundaries originated concomitantly with the recruitment of RNA Pol II (Figures 4A–4C). Quantification of the dynamics for the establishment of this insulation resembled that of RNA Pol II recruitment (Figure 4B). In addition, a small number of TAD boundaries (
5% of the total at nc14) do not associate with RNA Pol II, which hints at the presence of RNA Pol IIindependent TAD boundaries in the genome. Together these

A

C

B

D

E

F

G

Figure 5. Inhibition of Transcription Alters Chromatin Conformation (A) Schematic representation of the injection experiments. (B) Chromatin contact probability maps in untreated, water-treated, alpha-amanitin-treated, and triptolide-treated embryos each at nc14. Log2 fold changes in contact probability between the samples are displayed in between the Hi-C maps. (C) Insulation scores and difference in insulation relative to respective control conditions at TAD boundaries called for untreated nc14 embryos. (D) Insulation score at boundaries averaged across regions depicted in (C). (E and F) Distributions of contact probabilities of intra-TAD (E) and inter-TAD (F) interactions in control versus drug-treated embryos (n = 502). ***p value < 2.2 3 1016; Wilcoxon signed-rank test. (G) Aggregate Hi-C plots of boundary-boundary contacts at different distances for alpha-amanitin- and water-treated embryos as in Figure 3D. See also Figure S5 and Table S4.

results strongly associate the presence of RNA Pol II with the establishment of TAD boundary-like elements during development. The small number of early expressed loci prior to the major wave of ZGA that are rapidly switched off thereafter allowed us to test whether the loss of transcription would lead to changes in chromatin conformation. An example is the Bsg25A/Elba3 locus, which is expressed before nc14 only (Dai et al., 2015; Singer and Lengyel, 1997). This locus undergoes a significant remodelling of chromatin architecture when RNA Pol II binding is lost at nc14, resulting in the loss of its boundary-like structure (Figure 4D). This suggests that changes in RNA Pol II binding can be associated with a reorganization of chromatin conformation.

TAD Insulation and Boundary Co-localization Altered by Transcription Inhibition Next, we examined the role of RNA Pol II-mediated transcription in assembling chromatin organization. To do so, we injected embryos before the onset of zygotic transcription (< nc8) with the RNA Pol II pharmacological inhibitors alpha-amanitin or triptolide and with water as a control. We collected these embryos at nc14 and performed in situ Hi-C on them (Figures 5A and S5). RTqPCR quantification of embryonic transcripts showed that treated embryos were transcriptionally inactive at the time of collection (Figure S5B). Pan and phospho-Ser-5 RNA Pol II ChIP-seq in treated embryos showed that transcriptional elongation was severely reduced, while RNA Pol II association with promoters was slightly reduced upon treatment (Figure S5F

Cell 169, 216–228, April 6, 2017 221

and S5G). Consistent with these results, RNA Pol II ChIP qPCR quantification showed an overall decrease in the quantity of RNA Pol II bound to chromatin (Figure S5C). Remarkably, the resulting Hi-C maps from alpha-amanitinand triptolide-treated embryos displayed a significant degree of chromatin conformation, including TADs and extensive interTAD contacts (Figure 5B). In agreement with a recently proposed model (Li et al., 2015), because transcription is blocked in these embryos before they are transcriptionally active and prior to the establishment of TADs, our data demonstrate that the establishment of a TAD-based chromatin conformation is independent of transcription. However, after comparing the chromatin architecture at TAD boundaries, we found that the lack of transcription resulted in a significant loss of inter-TAD insulation (Figures 5C and 5D). Furthermore, the average contact density within TADs—a measure of their compactness—was significantly reduced, and inter-TAD interactions were stronger (Figures 5E and 5F). In addition, the strong enrichment of TAD boundary co-localization that we previously characterized in untreated nc14 embryos was severely affected in treated embryos (Figure 5G). These results imply that although higher-order chromatin conformation is still present in these embryos, the lack of transcription affects the properties of TAD organization. Chromatin Landscape at TAD Boundaries during Development We next investigated the chromatin landscape surrounding TAD boundary-like elements during development using stage-specific ChIP-seq data for 9 chromatin marks and 17 chromatinassociated proteins in Kc167 cells (see Table S2 for data sources). In agreement with previous results associating specific chromatin marks to early expressed genes (Li et al., 2014), we found that TAD boundaries at nc12 are enriched in H3K18ac, H3K27ac, and H4K8ac (Figures 6A and 6B). This hints at an association between these marks and early chromatin conformation events. In contrast, other marks of active chromatin such as H3K4me1, H3K4me3, H3K9ac, and H3K36me3 showed a marked enrichment only at later nuclear cycles, suggesting that they are not associated with the establishment of early topological nucleation sites (Figure 6B). TAD boundaries were also enriched in binding of Zelda—a zinc-finger protein associated with early zygotic transcription (Liang et al., 2008)—and architectural proteins such as BEAF-32, Cp190, Z4, Chromator, and dCTCF in Kc167 cells (Figure S6). Furthermore, TAD boundaries were highly enriched in open chromatin as measured by ATAC-seq data (Figure 6C). Similarly to that observed for RNA Pol II, a stage-specific chromatin opening occurs concomitantly with the establishment of TAD boundaries. On the contrary, marks of repressive chromatin, such as Polycomb-mediated H3K27me3, displayed a specific depletion in TAD boundary regions and enrichment toward the interior of TADs coinciding with the main wave of ZGA (Figure 6B). To more precisely understand the molecular mechanisms surrounding TAD boundary formation, we performed de novo motif enrichment analysis on the set of open chromatin regions at TAD boundaries. This analysis resulted in the identification of significantly enriched motifs in boundaries established at each nuclear

222 Cell 169, 216–228, April 6, 2017

cycle including those of the known architectural proteins BEAF32 and GAF (Figures 6D and 6E), as well as motifs with unknown binding partners. Consistent with the enrichment of housekeeping genes at TAD boundaries, we also found a significant enrichment of Motif 1 (bound by M1BP), which has been previously implicated in housekeeping gene regulation (Lam et al., 2012). Strikingly, the DNA binding motif of Zelda (Harrison et al., 2011) was the highest enriched motif at nuclear cycle 12, pinpointing a potential implication of Zelda in establishing TAD boundaries (Figure 6D). Zelda-Dependent Establishment of TAD Boundary Insulation To investigate the function of Zelda in the establishment of chromatin conformation, we performed in situ Hi-C on nc14 embryos obtained from a MTD-Gal4/UAS-shRNA-zld fly line, which has been previously shown to efficiently deplete Zelda during oogenesis and early embryonic development (Sun et al., 2015). Overall, Zelda-depleted embryos displayed TADs as well as TAD boundaries similar to those of nc14 embryos (Figure S7A). In contrast to our pharmacological transcription inhibition experiments, where transcription is altered globally, Zelda depletion affected RNA Pol II binding only at highly bound loci (Blythe and Wieschaus, 2015). Accordingly, Zelda-depleted embryos—despite undergoing developmental arrest before cellularization at nc14—displayed normal levels of insulation at most TAD boundaries and spatial co-localization of TAD boundaries (Figure S7B and S7G). However, Zelda-depleted embryos displayed distinct locus-specific defects, characterized by a loss of insulation at strong Zelda-bound loci genome-wide and an associated locus-specific loss of insulation at TAD boundaries (Figures 7A and 7B), which was more pronounced at early established boundaries (Figure S7F). Intriguingly, loci with mid to low levels of Zelda binding did not display a change in insulation properties when compared with control embryos, and regions with high Zelda binding lose insulation only partially, equivalent to levels observed for mid- to low-bound loci (Figure S7E), which suggest the presence of additional mechanisms that contribute to this insulation. Taken together, our results identify a role for Zelda in establishing insulation at TAD boundaries. DISCUSSION Our results provide strong evidence that chromatin conformation in Drosophila is established during zygotic genome activation in a transcription-independent manner. We find that, before the main wave of ZGA, the genome is mostly unstructured, except for
180 regions enriched in RNA Pol II binding and housekeeping genes that display TAD boundary-like insulation. These early transcribed loci serve as nucleation sites for the establishment of TAD boundaries, and inter-TAD insulation is progressively gained at loci marked for transcription (Figure 7C). These observations have implications for our understanding of how chromatin conformation is regulated that we discuss below. Establishment and Maintenance of Chromatin Conformation Previous studies have examined chromatin conformation in a range of cell lines and tissues across different species leading

A

B

Figure 6. TAD Boundaries Are Enriched in Active Chromatin Marks, Open Chromatin, and Binding Sites for Architectural Proteins and Zelda (A) ChIP-seq enrichment of histone modifications at boundary elements in nc12, nc13, and nc14. (B) Enrichment of chromatin marks at boundary elements. NS: not significant; *adjusted p value < 1 3 102; **adjusted p value < 1 3 103; ***adjusted p value < 1 3 104; permutation test; n = 10,000; dashed lines, mean of ChIP-seq signal in permutation test; gray shaded areas, 95% intervals. (C) Open chromatin at newly established boundaries measured using ATAC-seq signal from Blythe and Wieschaus (2016). (D) Sequence enrichment analysis at open chromatin regions at boundaries shown in (C). Displayed are the motif logo (left column), enrichment adjusted p value (middle column), and best transcription factor match, corresponding p value and expression level at 0–2 hpf according to Graveley et al. (2011) (right column). (E) Known binding sequence motifs of Zelda and architectural proteins BEAF-32 and GAF. See also Figure S6 and Table S2.

C

D

to the observation of TADs, which seem to be mostly invariant between cell lines and developmental stages and across evolutionary conserved loci (Dixon et al., 2012; Eagen et al., 2015; Ulianov et al., 2016; Vietri Rudan et al., 2015). This organization remains prevalent even when the presence of proteins associated with chromatin organization such as CTCF or cohesin is perturbed (Li et al., 2015; Seitan et al., 2013; Sofueva et al., 2013; Zuin et al., 2014). However, a common feature of these studies is that they were performed in cells actively engaged in transcribing their genomes.

In contrast, using Drosophila embryos allowed us to probe chromatin conformation in vivo in a unique nuclear setting, before the activation of the zygotic genome. The resulting in situ HiE C maps display a dramatic remodelling of the three-dimensional organization of the chromatin, coinciding with the main activation of the zygotic genome. Furthermore, we demonstrate that TAD boundaries throughout development are enriched in housekeeping genes, which have been previously found to be arranged in clusters in the genome (Weber and Hurst, 2011). Once established, TAD boundaries are maintained in vivo throughout development, in agreement with previous observations in cell lines (Ulianov et al., 2016). Given that these boundaries are enriched in housekeeping genes, which are actively transcribed in most cellular types in an organism, our results provide further evidence that the genomic arrangement of housekeeping genes in the genome is a major determinant for the conserved TAD organization found in most tissues (Ramı´rez et al., 2015; Ulianov et al., 2016). In line with this suggestion, previous studies have found an association between TAD boundaries and constitutively expressed genes that seems to be conserved throughout eukaryotic genomes ranging from yeast to mammals (Dixon et al., 2012; Hsieh et al., 2015). Interestingly, a recent examination of the

Cell 169, 216–228, April 6, 2017 223

A RNA Pol II nc13 sh zld

Zelda occupancy

fold-enrichment

RNA Pol II nc13

Figure 7. Loss of Insulation and TAD Boundaries at Zelda-Bound Loci

nc14

B

insulation score Zelda sh zld nc14 vs nc14 nc14

35

RNA Pol II nc14

0 50

RNA Pol II nc13

0 120

RNA Pol II nc13 sh zld

0 400

Zelda nc14

0 35

ATAC-seq nc13 18 min

0 sh zld

-100kb

0 +100kb

difference insulation score

ChIP ChIP fold-enrichment fold-enrichment Zelda RNA Pol II

-0.5 0 +0.5 0 more less insulated insulated

20

40 0

5

sh zld vs. nc14

10

centred on Zelda peaks

3R:14550kb 14700kb 14850kb 15000kb 15050kb contact probability -4

-3

-2

10 10 10 10 nc12

C

3-4h genome activation

?

?

? transcription inhibition

Boundaries located at early expressed genes RNA Pol II

?

?

?

Zelda depletion

Architectural proteins (BEAF-32, CP190, GAF, etc.) Zelda

?

relationship between dosage compensation and chromatin conformation in mice has also highlighted a strong functional link between transcription and the establishment of local chromatin architecture, including TAD-based conformation at escapee loci on the inactive X chromosome (Giorgetti et al., 2016). Together, these results suggest that while genomes of different sizes have evolved specific mechanisms to ensure the right three-dimensional organization of chromatin in the nuclear space (Crane et al., 2015; Hsieh et al., 2015; Rao et al., 2014; Van Bortle et al., 2014), the insulating properties found at transcribed loci are a conserved feature across evolution.

224 Cell 169, 216–228, April 6, 2017

(A) Insulation score difference between Zeldadepleted and control embryos at nc14 centered on Zelda peaks. RNA Pol II binding is shown for the same regions for wild-type and Zelda-depleted embryos at nc13 from Blythe and Wieschaus (2015). (B) Example showing the loss of a TAD boundary after Zelda depletion. Log2 fold changes in contact probability between the samples are displayed underneath the Hi-C maps. ChIP-seq fold enrichments for RNA Pol II and Zelda are also shown. (C) Model for the establishment of TADs and TAD boundaries. Prior to ZGA, the genome is mostly unstructured except for regions enriched in housekeeping genes that resemble TAD boundaries. During ZGA, TADs and inter-TAD contacts are established genome-wide in association with RNA Pol II and architectural proteins. TADs and TAD boundary formation is independent of transcription, although inhibition of transcription leads to a loss of insulation and spatial clustering at TAD boundaries. A subset of highly Zelda-bound early boundaries fails to establish proper insulation upon Zelda depletion. See also Figure S7 and Tables S2 and S4.

log2 fold change

Spatial Clustering of Housekeeping Genes An observation of our chromatin organization maps is the presence of contacts between TAD boundary regions at large genomic distances, which resemble the interactions of ‘‘A’’ compartments previously described for mammals (Lieberman-Aiden et al., 2009). These enriched contacts suggest that housekeeping, gene-rich boundary regions cluster together in the nucleus. These observations are reminiscent of transcription foci that have been observed previously by microscopy (Buckley and Lis, 2014; Cisse et al., 2013; Ghamari et al., 2013; Zhou et al., 2006) and are in agreement with the extensive promoter-promoter contacts between housekeeping genes found in mammals using ChIA-PET (Li et al., 2012). The long-range contacts between housekeeping genes are weaker in embryos before nc14. It is worth noting that, prior to nc14, nuclear divisions take place consecutively in a short period of time (
10–19 min each) (Blythe and Wieschaus, 2015), which might preclude the formation of those long-range interactions. Examinations of interaction dynamics in embryos with extended nuclear division timing (Ruden and Ja¨ckle, 1995) will help to determine the impact of short nuclear cycle divisions.

-1

-1.6

0

+1.6

Cause or Consequence: Mechanisms of TAD Formation The precise mechanisms by which chromatin conformation emerges during development remain to be fully unveiled. The

relationship between transcription and chromatin organization has been examined before, suggesting different levels of interdependency between the two (Branco and Pombo, 2006; Jin et al., 2013; Li et al., 2015; Mitchell and Fraser, 2008; Palstra et al., 2008). These studies, however, were unable to examine the role of transcription in establishing TADs and TAD boundaries since they were performed in cells with an already established chromatin conformation. By inhibiting transcription in embryos without prior apparent chromatin organization, we demonstrate that the establishment of TADs and TAD boundaries is independent of transcription and does not require zygotic factors. However, inhibition of RNA Pol II function reduced inter-TAD insulation and spatial clustering of housekeeping genes. Our ChIP data for RNA Pol II showed an almost complete absence of binding on gene bodies in both transcription inhibition experiments and reduced binding at promoters after treatment, which argues for a quantitative lack of RNA Pol II binding in these embryos. Given that both treatments can lead to degradation of RNA Pol II (Chen et al., 2015; Nguyen et al., 1996), the loss of insulation therefore might arise from a dosage-dependent lack of RNA Pol II occupancy at TAD boundaries. Together these data suggest that while the emergence of chromatin architecture is independent of transcription, mechanisms associated with this process, such as chromatin opening, the assembly of the preinitiation complex (PIC), or recruitment of RNA Pol II or other basic transcriptional machinery, might play a role in establishing and maintaining chromatin conformation. We also observe a small number of boundaries that are not associated with RNA Pol II, which indicates that other mechanisms might be involved in the formation of those boundaries. It should be noted that global inhibition of transcription might lead to secondary effects that affect the organization of chromatin. Intriguingly, binding of dCTCF in Drosophila does not seem to be as predictive of TAD boundaries as in mammals (Mourad and Cuvier, 2016; Ulianov et al., 2016). Therefore, in agreement with previous observations identifying insulation properties for RNA Pol II binding (Chopra et al., 2009), our results suggest RNA Pol II as a possible major determinant of inter-TAD insulation in Drosophila. Our analyses of the temporal dynamics of TAD boundary formation and the open chromatin state of TAD boundaries suggest that events preceding the recruitment of RNA Pol II might be important for the establishment of TAD boundaries. Of interest is the NSL complex, which has been previously associated with the eviction of nucleosomes and assembly of the PIC at promoters of housekeeping genes (Lam et al., 2012). Among the marks for active chromatin that are enriched at TAD boundaries, we find distinct temporal dynamics, suggesting that different mechanisms might be associated with the establishment and maintenance of these boundaries. The association with some of these marks might only be consequential and not causal—for example H3K36me3—given that chromatin organization still occurs upon inhibition of transcription, which impairs the establishment of this chromatin mark (Gerber and Shilatifard, 2003). Little is known about the role of two of the associated histone marks, namely H3K18ac and H4K8ac. Interestingly, H3K18ac is also enriched in yeast micro TAD boundaries (Hsieh et al., 2015), suggesting an evolutionarily conserved role. H3K18ac has been associated with the regulation of double-strand breaks and genome integrity

in mammals (Tasselli et al., 2016; Vazquez et al., 2016), which together with the recent association of topoisomerase (DNA) II beta (TOP2B) at TAD boundaries (Uusku¨la-Reimand et al., 2016) provide further evidence that TAD boundaries might be involved in the regulation of replication timing (Pope et al., 2014). Current models of TAD formation include the ‘‘loop-extrusion’’ model, based on the presence of cohesin-mediated loops and CTCF anchors (Fudenberg et al., 2016; Nichols and Corces, 2015; Sanborn et al., 2015), and the ‘‘self-assembly’’ model, based on the aggregation ability of nucleosomes from inactive chromatin (Ulianov et al., 2016), among others (Dixon et al., 2016). So far our analyses do not allow determining which model would be more plausible. The architectural proteins involved in the establishment of such conformation in the loop extrusion model are maternally provided and ubiquitously expressed thereafter (Graveley et al., 2011). Furthermore, according to computational simulations (Fudenberg et al., 2016), the rapid succession of nuclei divisions might not allow for the necessary time to form TAD-like structures at early stages of development. Irrespective of how chromatin is compacted into TADs, our analyses suggest that the genome is partitioned in the three-dimensional nuclear space into specific functional units consisting of TADs and inter-TAD regions, which might allow for fine-tuned regulation of developmentally regulated and housekeeping genes separately. Finally, our genetic analyses on the role of Zelda in establishing chromatin conformation demonstrate that this transcription factor is necessary for the formation of Zelda-specific TAD boundaries. Given the ability of Zelda to open chromatin, it is plausible that generating a relaxed local chromatin environment could be one of the initial steps in the establishment of chromatin conformation (Sun et al., 2015). Zelda binding is necessary for chromatin accessibility and transcription factor binding during ZGA (Schulz et al., 2015), and hence boundaries that lose insulation upon Zelda depletion might do so because of a combination of lack of architectural protein binding, binding of other transcription factors, or binding of RNA Pol II. Our analyses of open chromatin at TAD boundaries identify other enriched DNA motifs such as BEAF-32 and GAF, which, together with the observation of Zelda-independent TAD boundaries, suggest that other DNA binding proteins might act in a coordinated fashion to orchestrate the emergence of chromatin architecture. Of interest will be the novel DNA binding motifs that we identify here, as well as factors such as GAF, which has been shown to be able to open chromatin in a Zelda-independent manner (Schulz et al., 2015). Overall, our data provide new insight regarding when and how chromatin conformation is established during development. The powerful tools of Drosophila genetics such as maternal knockdowns and forward genetic screens combined with the recent development of genome editing tools such as CRISPR-Cas9 and targeted approaches to modify chromatin environments (Therizols et al., 2014; Wijchers et al., 2016) will allow for precise examinations of other factors implicated in the establishment of chromatin conformation. Therefore, explorations of the ZGA in Drosophila will constitute an ideal platform for examining the molecular mechanisms that determine how chromatin organization is regulated and how it affects gene expression during development and disease.

Cell 169, 216–228, April 6, 2017 225

STAR+METHODS

REFERENCES

Detailed methods are provided in the online version of this paper and include the following:

Bailey, T.L., Johnson, J., Grant, C.E., and Noble, W.S. (2015). The MEME Suite. Nucleic Acids Res. 43(W1), W39–W49.

d d d d

d

d

KEY RESOURCES TABLE CONTACT FOR REAGENT AND RESOURCE SHARING EXPERIMENTAL MODEL AND SUBJECT DETAILS METHOD DETAILS B Embryo collection and sorting B Embryo injections B Immunohistochemistry B Chromatin immunoprecipitation B RNA extraction and RT-qPCR B In situ Hi-C B Hi-C sequencing library preparation QUANTIFICATION AND STATISTICAL ANALYSIS B Hi-C data processing B Insulation score and boundary calling B Boundary overlap analysis and selection of newly established boundaries B Housekeeping gene enrichment analysis B Boundary-boundary contact analysis B ChIP-seq data analysis B ChIP-seq metaregion analysis B ATAC-seq analysis B Sequence motif enrichment analysis B Intra- and inter-TAD interaction analysis B Analysis of boundaries in Zelda depleted embryos DATA AND SOFTWARE AVAILABILITY

SUPPLEMENTAL INFORMATION Supplemental Information includes seven figures and four tables and can be found with this article online at http://dx.doi.org/10.1016/j.cell.2017.03.024. AUTHOR CONTRIBUTIONS J.M.V. conceived and supervised the study. J.M.V. and C.B.H. designed the experiments and analyses. A.G.G. performed RT-qPCR, ChIP-seq, and immunohistochemistry experiments. K.K. provided methods to analyze Hi-C data. C.B.H. performed all other experiments and computational analyses. J.M.V. and C.B.H. analyzed and interpreted the data. All authors participated in the discussion of the results. J.M.V. and C.B.H. wrote the manuscript. ACKNOWLEDGMENTS

Bickmore, W.A., and van Steensel, B. (2013). Genome architecture: domain organization of interphase chromosomes. Cell 152, 1270–1284. Blythe, S.A., and Wieschaus, E.F. (2015). Zygotic genome activation triggers the DNA replication checkpoint at the midblastula transition. Cell 160, 1169–1181. Blythe, S.A., and Wieschaus, E.F. (2016). Establishment and maintenance of heritable chromatin structure during early Drosophila embryogenesis. eLife 5, 5. Branco, M.R., and Pombo, A. (2006). Intermingling of chromosome territories in interphase suggests role in translocations and transcription-dependent associations. PLoS Biol. 4, e138. Buckley, M.S., and Lis, J.T. (2014). Imaging RNA Polymerase II transcription sites in living cells. Curr. Opin. Genet. Dev. 25, 126–130. Chen, F., Gao, X., and Shilatifard, A. (2015). Stably paused genes revealed through inhibition of transcription initiation by the TFIIH inhibitor triptolide. Genes Dev. 29, 39–47. Chopra, V.S., Cande, J., Hong, J.-W., and Levine, M. (2009). Stalled Hox promoters as chromosomal boundaries. Genes Dev. 23, 1505–1509. Cisse, I.I., Izeddin, I., Causse, S.Z., Boudarene, L., Senecal, A., Muresan, L., Dugast-Darzacq, C., Hajj, B., Dahan, M., and Darzacq, X. (2013). Real-time dynamics of RNA polymerase II clustering in live human cells. Science 341, 664–667. Crane, E., Bian, Q., McCord, R.P., Lajoie, B.R., Wheeler, B.S., Ralston, E.J., Uzawa, S., Dekker, J., and Meyer, B.J. (2015). Condensin-driven remodelling of X chromosome topology during dosage compensation. Nature 523, 240–244. Dai, Q., Ren, A., Westholm, J.O., Duan, H., Patel, D.J., and Lai, E.C. (2015). Common and distinct DNA-binding and regulatory activities of the BEN-solo transcription factor family. Genes Dev. 29, 48–62. De Renzis, S., Elemento, O., Tavazoie, S., and Wieschaus, E.F. (2007). Unmasking activation of the zygotic genome using chromosomal deletions in the Drosophila embryo. PLoS Biol. 5, e117. Dixon, J.R., Selvaraj, S., Yue, F., Kim, A., Li, Y., Shen, Y., Hu, M., Liu, J.S., and Ren, B. (2012). Topological domains in mammalian genomes identified by analysis of chromatin interactions. Nature 485, 376–380. Dixon, J.R., Jung, I., Selvaraj, S., Shen, Y., Antosiewicz-Bourget, J.E., Lee, A.Y., Ye, Z., Kim, A., Rajagopal, N., Xie, W., et al. (2015). Chromatin architecture reorganization during stem cell differentiation. Nature 518, 331–336. Dixon, J.R., Gorkin, D.U., and Ren, B. (2016). Chromatin domains: the unit of chromosome organization. Mol. Cell 62, 668–680. Eagen, K.P., Hartl, T.A., and Kornberg, R.D. (2015). Stable chromosome condensation revealed by chromosome conformation capture. Cell 163, 934–946.

This research was funded by the Max Planck Society. C.B.H. was supported by a fellowship from the International Max Planck Research School – Molecular Biomedicine. We thank Shelby Blythe and Eric Wieschaus for kindly providing the eGFP-PCNA Drosophila melanogaster line and Christine Rushlow for providing the UAS-shRNA-zld line; Georg Steffes and Christian Klaembt for support with Drosophila culturing; Celeste Brennecka for editing support; Maria Elena Torres-Padilla for critical reading of the manuscript; and the Genomics core facility of the Medical Faculty Muenster for sequencing. Stocks obtained from the Bloomington Drosophila Stock Center (NIH P40OD018537) were used in this study.

Flavahan, W.A., Drier, Y., Liau, B.B., Gillespie, S.M., Venteicher, A.S., Stemmer-Rachamimov, A.O., Suva`, M.L., and Bernstein, B.E. (2016). Insulator dysfunction and oncogene activation in IDH mutant gliomas. Nature 529, 110–114.

Received: September 21, 2016 Revised: February 21, 2017 Accepted: March 16, 2017 Published: April 6, 2017

Gerber, M., and Shilatifard, A. (2003). Transcriptional elongation by RNA polymerase II and histone methylation. J. Biol. Chem. 278, 26303–26306.

226 Cell 169, 216–228, April 6, 2017

Foe, V.E., and Alberts, B.M. (1983). Studies of nuclear and cytoplasmic behaviour during the five mitotic cycles that precede gastrulation in Drosophila embryogenesis. J. Cell Sci. 61, 31–70. Fudenberg, G., Imakaev, M., Lu, C., Goloborodko, A., Abdennur, N., and Mirny, L.A. (2016). Formation of chromosomal domains by loop extrusion. Cell Rep. 15, 2038–2049.

Ghamari, A., van de Corput, M.P.C., Thongjuea, S., van Cappellen, W.A., van Ijcken, W., van Haren, J., Soler, E., Eick, D., Lenhard, B., and Grosveld, F.G.

(2013). In vivo live imaging of RNA polymerase II transcription factories in primary cells. Genes Dev. 27, 767–777.

(2009). Comprehensive mapping of long-range interactions reveals folding principles of the human genome. Science 326, 289–293.

Ghavi-Helm, Y., Klein, F.A., Pakozdi, T., Ciglar, L., Noordermeer, D., Huber, W., and Furlong, E.E.M. (2014). Enhancer loops appear stable during development and are associated with paused polymerase. Nature 512, 96–100.

Lott, S.E., Villalta, J.E., Schroth, G.P., Luo, S., Tonkin, L.A., and Eisen, M.B. (2011). Noncanonical compensation of zygotic X transcription in early Drosophila melanogaster development revealed through single-embryo RNA-seq. PLoS Biol. 9, e1000590.

Giorgetti, L., Lajoie, B.R., Carter, A.C., Attia, M., Zhan, Y., Xu, J., Chen, C.J., Kaplan, N., Chang, H.Y., Heard, E., and Dekker, J. (2016). Structural organization of the inactive X chromosome in the mouse. Nature 535, 575–579. Graveley, B.R., Brooks, A.N., Carlson, J.W., Duff, M.O., Landolin, J.M., Yang, L., Artieri, C.G., van Baren, M.J., Boley, N., Booth, B.W., et al. (2011). The developmental transcriptome of Drosophila melanogaster. Nature 471, 473–479. Harrison, M.M., Li, X.-Y., Kaplan, T., Botchan, M.R., and Eisen, M.B. (2011). Zelda binding in the early Drosophila melanogaster embryo marks regions subsequently activated at the maternal-to-zygotic transition. PLoS Genet. 7, e1002266. Hnisz, D., Weintraub, A.S., Day, D.S., Valton, A.-L., Bak, R.O., Li, C.H., Goldmann, J., Lajoie, B.R., Fan, Z.P., Sigova, A.A., et al. (2016). Activation of proto-oncogenes by disruption of chromosome neighborhoods. Science 351, 1454–1458. Hou, C., Li, L., Qin, Z.S., and Corces, V.G. (2012). Gene density, transcription, and insulators contribute to the partition of the Drosophila genome into physical domains. Mol. Cell 48, 471–484. Hsieh, T.-H.S., Weiner, A., Lajoie, B., Dekker, J., Friedman, N., and Rando, O.J. (2015). Mapping nucleosome resolution chromosome folding in yeast by micro-C. Cell 162, 108–119. Jin, F., Li, Y., Dixon, J.R., Selvaraj, S., Ye, Z., Lee, A.Y., Yen, C.-A., Schmitt, A.D., Espinoza, C.A., and Ren, B. (2013). A high-resolution map of the threedimensional chromatin interactome in human cells. Nature 503, 290–294. Jullien, J., Astrand, C., Szenker, E., Garrett, N., Almouzni, G., and Gurdon, J.B. (2012). HIRA dependent H3.3 deposition is required for transcriptional reprogramming following nuclear transfer to Xenopus oocytes. Epigenetics Chromatin 5, 17. Knight, P.A., and Ruiz, D. (2012). A fast algorithm for matrix balancing. IMA J. Numer. Anal. 33, 1029–1047. Kruse, K., Hug, C.B., Herna´ndez-Rodrı´guez, B., and Vaquerizas, J.M. (2016). TADtool: visual parameter identification for TAD-calling algorithms. Bioinformatics 32, 3190–3192. Lam, K.C., Mu¨hlpfordt, F., Vaquerizas, J.M., Raja, S.J., Holz, H., Luscombe, N.M., Manke, T., and Akhtar, A. (2012). The NSL complex regulates housekeeping genes in Drosophila. PLoS Genet. 8, e1002736. Langmead, B., and Salzberg, S.L. (2012). Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359. Li, J., and Gilmour, D.S. (2013). Distinct mechanisms of transcriptional pausing orchestrated by GAGA factor and M1BP, a novel transcription factor. EMBO J. 32, 1829–1841. Li, G., Ruan, X., Auerbach, R.K., Sandhu, K.S., Zheng, M., Wang, P., Poh, H.M., Goh, Y., Lim, J., Zhang, J., et al. (2012). Extensive promoter-centered chromatin interactions provide a topological basis for transcription regulation. Cell 148, 84–98. Li, X.-Y., Harrison, M.M., Villalta, J.E., Kaplan, T., and Eisen, M.B. (2014). Establishment of regions of genomic activity during the Drosophila maternal to zygotic transition. eLife 3, 3. Li, L., Lyu, X., Hou, C., Takenaka, N., Nguyen, H.Q., Ong, C.-T., Cuben˜asPotts, C., Hu, M., Lei, E.P., Bosco, G., et al. (2015). Widespread rearrangement of 3D chromatin organization underlies polycomb-mediated stress-induced silencing. Mol. Cell 58, 216–231.

Lupia´n˜ez, D.G., Kraft, K., Heinrich, V., Krawitz, P., Brancati, F., Klopocki, E., Horn, D., Kayserili, H., Opitz, J.M., Laxova, R., et al. (2015). Disruptions of topological chromatin domains cause pathogenic rewiring of geneenhancer interactions. Cell 161, 1012–1025. Mitchell, J.A., and Fraser, P. (2008). Transcription factories are nuclear subcompartments that remain in the absence of transcription. Genes Dev. 22, 20–25. Mourad, R., and Cuvier, O. (2016). Computational identification of genomic features that influence 3D chromatin domain formation. PLoS Comput. Biol. 12, e1004908. Narendra, V., Rocha, P.P., An, D., Raviram, R., Skok, J.A., Mazzoni, E.O., and Reinberg, D. (2015). CTCF establishes discrete functional chromatin domains at the Hox clusters during differentiation. Science 347, 1017–1021. Naumova, N., Imakaev, M., Fudenberg, G., Zhan, Y., Lajoie, B.R., Mirny, L.A., and Dekker, J. (2013). Organization of the mitotic chromosome. Science 342, 948–953. Nguyen, V.T., Giannoni, F., Dubois, M.F., Seo, S.J., Vigneron, M., Ke´dinger, C., and Bensaude, O. (1996). In vivo degradation of RNA polymerase II largest subunit triggered by alpha-amanitin. Nucleic Acids Res. 24, 2924–2929. Nichols, M.H., and Corces, V.G. (2015). A CTCF code for 3D genome architecture. Cell 162, 703–705. Nora, E.P., Lajoie, B.R., Schulz, E.G., Giorgetti, L., Okamoto, I., Servant, N., Piolot, T., van Berkum, N.L., Meisig, J., Sedat, J., et al. (2012). Spatial partitioning of the regulatory landscape of the X-inactivation centre. Nature 485, 381–385. Ohler, U., Liao, G.C., Niemann, H., and Rubin, G.M. (2002). Computational analysis of core promoters in the Drosophila genome. Genome Biol. 3, H0087. Palstra, R.-J., Simonis, M., Klous, P., Brasset, E., Eijkelkamp, B., and de Laat, W. (2008). Maintenance of long-range DNA interactions after inhibition of ongoing RNA polymerase II transcription. PLoS ONE 3, e1661. Patel, N.H. (1994). Imaging neuronal subsets and other cell types in wholemount Drosophila embryos and larvae using antibody probes. Methods Cell Biol. 44, 445–487. Pope, B.D., Ryba, T., Dileep, V., Yue, F., Wu, W., Denas, O., Vera, D.L., Wang, Y., Hansen, R.S., Canfield, T.K., et al. (2014). Topologically associating domains are stable units of replication-timing regulation. Nature 515, 402–405. Ramı´rez, F., Lingg, T., Toscano, S., Lam, K.C., Georgiev, P., Chung, H.-R., Lajoie, B.R., de Wit, E., Zhan, Y., de Laat, W., et al. (2015). High-affinity sites form an interaction network to facilitate spreading of the MSL complex across the X chromosome in Drosophila. Mol. Cell 60, 146–162. Ramı´rez, F., Ryan, D.P., Gru¨ning, B., Bhardwaj, V., Kilpert, F., Richter, A.S., Heyne, S., Du¨ndar, F., and Manke, T. (2016). deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res. 44(W1), W160–W165. Rao, S.S.P., Huntley, M.H., Durand, N.C., Stamenova, E.K., Bochkov, I.D., Robinson, J.T., Sanborn, A.L., Machol, I., Omer, A.D., Lander, E.S., and Aiden, E.L. (2014). A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping. Cell 159, 1665–1680. Ruden, D.M., and Ja¨ckle, H. (1995). Mitotic delay dependent survival identifies components of cell cycle control in the Drosophila blastoderm. Development 121, 63–73.

Liang, H.-L., Nien, C.-Y., Liu, H.-Y., Metzstein, M.M., Kirov, N., and Rushlow, C. (2008). The zinc-finger protein Zelda is a key activator of the early zygotic genome in Drosophila. Nature 456, 400–403.

Sanborn, A.L., Rao, S.S.P., Huang, S.-C., Durand, N.C., Huntley, M.H., Jewett, A.I., Bochkov, I.D., Chinnappan, D., Cutkosky, A., Li, J., et al. (2015). Chromatin extrusion explains key features of loop and domain formation in wildtype and engineered genomes. Proc. Natl. Acad. Sci. USA 112, E6456–E6465.

Lieberman-Aiden, E., van Berkum, N.L., Williams, L., Imakaev, M., Ragoczy, T., Telling, A., Amit, I., Lajoie, B.R., Sabo, P.J., Dorschner, M.O., et al.

Schulz, K.N., Bondra, E.R., Moshe, A., Villalta, J.E., Lieb, J.D., Kaplan, T., McKay, D.J., and Harrison, M.M. (2015). Zelda is differentially required for

Cell 169, 216–228, April 6, 2017 227

chromatin accessibility, transcription factor binding, and gene expression in the early Drosophila embryo. Genome Res. 25, 1715–1726. Seitan, V.C., Faure, A.J., Zhan, Y., McCord, R.P., Lajoie, B.R., Ing-Simmons, E., Lenhard, B., Giorgetti, L., Heard, E., Fisher, A.G., et al. (2013). Cohesinbased chromatin interactions enable regulated gene expression within preexisting architectural compartments. Genome Res. 23, 2066–2077. Sexton, T., and Cavalli, G. (2015). The role of chromosome domains in shaping the functional genome. Cell 160, 1049–1059. Sexton, T., Yaffe, E., Kenigsberg, E., Bantignies, F., Leblanc, B., Hoichman, M., Parrinello, H., Tanay, A., and Cavalli, G. (2012). Three-dimensional folding and functional organization principles of the Drosophila genome. Cell 148, 458–472. Singer, J.B., and Lengyel, J.A. (1997). Expression and sequence analysis of the Drosophila blastoderm-specific gene bsg25A. Gene 197, 379–382. Sofueva, S., Yaffe, E., Chan, W.-C., Georgopoulou, D., Vietri Rudan, M., MiraBontenbal, H., Pollard, S.M., Schroth, G.P., Tanay, A., and Hadjur, S. (2013). Cohesin-mediated interactions organize chromosomal domain architecture. EMBO J. 32, 3119–3129. Staller, M.V., Yan, D., Randklev, S., Bragdon, M.D., Wunderlich, Z.B., Tao, R., Perkins, L.A., Depace, A.H., and Perrimon, N. (2013). Depleting gene activities in early Drosophila embryos with the ‘‘maternal-Gal4-shRNA’’ system. Genetics 193, 51–61. Sun, Y., Nien, C.-Y., Chen, K., Liu, H.-Y., Johnston, J., Zeitlinger, J., and Rushlow, C. (2015). Zelda overcomes the high intrinsic nucleosome barrier at enhancers during Drosophila zygotic genome activation. Genome Res. 25, 1703–1714. Tasselli, L., Xi, Y., Zheng, W., Tennen, R.I., Odrowaz, Z., Simeoni, F., Li, W., and Chua, K.F. (2016). SIRT6 deacetylates H3K18ac at pericentric chromatin to prevent mitotic errors and cellular senescence. Nat. Struct. Mol. Biol. 23, 434–440. Therizols, P., Illingworth, R.S., Courilleau, C., Boyle, S., Wood, A.J., and Bickmore, W.A. (2014). Chromatin decondensation is sufficient to alter nuclear organization in embryonic stem cells. Science 346, 1238–1242. Ulianov, S.V., Khrameeva, E.E., Gavrilov, A.A., Flyamer, I.M., Kos, P., Mikhaleva, E.A., Penin, A.A., Logacheva, M.D., Imakaev, M.V., Chertovich, A., et al. (2016). Active chromatin and transcription play a key role in chromosome partitioning into topologically associating domains. Genome Res. 26, 70–84. Uusku¨la-Reimand, L., Hou, H., Samavarchi-Tehrani, P., Rudan, M.V., Liang, M., Medina-Rivera, A., Mohammed, H., Schmidt, D., Schwalie, P., Young, E.J., et al. (2016). Topoisomerase II beta interacts with cohesin and CTCF at topological domain borders. Genome Biol. 17, 182.

228 Cell 169, 216–228, April 6, 2017

Van Bortle, K., Nichols, M.H., Li, L., Ong, C.-T., Takenaka, N., Qin, Z.S., and Corces, V.G. (2014). Insulator function and topological domain border strength scale with architectural protein occupancy. Genome Biol. 15, R82. Vazquez, B.N., Thackray, J.K., Simonet, N.G., Kane-Goldsmith, N., MartinezRedondo, P., Nguyen, T., Bunting, S., Vaquero, A., Tischfield, J.A., and Serrano, L. (2016). SIRT7 promotes genome integrity and modulates nonhomologous end joining DNA repair. EMBO J. 35, 1488–1503. Vietri Rudan, M., Barrington, C., Henderson, S., Ernst, C., Odom, D.T., Tanay, A., and Hadjur, S. (2015). Comparative Hi-C reveals that CTCF underlies evolution of chromosomal domain architecture. Cell Rep. 10, 1297–1309. Weber, C.C., and Hurst, L.D. (2011). Support for multiple classes of local expression clusters in Drosophila melanogaster, but no evidence for gene order conservation. Genome Biol. 12, R23. Wijchers, P.J., Krijger, P.H.L., Geeven, G., Zhu, Y., Denker, A., Verstegen, M.J.A.M., Valdes-Quezada, C., Vermeulen, C., Janssen, M., Teunissen, H., et al. (2016). Cause and consequence of tethering a SubTAD to different nuclear compartments. Mol. Cell 61, 461–473. Wood, A.M., Van Bortle, K., Ramos, E., Takenaka, N., Rohrbaugh, M., Jones, B.C., Jones, K.C., and Corces, V.G. (2011). Regulation of chromatin organization and inducible gene expression by a Drosophila insulator. Mol. Cell 44, 29–38. Zabidi, M.A., Arnold, C.D., Schernhuber, K., Pagani, M., Rath, M., Frank, O., and Stark, A. (2015). Enhancer-core-promoter specificity separates developmental and housekeeping gene regulation. Nature 518, 556–559. Zeitlinger, J., Stark, A., Kellis, M., Hong, J.-W., Nechaev, S., Adelman, K., Levine, M., and Young, R.A. (2007). RNA polymerase stalling at developmental control genes in the Drosophila melanogaster embryo. Nat. Genet. 39, 1512–1516. Zhang, Y., Liu, T., Meyer, C.A., Eeckhoute, J., Johnson, D.S., Bernstein, B.E., Nusbaum, C., Myers, R.M., Brown, M., Li, W., and Liu, X.S. (2008). Modelbased analysis of ChIP-Seq (MACS). Genome Biol. 9, R137. Zhou, G.-L., Xin, L., Song, W., Di, L.-J., Liu, G., Wu, X.-S., Liu, D.-P., and Liang, C.-C. (2006). Active chromatin hub of the mouse alpha-globin locus forms in a transcription factory of clustered housekeeping genes. Mol. Cell. Biol. 26, 5096–5105. Zuin, J., Dixon, J.R., van der Reijden, M.I.J.A., Ye, Z., Kolovos, P., Brouwer, R.W.W., van de Corput, M.P.C., van de Werken, H.J.G., Knoch, T.A., van IJcken, W.F.J., et al. (2014). Cohesin and CTCF differentially affect chromatin architecture and gene expression in human cells. Proc. Natl. Acad. Sci. USA 111, 996–1001.

STAR+METHODS KEY RESOURCES TABLE

REAGENT or RESOURCE

SOURCE

IDENTIFIER

Antibodies Anti-RNA Pol II pSer5

Abcam

cat#ab5408; RRID: AB_304868

Anti-pan RNA Pol II

Abcam

cat#ab103968

Anti-H3

Abcam

cat#ab1791; RRID: AB_302613

Anti-Rabbit IgG Alexa Fluor 488

Invitrogen

cat#A11029; RRID: AB_2534088

Anti-Mouse IgG Alexa Fluor 568

Invitrogen

cat#A11004; RRID: AB_2534072

Chemicals, Peptides, and Recombinant Proteins Alpha-amanitin

Sigma-Aldrich

cat#A2263

Triptolide

Selleck Chemicals

cat#S3604

Biotin-14-dATP

Life Technologies

cat#19524016

MboI

New England Biolabs

cat#R0147L

DNA Polymerase I Klenow Fragment

New England Biolabs

cat#M0210L

T4 DNA Ligase

Thermo Fisher

cat#EL0012

T4 DNA Polymerase

New England Biolabs

cat#M0203L

Proteinase K

AppliChem

cat#A4392

GlycoBlue

Life Technologies

cat#AM9516

Complete Ultra EDTA-free protease inhibitors

Roche

cat#05892791001

TRIzol Reagent

Thermo Fisher

cat#15596018

Normal Goat Serum

Sigma-Aldrich

cat#G9023

Dynabeads MyOne Streptavidin C1

Life Technologies

cat#65002

AMPure XP beads

Beckman Coulter

cat#A63881

RNase-free DNase I

Qiagen

cat#79254

Critical Commercial Assays Magna ChIP A/G Chromatin Immunoprecipitation Kit

Millipore

cat#17-10085

NEBNext Ultra II DNA Library Prep Kit

New England Biolabs

cat#E7645

SuperScript III First-Strand Synthesis SuperMix for qRT-PCR kit

Thermo Fisher

cat#11-752-050

MESA GREEN qPCR MasterMix Plus for SYBR

Eurogentec

cat#RT-SY2X-20+WOU

RNeasy Mini

QIAGEN

cat#74104

Deposited Data Hi-C from staged embryos

this paper

ArrayExpress: E-MTAB-4918

Hi-C from Kc cells

Li et al. (2015)

GEO: GSE63515

Hi-C from 16-18hpf embryos

Sexton et al. (2012)

GEO: GSE34453

RNA Pol II ChIP-seq reads from injected staged embryos

this paper

ArrayExpress: E-MTAB-4918

RNA Pol II ChIP-seq from staged embryos

Blythe and Wieschaus (2015)

GEO: GSE62925

Histone ChIP-seq from staged embryos

Li et al. (2014)

GEO: GSE58935

Zld ChIP-seq from 2-3hpf embryos

Sun et al. (2015)

GEO: GSE65441

Zld ChIP-seq from staged embryos

Harrison et al. (2011)

GEO: GSE30757

ATAC-seq from staged embryos

Blythe and Wieschaus (2016)

GEO: GSE83851

BEAF ChIP-seq from Kc cells

Li et al. (2015)

GEO: GSE62904

CapH2 ChIP-seq from Kc cells

Van Bortle et al. (2014)

GEO: GSE54529

CBP ChIP-seq from Kc cells

Li et al. (2015)

GEO: GSE62904

Chromator ChIP-seq from Kc cells

Li et al. (2015)

GEO: GSE62904

Barren ChIP-seq from Kc cells

Li et al. (2015)

GEO: GSE62904 (Continued on next page)

Cell 169, 216–228.e1–e9, April 6, 2017 e1

Continued REAGENT or RESOURCE

SOURCE

IDENTIFIER

CP190 ChIP-seq from Kc cells

Li et al. (2015)

GEO: GSE62904

CTCF ChIP-seq from Kc cells

Van Bortle et al. (2014)

GEO: GSE54529

DREF ChIP-seq from Kc cells

Li et al. (2015)

GEO: GSE62904

GAF ChIP-seq from Kc cells

Van Bortle et al. (2014)

GEO: GSE54529

IIC220 ChIP-seq from Kc cells

Van Bortle et al. (2014)

GEO: GSE54529

L3mbt ChIP-seq from Kc cells

Li et al. (2015)

GEO: GSE62904

Modmdg4 ChIP-seq from Kc cells

Li et al. (2015)

GEO: GSE62904

Rad21 ChIP-seq from Kc cells

Li et al. (2015)

GEO: GSE62904

Su(Hw) ChIP-seq from Kc cells

Wood et al. (2011)

GEO: GSE30740

TFIIIC ChIP-seq from Kc cells

Li et al. (2015)

GEO: GSE62904

Z4 ChIP-seq from Kc cells

Li et al. (2015)

GEO: GSE62904

FlyBase RNA-seq profile

Graveley et al. (2011)

http://flybase.org/static_pages/rna-seq/ rna-seq_profile_search.html

Nuclear marker strain w[*]; P{EGFP-PCNA}attP2

Laboratory of Eric Wieschaus (Blythe and Wieschaus, 2016)

N/A

RNAi of Zelda w[*]; P{UAS-shRNA-zld}

Laboratory of Christine Rushlow (Sun et al., 2015)

N/A

Maternal triple driver P{otu-GAL4::VP16.R}1, w[*]; P{GAL4-nos.NGT}40; P{ GAL4::VP16-nos.UTR}CG6325

Bloomington Drosophila Stock Center

Stock number 31777; RRID: BDSC_31777

shRNA-zld passenger strand sequence CGGATGCAAGTTGCAGTGCAA

Sun et al. (2015)

N/A

Primers for RT-qPCR see Table S4

this paper; Synthesized by Eurofins Genomics

N/A

Langmead and Salzberg (2012)

http://bowtie-bio.sourceforge.net/ bowtie2/index.shtml

BBMap (v35.85)

N/A

https://sourceforge.net/projects/bbmap/

MACS2 (v2.1.0)

Zhang et al. (2008)

https://github.com/taoliu/MACS/

deeptools (v2.4.0)

Ramı´rez et al. (2016)

https://github.com/fidelram/deepTools/

MEME Suite (v4.11.2.2)

Bailey et al. (2015)

http://meme-suite.org

this paper

https://github.com/vaquerizaslab/ Hug-et-al-Cell-2017-Supp-Site

Experimental Models: Organisms/Strains

Oligonucleotides

Software and Algorithms Bowtie2 (v2.2.6)

Other Resource website for this publication

CONTACT FOR REAGENT AND RESOURCE SHARING Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Juan M. Vaquerizas ([email protected]). EXPERIMENTAL MODEL AND SUBJECT DETAILS Drosophila melanogaster fly stocks were maintained at room temperature with natural light/dark cycle. Flies were kept in roundbottom vials (4 3 11.5 cm) filled with standard corn-meal agar food. At least one day before embryo collections, three vials with flies not older than one week were moved to a collection cage (5.6 3 7.6 cm; Flystuff.com) and maintained at 25 C in incubators with 12 hr light/dark cycle. yw; eGFP-PCNA flies were kindly provided by Blythe S.A. and Wieschaus E.F (Blythe and Wieschaus, 2016). Briefly, the proliferating cell nuclear antigen (PCNA) locus was amplified by PCR from genomic DNA and fused with an EGFP cassette. Fly lines containing the fusion construct were established by phiC31-mediated insertion into the attP2 landing site.

e2 Cell 169, 216–228.e1–e9, April 6, 2017

To knock down zld in oocytes, we used the Maternal-Gal4-shRNA system (Staller et al., 2013). The UAS-shRNA-zld line we obtained from the Rushlow laboratory and the MTD-Gal4 line from the Bloomington Stock Center (stock number 31777). The zld deficient embryos were generated as described before (Sun et al., 2015). Briefly, we crossed UAS-shRNA-zld virgins with MTD-Gal4 males to produce heterozygous UAS-shRNA-zld/MTD-Gal4 offspring. These were mated with their siblings to generate embryos deficient in zld. METHOD DETAILS Embryo collection and sorting We adapted the fixation and sorting procedure described in (Blythe and Wieschaus, 2015) for in situ Hi-C (Rao et al., 2014). Following a pre-collection period of at least one hour, fly embryos were collected on yeasted 0.4% acetic acid agar plates at 25 C. After 30 min of collection, the embryos on the plate were incubated at 25 C until they reached the desired developmental stage. Embryos were dechorionated for 2 min in 2.5% sodium hypochlorite. After rinsing with water, embryos were suspended in 2 mL of PBS, 0.5% Triton X-100 and 6 mL of heptane. Cross-linking was initiated by adding 100 mL of 37% formaldehyde, followed by vigorous shaking. 10 min later, samples were spun at 500 g for 1 min and embryos were collected from the lower aqueous layer. 15 min after the start of fixation, 5 mL of PBS, 0.5% Triton X-100, 125 mM glycine was added to the embryos, followed by vigorous shaking for 1 min. The embryos were rinsed 3 times with PBS, 0.5% Triton X-100. Subsequent sorting was carried out on ice in the same solution including cOmplete Ultra protease inhibitors (Roche). We sorted embryos in small batches of approximately 100 embryos at a time under a Leica fluorescent binocular equipped with a GFP filter set. We used the fluorescent signal from the eGFP-PCNA fusion protein to determine the developmental stage and cell cycle status of each embryo. The developmental stage was readily established by examining the nuclear density, which roughly doubles with each nuclear division cycle. Embryos undergoing a mitotic division cycle were easily identifiable by the lack of nuclear eGFP-PCNA signal. Embryos for 3-4 hpf were sorted solely based on their morphology. We removed all embryos with inappropriate developmental stage or mitotic nuclei as well as damaged embryos or such with abnormal morphology. The remaining sorted population of embryos was imaged for quality control. Using this sorting procedure, we obtained > 99% nuclear-cycle specific populations of embryos for our experiments. After sorting, embryos were pooled so that a single tube contained enough embryos for one experiment. The whole supernatant was removed, the embryos were flash-frozen in liquid nitrogen and stored at 80 C. For each in situ Hi-C experiment we used a range between 20 embryos (for stages after cellularisation) and 80 embryos (for nuclear cycle 12 embryos). Embryo injections Embryos expressing the eGFP-PCNA transgene were collected and dechorionated as described above. Immediately following collection, embryos were inspected under a fluorescent binocular and embryos older than nuclear cycle 6 were discarded. For each batch, approximately 30-40 embryos were lined up individually on an agar block, transferred to a coverslip covered with heptane glue, air-dried for 3 min and covered with a drop of Voltalef 10S oil (VWR International). Injections were carried out using an Eppendorf Femtojet microinjection device. Alpha-amanitin (Sigma Aldrich) was dissolved in water at a concentration of 0.5 mg/ml. Triptolide (Biozol GmbH) was dissolved in DMSO at 1 mg/ml and diluted with water to 0.05 mg/ml for the injections. Each embryo was injected with approximately 0.2 nL of alpha-amanitin solution, triptolide solution or water as a control. Given that the average volume of a Drosophila embryo is 15 nL (Foe and Alberts, 1983) and assuming that the injected material disperses evenly within the embryo, we expect a final concentration of approximately 6.6 mg/ml alpha-amanitin or 1.8 mM triptolide in injected embryos. If embryos were older than nuclear cycle 8 at the time of injection they were discarded. After injection, embryos were incubated at 25 C and regularly monitored. Once the first embryos reached nuclear cycle 14 (which represent the oldest embryos at the time of injection), five embryos of each batch were removed from the coverslip and solubilized in 300 mL of Trizol reagent (Life Technologies) for RNA extraction and quality control. Once the remaining embryos reached nuclear cycle 14, they were flushed off the coverslip using heptane and formaldehyde-fixed using the same procedure as described above. Immunohistochemistry Immunohistochemistry (IHC) was performed as described in (Patel, 1994). Briefly, embryos were fixed for 20 min in a 1:1 solution of heptane and PBST (PBS-triton 0.1%) formaldehyde 4%, devitelinized by removing the aqueous phase, and adding an equal volume of methanol and strong agitation. Embryos were finally washed 3 times in fresh methanol and stored in methanol at 20 C. When used for IHC experiments, embryos were progressively rehydrated in PBST and subsequently washed for 30 min using the same solution. To block non-specific antibody binding, embryos were incubated for 30 min in PBS-NGS (PBS-triton 0.1%, BSA 0.1%, normal goat serum (NGS) 5%; G9023 Sigma). Primary antibodies (rabbit anti-H3: ab1791 Abcam; mouse anti-RNA Pol II pSer5: ab5408 Abcam) were then directly added to a final dilution of 1:1000, and embryos were incubated overnight at 4 C with light agitation. The following day, embryos were washed 4 times for 30 min each time with PBST, 30 min with PBS-NGS, and the secondary antibodies (Alexa Fluor 488 Goat Anti-Rabbit IgG; Alexa Fluor 568 Goat Anti-Mouse IgG; Invitrogen) were added to a final dilution of 1:100. After a 2 hr incubation at RT, embryos were washed 4 times 30 min with PBST, 1 hr with PBS-glycerol 50% DAPI 1 mg/ml, resuspended in DABCO solution (glycerol 70%; DABCO 2.5%), and mounted on slides. Image acquisition was performed using a LSM880 confocal microscope (Zeiss) equipped with a plan-apochrom 63x/1.40 oil DIC M27 objective. Laser lines at 488 and

Cell 169, 216–228.e1–e9, April 6, 2017 e3

561 nm were used, and the excitation light was reflected with an MBS 488/561 beam splitter. Filters allowed detection of the emission light at ranges of 508-535 and 589-618 nm respectively from the 488 nm and 561 excitation wavelengths. The microscope hardware was controlled by the Zen2 software (Zeiss). Chromatin immunoprecipitation Crosslinking and chromatin preparation from nuclear cycle 14 embryos was carried out as described in (Blythe and Wieschaus, 2015). After injection and storage at 80 C, 150 to 200 embryos where mechanically homogenized in 700 mL of RIPA (150 mM NaCl, 1.0% IGEPAL CA-630, 0.5% sodium deoxycholate, 0.1% SDS, 50 mM Tris-Cl, pH 8.0) + protease inhibitors (Sigma 04693159001) using a pellet pestle. After centrifugation for 5 min at 16000 g and 4 C, the supernatant was discarded taking care to pipette out lipids at the surface and not touch the pellet. The pellet was then resuspended in 130 mL of fresh RIPA + protease inhibitors and sonicated on a Covaris S2 sonicator at 5% duty cycle for 10 min. Following sonication, the samples were centrifuged at 16000 g for 5 min at 4 C. The supernatant (
120 ml), containing the chromatin, was finally carefully removed and split in 2x 50 mL aliquots. The remaining 20 mL were decrosslinked overnight at 61 C with proteinase K 1 mg/ml final and loaded on a 1.5% agarose gel with 1X GelRed (Biotium #41003), migrated for 1.5 hr at 65V, and imaged to control that the size of the fragments was between 150 and 500 bp. During the whole procedure, the samples were constantly kept on ice to minimize protein degradation. Chromatin immunoprecipitation was performed using commercial kits (17-10085 Millipore) in 1.5 mL DNA LoBind Eppendorf tubes (Eppendorf) as follows. For each ChIP, 450 mL of dilution buffer was supplemented with 2.25 mL of protease inhibitor cocktail, and added to a 50 mL chromatin aliquot. For each sample, 5 mL (corresponding to 1% of the total volume) is set apart and stored at 4 C for later processing. These aliquot contain non-immunoprecipitated chromatin, which will serve as input during the processing of the ChIP-Seq and the analysis of the qPCR. To each 500 mL reaction volume, 20 mL of magnetic beads were added. The antibodies (anti-RNA Pol II pSer5: ab5408, Abcam; anti-pan RNA Pol II: ab103968, Abcam) were respectively used at a final dilution of 6:1000 and 18:1000 (3 mL or 9 mL in a 500 mL volume, respectively). The samples were then incubated overnight at 4 C on a rotating wheel. The next day, the beads were washed twice with the low salt buffer, twice with the high salt buffer, once with the LiCl buffer, and once with TE. These more stringent washes allow for better ChIP enrichment. Both the beads and the 5 mL aliquots previously set aside were finally resuspended in 100 mL of elution buffer + proteinase K 1 mg / ml final and incubated overnight at 61 C. The following day, the volume of the samples were adjusted to 400 mL with TE and DNA purified with an equal volume of phenol:chloroform:isoamyl alcohol 25:24:1 (Sigma) in a phase-lock-gel tube (5prime). After centrifugation at 12000 g for 5 min at RT, the supernatant (400 ml) was transferred to a new tube. The DNA was precipitated by addition of 40 mL of sodium acetate 3M pH 5, 2 mL of glycogen 20 mg/ml and 640 mL of 20 C ethanol 100%. The samples were then incubated for 1h at 20 C and centrifuged at 16000 g for 5 min at 4 C. Finally, the DNA pellets were washed twice with 800 mL 20 C ethanol 70% and resuspended in 20 mL TE. For ChIP-qPCR, samples were diluted to 1:10. The sequence of primers used for these experiments are presented in Table S4. For qPCR analysis, abundance of immunoprecipitated material is expressed as percentage of input recovered. ChIP-Seq libraries were prepared directly from immunoprecipitated material using NEBNext Ultra II DNA Library Prep Kits (E7645, NEB) as follows. The immunoprecipitated material volume was adjusted to 50 mL with Tris-Cl pH 8.0 and transferred to a 0.2 mL PCR strip tube. To the reaction volume, 3 mL of End Prep Enzyme Mix and 7 mL of End Prep Reaction buffer were added. The mixture was incubated in a thermocycler running the following program: 20 C for 30 min, 65 C for 30 min, 4 C hold. 30 mL of Ligation Master Mix, 2.5 mL Adaptor for Illumina (1:15 diluted from stock) and 1 mL of Ligation enhancer were added to the sample and the mixture was incubated in a thermocycler at 20 C for 15 min. Next, 3 mL of USER enzyme was added and the mixture incubated at 37 C for 15 min in a thermocycler. Because of the low amount of initial material (2-5 ng), as recommended by the manufacturer, a size-selection step was not performed and fragments were cleaned up to eliminate unligated adaptors by adding 87 mL of Ampure XP beads (Beckman Coulter) and mixing it with the reaction volume. After a 5 min incubation at RT, the beads were separated on a magnetic stand and the supernatant containing the short fragments was discarded. The beads were then washed twice with 200 mL of ethanol 80%, air-dried for 5 min and resuspended in 15 mL of 10 mM Tris-Cl pH 8.0. After separating the beads on a magnetic stand, the supernatant was transferred to a new 0.2 mL PCR. For each library, a PCR reaction was prepared as follows: 25 mL of 2x NEBNext Ultra II Q5 Master Mix, 5 mL of 25 mM Universal PCR primer, 5 mL of 25 mM Indexed PCR primer and finally 15 mL of the ChIP-Seq library. The PCR reactions were run using the following program: 98 C for 30 s, (98 C for 10 s, 65 C for 75 s, ramping 1.5 C/s) repeated 12 times, 65 C for 5 min, 4 C hold. The PCR product was size-selected using Ampure XP beads. The PCR reactions were combined and adjusted to exactly 110 mL volume with water. 110 mL of beads were mixed with the sample and incubated at room temperature for 5 min. After separation on a magnetic stand, the supernatant was transferred to a new tube and the beads -containing fragments too large for sequencing- were discarded. Another 40 mL of beads was added to the sample and incubated at room temperature for 5 min. After separation, the supernatant was discarded and the beads washed twice with 700 mL of 80% ethanol, each time incubating for 30 s. Following complete removal of the ethanol, the beads were resuspended in 100 mL of 10 mM Tris-Cl pH 8.0 and incubated at room temperature for 1 min. The following second round of selection removes all traces of primers and possible primer dimers. Without separating, another 80 mL of fresh beads was added and the suspension was incubated at room temperature for 5 min. After separation of the beads, the supernatant was discarded and the beads were washed twice with 80% ethanol, as above. After all traces of ethanol were removed, the beads were resuspended in 20 mL of 10 mM Tris-Cl pH 8.0 and incubated at RT for 5 min. Finally, the supernatant, containing the final library was transferred to a 1.5 mL DNA LoBind Eppendorf tube. The libraries were sequenced on Illumina MiSeq (2x84 bp paired-end) instruments. Each ChIP-Seq was performed using two biological replicates and an average sequencing depth of 4 million reads per sample.

e4 Cell 169, 216–228.e1–e9, April 6, 2017

RNA extraction and RT-qPCR To validate the effect of alpha-amanitin and triptolide on RNA Pol II elongation, the relative amounts of transcripts for the early zygotic genes bnb, bsg25a, CG10035 and term, and the maternally provided gene rp49 were measured by RT-qPCR. Total RNA was extracted from frozen injected samples as follows. Embryos stored in 500 mL Trizol were pestle homogenized. 1/5 volumes (100 ml) of chloroform (Fluka) was added. After vigorous shaking and a brief incubation, the samples were centrifuged at 12,000 g for 15 min at 4 C. To purify the RNA, the supernatant was collected in a new tube and an equal volume of 100% ethanol was added. After vortexing, the whole volume was transferred to an RNeasy Mini (QIAGEN) column. After a brief 8000 g centrifugation, the RNA was washed once with 700 mL of RW1 buffer. To remove traces of contaminating genomic DNA, the columns were treated with 10 mL of DNaseI in 80 mL of Buffer RDD for 15 min at RT. 350 mL of buffer RW1 was then directly added and the columns were centrifuged for 30 s at 8000 g. Next, the columns were washed twice with 500 mL of buffer RPE, and RNA was finally eluted from the column in 14 mL of RNase free water. The quality and integrity of the RNA was assayed by loading 1 mL of each sample on a RNA 6000 Nano chip (Agilent) and visual inspection of the obtained electropherograms, which did not reveal any apparent degradation (presence of two strongly resolved peaks at around 2000 nt with no smear at lower molecular weights). Because the Drosophila 28S rRNA is processed into two fragments that individually migrate similarly to the 18S rRNA, RIN values could not be obtained. The concentrations of the samples were measured using Qubit RNA HS (ThermoFisher), and ranged from 1.23 to 2.92 ng/ml. Reverse transcription was carried out using the SuperScript III First-Strand Synthesis SuperMix for qRT-PCR kit (ThermoFisher) on the samples adjusted to 1.2 ng/ml as follows. A reaction volume containing 10 mL of 2x RT Reaction Mix, 2 mL Enzyme Mix and 8 mL of purified RNA was prepared on a 0.2ml PCR strip tube and incubated on a thermocycler at 25 C for 10 min, 50 C for 30 min and 85 C for 5 min. Finally, cDNA were diluted 1:10 with RNase free water. The RT-qPCR itself was carried out in triplicate for each sample in 96 well plates. Each reaction volume contained 1 mL of 10 mM primers (forward and reverse; Table S4), 5.5 mL of nuclease free water, 12.5 mL of MESA GREEN qPCR MasterMix Plus for SYBR (Eurogentec) and 5 mL of 1:10 template cDNA, in a total final volume of 25 ml. Cycling conditions were 95 C for 10 min followed by 40 cycles of 95 C for 30 s and 60 C for 60 s. rp49 was used as an endogenous control. Ct values were obtained from the raw data by using the auto-thresholding function of the SDS 1.3 software. Data analysis was performed by computing the delta-Ct using rp49 as reference. For each primer set, the delta-Ct expression levels were normalized by dividing by the mean of the expression level in the water injected samples. Thus the normalized expression levels average to 1 in the water-injected samples, and the values in the alpha-amanitin and triptolide treated samples can be interpreted as fraction of expression relative to the water-injected embryos. In situ Hi-C We followed the general procedure for in situ Hi-C as described in (Rao et al., 2014), which we adapted for use with fly embryos. Frozen embryos were placed on ice and suspended in 500 mL of ice-cold lysis buffer (10 mM Tris-Cl pH 8.0, 10 mM NaCl, 0.2% IGEPAL CA-630, cOmplete Ultra protease inhibitors). Using a tight-fitting metal pestle (Carl-Roth P985.1), the embryos were ground until a homogeneous suspension was achieved. After 15 min incubation on ice, the nuclei were spun down at 1000 g, 4 C for 5 min. The supernatant was removed and the nuclei pellet resuspended in 500 mL of lysis buffer. Following another spin, the supernatant was discarded and the pellet resuspended in 100 mL 0.5% SDS. The nuclei were permeabilised by incubating at 65 C for 10 min. The SDS was then quenched by adding 50 mL of 10% Triton X-100 and incubating at 37 C for 15 min. Chromatin was digested by adding 25 mL of CutSmart buffer and 50 U of MboI (New England Biolabs), followed by three hours of incubation at 37 C. The enzyme was added in two instalments of 25 U in the beginning and after 90 min of incubation. MboI was heat-inactivated by incubating at 62 C for 20 min. The overhangs left by MboI were filled in by adding 18.75 mL of 0.4 mM biotin-14-dATP (Life Technologies), 2.25 mL of a 3.3 mM dCTP/dGTP/dTTP mixture and 40 U of DNA Polymerase I Klenow Fragment (New England Biolabs), followed by incubation at 37 C for 90 min. The DNA fragments were then ligated by adding 657 mL of water, 120 mL of 10x T4 DNA Ligase buffer (Thermo Fisher Scientific), 100 mL of 10% Triton X-100, 6 mL of 20 mg/ml BSA and 25 Weiss U of T4 DNA Ligase (Thermo Fisher Scientific). The mixture was rotated slowly at room temperature for four hours. A second instalment of 25 Weiss U of T4 DNA Ligase was added after two hours. The nuclei were then spun down at 2500 g for 5 min, the supernatant discarded and the pellet resuspended in 500 mL of extraction buffer (50 mM Tris-Cl pH 8.0, 50 mM NaCl, 1 mM EDTA, 1% SDS). To digest the protein, 20 mL of 20 mg/ml Proteinase K (Applichem) was added and the mixture incubated at 55 C, 1000 rpm shaking for 30 min. 130 mL of 5 M sodium chloride was added and the mixture was incubated at 68 C, 1000 rpm shaking for at least 8 hr. For DNA precipitation, 63 mL of 3 M sodium acetate pH 5.2, 2 mL of 15 mg/ml GlycoBlue (Life Technologies) and finally 1008 mL of absolute ethanol were mixed with the sample, followed by incubation at 80 C for 15 min. The sample was then spun at 20,000 g, 4 C for one hour. The supernatant was removed and the DNA pellet washed two times using 800 mL of 70% ethanol each. All traces of remaining supernatant were removed; the pellet air-dried for 2 min and then solubilized in 50 mL of 10 mM Tris-Cl pH 8.0. RNA was digested by adding 1 mL of 20 mg/ml RNase A and incubating at 37 C for 15 min. Hi-C sequencing library preparation Biotin-14-dATP was removed from the ends of unligated fragments by mixing 50 mL of the Hi-C library with 12 mL of 10x NEB2 buffer (New England Biolabs), 3 mL of 1 mM dATP, 3 mL of 1 mM dGTP, 0.6 mL of 20 mg/ml BSA, 15 U of T4 DNA Polymerase (New England Biolabs) and water to 120 mL total. The mixture was incubated at 20 C for 30 min and then stopped by adding 3 mL of 0.5 M EDTA. The DNA was sheared to
200-400 bp with a Covaris S2 sonicator using the following program: 2 cycles of 50 s, 10% duty, intensity 5, Cell 169, 216–228.e1–e9, April 6, 2017 e5

200 cycles/burst. For the preparation of magnetic beads for biotin pull-down, 30 mL of Dynabeads MyOne Streptavidin C1 (Life Technologies) were separated on a magnetic stand. The supernatant was discarded. The beads were washed once and then resuspended in 120 mL of 1x B&W buffer (2x stock: 10 mM Tris-Cl pH 7.4, 1 mM EDTA, 2 M NaCl). 120 mL of the sheared Hi-C sample was mixed with 120 mL of the prepared magnetic beads and the suspension was rotated at room temperature for 15 min. The beads were washed twice with 600 mL each of 1x B&W buffer + 0.1% Triton X-100 by incubating at 55 C, 1000 rpm shaking for 2 min. This was followed by a wash with 10 mM Tris-Cl pH 8.0. The beads were then resuspended in 50 mL of the same solution. Next, library preparation for Illumina sequencing was performed using components from the NEBNext Ultra II DNA Library Prep Kit for Illumina (New England Biolabs). 3 mL of End Prep Enzyme Mix and 7 mL of End Prep Reaction buffer were combined in a 0.2 mL PCR strip tube. The bead suspension in 50 mL Tris-Cl pH 8.0 was mixed into the End Prep mixture and incubated in a thermocycler running the following program: 20 C for 30 min, 65 C for 30 min, 4 C hold. 30 mL of Ligation Master Mix, 2.5 mL Adaptor for Illumina (1:10 diluted from stock) and 1 mL of Ligation enhancer were added to the sample and the mixture was incubated in a thermocycler at 20 C for 15 min. Next, 3 mL of USER enzyme was added and the mixture incubated at 37 C for 15 min in a thermocycler. The beads were separated on a magnetic stand and the supernatant was discarded. The beads were washed twice using 1x B&W + 0.1% Triton X-100 and transferred to a 1.5 mL tube. The beads were washed once with 10 mM Tris-Cl pH 8.0 and then resuspended in 50 mL of the same solution. Next, the final Hi-C library was amplified using PCR. Determining the optimal number of amplification cycles requires empirical testing. In our experience 9-12 cycles are sufficient. If libraries required more amplification, they were usually of low quality and did not yield satisfactory results due to excessive duplication of reads. For each library, two PCR reactions were prepared as follows: 25 mL of 2x NEBNext Ultra II Q5 Master Mix, 1.5 mL of 10 mM Universal PCR primer, 1.5 mL of 10 mM Indexed PCR primer and finally 22 mL of the Hi-C library bead suspension. The PCR reactions were run using the following program: 98 C for 1 min, (98 C for 15 s, 65 C for 75 s, ramping 1.5 C/s) repeated 9-12 times, 65 C for 5 min, 4 C hold. The PCR product was size-selected using Ampure XP beads (Beckman Coulter). The PCR reactions were combined and adjusted to exactly 110 mL volume with water. 110 mL of beads were mixed with the sample and incubated at room temperature for 5 min. After separation on a magnetic stand, the supernatant was transferred to a new tube and the beads -containing fragments too large for sequencing- were discarded. Another 40 mL of beads was added to the sample and incubated at room temperature for 5 min. After separation, the supernatant was discarded and the beads washed twice with 700 mL of 80% ethanol, each time incubating for 30 s. Following complete removal of the ethanol, the beads were resuspended in 100 mL of 10 mM Tris-Cl pH 8.0 and incubated at room temperature for 1 min. The following second round of selection removes all traces of primers and possible primer dimers. Without separating, another 80 mL of fresh beads was added and the suspension was incubated at room temperature for 5 min. After separation of the beads, the supernatant was discarded and the beads were washed twice with 80% ethanol, as above. After all traces of ethanol were removed, the beads were resuspended in 50 mL of 10 mM Tris-Cl pH 8.0 and incubated at 37 C for 5 min. Following separation, the supernatant contains the final in situ Hi-C library including Illumina sequencing adapters. Libraries were sequenced on Illumina MiSeq (2x84 bp paired-end) and NextSeq 500 (2x80 bp paired-end) instruments. All experiments were performed using at least two biological replicates and an average sequencing depth of > 150 million reads per sample (Table S1). QUANTIFICATION AND STATISTICAL ANALYSIS Hi-C data processing First, each end of the paired-end reads was scanned separately for possible MboI ligation junctions (GATCGATC) and trimmed at these sites to improve mappability. Each half of the trimmed mate pairs were mapped individually onto the Drosophila melanogaster reference genome (FlyBase r6.07) using Bowtie 2 (version 2.2.6) (Langmead and Salzberg, 2012) with the following parameters: ‘-–very-sensitive -–no-unal’. To increase the number of mappable reads, we follow an iterative mapping strategy: First, reads were trimmed to 20 base pairs (bp) and mapped to the reference genome. If the read did not align uniquely, or its mapping quality was lower than 30, the read was extended by 5 bp in the 30 direction and sent for mapping again. This extension process was repeated until reads mapped uniquely with a mapping quality of at least 30 or until the full length of the read was reached, in which case the best unique alignment was used. Reads that mapped to multiple loci with the same score (multimapping reads) were removed. Additionally, reads with a mapping quality of less than 10 were filtered. The filtered set of aligned read pairs was then assigned to MboI restriction fragments. Inward facing pairs, i.e., pairs that aligned head-to-head on opposite strands, were filtered if they were less than 500 bp away from each other, since they likely represent undigested, unligated fragments. Outward-facing pairs, i.e., pairs oriented tail-to-tail on opposite strands, were likewise filtered if they mapped less than 500 bp apart, as they likely represent undigested, self-ligated fragments (Jin et al., 2013). We considered two read pairs to be PCR duplicates if both pairs’ forward and both pairs’ reverse reads aligned within 1 bp of each other. If multiple read pairs were found to be duplicates only one pair was retained. To generate contact maps, we divided the reference genome into equally sized, contiguous bins and assigned the aligned read pairs to these bins. This resulted in a square, symmetric contact matrix. Bins with less than 10% of median number of assigned read pairs per bin were masked and not considered in the downstream analysis, since they represented regions with low mappability, which often resulted in artifacts during matrix normalization. These are displayed as gray lines in our Hi-C maps. In order to reduce biases arising from GC content, fragment size, mappability, or other unspecified sources, the contact matrices were balanced using the Knight-Ruiz matrix balancing algorithm (Knight and Ruiz, 2012), as described for Hi-C data in (Rao et al., 2014). Using this approach the Hi-C maps are normalized so that interactions for each bin sum up to one, and therefore the values

e6 Cell 169, 216–228.e1–e9, April 6, 2017

can be interpreted as contact probabilities. Principal component analysis (Figure S2A) and comparison of boundary calls between biological replicates (Figure S3F) revealed a high similarity between them, which prompted us to merge replicates for subsequent analyses to maximize resolution, unless stated otherwise. Hi-C datasets were merged by combining the sets of aligned and filtered reads before generating contact matrices. The representation of Hi-C maps is shown at 10kb resolution unless otherwise stated. Insulation score and boundary calling We used the insulation score metric introduced by (Crane et al., 2015) to assess the extent of insulation genome-wide. Briefly, a square window is moved along the diagonal of the contact matrix and the average probability of contacts within the window is recorded (Figure S3A). If C is the binned contact matrix and w is the window size in bins, then the insulation score I for bin number i can be calculated as follows: P Cðj; kÞ i+1 % j < i+w+1 Iði; wÞ =

iw % k < i

w2

For boundary calling, we used balanced contact matrices at 5kb resolution and a window size of 8 bins throughout this work. To normalize the insulation score, it was divided by the central moving average (300 bins window) of the insulation score and converted to log2 scale. Boundaries were called using the delta method, as described in (Crane et al., 2015). Briefly, the delta vector D was calculated from the insulation score vector I for position i and the window size w using the following formula: P P IðkÞ  iw % k < i IðkÞ Dði; wÞ = i + 1 % k < i + w + 1 w The zero crossings of the delta vector represent local minima (putative boundaries) or local maxima (putative TAD centers). For each putative boundary, we determined the local minimum of the delta vector on the right and the local maximum to the left of the boundary bin. These were calculated by analyzing the first and second derivative of the Savitzky–Golay filtered delta vector (window size 7 bins, polynomial degree 2). The difference in delta vector at these two local extrema surrounding a boundary was used as a boundary score. The stepwise increase in TAD boundary calls during development is irrespective of the threshold used or the resolution of the Hi-C map employed (Figure S3G). Throughout this work we call boundaries using Hi-C contact maps binned at 5 kb resolution and requiring a minimum boundary score of 0.5 for calling boundaries. This threshold was selected to provide a tradeoff between sensitivity for detecting weak TAD-boundaries and specificity. We noticed that when we decreased the boundary threshold the number of boundaries increases rapidly (Figure S3G). Manual review of these boundary calls, revealed a large proportion of false-positives (as determined by visual inspection) when a threshold lower than 0.5 was used. When the threshold approaches 0, even minute dips in the insulation score, which are not readily detectable visually, represent minima that will pass the boundary threshold. While these might correspond to interesting unidentified features, we decided to use a more conservative threshold that resulted in robust TAD boundary calls. The sensitivity of our approach is highlighted by the fact that we can reproduce a large proportion of the boundary calls for c12 in Hi-C maps for mitotic embryos (Figure S3I), even though TAD boundaries are not apparent by eye in the contact maps. Similarly, we decided to use boundary calls made using Hi-C maps binned at 5kb resolution for our analysis. We reasoned that the gain in resolution by using 2kb bins was outweighed by the robustness of the boundary calls at 5kb resolution. However, to ensure the reproducibility of our main findings at different resolutions, we computed the enrichment of RNA Pol II at boundaries called using 2kb Hi-C maps and found very similar or higher enrichments (Table S3). For each developmental stage, a consensus set of boundaries was created using the merged Hi-C maps, combining all biological replicates for each stage. Insulation plots for different window sizes were computed as described in (Kruse et al., 2016). Wilcoxon signed-rank tests were used to compare changes in insulation at TAD boundaries between control and transcription inhibited embryos (p value < 2.2e-16 for both treatments; Figure 5D). A list of boundary calls is available at https://github.com/vaquerizaslab/Hug-et-al-Cell-2017-Supp-Site. Boundary overlap analysis and selection of newly established boundaries To analyze the overlap of boundaries between biological replicates or between developmental stages, we considered a sliding window of 20kb, since we estimated that the precision of boundary calls was limited to approximately ± 10kb. Boundaries that intersected a common 20kb sliding window were considered to overlap. To test the statistical significance of the overlaps we assume a universe of 3,347 possible boundaries (26,779 Hi-C bins divided by 4 bins overlap window at 5kb resolution). Fisher’s Exact tests were used to assess statistical significance (adjusted p value < 2.2e-16 between all consecutive stages). The total number of TAD boundary calls per sample can change slightly in each comparison since boundaries overlapping a common window surrounding 20kb are merged by default among all the samples being compared. This results in 3-4 hpf embryos displaying different total number of TAD boundaries in Figures 2B (850) and 2C (851). For each nuclear cycle a set of newly established boundaries was determined by removing all boundaries called in the preceding nuclear cycles. Housekeeping gene enrichment analysis Positions of housekeeping genes were obtained from Flybase (Graveley et al., 2011), including genes with at least moderate (RPKM > 11) expression level at all assayed developmental times and tissues. Housekeeping enhancers were defined using

Cell 169, 216–228.e1–e9, April 6, 2017 e7

processed STARR-seq data from (Zabidi et al., 2015). Odds ratios and p values were calculated using Fisher’s exact test on a contingency table of housekeeping versus not housekeeping genes and boundary versus no boundary with Bonferroni correction. Boundary-boundary contact analysis In order to assess the extent of contacts between boundaries, we used aggregate peak analysis as described in (Rao et al., 2014). We used boundary calls from the insulation score method. To do so, for each particular reference boundary (red arrow; Figure 3C), we investigated the contact density between neighboring boundaries (labeled as ‘1’; Figure 3C) and contacts between boundaries separated by one or more intermediate boundaries (labeled as ‘2’; Figure 3C). The analysis is extended for subsequent boundaries up to a distance of 1Mb from the reference boundary. All pairwise combinations of boundaries that were within 1Mb of each other were considered for these analyses. We did not extend the analysis further since beyond this distance the contacts were not significant above background level. For each pair, we extracted a sub-matrix from the contact map corresponding to the crossing point of the two boundaries including an extra 50kb on each side of the boundary. Sub-matrices were classified according to the number of intermediate boundaries separating the boundary pair (ordinal boundary-boundary distance). For all sub-matrices at a given distance, sub-matrices were aligned and the average value at each pixel calculated. As a control, we used sub-matrices mirroring boundary pairs (Figure 3C). The mirror-contact corresponds to the region between the first boundary in the pair and a control region upstream, which is the same distance to the boundary as in the pair being examined. This strategy minimizes biases based on the distance decay of contact probability as well as biases arising from the insulating properties of boundary elements. The signal for the first two columns is dominated by the strong contact frequency observed in TAD corners, which masks signal from boundary-boundary interactions at these distances despite the contacts being visible in the Hi-C maps. Mann Whitney U tests were used to compare changes in contacts between boundary and mirror-regions (Figure 3E) or changes in boundaries-boundary interactions between control and transcription inhibited embryos (p value < 2.2e-16 for both treatments; Figure 5G). In the box plots, the lower and upper limits of the box represent the 25th and 75th percentiles. Notches correspond to the 95% confidence interval of the median. The bar in the middle of the box corresponds to the median. The lower and upper whiskers extend to the highest value that is within ± 1.5 3 IQR of the limits of the boxes, respectively. ChIP-seq data analysis ChIP-seq sequencing reads obtained from SRA were mapped to the Flybase r6.07 Drosophila reference genome using BBMap (version 35.85, https://sourceforge.net/projects/bbmap/) with the following parameters: ‘slow –k = 12 sam = 1.3 mappedonly’. Mapped reads with a mapping quality below 10 were discarded. Using the aligned ChIP and matched input reads, peaks were then called using MACS2 (version 2.1.0; default parameters) (Zhang et al., 2008) in the narrow peak mode for RNA Pol II and architectural protein peaks. For the ChIP-seq analysis in Figure 2E, only the 1,000 strongest RNA Pol II peaks by the adjusted p value of the peak calls were included. For nuclear cycle 12 only 287 peaks were called. RNA Pol II ChIP-seq data for 3-4hpf embryos were plotted on a different scale in order to compensate for the overall lower ChIP enrichment in this dataset. For all histone modifications the broad peak mode was used. Unless otherwise stated, we used a q-value cutoff of 10 for calling significant peaks. Fold enrichment and log2 fold enrichment of ChIP versus input was calculated using the MACS2 bdgcmp command. The stalling index of RNA Pol II bound transcripts was calculated as described in (Zeitlinger et al., 2007). Briefly, the average RNA Pol II occupancy on promoter regions (200 bp upstream until 200 bp downstream of the TSS) and on gene bodies (200 bp downstream to 600 bp downstream of the TSS) was calculated for each transcript. The stalling index is then calculated as the ratio between promoter signal and gene body signal. Higher values representing stalled polymerase and lower values elongating polymerase. For this study we considered genes elongating if the stalling index was lower than 1, the promoter overlapped with a significant RNA Pol II peak and the entire gene body has full ChIP-seq coverage. A list of all ChIP-seq datasets used in this study and their references can be found in Table S2. Throughout the manuscript ChIP-seq data presented in heatmaps are centered according to the indicated feature in the panel. ChIP-seq metaregion analysis To generate metaregions plots of histone modification ChIP-seq signals surrounding boundary elements, the fold enrichment signal was averaged over 1kb bins. The signal vectors were aligned at the center of the boundary elements and averaged across boundaries. For factors with a smaller footprint, such as RNA Polymerase II and architectural proteins, we instead divided the genome into contiguous 1kb bins, and counted for each bin whether one or more significant peaks were overlapping that bin. In order to generate a null hypothesis for the expected distribution of signal, the positions of the boundaries were randomized 10,000 times and the ChIPseq signal was computed for each randomized set of boundaries, as described above. In each round of randomizations, in order to preserve the relative spatial distribution of boundaries, the positions of all boundaries were shifted by the same random amount while treating the chromosomes as circular. The distribution of signal from this set of randomized boundaries was then used as a null-hypothesis for statistical testing. To test the enrichment of features at boundaries, we computed the ratio between the ChIP-seq signal in a 7.5kb window around the center of the boundary and two 10kb flanking regions 30kb upstream and downstream of the center of the boundary, which was used as a proxy for the local background signal. This ratio was also calculated for the set of randomized boundaries. Permutation tests were then used to calculate p values for the significance of enrichments at boundaries. The p values reported were corrected for multiple comparisons across all enrichment tests using the Benjamini-Hochberg procedure. Metaregion plots summarizing RNA Pol II ChIP-seq signal in Figure 4C were produced using a 5% trimmed mean of the ChIP-seq enrichment

e8 Cell 169, 216–228.e1–e9, April 6, 2017

across all boundaries in the set to reduce the impact of extreme outliers. The metagene plots in Figures S5F and S5G were generated using deeptools (Ramı´rez et al., 2016) using all genes whose promoter region overlapped a significant RNA Pol II peak at nuclear cycle 14. Similarity between the enrichment of chromatin marks at boundaries was computed by applying hierarchical clustering using Euclidian distance and average linkage (Figure 6A). ATAC-seq analysis Processed genome tracks for open chromatin and the respective peak calls for each nuclear cycle were downloaded from GEO (GSE83851; Table S2) (Blythe and Wieschaus, 2016). The ATAC-seq data are available for multiple tightly staged time points at nc11, nc12 and nc13 subdivided in three-minute steps after each round of mitosis. As for ChIP-data, the ATAC-seq open chromatin tracks were plotted for a region of 80kb centered on TAD boundary calls newly established at each nuclear cycle using a 5kb bin size to smooth the data. Metaregion plots where generated using a 5% trimmed mean, as with the ChIP metaregion plots. Sequence motif enrichment analysis Sequence motif analysis was carried out on the open chromatin peak calls from the ATAC-seq data (Blythe and Wieschaus, 2016) using tools from the MEME suite (Bailey et al., 2015). For each nuclear cycle, the set of newly established boundaries was intersected with all ATAC-seq peaks at the same nuclear cycle (for boundaries at nuclear cycle 14, ATAC-peak calls from nuclear cycle 13 were used as a proxy). Only ATAC-seq peaks that overlapped at region of 20kb around TAD boundary calls were used for the analysis. The genomic sequences underlying the set of boundary associated ATAC-peaks were retrieved and de-novo sequence motif enrichment analysis was carried out using DREME with default settings and using dinucleotide scrambled versions of the input sequences as control. To compare the de novo motifs with known motifs, we used Tomtom on the set of the 8 most significantly enriched motifs for each nuclear cycle using the following databases included in MEME: OnTheFly_2014; Fly Factor Survey, FLYREG v2; and iDMMPMM. Additionally, we manually compared the core promoter consensus motifs described in (Ohler et al., 2002) with our de novo motifs using Tomtom, since these are not included in the databases used by MEME. This resulted in significant matches with Ohler’s motif 1 that has been previously reported to be targeted by the Motif 1 binding protein (M1BP) (Li and Gilmour, 2013). We report the most significant hit found in these databases for each of the de novo motifs that we called using DREME, unless the top hit was not expressed (measured as an RPKM value of zero at 0-2 hpf according to Flybase (Graveley et al., 2011)) in which case the search was extended for the top 3 hits. This occurred for only one motif, namely the sixth de novo motif hit found for nuclear cycle 13 embryos, which we report as BEAF-32. The top hit in this case was a motif reportedly bound by the unexpressed transcription factor pnr. For each de novo motif found for nc12, nc13 and nc14 embryos, we report the motif logo (left column), the Bonferroni adjusted p value of the motif enrichment (middle column), and the name and p value of the best match of the motif with entries in the transcription factor databases and expression state at 0-2 hpf according to (Graveley et al., 2011) (right column) (Figure 6D). Intra- and inter-TAD interaction analysis For intra- and inter-TAD interaction analysis, TADs with sizes between 50kb and 1Mb were considered (n = 502). For intra-TAD interactions, the mean probability of all contacts within TADs was used. For inter-TAD interactions, the mean probability of contacts between neighboring TADs was employed. In the box plots, the lower and upper limits of the box represent the 25th and 75th percentiles. Notches correspond to the 95% confidence interval of the median. The bar in the middle of the box corresponds to the median. The lower and upper whiskers extend to the highest value that is within ± 1.5 3 IQR of the limits of the boxes, respectively. Analysis of boundaries in Zelda depleted embryos In order to characterize the changes in chromatin architecture in Zelda depleted embryos, we divided the genome into 5kb bins matching the bins of the Hi-C matrix and aggregated ChIP-seq data for Zelda at nuclear cycle 14 (Harrison et al., 2011) for these bins. Specifically, we calculated the sum of the -log10 adjusted p values of all peaks that were contained in each bin. These aggregated data were used to generate a ranking of the highest Zelda bound regions, which we then used to plot the insulation score of the Hi-C maps around all Zelda bound sites, ranked by Zelda binding. Each 5kb bin is depicted only once in Figure 7A, even if it contained multiple Zelda peaks. To get an estimate of the number of boundaries that lose insulation upon Zelda depletion, we calculated the difference in insulation score at the boundary bins for wild-type and Zelda depleted embryos. The standard deviation of this difference was 0.2. Any boundaries whose insulation score differed by more than 1.5 standard deviations (0.3) from its wild-type value were considered to be affected by the Zelda depletion. DATA AND SOFTWARE AVAILABILITY The accession number for the in situ Hi-C data and RNA Pol II ChIP-seq data from injected embryos reported in this paper is ArrayExpress: E-MTAB-4918. The accession numbers of additional datasets used in this study are listed in Table S2.

Cell 169, 216–228.e1–e9, April 6, 2017 e9

Supplemental Figures

Figure S1. Hand Sorting of Embryos using the eGFP-PCNA Fusion Protein, Related to Figure 1 (A) Fluorescent microscopy images of DAPI-counterstained embryos from a fly line expressing the eGFP-PCNA transgene. S-phase embryo (top row). Embryo undergoing synchronous nuclear division (bottom row). Signal from the DAPI reveals the location of nuclei. eGFP-PCNA is present in the nuclei during S-phase but dispersed throughout the syncytial blastoderm during M-phase. (B to D) Representative images of embryo populations after hand sorting for nuclear cycle 12 (B), 13 (C) and 14 embryos (D). Bars represent 250 mm. Insets show magnified views of highlighted regions.

Figure S2. Global Characterization and Quality Control of Contact Maps, Related to Figure 1 (A) Principal component analysis (PCA) of all Hi-C datasets generated in this study. All intra-chromosomal contact bins of loci between 50kb and 1Mb from each other were used as input for the PCA. We used corrected contact maps binned at 10kb resolution. Each point represents a single replicate in the two-dimensional projection of the first and second principal components. The proportion of variance explained by each component is displayed underneath. (B) Contact probabilities plotted against genomic distance. The mean contact probability at each genomic distance is shown for the combined contact maps from each developmental time point. This gives an overview of how likely contacts at a certain distance are in each sample. Note the markedly different distribution of contact probabilities between mitotic and non-mitotic samples. In mitotic nuclei, chromatin is much more likely to engage in long-range contacts over 150kb. (C) Contact maps for a 3Mb region on chromosome arm 3R binned at 10kb resolution. Shown are our own data for 3-4h embryos, previously published Hi-C data from Kc167 cells (Li et al., 2015) and Hi-C data from 16-18h embryos (Sexton et al., 2012). (D) Scatterplots comparing contact probabilities derived from our own in situ Hi-C data (x axis) with previously published Hi-C data (y axis). We used corrected maps binned at 10kb resolution. Only contacts at distances between 50kb and 1Mb were considered. The Pearson correlation coefficient is displayed at the top left of each panel. All correlations were significant at p < 2.2e-16. (E) Corrected Hi-C contact maps (2kb resolution) at 3-4 hpf. Shown are three loci containing mesoderm specific enhancers (center), which were previously shown to engage in long-range genomic contacts in 3-4 and 6-8 hpf embryos using 4C-seq (Ghavi-Helm et al., 2014) (data were obtained from http://furlonglab.embl.de/ 4CBrowser). The 4C signals for each stage are plotted in the bottom panels with the position of the bait sequence indicated by the black box. The 4C-seq signal is defined as the normalized number of reads in the 4C-seq library and is roughly proportional to the contact probability with the given bait locus. The baits are CRM_6190 (left), CRM_en (center) and CRM_4319 (right), as described in (Ghavi-Helm et al., 2014). Since the 4C-seq signal is disproportionally high directly at the bait locus, the signal is truncated at the value indicated by the zigzag line.

Figure S3. Calculation of Insulation Scores and TAD Boundary Calls, Related to Figure 2 (A) Corrected contact map for a 2Mb region on chromosome 3R. The black outlined squares represent the area of the contact map that is used for the calculation of the insulation score for the bins indicated by the dashed line. The square is slid along the genome to calculate scores for every bin. (B) Insulation score values for the displayed region using a window size of 8 bins (40kb). Higher insulation scores reflect areas of highly interacting loci, such as intra-TAD regions, whereas lower insulation scores reflect areas with strong insulating properties, such as TAD boundary regions or inter-TADs. The local minima of the insulation score were used to call boundary elements genome-wide. (C) Heatmap of insulation scores for different window sizes (y axis units correspond to number of bins in the Hi-C map). (D) Violin plots representing the insulation scores genome-wide for each merged sample. Black bars represent the values of the quartiles for each distribution. Levene’s tests were used to evaluate whether variances were significantly different between conditions. In order to avoid overestimating differences due to the

(legend continued on next page)

interdependence of neighboring insulation scores, we sampled the insulation score every 16 Hi-C contact matrix bins (twice the window size of 8 bins). ***p value < 2.2e-16; Levene’s test for equality of variances. Bars within the plot correspond to the 25th, 50th, and 75th percentiles. (E) Interquartile ranges (IQR) of genome-wide insulation scores for each Hi-C replicate (points). Black bars indicate the mean. The insulation IQR serves as a measure of overall level of structural variation in the Hi-C contact maps. (F) Venn diagrams of overlap between boundary calls for all replicates in each experimental condition. Boundaries were considered to overlap if their positions varied by no more than 20kb in linear distance. (G) Number of TAD-boundary calls for different Hi-C map resolutions and boundary call thresholds. (H) Overlap of TAD boundary calls for Hi-C maps at different resolutions. (I) Overlap of TAD boundary calls between nuclear cycle 12, 13 and 14, and mitotic embryos. (J) Comparison of insulation scores at stalled and elongating RNA Pol II loci (Wilcoxon rank-sum test). (K) Distribution of RNA Pol II binding peaks on normalized TAD length between boundary calls. The histogram represents the normalized density of RNA Pol II binding at TAD boundaries and the regions in between.

distance

1

2

3

4

5

7

8

9

ND

ND

ND

ND

boundarymirror

ND

ND

ND

ND

c12

6

boundaryboundary

c13

boundaryboundary boundarymirror

c14

boundaryboundary boundarymirror

3-4h

boundaryboundary

16-18h Sexton et al.

Kc167 Li et al.

boundarymirror

boundaryboundary boundarymirror

boundaryboundary +50 kb boundarymirror

0 -50 kb

-50 kb 0 +50 kb

Figure S4. Comparison of Boundary-Boundary and Boundary-Mirror Contacts across Stages, Related to Figure 3 Aggregate Hi-C plots of boundary-boundary contact regions, as in Figure 3. The boundary-boundary contacts as well as the matched boundary-mirror control contacts are shown (see Figure 3C for the explanatory schematic and Methods). As comparison, the same analysis for previously published Hi-C data from Kc167 cells and 16-18h embryos is shown. ND: not detected.

Figure S5. Validation of Transcription Inhibition Treatments, Related to Figure 5 (A) Representative fluorescent microscopy images of treated eGFP-PCNA embryos at 3:30 hpf. Embryos were injected with water, alpha-amanitin or triptolide. Beyond nuclear cycle 14 (2:30 hpf) alpha-amanitin and triptolide treated embryos showed strong morphological abnormalities with an uneven distribution of nuclei in the blastoderm. At the time of fixation (2:30 hpf), embryos were morphologically typical (see H). We observed no cases where transcription inhibited embryos proceeded to gastrulation, suggesting a 100% penetrant phenotype. White bars represent 250 mm. (B) RT-qPCR analysis of exclusively zygotically expressed transcripts. Shown are expression levels (delta-Ct) calculated using the maternally provided rp49 transcript as reference and then normalized to the mean expression level of uninjected samples. For each of the two biological replicates, embryos were injected prior to nuclear cycle 8. The RNA extraction was performed from two independent batches of five embryos. RNA was collected at nuclear cycle 14. Means are indicated using horizontal bars. (C) qPCR Quantification of ChIP experiments using pSer5 RNA-Pol II antibody on uninjected embryos (control) or inhibitor injected embryos (alpha-amanitin and triptolide). The RNA Pol II occupancy was separately quantified at the promoters and on gene bodies of three zygotically transcribed genes. Values are represented as percent ChIP signal normalized to corresponding sample-specific ChIP inputs. (D) Fold-change Hi-C maps of transcription inhibited embryos at nuclear cycle 14 compared to untreated embryos at earlier stages and between the treatments. Depicted is the same region as in Figure 5B. (E) Hi-C contact maps of control and treated Kc167 cells with the transcription inhibitors triptolide or flavopiridol from (Li et al., 2015). (F) Metaregion plots showing the coverage of expressed genes at nuclear cycle 14 by RNA Pol II ChIP in uninjected embryos compared to embryos injected with either transcription inhibitor using pSer5- and pan-RNA Pol II antibodies. Data are plotted as ChIP-seq fold-enrichment over input on the normalized gene length between transcription start (TSS) and transcription end site (TES). Additionally, a region of 500bp is plotted before and after each gene. (G) Data as in (F) plotted as heatmaps with one gene in every row. Heatmaps are sorted by the row mean fold-enrichment.

(legend continued on next page)

(H) Confocal images of embryos injected with water (control) or with triptolide showing staining of pSer5 RNA Pol II, histone H3 and DAPI. Transcription inhibited embryos are largely normal in their morphology at the time of fixation for Hi-C at nuclear cycle 14. The distribution of histone H3 is slightly altered in treated samples, which might reflect the transcription interdependence of H3.3 deposition (Jullien et al., 2012) or the lack of DNA replication after the treatment (Blythe and Wieschaus, 2015). Bars represent 20 mm.

Fraction

Figure S6. Enrichment of Architectural Proteins, RNA Polymerase II, and Zelda at TAD Boundaries, Related to Figure 6 Enrichment of architectural protein binding (Li et al., 2015; Van Bortle et al., 2014; Wood et al., 2011), RNA Pol II (Blythe and Wieschaus, 2015) and Zelda (Sun et al., 2015) binding averaged across all boundaries identified at nuclear cycle 12, 13 and 14. With the exception of the RNA Pol II data, which is nuclear cycle

(legend continued on next page)

specific, and the Zelda data, which was obtained from 2-3h old embryos, all ChIP-seq data were produced from Kc167 cells and are not stage specific (the complete list of datasets is included in Table S2). The genome was binned in 500 bp intervals for a 80kb region centered at the boundary calls. For each 500 bp bin we considered whether a MACS peak call overlapped and plotted the proportion of boundaries that had at least one overlapping peak within each bin. The gray shading and dashed line correspond to the 95% confidence interval and mean of ChIP signal at the randomized boundaries, respectively. NS: not significant; ** adjusted p value < 1e-3; *** adjusted p value < 1e-4; permutation test; n = 10,000.

Figure S7. Validation of Zelda Depletion and Its Effect on TAD Boundaries, Related to Figure 7 (A) Hi-C contact maps of wild-type and Zelda depleted embryos at nuclear cycle 14. (B) Heatmaps display insulation scores in wild-type and Zelda depleted embryos and the difference in insulation between the two samples. Shown are the same 200 kb regions centered on boundaries in wild-type nuclear cycle 14 embryos. Positive values (purple) indicate a loss of insulation in Zelda depleted embryos; negative values (orange) indicate a gain of insulation. (C) Metaregion plots showing the distribution of Zelda binding centered on boundaries established at different stages. The cycle-specifc Zelda ChIP-data are plotted as ChIP fold-enrichment over input (Harrison et al., 2011). Zelda is largely concentrated on boundaries established at nuclear cycles 12 and 13. (D) RT-qPCR validation of the Zelda depletion phenotype. Shown are expression levels (delta-Ct) calculated using the maternally provided rp49 transcript as reference normalized to the mean expression level of the wild-type samples. Dots represent single replicates (n = 3); dashes represent the mean of all replicates. Term, bnb, Bsg25A, CG10035, and CG14317 are known to depend on Zelda for their activation (Blythe and Wieschaus, 2015). Mdr49 is used as control since it is a strictly zygotic gene that has been shown to be transcribed in the absence of Zelda (Blythe and Wieschaus, 2015). (E) Metaregion plots showing the insulation score in wild-type and Zelda depleted nuclear cycle 14 embryos centered on Zelda peaks at the same stage as in Figure 7A. The insulation signal is stratified by the strength of Zelda peaks. For this purpose, the genome was divided into 5kb bins and Zelda binding was aggregated for each bin. The bins were then divided evenly into five quintiles from strongest to weakest Zelda binding. (F) Barplot representing the number of boundaries that lose insulation following Zelda depletion. **p < 1e-2; ***p < 1e-3; NS, not significant; Fisher’s exact test. (G) Aggregate Hi-C plots of boundary-boundary contact regions, as in Figure 3D for wild-type and Zelda depleted nuclear cycle 14 embryos.