1. Genes and Chromosomes
  2. Genomics and Evolutionary Biology
Download icon

Convergence of topological domain boundaries, insulators, and polytene interbands revealed by high-resolution mapping of chromatin contacts in the early Drosophila melanogaster embryo

  1. Michael R Stadler
  2. Jenna E Haines
  3. Michael B Eisen Is a corresponding author
  1. University of California, United States
  2. Howard Hughes Medical Institute, United States
Research Article
  • Cited 0
  • Views 2,189
  • Annotations
Cite as: eLife 2017;6:e29550 doi: 10.7554/eLife.29550

Abstract

High-throughput assays of three-dimensional interactions of chromosomes have shed considerable light on the structure of animal chromatin. Despite this progress, the precise physical nature of observed structures and the forces that govern their establishment remain poorly understood. Here we present high resolution Hi-C data from early Drosophila embryos. We demonstrate that boundaries between topological domains of various sizes map to DNA elements that resemble classical insulator elements: short genomic regions sensitive to DNase digestion that are strongly bound by known insulator proteins and are frequently located between divergent promoters. Further, we show a striking correspondence between these elements and the locations of mapped polytene interband regions. We believe it is likely this relationship between insulators, topological boundaries, and polytene interbands extends across the genome, and we therefore propose a model in which decompaction of boundary-insulator-interband regions drives the organization of interphase chromosomes by creating stable physical separation between adjacent domains.

https://doi.org/10.7554/eLife.29550.001

eLife digest

The DNA inside a cell is packaged into threaded structures called chromosomes. Early studies of chromosomes using insect larvae revealed a pattern of dark- and light-colored bands on the chromosomes that was unique to every region. For decades, it remained unclear if the bands had a specific role.

More advanced techniques have shown that chromosomes are organized into a series of compact domains that contain separate regions of genes. Each region can be turned on or off at different times, depending on the needs of different cells. This allows cells to specialize into different cell types and tissues. Until now, it was unclear how these different regions are formed and controlled, and how they relate to the chromosome bands.

Here, Stadler, Haines and Eisen used a specific technique to map the structure of chromosomes in early fly embryos, by using a chemical trap to capture closely located DNA pieces. The results showed that the chromosome domains corresponded to the banding patterns seen in the early studies. This suggest that light bands represent extended DNA that act as spacers between the dark gene regions.

This study adds to the view that the way the DNA is organized influences gene activity. Creating a high-resolution model of the chromosomes will help us to better understand how their structure can influence the activity of genes. In future, scientists may be able to identify diseases that are caused by errors in the chromosome organization.

https://doi.org/10.7554/eLife.29550.002

Introduction

Beginning in the late 19th century, cytological investigations of the polytene chromosomes of insect salivary glands implicated the physical structure of interphase chromosomes in their cellular functions (Balbiani, 1881; Balbiani, 1890; Heitz and Bauer, 1933; King and Beams, 1934; Painter, 1935). Over the next century plus, studies in the model insect species Drosophila melanogaster were instrumental in defining structural features of animal chromatin. Optical and electron microscopic analysis of fly chromosomes produced groundbreaking insights into the physical nature of genes, transcription and DNA replication (Benyajati and Worcel, 1976; Laird and Chooi, 1976a; Laird et al., 1976b; McKnight and Miller, 1976; McKnight and Miller, 1977; Vlassova et al., 1985; Belyaeva and Zhimulev, 1994).

Detailed examination of polytene chromosomes in Drosophila melanogaster revealed a stereotyped organization, with compacted, DNA-rich ‘bands’ alternating with extended, DNA-poor ‘interband’ regions (Bridges, 1934; Rabinowitz, 1941; Lefevre, 1976; Benyajati and Worcel, 1976; Laird and Chooi, 1976a), and it appears likely that this structure reflects general features of chromatin organization shared by non-polytene chromosomes. While these classical studies offered extensive structural and molecular characterization of chromosomes in vivo, the question of what was responsible for organizing chromosome structure remained unanswered.

A critical clue came with the discovery of insulators, DNA elements initially identified based on their ability to block the activity of transcriptional enhancers when located between an enhancer and its targeted promoters (Kellum and Schedl, 1991; Holdridge and Dorsett, 1991; Geyer and Corces, 1992; Kellum and Schedl, 1992). Subsequent work showed that these elements could also block the spread of silenced chromatin states (Roseman et al., 1993; Sigrist and Pirrotta, 1997; Mallin et al., 1998; Recillas-Targa et al., 2002; Kahn et al., 2006) and influence the structure of chromatin. Through a combination of genetic screens and biochemical purification, a number of protein factors have been identified that bind to Drosophila insulators and modulate their function, including Su(Hw), BEAF-32, mod(mdg4), CP190, dCTCF, GAF, Zw5, and others (Lindsley and Grell, 1968; Lewis, 1981; Parkhurst and Corces, 1985; Parkhurst and Corces, 1986; Spana et al., 1988; Parkhurst et al., 1988; Zhao et al., 1995; Gerasimova et al., 1995; Bell et al., 1999; Gaszner et al., 1999; Scott et al., 1999; Büchner et al., 2000; Pai et al., 2004; Melnikova et al., 2004; Moon et al., 2005). Except for CTCF, which is found throughout bilateria, all of these proteins appear to be specific to arthropods (Heger et al., 2013).

Staining of polytene chromosomes with antibodies against such insulator proteins showed that many of them localize to polytene interbands (Belyaeva and Zhimulev, 1994; Zhao et al., 1995; Byrd and Corces, 2003; Eggert et al., 2004; Pai et al., 2004; Gortchakov et al., 2005; Gilbert et al., 2006; Berkaeva et al., 2009; Vatolina et al., 2011b), with some enriched at interband borders. Further, some, though not all, insulator protein mutants disrupt polytene chromosome structure (Roy et al., 2007). Together, these data implicate insulator proteins, and the elements they bind, in the organization of the three-dimensional structure of fly chromosomes.

Several high-throughput methods to probe three-dimensional structure of chromatin have been developed in the last decade (Lieberman-Aiden et al., 2009; Fullwood et al., 2009; Rao et al., 2014; Beagrie et al., 2017). Principle among these are derivatives of the chromosome conformation capture (3C) assay (Dekker et al., 2002), including the genome-wide ‘Hi-C’ (Lieberman-Aiden et al., 2009). Several groups have performed Hi-C on Drosophila tissues or cells and have shown that fly chromosomes, like those of other species, are organized into topologically associated domains (TADs), regions within which loci show enriched 3C linkages with each other but depleted linkages with loci outside the domain. Disruption of TAD structures by gene editing in mammalian cells has been shown to disrupt enhancer-promoter interactions and significantly alter transcriptional activity (Guo et al., 2015; Lupiáñez et al., 2015).

Although TADs appear to be a common feature of animal genomes, the extent to which TAD structures are a general property of a genome or if they are regulated as a means to control genome function remains unclear, and the question of how TAD structures are established remains largely open. Previous studies have implicated a number of features in the formation of Drosophila TAD boundaries, including transcriptional activity and gene density, and have reached differing conclusions about the role played by insulator protein binding (Sexton et al., 2012; Hou et al., 2012; Van Bortle et al., 2014; Ulianov et al., 2016; Li et al., 2015). Tantalizingly, Eagen et al., using 15 kb resolution Hi-C data from D. melanogaster have shown that there is a correspondence between the distribution of large TADs and polytene bands (Eagen et al., 2015).

We have been studying the formation of chromatin structure in the early D. melanogaster embryo because of its potential impact on the establishment of patterned transcription during the initial stages of development. We have previously has shown that regions of ‘open’ chromatin are substantially remodeled at enhancers and promoters during early development (Harrison et al., 2011; Li et al., 2014) and were interested in the role three-dimensional chromatin structure plays in spatial patterning.

We therefore generated high-resolution Hi-C datasets derived from nuclear cycle 14 Drosophila melanogaster embryos (Foe and Alberts, 1983), and from the anterior and posterior halves of hand-dissected embryos at the same developmental stage. We show that high-resolution chromatin maps of anterior and posterior halves are nearly identical, suggesting that chromatin structure neither drives nor directly reflects spatially patterned transcriptional activity. However, we show that stable long-range contacts evident in our chromatin maps generally involve known patterning genes, implicating chromatin conformation in transcriptional regulation.

To investigate the origins of three-dimensional chromatin structure, we carefully map the locations of the boundaries between topological domains using a combination of manual and computational annotation. We demonstrate that these boundaries resemble classical insulators: short (500–2000 bp) genomic regions that are strongly bound by (usually multiple) insulator proteins and are sensitive to DNase digestion. Additionally, we find that boundaries share the molecular features of polytene interband regions. Finally, we show that for a region in which the fine polytene banding pattern has been mapped to genomic positions, boundaries show precise colocalization with interband regions that separate compacted bands corresponding to TADs. We propose that this relationship between insulators, TADs and polytene interbands extends across the genome, and suggest a model in which the decompaction of these regions drives the organization of interphase fly chromosomes by creating stable physical separation between adjacent domains.

Results

Data quality and general features

We prepared and sequenced in situ Hi-C libraries from two biological replicates of hand-sorted cellular blastoderm (mitotic cycle 14; mid-stage 5) embryos using a modestly adapted version of the protocol described in Rao et al., 2014. To examine possible links between chromatin maps and transcription, we sectioned hand-sorted mitotic cycle 14 embryos along the anteroposterior midline, and generated Hi-C data from the anterior and posterior halves separately, also in duplicate. In total, we produced ~452 million informative read pairs (see Supplementary file 1).

We assessed the quality of these data using metrics similar to those described by (Lieberman-Aiden et al., 2009; Rao et al., 2014). Specifically, the strand orientations of our reads were approximately equal in each sample (as expected from correct Hi-C libraries but not background genomic sequence; see Supplementary file 1), the signal decay with genomic distance was similar across samples, and, critically, visual inspection of heat maps prepared at a variety of resolutions showed these samples to be very similar both to each other and to previously published data prepared using similar methods (Sexton et al., 2012). We conclude that these Hi-C are of high quality and reproducibility.

We next sought to ascertain the general features of the data at low resolution. We examined heatmaps for all D. melanogaster chromosomes together using 100 kb bins, as shown in Figure 1. Several features of the data are immediately apparent. The prominent ‘X’ patterns for chromosomes 2 and 3, which indicate an enrichment of linkages between chromosome arms, reflects the known organization of fly chromosomes during early development known as the Rabl configuration (Csink and Henikoff, 1998; Wilkie et al., 1999; Duan et al., 2010): telomeres are located on one side of the nucleus, centromeres are located on the opposite side, and chromosome arms are arranged roughly linearly between them. Centromeres and the predominantly heterochromatic chromosome 4 cluster together, as, to a lesser extent, do telomeres, reflecting established cytological features that have been detected by prior Hi-C analysis (Sexton et al., 2012) and fluorescence in situ hybridization (FISH) (Lowenstein et al., 2004). These features were evident in all replicates, further confirming both that these datasets are reproducible and that they capture known features of chromatin topology and nuclear arrangement.

Hi-C map of the stage 5 Drosophila melanogaster genome at 100 kb resolution.

Data from all nc14 datasets was aggregated and normalized by the ‘vanilla coverage’ method. To enhance contrast, the logarithm values of the normalized counts were histogram equalized, and maximum and minimum values were adjusted for optimal display.

https://doi.org/10.7554/eLife.29550.003

TAD boundaries are short elements bound by insulator proteins

Because we used a 4-cutter restriction enzyme and deep sequencing, and because the fly genome is comparatively small, we were able to resolve features at high resolution. We visually inspected genome-wide maps of a number of genomic regions constructed using bins of 500 bp, and were able to see a conspicuous pattern of TADs across a wide range of sizes, some smaller than 5 kb (Figure 2, Figure 2—figure supplements 15). When we compared maps for several of these regions with available functional genomic data from embryos, we observed that the boundaries between these domains showed a remarkably consistent pattern: they were formed by short regions of DNA (500–2000 bp) that are nearly always associated with high chromatin accessibility, measured by DNase-seq (Li et al., 2011), strong occupancy by known insulator proteins as measured by chromatin immunoprecipitation (ChIP) (Nègre et al., 2010) (Figure 2, Figure 2—figure supplements 15) properties characteristic of classical Drosophila insulator elements.

Figure 2 with 5 supplements see all
Example region of Hi-C data at 500 bp resolution.

Heat map of aggregate Hi-C data for all nc14 datasets binned at 500 bp is shown for the region located at 3R:24924500–25174500 (dm3: 3R:20750000–20999999). Raw counts were normalized by the vanilla coverage method, the logarithm was taken, and minimum and maximum values were adjusted for visual contrast. A UCSC browser (Kent et al., 2002) window for the corresponding coordinates is shown with tracks for Hi-C directionality (calculated from the Hi-C data shown in the heatmap), DNase accessibility (X.-Y. Li et al., 2011), RNA polII and TFIIB (Li et al., 2008), and the insulator proteins CP190, BEAF-32, dCTCF, GAF, mod(mdg4), and Su(Hw) from (Nègre et al., 2010). Dashed red lines are visual guides and are manually drawn at locations of apparent boundaries; they do not reflect algorithmically or unbiased hand-curated boundary calls.

https://doi.org/10.7554/eLife.29550.004

To confirm this visually striking association, we systematically called TAD boundaries by visual inspection of panels of raw Hi-C data covering the entire genome. Critically, these boundary calls were made from Hi-C data alone, and the human caller lacked any information about the regions being examined, including which region (or chromosome) was represented by a given panel. In total, we manually called 3122 boundaries in the genome for nc14 embryos. Taking into account the ambiguity associated with intrinsically noisy data, the difficulty of resolving small domains, and the invisibility of sections of the genome due to repeat content or a lack of MboI cut sites, we consider 4000–5500 to be a reasonable estimate for the number of boundaries in the genome.

To complement these manual calls, we developed a computational approach for calling boundaries that is similar to methods used by other groups (Lieberman-Aiden et al., 2009; Sexton et al., 2012; Rao et al., 2014; Crane et al., 2015). In brief, we assigned a directionality score to each genomic bin based on the number of Hi-C reads linking the bin to upstream versus downstream regions, and then used a set of heuristics to identify points of transition between regions of upstream and downstream bias. We adjusted the parameters of the directionality score and the boundary calling to account for features of the fly genome, specifically the relatively small size of many topological domains.

Attempts to exhaustively and definitively identify features within genomic data are necessarily variable due to differences in the choice of algorithm, parameters, cutoffs, and unavoidable tradeoffs between sensitivity and accuracy. We therefore sought a representative set of TAD boundaries with which to analyze features of these elements. Of our top 1000 computationally-identified domain boundaries, we found that 952 were matched by a manually-called boundary within 1 kb. This high level of agreement suggested that the computational approach robustly identified the domain features that are apparent by eye. By taking the union of our computational calls, applied with a stringent cutoff, and our manual calls, we developed a very conservative set of exceptionally high confidence boundaries. We emphasize that this set represents only a subset of the boundaries identified by manual and computational approaches, the complete lists of which are provided in Supplementary file 1.

Comparing these 952 boundaries to other genomic datasets confirms our initial observations and reveals a highly stereotyped pattern of associated genomic features. Most strikingly, boundaries are enriched for the binding of the known insulator proteins CP190, BEAF-32, mod(mdg4), dCTCF, and to a lesser extent GAF and Su(Hw) (Figure 3). CP190 and BEAF-32 show the strongest enrichment, and indeed, virtually all (95.1%) of the examined boundaries appear to be associated with CP190 binding (Figure 3—figure supplement 1). Domains of H3K27 trimethylation, a marker of polycomb silencing, showed a strong tendency to terminate at boundaries, and the enhancer mark H3K4me1 showed an interesting pattern of depletion at boundaries but enrichment immediately adjacent to boundary locations (Figure 3). Boundaries also exhibit peaks of DNase accessibility and nucleosome depletion (Figure 3), as well as marks associated with promoters, including the general transcription factors TFIIB and the histone tail modification H3K4me3. Despite the presence of promoter marks, we find that RNA polII is present at only a subset (45.1%) of stage 5 boundaries (Figure 3, Figure 3—figure supplement 1).

Figure 3 with 4 supplements see all
Topological domain boundaries show distinct patterns of associated proteins and genomic features.

Heatmaps showing the distribution of signals from embryonic ChIP and DNase-seq datasets around 952 topological boundaries identified jointly by computational and manual curation. All plots show 500 bp genomic bins in 100 kb windows around boundaries. All plots in blue are sorted by boundary strength, calculated from the difference in upstream and downstream Hi-C directionality scores. The plot for H3K27me3 (in red) is specially sorted to highlight the tendency for enriched domains to terminate at boundaries. Rows for this plot were sorted by calculating the total H3K27me3 signal in the 50 kb windows upstream and downstream of the boundary and then sorting, top to bottom: upstream signal above median and downstream signal below the median, upstream below and downstream above, upstream and downstream both above, upstream and downstream both below the median. For comparison, identically prepared and sorted plots around H3K4me3 peaks are shown in Figure 3—figure supplement 2. Percentages are calculated as the percentage of boundaries with a >2 fold enrichment for the given signal within a 3 kb window centered on the boundary (±1.5 kb). Data for insulator proteins, DNase accessibility, RNA polII and TFIIB are from the same sources indicated in Figure 2. ChIP for H3, H3K4me1 are taken from (Li et al., 2014), and H3K27me3 are from modEncode (Contrino et al., 2012).

https://doi.org/10.7554/eLife.29550.010

It is striking that we observe that not only are sites of combinatorial insulator protein binding enriched at TAD boundaries, but they are highly predictive. Of our representative set of boundaries, 95.1% are are enriched >2 fold for CP190 binding within a 1.5 kb window. Conversely, of the strongest 1000 CP190 peaks, 75.2% are within 2 kb of a manual or computationally-called boundary (compared to 37.4% of the top 1000 RNAPII peaks). It is important to note that we do identify a small subset of boundaries that are not apparently associated with sites of insulator binding (~1–2% show no enrichment for CP190, BEAF-32, or dCTCF, depending on thresholds used), suggesting that there are multiple phenomena that can create topological boundaries in flies (e.g., see Figure 6). However, the overwhelming majority of topological boundaries identified in this study coincide with elements that match the properties of CP190-associated insulators.

An important confounding factor in sorting out the nature of topological boundaries is the strong tendency, observed by multiple authors, of insulator proteins to bind near promoters specifically between divergent promoters (Nègre et al., 2010; Ramirez et al., 2017; Schwartz et al., 2012). Indeed, we find that boundary elements, as identified from Hi-C, are often found proximal to promoters and show a general enrichment of promoter-associated marks (Figure 3), raising the possibility that transcriptional activity at promoters may drive topological boundary formation. However, a number of features of the data argue against this possibility. First and most critically, many of the topological boundaries (54.9%) we identify are not associated with RNAPII binding in nc14 embryos. Similarly, there are many active promoters that do not appear to form topological boundaries (e.g., Figure 8 and supplements). Hug et al., 2017 pharmacologically inhibited transcription in early embryos and observe that TADs remain intact. Finally, topological boundaries are invariant between anterior and posterior sections of embryos despite substantial differences in the transcriptional profiles of these regions (see below). We further examined the distributions of the same genomic features around the top 1000 peaks of H3K4me3, a marker of active promoters, in data from stage 5 embryos (Figure 3—figure supplement 2) (Li et al., 2014). While these sites show enrichments for insulator proteins, these enrichments are substantially weaker than those observed at topological boundaries, while RNA polII enrichment is much stronger at promoters than boundaries. The tendency for polycomb domains to terminate at promoters is also much less pronounced at promoters than boundaries. Together, these data argue that boundaries constitute a distinct class of genetic elements that are not formed by promoter transcription, but are instead frequently located near promoters, possibly as a result of selective pressure to insulate these proximal promoters from distal regulatory elements. While we cannot rule out any role for promoter-bound transcription machinery in the formation of topological boundaries (notably, TFIIB is enriched at 69.1% of boundaries), we think it is unlikely that transcriptional activity plays a major role in establishing the topological domains of interphase fly chromosomes.

Finally, we examined the sequence composition of boundary elements by comparing the frequency of DNA words of up to seven base pairs in the set of high confidence boundaries to flanking sequence. The most enriched sequences correspond to the known binding site of BEAF-32 and to a CACA-rich motif previously identified as enriched in regions bound by CP190 (Nègre et al., 2010; Yang and Corces, 2012), both of which show strong association with the set of boundary sequences as a whole (Figure 4).

Sequence features of TAD boundary elements.

(A) Histograms showing the frequency of enriched 7-mers in 5 kb windows around 952 high-confidence TAD boundaries. (B) Scatter plots of occurrences of words matching known BEAF-32 binding motifs (left) and CACA motif (right) in 10 kb windows around high-confidence TAD boundaries. Points are plotted with low opacity, such that darker points correspond to positions where multiple words occurs close together in sequence.

https://doi.org/10.7554/eLife.29550.015

Boundary elements correspond to polytene interbands

The examination of these boundary elements led us to consider the physical basis of topological domain separation. Chromosome conformation capture is a complex assay (Gavrilov et al., 2013; Gavrilov et al., 2015), and inferring discrete physical states of the chromatin fiber from Hi-C signals generally requires orthogonal experimental data. To address this problem, we sought to leverage information from polytene chromosomes to draw associations between features of Hi-C data and physical features of chromosomes.

The Zhimulev laboratory has extensively studied the nature and composition of polytene bands and interbands for decades. Using a combination of approaches, they have identified interbands as a set of ~5700 short decompacted regions that tend to be located near divergent promoters and are characterized by DNase hypersensitivity and the binding of characteristic proteins, including insulator proteins (Zhimulev et al., 2014). It was immediately apparent to us that these elements bore significant similarity to the topological boundary elements we identified. We thus sought to compare our Hi-C data to known polytene chromosome structures.

There is surprisingly little data mapping features of polytene structure to specific genomic coordinates at high resolution. Vatolina et al., 2011a used exquisitely careful electron microscopy to identify the fine banding pattern of the 65 kb region between polytene bands 10A1-2 and 10B1-2, revealing that this region, which appears as a single interband under a light microscope, actually contains six discrete, faint bands and seven interbands. The region is flanked by two large bands, whose genomic locations have been previously mapped and refined by FISH (Vatolina et al., 2011a). Vatolina et al. then used available molecular genomic data to propose a fine mapping of these bands and interbands to genomic coordinates.

Figure 5 shows the correspondence between Vatolina et al.’s proposed polytene map from this region and our high-resolution Hi-C data, along with measures of early embryonic DNase hypersensitivity from (Li et al., 2011) and the binding of six insulator proteins (Nègre et al., 2010). There is a striking correspondence between the assignments of Vatolina et al. and our Hi-C data: faint polytene bands correspond to TADs, and interbands correspond to the boundary elements that separate the TADs.

Figure 5 with 2 supplements see all
Topological boundary elements correspond to polytene interbands.

Heat map of aggregate Hi-C data for all nc14 datasets binned at 500 bp and UCSC browser data shown for the region X:11077500–11181000 (dm3: X:10971500–11075000) for which Vatolina et al. provided fine-mapping of polytene banding structure. Hi-C and browser data were prepared and sourced as indicated in Figure 2. Dashed red lines are visual guides drawn from the interband assignments of Vatolina et al. Top: accurately-scaled representations of the size of the mapped bands and interbands in base pairs (‘Genomic’) and the corresponding physical distances in polytene chromosomes derived from electron microscopic analysis of polytene chromosomes by Vatolina et al. Increased relative physical size of interband regions demonstrates their lower compaction ratios.

https://doi.org/10.7554/eLife.29550.016

This correspondence is not perfect. Specifically, the evidence in our Hi-C data for the separation between the major band 10A1-2 and the minor band 10A3 is weak, though that may be partly explained by the absence of MboI cut sites obscuring much of this region. This minor band is barely detectable in polytene spreads, and the combination of this and weak support in Hi-C data may suggest that this band is not real or is perhaps only present in a minority of nuclei. Similarly, the Vatolina et al. report that they only rarely observe the interband between bands 10A6 and 10A7, and we indeed observe substantial contact between these two putative bands in Hi-C maps (the light orange region near the peak of the ‘pyramid’ formed by 10A6 and 10A7 in Figure 5), though each shows stronger intra- than inter-domain interactions. One possible explanation for this observation is that the interband separating these two domains is not constitutive but rather is formed in only a fraction of nuclei. The pattern exhibited by these two domains--adjacent domains that show a clear separation but also a substantial interaction signal--is one we observe frequently in our early embryonic Hi-C data, suggesting that variable boundaries may be common features of the fly genome.

Overall, the alignment between polytene band mapping and Hi-C data in this region supports a strong correspondence between these two types of data. For five interbands which were easily visible in polytene spreads (10A3/4-5, 10A4-5/10A6, 10A7/10A8-9, 10A8-9/10A10, 10A10-11/10B1), we observe strong domain boundaries in Hi-C data. For two interbands supported by weaker evidence in polytene analysis, we observe in Hi-C maps a weak or non-existent boundary (10A1-2/10A3) and a boundary with significant interaction across it, possibly representing heterogeneity between nuclei matching heterogeneity in polytenes (10A6/10A7).

The 5’ region of the Notch gene has also been carefully mapped. Rykowski et al. used high-resolution in situ hybridization to determine that the coding sequences of Notch lies within polytene band 3C7, while the sequences upstream of the transcription start site (TSS) lie in the 3C6-7 interband. Examining the Notch locus in our Hi-C data, we see that the gene body is located within an ~20 kb TAD, and the TSS directly abuts a TAD boundary that is strongly bound by CP190 and dCTCF (Figure 5—figure supplement 1), an arrangement consistent with the correspondence of boundaries and interbands.

The chromodomain-containing protein Chriz has been suggested as the strongest diagnostic feature of polytene interbands (Zhimulev et al., 2014). Using publicly available ChIP datasets from Kc167 cells (derived from late embryonic tissue), we observed a strong enrichment of Chriz binding at our representative boundaries (87.9% >2 fold enriched within 1.5 kb, Figure 5—figure supplement 2A). Further, Hi-C directionality around Chriz peaks shows the characteristic pattern of boundary formation, and Chriz profiles at representative loci show substantial correspondence between boundary regions and Chriz binding (Figure 5—figure supplement 2B and C), offering further support for the association between boundaries and interbands.

Eagen et al. previously identified a broad correspondence between polytene interbands and inter-TAD regions from Hi-C data at 15 kb resolution (Eagen et al., 2015). Our Hi-C data allows the detection of fine structure within these inter-TAD regions, down to individual boundary elements. Owing to the dearth of finely mapped polytene regions, the association between topological boundaries and interband regions is necessarily based on a limited number of example loci. However, the combination of data from these loci with the close agreement of the molecular composition of these regions, specifically the strong localization of the interband marker Chriz to topological boundaries, leads us to propose that the precise relationship between topological boundaries, insulator elements, and decompacted interband regions we observe is a general one, and that it extends across the genome.

The association between boundary elements and interbands suggests a simple model for insulator function. A key feature that distinguishes polytene interbands from bands is their low compaction ratio: they span a larger physical distance per base pair. The association between insulator binding and genomic regions with low compaction ratios suggests insulators may function by simply increasing the physical distance between adjacent domains via the unpacking and extension of intervening chromatin. Figure 5 (top) shows a representation of the conversion of genomic distance to physical distance for the 10A1-B1 region, as measured by Vatolina et al. Any model for insulator function must explain several features of insulator function, including the ability to organize chromatin into physical domains, block interactions between enhancers and promoters exclusively when inserted between them, protect transgenes from position effect variegation and block the spread of chromatin silencing states. This chromatin extension model for fly insulator function can potentially explain these defining characteristics via simple physical separation.

Hi-C data can elaborate fine polytene structures

We reasoned that if our Hi-C data is capable of resolving fine banding patterns such as that at the 10A1-B1 locus, we should be able to resolve the borders of major bands with precision. We focused on a region of chromosome 2L that had previously been shown by Eagen et al. to appear as a single ~500 kb TAD using Hi-C at 15 kb resolution, but contains a faint interband in Bridge’s map. Our Hi-C data reveal an intricate structure at this locus (Figure 6A). There are two large TADs on the left and right, divided by a series of smaller domains in the center. We suspect that this middle region accounts for the interband in Bridge’s map, in a manner similar to the 10A1-1/10B1-2 region: a complex region consisting of several minor bands bounded by decompacted boundary regions appears as a single interband region under optical microscopy.

Complex topological structure of a region of chromosome 2L.

Hi-C maps using 500 bp bins of the region of chromosome 2L corresponding to polytene band 22A1-2. This regions was shown by Eagen et al. to comprise a single TAD in Hi-C data viewed at 15 kb resolution, and is occasionally observed to contain an interband in polytene spreads. (A) View of the entire region, revealing complex internal structure. (B–D) Zoomed-in views of three regions comprising the left border (B), complex middle section (C), and right border of the larger region corresponding to the band/TAD investigated by Eagen et al., with associated stage 5 DNase accessibility, CP190, and Chriz (kc167 cells) profiles. Coordinates for this region are identical in dm3 and dm6.

https://doi.org/10.7554/eLife.29550.019

This region provides examples of a number of interesting features that we observe in our Hi-C data. First, the large TADs are bounded on both sides by gene-rich regions consisting of a number of smaller topological domains (Figure 6B,D). The boundaries of large and small domains in this region nearly all share the common features of boundary elements: DNase hypersensitivity and binding of diagnostic insulator (e.g. CP190) and interband (CHRIZ) proteins. This region also contains a prominent example of an exception to this pattern: a loop is formed that appears to generate boundaries not associated with these characteristic protein binding events (Figure 6C, indicated by dotted yellow lines and loop). This example highlights a critical point: while the description we provide of the association between TAD boundaries, insulator elements, and decompacted interbands appears to describe the overwhelming majority of cases, there are counter-examples. Indeed given the extraordinary capacity of nature to innovate with respect to gene regulation and structures, we expect that animal genomes will provide no shortage novel chromosome topological and structural features for future investigations.

Topological boundaries are nearly identical in anterior and posterior sections of the embryo

We next asked whether the boundaries we identified as boundary elements represent constitutive features of chromatin organization or whether their function might be regulated in a cell-type specific or developmental manner. We reasoned that, since different sets of patterning genes are transcribed in the anterior and posterior portions of the pre-gastrula D. melanogaster embryos, a comparison of chromatin interaction maps between anterior and posterior regions would reveal whether context, especially transcriptional state, affects the TAD/boundary structure of the genome. To this end, we performed two separate biological replicates of an experiment in which we sectioned several hundred mid stage 5 embryos along the anteroposterior midline, and produced deep-sequenced Hi-C libraries from the anterior and posterior halves in parallel.

Resulting Hi-C signals at boundaries are virtually identical in the two halves, despite substantially different gene expression profiles in these two embryonic regions (Figure 7A). Indeed, overall Hi-C signals are remarkably similar, with anterior and posterior samples correlating as strongly as replicates. Examination of individual loci at high resolution reveal consistent profiles and boundaries, notably including genes expressed differentially in the anterior or posterior (Figure 7B).

Figure 7 with 1 supplement see all
Hi-C signals from anterior and posterior halves of stage 5 embryos reveal highly similar chromatin topologies.

(A) The distribution of Hi-C directionality scores in whole embryos, anterior, and posterior halves is shown around 952 topological boundaries identified jointly by computational and manual curation. (B) Heat maps of Hi-C data at 500 bp resolution at four example regions in anterior and posterior embryo halves. Plots represent the aggregate data of two biological and technical replicates each for anterior and posterior samples, and were prepared as in Figure 2. The regions shown are the region mapped by Vatolina et al. (dm6: X:11077500–11181000, dm3: X: 10971500–11075000), the example region from Figure 2—figure supplement 4 (dm6: 2R:9124500–9223500, dm3: 2R:5012000–5111000)), and the genomic regions surrounding the eve (dm6: 2R:9903060–10056959, dm3: 2R:5790565–5944464) and ftz (dm6: 3R:6769234–6961333, dm3: 3R:2594956–2787055) loci. (C) Chromatin accessibility around topological boundaries as measured by ATAC-seq in anterior and posterior nc14 (S5) embryos and by DNase-seq on stage 11 and 14 embryos (X.-Y. Li et al., 2011).

https://doi.org/10.7554/eLife.29550.020

The correspondence of insulator boundary elements and interbands, and the chromatin extension model, implies that the chromatin accessibility of insulator regions will be a useful proxy for their functionality in structurally organizing the genome. Intriguingly, (Van Bortle et al., 2014) found that DNase accessibility of insulator protein-bound regions tracked with the ability of these sequences to block enhancer-promoter interactions in a cell-culture assay. We again sectioned embryos into anterior and posterior halves and performed ATAC-seq (Buenrostro et al., 2013) on pools of 20 embryo halves. ATAC-seq is a technique in which intact chromatin is treated with Tn5 transposase loaded with designed DNA sequences which are preferentially inserted into open, accessible chromatin regions. These insertions can be used to generate high-throughput sequencing libraries, producing data that is largely analogous to DNase-seq data.

Analysis of ATAC-seq signal at insulator boundary elements in anterior and posterior halves showed that these elements have nearly identical accessibility in these two samples (Figure 7C). Additionally, DNase-seq data from later embryonic stages that feature substantial tissue differentiation, transcription, and chromatin changes show highly consistent profiles at boundaries (Figure 7C, Figure 7—figure supplement 1). It is also striking that we observe significant enrichment of insulator proteins and Chriz at boundaries, despite the fact that boundaries were identified from Hi-C data from carefully-staged nc14 embryos (2–3 hr), whereas these ChIP datasets are derived from 0 to 12 hr old embryos or late embryonic cultured cells (Chriz). Together, these results are consistent with a model in which insulator-mediated chromatin organization is a constitutive feature of interphase chromatin of D. melanogaster embryos.

Distal chromatin contacts in the early fly embryo

Many models of insulator function invoke physical contact between insulators to form ‘looped’ chromatin domains (Fujioka et al., 2009; Yang and Corces, 2012; Kyrchanova and Georgiev, 2014; Kravchenko et al., 2005), and a substantial literature exists demonstrating that many insulator proteins are able to interact with each other and to self-associate (Büchner et al., 2000; Gause et al., 2001; Ghosh et al., 2001; Blanton et al., 2003; Pai et al., 2004; Mohan et al., 2007; Golovnin et al., 2007; Vogelmann et al., 2014). In general, we do not observe looping interactions between domain boundaries in our Hi-C data. However, during manual calling of topological boundaries for the entire genome, we noted 46 prominent examples of interactions between non-adjacent domains (Figure 8 and Figure 8—figure supplements 110, Supplementary file 1), in addition to the previously noted clustering of PcG-regulated Hox gene clusters (Sexton et al., 2012). Because the interactions we observed were not of a uniform character, we did not attempt to computationally search for all such phenomena in our data, nor do we claim that this list is necessarily complete. It is merely the union of two sets of ‘interesting’ loci identified in two independent rounds of visual inspection Hi-C maps for the entire genome, and we feel it is informative with respect to the nature and significance of distal interaction in the fly embryo.

Figure 8 with 10 supplements see all
Looping and domain-skipping activity observed in nc14 chromatin.

(A) An example of domain-skipping and looping at the Scr-ftz-Antp locus. ftz is contained within a domain that shows enriched Hi-C interactions between its boundaries, indicative of the formation of a looped domain. Adjacent domains show depleted interaction with the ftz domain and enriched interaction with each other, with especially strong contacts between the region containing the Scr promoter and characterized Scr regulatory elements 3’ of the Antp locus (Calhoun et al., 2002; Calhoun and Levine, 2003). Dotted lines connect features in the Hi-C map to the genomic locations of genes in this region. (B) A strong looping interaction between the kni locus and the 5’ end of the related knrl (kni-like) gene. kni and knrl are known to have identical expression patterns and partially redundant, though distinct, domains of biochemical activity (González-Gaitán et al., 1994).

https://doi.org/10.7554/eLife.29550.022

The most visually striking locus, which we emphasize was identified in an unbiased manner without knowing its identify, is the locus containing the Scr, ftz, and Antp genes (Figure 8A). This locus has been extensively studied, and a number of regulatory elements have been identified that reside between the ftz and Antp genes but ‘skip’ the ftz promoter to regulate Scr (Calhoun et al., 2002; Calhoun and Levine, 2003). Consistent with this, we observe enriched contacts between the region containing the Scr promoter and a domain on the other side of ftz that contains the known Scr-targeting cis regulatory elements, while the ftz-containing domain makes minimal contact with its neighboring domains. Critically, we observe hot spots of apparent interaction between two sets of boundary elements (Figure 8A: 1 and 4, 2 and 3), suggesting that physical association of boundary elements (or their associated proteins) may play a role in this interaction.

Curiously, we detected a similar situation on the other side of Scr, where a domain containing the hox gene Dfd is ‘skipped’ over by the Ama locus to interact with a short element 3’ of the Scr transcription unit (Figure 8—figure supplement 1). We also observe a similar arrangement near the eve locus (Figure 8—figure supplement 2). In these cases, a plausible topology is that the skipped domain is ‘looped out’, preventing interaction with neighbors, while the adjacent domains are brought into proximity.

In addition to these domain-skipping events, we observe a small number of looping interactions, where two distal loci show high levels of interaction, without the associated enriched interactions between the domains flanking the loop. In every case we observe, the loop forms between two domain boundaries. As shown in Figure 8B, one of these loops brings together the promoters of knirps and the related knrl (knirps-like) genes. Other loops connect the achaete and scute genes (Figure 8—figure supplement 3), sloppy paired 1 and sloppy paired 2 (Figure 8—figure supplement 4), and the promoter of Ultrabithorax  with an element in its first intron (Figure 8—figure supplement 5).

These loci demonstrate that looping and domain-skipping events can be detected in our Hi-C data, but it appears that such interactions are rare and that looping does not occur between the overwhelming majority of insulator boundary elements. Nevertheless, it is striking that of the limited number of distal interactions we observed, many of them involve genes that are transcriptionally active during stage 5 of embryogenesis. This raises the possibility that these interactions may be stage or tissue-specific regulatory phenomena, and that more may be present in other tissues, developmental time points, or conditions.

Discussion

Several Hi-C studies in flies have identified enrichments of insulator proteins at TAD boundaries (Sexton et al., 2012; Ulianov et al., 2016; Eagen et al., 2015; Mourad and Cuvier, 2016) These studies varied in their resolution (due to use of 4- vs. 6-cutter enzymes and sequencing depth), methods (solution vs. in situ Hi-C), and, critically, in the methods used to identify TAD boundaries. As a result, each study relied on distinct sets of boundaries for analyses of the molecular features of these structures. We explored several methods to identify topological domains and associated boundaries and found that no single approach was sufficient to exhaustively identify all of these features in the genome. Rather, by using a combination of visual inspection of Hi-C maps at a large number of loci, unbiased hand-calling, and computational searches, we consistently observed a very close, two-way association between sites of combinatorial insulator protein binding (insulators) and the boundaries between topological domains. This result supports prior studies which found enriched insulator protein binding at topological boundaries, and extends this finding by localizing boundaries to discrete insulator elements. Hi-C data are exceptionally complex and reveal many layers of genomic organization, and we suspect that many questions in this field will only be resolved by the combined work of multiple groups using distinct analysis strategies and techniques.

Our most intriguing finding is the association of TAD boundaries with polytene interbands. The implication that these elements are decompacted, extended chromatin regions provides an attractive model in which simple physical separation explains multiple activities associated with insulators, including the ability to block enhancer-promoter interactions, prevent the spread of silenced chromatin, and organize chromatin structure.

A number of prior observations are consistent with the identity of insulators/boundaries as interbands. First, estimates suggest that there are ~5000 interbands constituting 5% of genomic DNA, with an average length of 2 kb (Zhimulev, 1996; Vatolina et al., 2011a), numbers that are in line with our estimates of boundary element length and number. Second, interbands are associated with insulator proteins, with CP190 appearing to be a constitutive feature of all or nearly all interbands (Pai et al., 2004; Gerasimova et al., 2007), which is precisely what we observe for boundary elements. Third, interbands and boundary elements are highly sensitive to DNase digestion (Vatolina et al., 2011b). Fourth, interbands have been shown to contain the promoters and 5’ ends of genes (Jamrich et al., 1977; Sass, 1982; Sass and Bautz, 1982a, 1982b; Rykowski et al., 1988), and we see a strong enrichment for promoters oriented to transcribe away from boundaries, which would place upstream regulatory elements within or near the interband. Finally, deletion of both isoforms of BEAF-32, the second-most highly enriched insulator protein at boundary elements, results in polytene X chromosomes that exhibit loss of banding and are wider and shorter than wild type, consistent with a loss of decompacted BEAF-32-bound regions (Roy et al., 2007). It is possible that interbands in polytene chromosomes result from multiple underlying molecular phenomena, but we believe it is likely that decompacted insulator elements constitute a significant fraction of these structures.

While we and others have not observed frequent looping of insulators in Hi-C data from fly tissue, our model of chromatin compaction at insulators is not mutually exclusive with a role for looping in the function of some insulators. Indeed, we have observed a limited set of cases in which interactions between boundaries seem to organize special genome structures with, at least in the case of the Scr locus, clear functional implications. It is likely that additional boundary-associated distal interactions will be found in other tissues and stages of fly development. However, we emphasize that these interactions are rare and do not appear to be general features of the function of boundary elements.

Conclusions

The data presented here offer a picture of the structure of the interphase chromatin of Drosophila that attempts to unify years of studies of polytene chromosomes with modern genomic methods (Figure 9). In this picture, interphase chromatin consists of alternating stretches of compacted, folded chromatin domains separated by regions of decompacted, stretched regions. The compacted regions vary in size from a few to hundreds of kilobases and correspond to both polytene band regions and TADs in Hi-C data. Decompacted regions that separate these domains are short DNA elements that are defined by the strong binding of insulator proteins and correspond to polytene interbands and TAD boundaries (insulators). An intuitive view of this structure in a non-polytene context might resemble the well-worn ‘beads on a string’, in which insulator/interband regions are the string and bands/TADs form beads of various sizes. Future work, including experimental manipulation of the sequences underlying these structures, will focus on validating and refining this model, exploring how it fits into hierarchical levels of genome organization, and understanding its implications for genome function.

A chromatin extension model of insulator function.

We propose a model in which insulators achieve domain separation by lowering the compaction ratio of bound chromatin, thereby converting the short lengths of insulator DNA (measured in base pairs) into large relative physical distances. By increasing the distance between domains, this model plausibly explains how insulators can achieve their diverse effects, including organizing chromatin structure, blocking enhancer-promoter interactions, and limiting the spread of chromatin silencing states.

https://doi.org/10.7554/eLife.29550.033

Materials and methods

Embryo collection, sorting, and sectioning

OregonR strain D. melanogaster (RRID:FlyBase_FBst1000080) embryos were collected on molasses plates seeded with fresh yeast paste from a population cage and aged to appropriate developmental stages, all at 25°C. Embryos were washed into nitex meshes, dechorionated by treatment with dilute bleach for 2 min, dipped briefly (15–20 s) in isopropanol, and gently rocked in fixative solution of (76.5% hexanes, 5% formaldehyde in 1x PBS) for 28–30 min. Embryos were then thoroughly washed in PBS with 0.5% triton and stored for no more than fivethree days at 4°C. For sample HiC-2/4, embryos were inspected under a light microscope to confirm that the vast majority corresponded to early cellularized blastoderm, and approximately 4000 embryos were used in the Hi-C protocol. For samples HiC-10, 12, 13–16, fixed embryos were hand-sorted under a light microscope as described in (Harrison et al., 2011), using morphological markers to identify early cellularized embryos (nc14, stage 5). For whole embryo experiments, sorted embryos were placed directly into the Hi-C protocol, with no more than 3 days having elapsed since fixation.

For sectioned embryos, hand-sorted embryos of precise developmental stages were first arranged in rows on a block of 1% agarose with bromophenol blue in a shared anterior-posterior orientation, with between 20–40 embryos per block. Aligned embryos were then transferred to the bottom of a plastic embedding mold (Sigma Aldrich E6032), the bottom of which had previously been coated with hexane glue, carefully keeping track of the anterior-posterior orientation of embryos by marking the cup with marker. Embryos were covered with clear frozen section compound (VWR 95057–838) and frozen at −80°C for up to two months. Frozen blocks wer4)e retrieved from the freezer and embryos rapidly sliced at approximately the mid-point by hand using a standard razor blade under a dissecting microscope. Anterior and posterior halves were separately transferred to microcentrifuge tubes containing ~200 µL PBS with 0.5% triton using an embryo pick (a tool of mysterious provenance that appears to be a clay sculpting tool). Successful transfer was confirmed visually by the presence of blue embryos which had absorbed bromophenol blue from the agarose block. Between transferring anterior and posterior halves, the pick was washed thoroughly with water and ethanol, and rubbed vigorously with kimwipes. We note that anterior and posterior half samples are precisely matched: samples HiC-13 and 14 contain the anterior and posterior halves (respectively) of the same embryos, and the same is true for HiC-15 and 16.

Hi-C

Experimental procedure

Hi-C experiments were conducted as described in Rao (Rao et al., 2014), with slight modifications. For completeness, we describe the detailed protocol: Embryos (or halves) were suspended in 1X NEB2 buffer (NEB B7002) and homogenized on ice by douncing for several minutes each with the loose and tight dounces. Insoluble material (including nuclei) was pelleted by spinning for 5 min at 4500 x g in microcentrifuge cooled to 4°C (all wash steps used these conditions for pelleting). Nuclei were washed twice with 500 µL of 1x NEB2 buffer and then suspended in 125 µL of the same. 42.5 µL of 2% SDS was added and tubes placed at 65°C for 10 min, then returned to ice, followed by addition of 275 µL of 1x NEB2 buffer and 22 µL of 20% Triton X-100, then incubated at room temperature for 5 min. Samples were digested overnight with 1500 units of MboI by shaking at 37°C. The next day, samples were washed twice with 1X NEB2, resuspended in 100 µL 1X NEB2, and 15 µL of fill-in mix (1.5 µL 10x NEB2, 0.4 µL each of 10 mM dATP, dGTP, dTTP, 9 µL 0.4 mM biotin-14-dCTP, 2.5 µL 5 U/µL Klenow (NEB M0210), 1 µL water) was added, followed by 1.5 hr at 37°C. Samples were then washed twice with 500 µL 1X ligation buffer (10X: 0.5 M Tris-HCl pH7.4, 0.1M MgCl2, 0.1M DTT), resuspended in 135 µL of the same, then supplemented with 250 µL of ligation mix (25 µL 10x ligation buffer, 25 µL 10% Triton X-100, 2.6 µL 10 mg/ml BSA, 2.6 µL 100 mM ATP, 196 µL water) and 2000 units of T4 DNA ligase (NEB M0202T) and incubated for 2 hr (or overnight) at room temperature. An additional 2000 units of ligase were added, followed by another 2 hr at room temperature. Cross-link reversal was carried out by adding 50 µL of 20 mg/mL proteinase K and incubating overnight at 65°C. An additional 50 µL proteinase K was then added followed by a 2 hr 65°C incubation. 0.1 volumes of 3M NaCl and 2 µL of glycoblue (Thermo Fisher AM9515) were added, then samples were extracted once with one volume of phenol pH 7.9, once with phenol-chloroform pH7.9, then precipitated with 3 volumes of EtOH. Washed pellets were resuspended in 130 µL water and treated with 1 µL of RNase A for 15 min at 37°C. DNA was fragmented using the Covaris instrument (Covaris, Woburn, MA) with peak power 140.0, duty factor 10.0, cycles/burst 200 for 80 s. Samples are brought to 300 µL total volume with water.

75 µL of Dynabeads MyOne Streptavidin C1 beads (Thermo Fisher 65001) were washed twice with 400 µL of tween wash buffer (TWB) (2X binding buffer [BB]: 100 µL of 1M Tris-HCl pH8, 20 µL 0.5 M EDTA, 4 mL of 5M NaCl, 5.88 mL water; TWB: 5 ml 2X binding buffer, 50 µL 10% Tween, 4.95 µL water), resuspended in 300 µL 2X BB, then added to 300 µL DNA. Samples were rocked at room temperature for 15 min, then washed once with TWB, twice with 1X BB, reclaimed on magnetic stand and resuspended in 100 µL 1X T4 DNA ligase buffer. Samples were then supplemented with end-repair mix (78 µL water, 10 µL 10X T4 DNA ligase buffer with ATP, 2 µL 25 mM dNTPs, 1 µL 10 U/µL T4 PNK (Thermo Fisher EK0031), 2 µL 5 U/µL Klenow, 3 µL 3 U/µL T4 DNA polymerase (Thermo Fisher EP0061)), incubated 30 min at room temp, washed as before, washed once with 100 µL 1X NEB2, and resuspended in 90 µL 1X NEB2. dA overhangs were added by adding 2 µL 10 mM dATP and 1 µL Klenow exo minus (NEB M0212S), incubating at 37°C for 30 min. Beads were washed as before, washed once with 100 µL 1X Quickligase (NEB M2200S) buffer, resuspended in 50 µL 1X Quickligase buffer, then supplemented with 3 µL Illumina adaptors and 1 µL Quickligase. Samples were incubated 15 min at room temperature, then washed twice with TWB, twice with 1X BB, twice with 200 µL TLE, and resuspend in 50 µL TLE. Beads are stable at 4°C, but were always amplified quickly. 100 µL (or more) of phusion PCR reaction was prepared (50 µL 2X Phusion master mix, 1 µL 100 µM forward primer [5-AATGATACGGCGACCACCGAG-3], 1 µL 100 µM reverse primer [5-CAAGCAGAAGACGGCATACGAG-3], 10 µL of beads with Hi-C library attached, 38 µL water). Reaction was mixed well and split into separate 12 µL reactions. Thermocycler conditions were 16 cycles of 98°C for 30 s, 63°C for 30 s, 72°C for 2 m. Reactions were pooled and loaded on a 2% agarose gel. Fragments corresponding to an insert size of ~300 bp (amplicon size of 421 bp) were excised from the gel, purified with the Zymo Gel DNA Recovery Kit (D4001T, Zymo), and submitted for sequencing at the Vincent J. Coates Genomic Sequencing Laboratory (Berkeley, CA).

Read processing and mapping

Our analysis routine was adapted by examining the approaches of multiple groups (Lieberman-Aiden et al., 2009; Sexton et al., 2012; Rao et al., 2014; Crane et al., 2015) in addition to procedures we developed independently. All analysis was performed with custom Python, R, and Perl scripts (Stadler, 2017; copy archived at https://github.com/elifesciences-publications/Stadlerlab-hi-c) except where noted. Single-ends of demultiplexed reads were separately mapped using Bowtie (Langmead et al., 2009) (parameters: -m1 --best --strata) to the D. melanogaster genome dm6 (R6.17). Due to the formation of chimeric reads intrinsic to the Hi-C procedure, reads can fail to properly map if the ligation junction lies within the 100 bp read. To address this, we used an iterative mapping procedure, in which we began by mapping the first 20 nt of the reads (using Bowtie’s --trim3 feature). Unique mappings were kept, reads that failed to map were stored, and the procedure was repeated on the multiply-mapping reads, incrementing the length of sequence to map by 7 nt each round (attempt to uniquely map using first 20, first 27, first 34...). We found that this method gave 5–10% increases in yield of mapped reads over a procedure in which we attempted to explicitly detect and trim ligation junctions from reads. Uniquely mapping reads from all iterations were collated as a single file.

Uniquely-mapping single-ends were paired based on read identity, and only pairs with two uniquely-mapping ends were retained. Duplicate reads that shared identical 5’ mapping positions were removed. Resulting paired, collapsed, uniquely mapping reads were then inspected for quality. Primary indicators of successful Hi-C libraries were the distance distribution of mapped pairs and the relative frequencies of reads in the four orientations described by Rao et al., 2014, in-in, in-out, out-in, and out-out. In all of our libraries, we detect some ~3–15% reads that appear to be simple genomic sequence, not the result of a Hi-C ligation event. These reads are readily detected by examining the size distributions of in-out reads (the orientation expected from standard genomic sequence) compared with the other three orientations. The in-out reads have a unique hump of reads showing a distance distribution of ~150–500 bp, varying slightly from sample to sample. In-out reads pairs spanning less than 500 bp were removed from further analysis.

Computational topological boundary detection

We explored a number of ways of identifying boundaries from directionality data. In the end, the most robust was to use a simple heuristic that at a boundary, by definition, regions to the left show left-bias and regions to the right show right bias. While attempts to derive a boundary score from a comparison of directionality scores upstream and downstream showed susceptibility to noise and artifacts, requiring expected upstream and downstream behavior allowed robust detection of sets of boundary elements. We describe the complete procedure below.

Read counts were assigned to 500 bp bins for all genomic bin combinations within 500 kb of the diagonal. Local directionality scores were calculated for each bin by summing the counts linking the bin to regions in a window encompassing the genomic regions between 1 and 15 kb from the bin (skipping the two proximal 500 bp bins, summing the next 28) upstream and downstream, then taking the log (10) ratio of downstream to upstream. These parameters were determined by visually comparing local directionality scores from a range of inputs to Hi-C heat maps for a number of genomic regions, identifying parameters in which directionality transitions reflected boundaries evident in the heat maps. We observed high levels of noise in the directionality metric in regions of low read coverage. To suppress these noisy signals, we devised a weighted local directionality score to weight these scores based on the total number of reads used in the calculation. We experimented with a variety of scaling factors a such that w = [read count]^a and found that a weighting of a = 0.5 worked well to reduce signal from low-read regions. From these directional scores, sites were first selected for which the mean directionality score of the 5 adjacent upstream bins was less than −2, and the mean for the 5 adjacent downstream bins was greater than 2. Boundary scores were assigned to resulting bins by subtracting the sum of the directionality scores for the 5 adjacent upstream bins from the 5 adjacent downstream bins. An issue with this scoring system is that bins that lack MboI sites can cause inflated directionality scores in adjacent regions. To address this, we simply assigned a boundary score of 0 to any bin with more than one such bin in its radius. The resulting distribution of boundary scores is dominated by series of consecutive bins with large boundary score maximums, which is uninformative since these scores are essentially derived from the same data (window shifted by one bin). We therefore merged adjacent bins that passed the cutoff and selected only the bin with the maximum boundary score within a contiguous block. By sorting the resulting table on the boundary score, we were able to select sets of candidate boundaries of various strengths for analysis.

In additional to these computationally-identified boundary locations, we manually called boundaries for the entire genome. An R script serially displayed Hi-C heat maps of 250 kb genomic windows and recorded the genomic coordinates of mouse clicks made at visually-identified boundaries. The human caller was unaware of any features of the regions examined other than the Hi-C maps, and was unaware of the locations being displayed in a given plot.

Sequence analysis

We used simple custom Python scripts to count the occurrences of all words of length 4, 5, 6 and 7 in 500 bp windows from 10,000 bp upstream to 10,000 bp downstream of the 500 bp window identified as a boundary. We then computed a simple enrichment score for each unique word equal to the counts of that word and its reverse complement in the boundary divided by the mean counts for the word and its reverse complement in the remaining windows. We noticed that many of the words identified as enriched in this analysis were also enriched in the 500 bp bins immediately flanking the boundary. We therefore updated our enrichment score for each word to be the mean of the counts of the word and its reverse complement in the boundary and the 500 bp bins immediately adjacent to it (three bins in total) divided by the mean counts of the word and its reverse complement in the remaining 38 bins. Counts and scores for all words are provided in the supplemental materials.

ATAC-seq

Experimental procedure

Early nc14 embryos were placed in ATAC-seq lysis buffer (Buenrostro et al., 2013) without detergent, with 5% glycerol added. Embryos were then taken out of the freezing solution and placed onto a glass slide which was then put on dry ice for 2 min. Once embryos were completely frozen, the glass slide was removed and embryos were sliced with a razor blade chilled in dry ice. Once sliced embryo halves were moved to tubes containing ATAC-seq lysis buffer with 0.15 mM spermine added to help stabilize chromatin. Embryo halves were then homogenized using single use plastic pestles. IGEPal CA-630 was added to a final concentration of 0.1%. After a 10 min incubation nuclei were spun down and resuspended in water. Twenty halves were added to the transposition reaction containing 25 µl of 2x TD buffer (Illumina), and 2.5 ul of Tn5 enzyme (Illumina) and the reaction was incubated at 37°C for 30 min as in Buenrostro et al. (2013). Transposed DNA was purified using Qiagen Minelute kit. Libraries were then amplified using phusion 2x master mix (NEB) and indexed primers from Illumina. Libraries were then purified with Ampure Beads and sequenced on the Hiseq4000 using 100 bp paired end reads.

Analysis

Fastq files were aligned to the Drosophila Dm3 genome with Bowtie2 (Langmead and Salzberg, 2012) using the following parameters: −5 5–3 5 N 1 -X 2000 --local --very-sensitive-local. Sam files were then sorted and converted to Bam files using Samtools (Li et al., 2009), only keeping mapped, properly paired reads with a MAPq score of 30 or higher using -q 30. Bams were then converted to Bed files with bedtools and shifted using a custom shell script to reflect a 4 bp increase on the plus strand and a 5 bp decrease on the minus strand as recommended by Buenrostro et al. (2013). Finally shifted bed files were converted into wig files using custom scripts and wig files which were uploaded to the genome browser. Wig files were normalized to reflect 10 million mapped reads.

Sample size determination

No explicit statistical method was used to compute sample size. All unique experiments were prepared in duplicate.

References

  1. 1
    Sur La Structure Du Noyau Des Cellules Salivares Chez Les Larves de Chironomus
    1. EG Balbiani
    (1881)
    Zoologischer Anzeiger 4:637662–641666.
  2. 2
    Sur La Structure Intime Du Noyau Du Loxophyllum Meleagris
    1. EG Balbiani
    (1890)
    Zoologischer Anzeiger 13:110–15.
  3. 3
  4. 4
  5. 5
    Molecular and cytogenetical characterization of the 1OA1-2 band and adjoining
    1. ES Belyaeva
    2. IF Zhimulev
    (1994)
    Genetics 136:1063–1073.
  6. 6
  7. 7
  8. 8
  9. 9
    Salivary chromosome maps with a key to the banding of the chromosomes of Drosophila melanogaster
    1. CB Bridges
    (1934)
    The Journal of Heredity pp. 60–64.
  10. 10
  11. 11
    Genetic and molecular complexity of the position effect variegation modifier mod(mdg4) in Drosophila
    1. K Büchner
    2. P Roth
    3. G Schotta
    4. V Krauss
    5. H Saumweber
    6. G Reuter
    7. R Dorn
    (2000)
    Genetics 155:141–157.
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
    Studies of nuclear and cytoplasmic behaviour during the five mitotic cycles that precede gastrulation in Drosophila embryogenesis
    1. VE Foe
    2. BM Alberts
    (1983)
    Journal of Cell Science 61:31–70.
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
    A photographic representation and interpretation of the polytene chromosomes of Drosophila melanogaster salivary glands
    1. G Lefevre
    (1976)
    The Genetics and Biology of Drosophila.
  57. 57
    Developmental Biology Using Purified Genes
    1. E Lewis
    (1981)
    189–208, Developmental genetics of the bithorax complex in Drosophila, Developmental Biology Using Purified Genes, 10.1016/B978-0-12-137420-4.50022-X.
  58. 58
  59. 59
  60. 60
  61. 61
  62. 62
  63. 63
  64. 64
    Genetic Variations of Drosophila Melanogaster
    1. DL Lindsley
    2. EH Grell
    (1968)
    Carnegie Inst. Wash. Publication.
  65. 65
  66. 66
  67. 67
    Polycomb Group Repression Is Blocked by the Drosophila Suppressor of Hairy-Wing [SU(HW)] Insulator
    1. DR Mallin
    2. JS Myung
    3. JS Patton
    4. PK Geyer
    (1998)
    Genetics 148:331–339.
  68. 68
  69. 69
  70. 70
  71. 71
  72. 72
  73. 73
  74. 74
  75. 75
  76. 76
    The Morphology of the Third Chromosome in the Salivary Gland of Rosophila Melanogaster and a New Cytological Map of This Element
    1. TS Painter
    (1935)
    Genetics 20:301–326.
  77. 77
  78. 78
  79. 79
  80. 80
  81. 81
    High-resolution TADs reveal DNA sequences underlying Genome organization in flies
    1. F Ramirez
    2. V Bhardwaj
    3. J Villaveces
    4. L Arrigoni
    5. BA Gruening
    6. KC Lam
    7. B Habermann
    8. A Akhtar
    9. T Manke
    (2017)
    bioRxiv, 10.1101/115063.
  82. 82
  83. 83
  84. 84
    The su(Hw) Protein Insulates Expression of the Drosophila Melanogaster White Gene from Chromosomal Position-Effects
    1. RR Roseman
    2. V Pirrotta
    3. PK Geyer
    (1993)
    The EMBO Journal 12:435–442.
  85. 85
  86. 86
  87. 87
  88. 88
  89. 89
  90. 90
  91. 91
    Enhancer blocking by the drosophila gypsy insulator depends upon insulator anatomy and enhancer strength
    1. KC Scott
    2. AD Taubman
    3. PK Geyer
    (1999)
    Genetics 153:787–798.
  92. 92
  93. 93
    Chromatin insulator elements block the silencing of a target gene by the Drosophila polycomb response element (PRE) but allow trans interactions between PREs on different chromosomes
    1. CJ Sigrist
    2. V Pirrotta
    (1997)
    Genetics 147:209–221.
  94. 94
  95. 95
    Hi-C, version 95e304b
    1. M Stadler
    (2017)
    Bitbucket.
  96. 96
  97. 97
  98. 98
    Identical functional organization of nonpolytene and polytene chromosomes in Drosophila melanogaster
    1. TY Vatolina
    2. LV Boldyreva
    3. OV Demakova
    4. SA Demakov
    5. EB Kokoza
    6. VF Semeshin
    7. VN Babenko
    8. FP Goncharov
    9. ES Belyaeva
    10. IF Zhimulev
    (2011)
    PLoS ONE, 6, 10.1371/annotation/45b44e2a-c751-418b-bbb7-7023998abdfc, 22022482.
  99. 99
  100. 100
  101. 101
  102. 102
  103. 103
  104. 104
  105. 105
  106. 106
    Morphology and structure of polytene chromosomes
    1. IF Zhimulev
    (1996)
    Advances in Genetics 34:1–490.

Decision letter

  1. Allan C Spradling
    Reviewing Editor; Howard Hughes Medical Institute, Carnegie Institution for Science, United States

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "Convergence of topological domain boundaries, insulators, and polytene interbands in the early Drosophila embryo" for consideration by eLife. Your article has been favorably evaluated by Jessica Tyler (Senior Editor) and three reviewers, one of whom is a member of our Board of Reviewing Editors. The reviewers have opted to remain anonymous.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

Stadler et al. present a HiC analysis of Drosophila chromatin from nuclear cycle 12 and cycle 14 embryos in order to address a variety of questions, with an emphasis on how insulators are related to chromatin topological domains (TADs). They conclude that TAD boundary regions are enriched for insulator binding sites, contain promoters, and are marked by accessible chromatin. The TAD boundary regions are found to be constant despite differences in gene expression between the anterior and posterior of the embryo, and they are detectable prior to the onset of zygotic gene expression. By comparing TAD boundaries to sites of insulator protein binding, and to very limited available information on polytene interbands, the authors propose a model of insulator function and suggest that their work "unifies years of polytene chromosome research.

Essential revisions:

Unfortunately, there are significant problems with the supporting data, requiring major revisions. The datasets upon which these conclusions are based (19 or 88 million contacts) appear to be significantly smaller than in recent publications (Hug et al.2017: – 400 million contacts; or even in older studies Corces 2012- 120 million contacts.) Barring some misunderstanding, neither 19 million or 88 million contacts would be adequate to achieve the claimed resolution of less than 1 kb. The questions addressed in the manuscript depend on achieving a high level of resolution, consequently, for this paper to be considered further, the authors need to address this problem. Additional sequencing of their libraries or combination of data with published datasets are possibilities to consider. If this problem can be resolved, and the resulting higher resolution analysis still supports the correlation between TAD boundaries and interbands, then a revised version of the paper will be considered for publication.

In a revised version, the claimed correspondence between TAD boundaries and polytene interbands needs to be solidified. The possible confounding effect of zones of active transcription should be addressed. They should acknowledge the highly tentative nature of a theory that has been tested at the very most in 5 out of > 5,000 cases.

Reviewer #1:

First, the claim that the data achieves a resolution of 1kb or less needs to be better justified. If the authors' HiC data is so accurate, why were only 947 high confidence junctions defined? Since the most robust domain junctions are those associated with large bands, why did the authors not identify all the junctions of those bands? How could the very fine junctions between almost invisible bands in region 10A be defined, when so many domains in the genome as a whole were not?

Second, it is not clear from the data presented in Figure 2 and supplements, that one can robustly define regions whose function is to separate domains from simple uncertainty in the domain edges. The surprising uniformity in size (500-1000bp) of the authors elements might have been an artifact of the requirement for a major directionality change. The "insulator" proteins used in the definition of these junctions do not appear to be very specific, as the authors own and third party data shows that these proteins bind at multiple sites in the genome, most of which are not at domain boundaries. The paper should calculate the numerical enrichment of each of these proteins at domain junctions vs. other sites in the genome.

The data presented in cytogenetic region 10A-10B, which has been intensively studied in the Zhimulev laboratory, was quite interesting, although it falls a little short of the "perfect" correspondence claimed by the authors. In particular, in Figure 5, the junction between the heavy band 10A1-2 and 10A3 (a virtually invisible band in polytenes) was not actually detected in the hand or high confidence boundary calls, but was arbitrarily added to the figure as a dashed line so it would fit expectations. As a result, the next junction, the strongest in the region, is scored as the junction between 10A3 and 10A4, i.e. between an invisible band and a weak band. It seems more likely that this strong boundary is actually the end of 10A1-2, and that there are only 4 rather than 5 domains between the two major bands in this region. This would be consistent with the authors own detection of only 4 boundaries in their unbiased and in the high confidence datasets. The authors separate band 10A6 from 10A7 following Zhimulev, even though he states that this separation is almost never observed, nor did Bridges report separate bands in his revised map. There is no DNAase or convincing directionality data for this separation. Since HiC would average conformational differences between cells, the data should reflect the common configuration, not a rare configuration. The authors should clearly state in the text that they are not using their own boundary measurements but making a special interpretation specifically for this experiment, rather than hiding that admission at the end of the figure legend.

Finally, the authors should comment on the correspondence between the boundaries and polytene bands. It stands to reason that a method capable of resolving interbands would need to have superb resolution of bands. How well do the boundaries defined in this manuscript fit with previous cytological and HiC reports for major bands? It is an interesting question, not only for its relevance to the question of HiC reproducibility between labs, but because different tissue sources, including salivary glands, tissue culture cells and embryos have been used. Is it possible that different tissues have the same general organization into TADs, but the detailed boundaries of the TADs differ in some cases?

Reviewer #2:

The data largely are clear and compelling except for a few issues:

1) The ChIP data for insulator binding sites are from 0-12 hour embryos, a much broader developmental window than the two time points examined for the TAD mapping. Thus authors should discuss how this may or may not affect their conclusions. Given that the insulator ChIP data are not from the same restricted developmental time as the TAD maps, it would be useful to discuss how the boundary elements mapped here compare to sites of localization of CHRIZ (Chromator), even though they have been determined in S2 and Kc cells. Zhimulev et al. concluded that CHRIZ binding was the most diagnostic feature of interbands, so this could serve as further confirmation of the coincidence of boundary regions with interbands.

2) Please provide color scale bars for the heatmaps of contact probabilities. This will aid readers in comparing Figures 1 and 7. Please provide these also for Figure 3 and the other protein binding figures.

3) The looped regions shown in Figure 8 and the corresponding supplemental figures are difficult for readers unfamiliar with Hi-C to visualize. The interaction in Figure 8B is clear, but could the authors provide a clearer graphic to explain the loops in Figure 8A? Could they also highlight the diagnostic loop signatures in the supplemental figures?

Reviewer #3:

The manuscript is interesting but many of the observations, including chromatin 3D organization in Drosophila embryos and association of boundaries with insulator proteins, have been published previously. The following are other more specific concerns:

1) Authors should explain how they calculate the resolution of Hi-C data. Because the frequency of interactions decays with distance, authors should keep in mind that resolution will depend on distance. If I understand Table S1 Sample Information in Supplementary file 1 correctly, the authors made 2 libraries from nc12 and 3 libraries from nc14. I used Juicer to process the data for these two libraries and the results are as follows:

Experiment description: nc12

Sequenced Read Pairs: 93,483,816

Normal Paired: 65,205,173 (69.75%)

Chimeric Paired: 16,627,312 (17.79%)

Chimeric Ambiguous: 6,549,875 (7.01%)

Unmapped: 5,101,456 (5.46%)

Ligation Motif Present: 45,319,979 (48.48%)

Alignable (Normal+Chimeric Paired): 81,832,485 (87.54%)

Unique Reads: 26,913,205 (28.79%)

PCR Duplicates: 54,884,068 (58.71%)

Optical Duplicates: 35,212 (0.04%)

Library Complexity Estimate: 28,537,160

Intra-fragment Reads: 760,888 (0.81% / 2.83%)

Below MAPQ Threshold: 6,649,224 (7.11% / 24.71%)

Hi-C Contacts: 19,503,093 (20.86% / 72.47%)

Ligation Motif Present: 7,811,696 (8.36% / 29.03%)

3' Bias (Long Range): 84% – 16%

Pair Type% (L-I-O-R): 25% – 25% – 25% – 25%

Inter-chromosomal: 1,545,226 (1.65% / 5.74%)

Intra-chromosomal: 17,957,867 (19.21% / 66.73%)

Short Range (<20Kb): 8,504,570 (9.10% / 31.60%)

Long Range (>20Kb): 9,451,838 (10.11% / 35.12%)

Experiment description: nc14

Sequenced Read Pairs: 393,928,165

Normal Paired: 144,374,310 (36.65%)

Chimeric Paired: 24,393,845 (6.19%)

Chimeric Ambiguous: 9,682,276 (2.46%)

Unmapped: 6,884,964 (1.75%)

Ligation Motif Present: 127,838,532 (32.45%)

Alignable (Normal+Chimeric Paired): 168,768,155 (42.84%)

Unique Reads: 119,180,726 (30.25%)

PCR Duplicates: 49,523,419 (12.57%)

Optical Duplicates: 64,010 (0.02%)

Library Complexity Estimate: 227,824,974

Intra-fragment Reads: 1,146,075 (0.29% / 0.96%)

Below MAPQ Threshold: 29,415,704 (7.47% / 24.68%)

Hi-C Contacts: 88,618,947 (22.50% / 74.36%)

Ligation Motif Present: 41,585,911 (10.56% / 34.89%)

3' Bias (Long Range): 83% – 17%

Pair Type% (L-I-O-R): 25% – 25% – 25% – 25%

Inter-chromosomal: 5,442,733 (1.38% / 4.57%)

Intra-chromosomal: 83,176,214 (21.11% / 69.79%)

Short Range (<20Kb): 49,153,740 (12.48% / 41.24%)

Long Range (>20Kb): 34,004,040 (8.63% / 28.53%)

This suggests that nc12 has 19 million HiC contacts and nc14 has 88.6 million HiC contacts. It is impossible to have a 500 bp resolution, as the authors claim, with these numbers of reads. The authors state that "We conclude that these Hi-C are of high quality and reproducibility". There are currently two HiC datasets from Drosophila available in Juicebox, and both of them have around 1 billion mapped reads. The HiC libraries made with Drosophila embryos by Hug et al. have around 400 million mapped reads. Therefore, the libraries presented in this study cannot be described as high quality compared to others already available. Furthermore, the nc12 library is of much lower quality than the nc14 library. This is a critical issue when interpreting the results.

2) Authors should use the dm6 genome instead of dm3.

3) It is strange that interbands in polytene chromosomes correspond to boundaries. These structures have always been thought to contain active genes. Hug et al. identify very small domain that they also claim are boundaries. In their Hi-C data (see for example Figure 1B), once can observed interactions between these small domains that resemble compartmental interactions in mammalian Hi-C data, suggesting that these small domains do not correspond to boundaries, they correspond to active domains. These features are not observed in the Hi-C data presented in this manuscript, which may be due to the fewer reads in the Hi-C data. Authors need to consider this and try to differentiate between effects of transcription and effects of insulator proteins.

4) Authors cannot directly compare features of nc12 and nc14 embryos unless they have the same number of Hi-C contacts in both samples. I think it is critical to obtain more reads for the nc12 libraries.

https://doi.org/10.7554/eLife.29550.040

Author response

Essential revisions:

Unfortunately, there are significant problems with the supporting data, requiring major revisions. The datasets upon which these conclusions are based (19 or 88 million contacts) appear to be significantly smaller than in recent publications (Hug et al.2017: – 400 million contacts; or even in older studies Corces 2012- 120 million contacts.) Barring some misunderstanding, neither 19 million or 88 million contacts would be adequate to achieve the claimed resolution of less than 1 kb. The questions addressed in the manuscript depend on achieving a high level of resolution, consequently, for this paper to be considered further, the authors need to address this problem. Additional sequencing of their libraries or combination of data with published datasets are possibilities to consider. If this problem can be resolved, and the resulting higher resolution analysis still supports the correlation between TAD boundaries and interbands, then a revised version of the paper will be considered for publication.

There does appear to be some misunderstanding. Our original manuscript was based on 263,868,080 Hi-C linkages from nc14 embryos. A complete description of the read totals for individual datasets, including a full accounting of reads removed at various processing steps, was supplied in Table S2 Sequencing Statistics in Supplementary file 1. Further, the processing pipeline used was described in detail in the Materials and methods section.

We appreciate the interest and efforts of reviewer 3 to analyze our raw data. We are unsure how they arrived at the total of 88 million nc14 contacts, but have simplified the naming of the original raw data files to make their relationships clearer in hopes that future users of the data will be able to work with them optimally.

Despite this confusion, we felt it was prudent to generate additional sequencing data from these libraries, and performed additional sequencing of four of our nc14 libraries and now have a total of 451,978,433 Hi-C contacts from nc14 embryos. The overall character of the data is unchanged, and none of our original observations are altered. We believe this is sufficient to achieve sufficient resolution to support the conclusions of the manuscript.

In a revised version, the claimed correspondence between TAD boundaries and polytene interbands needs to be solidified. The possible confounding effect of zones of active transcription should be addressed. They should acknowledge the highly tentative nature of a theory that has been tested at the very most in 5 out of > 5,000 cases.

We agree and have added additional discussion of the possible effects of active transcription on TAD boundary structure. The fact that TAD boundaries are invariant between anterior and posterior embryo sections, show only partial correlation with RNA polymerase occupancy (compared to exceptionally strong correspondence with insulator protein binding), and that TAD structures are not significantly altered in the presence of transcriptional inhibitors (Hug et al.) argue that transcriptional activity is not a strong driver of TAD boundary formation. With respect to the band/interband pattern of the 10A-B region, only a subset of the promoters near the identified TAD boundaries appear to be transcribed in nc14 embryos, based on RNAPII occupancy and RNAseq. We have added additional commentary on the tentative nature of the TAD boundary-interband association and the limited availability of high-resolution polytene mapping data.

Reviewer #1:

First, the claim that the data achieves a resolution of 1kb or less needs to be better justified. If the authors' HiC data is so accurate, why were only 947 high confidence junctions defined? Since the most robust domain junctions are those associated with large bands, why did the authors not identify all the junctions of those bands? How could the very fine junctions between almost invisible bands in region 10A be defined, when so many domains in the genome as a whole were not?

In retrospect we agree that claims about the resolution of the data are fraught, especially because there is no clear standard for what it means to assign a numerical value to HiC resolution. We have therefore removed these specific claims about resolution from the manuscript and focus instead on what the data actually show.

The list of 947 high confidence junctions was not meant to be a comprehensive list of junctions. This set was chosen as a representative set to analyze the features of boundaries, not as a comprehensive list, and is a very conservative list generated from the union of the highest-scoring computationally-identified boundaries and hand calls. We have attempted to make this distinction clear in the text.

To the broader issue of calling boundaries: as with all genomic datasets it is difficult to define an exact number of features of a certain class, since one can identify them with various algorithms and choices must be made about cutoffs, parameters, etc. with different justifiable choices leading to wildly different numbers of identifiable features with differing degrees of overlap. We identified a much larger set (~3,200) of hand-called boundaries. The number of computationally-identified boundaries varies depending on the thresholds used, but is likely over ~2,500. We have included the full list of hand-called and computationally-identified boundaries in supplemental data, as well as supplying discussion about the complexities of boundary calling procedures and estimates of the total number of boundaries present in the genome.

Second, it is not clear from the data presented in Figure 2 and supplements, that one can robustly define regions whose function is to separate domains from simple uncertainty in the domain edges. The surprising uniformity in size (500-1000bp) of the authors elements might have been an artifact of the requirement for a major directionality change. The "insulator" proteins used in the definition of these junctions do not appear to be very specific, as the authors own and third party data shows that these proteins bind at multiple sites in the genome, most of which are not at domain boundaries. The paper should calculate the numerical enrichment of each of these proteins at domain junctions vs. other sites in the genome.

We agree that there is uncertainty in our ability to define the precise location of individual boundaries, but believe that collectively the data strongly support our conclusions about the correspondence between these regions. We have, as suggested, calculated the numerical enrichment and have added this to the paper.

The data presented in cytogenetic region 10A-10B, which has been intensively studied in the Zhimulev laboratory, was quite interesting, although it falls a little short of the "perfect" correspondence claimed by the authors. In particular, in Figure 5, the junction between the heavy band 10A1-2 and 10A3 (a virtually invisible band in polytenes) was not actually detected in the hand or high confidence boundary calls, but was arbitrarily added to the figure as a dashed line so it would fit expectations. As a result, the next junction, the strongest in the region, is scored as the junction between 10A3 and 10A4, i.e. between an invisible band and a weak band. It seems more likely that this strong boundary is actually the end of 10A1-2, and that there are only 4 rather than 5 domains between the two major bands in this region. This would be consistent with the authors own detection of only 4 boundaries in their unbiased and in the high confidence datasets. The authors separate band 10A6 from 10A7 following Zhimulev, even though he states that this separation is almost never observed, nor did Bridges report separate bands in his revised map. There is no DNAase or convincing directionality data for this separation. Since HiC would average conformational differences between cells, the data should reflect the common configuration, not a rare configuration. The authors should clearly state in the text that they are not using their own boundary measurements but making a special interpretation specifically for this experiment, rather than hiding that admission at the end of the figure legend.

We thank the reviewer for his/her critical eye, and we agree that the correspondence between our Hi-C maps and Zhimulev’s banding analysis should be clarified. First, it is important to point out that the band/interband assignments presented in Figure 5, and their genomic locations, are assigned according to Vatolina et al., not from our measurements. As the reviewer points out, our data do not strongly support the existence of the minor band 10A-3, and combined with the fact that it is virtually invisible in polytenes, we agree that it is likely not real and have indicated this in the text.

The separation of 10A6 and 10A7 is also a fascinating case. Unfortunately, the exact junction between these domains is obscured by a lack of MboI cut sites, an unfortunate limitation of our Hi-C data. Nevertheless, we do seem to observe evidence for the existence of these bands (the stronger signal within them and weaker signal for contacts between them). Additionally, while there is not a strong signal of DNase digestion at this putative boundary, there is strong binding of CP190, BEAF-32, and dCTCF. However, the data also clearly show significant interactions between these two putative domains (the light orange regions near the top of the pyramid). Indeed, this situation resembles a pattern we observe frequently in the genome as a whole: adjacent domains that are clearly distinct, yet show significant interaction with each other. The fact that the 10A6-7 junction is observed infrequently in polytene band analysis suggests that this interband may be somewhat variable, present on some fraction of chromosomes and absent in others. This is a very interesting possibility that we discuss in the revised manuscript.

Overall, we feel that these issues are consistent with a correspondence between polytene banding patterns and our Hi-C data in this region. A band (10A3) that is virtually invisible in polytene banding analysis is not detected by Hi-C, and a region that is often seen as a single band, and occasionally as separate bands, shows a Hi-C signature that would be expected from a region that existed in both states within a population of nuclei. Meanwhile, all of the interbands that are clearly visible in polytene spreads are clearly visible in Hi-C data.

Finally, the authors should comment on the correspondence between the boundaries and polytene bands. It stands to reason that a method capable of resolving interbands would need to have superb resolution of bands. How well do the boundaries defined in this manuscript fit with previous cytological and HiC reports for major bands? It is an interesting question, not only for its relevance to the question of HiC reproducibility between labs, but because different tissue sources, including salivary glands, tissue culture cells and embryos have been used. Is it possible that different tissues have the same general organization into TADs, but the detailed boundaries of the TADs differ in some cases?

We thank the reviewer for this suggestion, and took a closer look several regions of associated polytene bands and TADs that were previously examined by Eagen et al. In the revised manuscript, we provide a detailed examination of the region of band 22A1-2, which appeared as a single large TAD in Eagen’s data, but which had occasionally been observed to contain an interband. Our data reveal this locus to have a more complex topology, with two large TADs surrounding a middle region containing smaller domains and a fascinating looped region, which together are likely responsible for the observed interband. This analysis is included in Figure 6 in the revised manuscript.

Reviewer #2:

The data largely are clear and compelling except for a few issues:

1) The ChIP data for insulator binding sites are from 0-12 hour embryos, a much broader developmental window than the two time points examined for the TAD mapping. Thus authors should discuss how this may or may not affect their conclusions. Given that the insulator ChIP data are not from the same restricted developmental time as the TAD maps, it would be useful to discuss how the boundary elements mapped here compare to sites of localization of CHRIZ (Chromator), even though they have been determined in S2 and Kc cells. Zhimulev et al. concluded that CHRIZ binding was the most diagnostic feature of interbands, so this could serve as further confirmation of the coincidence of boundary regions with interbands.

We agree with the reviewer that the differences in the tissues used to generate Hi-C and ChIP data needs to be clarified, and have done so in the revised manuscript. We have also included CHRIZ in our analyses of the molecular features of our TAD boundaries, and found this protein to be a highly diagnostic feature of our boundaries.

2) Please provide color scale bars for the heatmaps of contact probabilities. This will aid readers in comparing Figures 1 and 7. Please provide these also for Figure 3 and the other protein binding figures.

We have included color scale bars in all relevant figures.

3) The looped regions shown in Figure 8 and the corresponding supplemental figures are difficult for readers unfamiliar with Hi-C to visualize. The interaction in Figure 8B is clear, but could the authors provide a clearer graphic to explain the loops in Figure 8A? Could they also highlight the diagnostic loop signatures in the supplemental figures?

We have revised Figure 8 and its supplements, and hope that the looping interactions are clearer.

Reviewer #3:

The manuscript is interesting but many of the observations, including chromatin 3D organization in Drosophila embryos and association of boundaries with insulator proteins, have been published previously. The following are other more specific concerns:

1) Authors should explain how they calculate the resolution of Hi-C data. Because the frequency of interactions decays with distance, authors should keep in mind that resolution will depend on distance. If I understand Table S1 Sample Information in Supplementary file 1 correctly, the authors made 2 libraries from nc12 and 3 libraries from nc14. I used Juicer to process the data for these two libraries and the results are as follows:

Experiment description: nc12

Sequenced Read Pairs: 93,483,816

Normal Paired: 65,205,173 (69.75%)

Chimeric Paired: 16,627,312 (17.79%)

Chimeric Ambiguous: 6,549,875 (7.01%)

Unmapped: 5,101,456 (5.46%)

Ligation Motif Present: 45,319,979 (48.48%)

Alignable (Normal+Chimeric Paired): 81,832,485 (87.54%)

Unique Reads: 26,913,205 (28.79%)

PCR Duplicates: 54,884,068 (58.71%)

Optical Duplicates: 35,212 (0.04%)

Library Complexity Estimate: 28,537,160

Intra-fragment Reads: 760,888 (0.81% / 2.83%)

Below MAPQ Threshold: 6,649,224 (7.11% / 24.71%)

Hi-C Contacts: 19,503,093 (20.86% / 72.47%)

Ligation Motif Present: 7,811,696 (8.36% / 29.03%)

3' Bias (Long Range): 84% – 16%

Pair Type% (L-I-O-R): 25% – 25% – 25% – 25%

Inter-chromosomal: 1,545,226 (1.65% / 5.74%)

Intra-chromosomal: 17,957,867 (19.21% / 66.73%)

Short Range (<20Kb): 8,504,570 (9.10% / 31.60%)

Long Range (>20Kb): 9,451,838 (10.11% / 35.12%)

Experiment description: nc14

Sequenced Read Pairs: 393,928,165

Normal Paired: 144,374,310 (36.65%)

Chimeric Paired: 24,393,845 (6.19%)

Chimeric Ambiguous: 9,682,276 (2.46%)

Unmapped: 6,884,964 (1.75%)

Ligation Motif Present: 127,838,532 (32.45%)

Alignable (Normal+Chimeric Paired): 168,768,155 (42.84%)

Unique Reads: 119,180,726 (30.25%)

PCR Duplicates: 49,523,419 (12.57%)

Optical Duplicates: 64,010 (0.02%)

Library Complexity Estimate: 227,824,974

Intra-fragment Reads: 1,146,075 (0.29% / 0.96%)

Below MAPQ Threshold: 29,415,704 (7.47% / 24.68%)

Hi-C Contacts: 88,618,947 (22.50% / 74.36%)

Ligation Motif Present: 41,585,911 (10.56% / 34.89%)

3' Bias (Long Range): 83% – 17%

Pair Type% (L-I-O-R): 25% – 25% – 25% – 25%

Inter-chromosomal: 5,442,733 (1.38% / 4.57%)

Intra-chromosomal: 83,176,214 (21.11% / 69.79%)

Short Range (<20Kb): 49,153,740 (12.48% / 41.24%)

Long Range (>20Kb): 34,004,040 (8.63% / 28.53%)

This suggests that nc12 has 19 million HiC contacts and nc14 has 88.6 million HiC contacts. It is impossible to have a 500 bp resolution, as the authors claim, with these numbers of reads. The authors state that "We conclude that these Hi-C are of high quality and reproducibility". There are currently two HiC datasets from Drosophila available in Juicebox, and both of them have around 1 billion mapped reads. The HiC libraries made with Drosophila embryos by Hug et al. have around 400 million mapped reads. Therefore, the libraries presented in this study cannot be described as high quality compared to others already available. Furthermore, the nc12 library is of much lower quality than the nc14 library. This is a critical issue when interpreting the results.

As discussed above there is some misunderstanding about our data. Data with our original submission contained > 250,000,000 contacts at nc14 and we have supplemented it bring the total to > 450,000,000. We have removed specific claims about the resolution, since this is a poorly defined concept, and instead let the data speak for itself.

2) Authors should use the dm6 genome instead of dm3.

We agree with the reviewer, and have re-analyzed all sequencing data by mapping to dm6 (and lifting over to dm3 where needed to correspond with published datasets). All relevant genomic coordinates are provided primarily in dm6 (in figures). Because so many existing datasets are in dm3, and because the UCSC browser has substantially richer data for the dm3 genome, we have also provided dm3 coordinates where possible to aid others in using and expanding our work.

3) It is strange that interbands in polytene chromosomes correspond to boundaries. These structures have always been thought to contain active genes. Hug et al. identify very small domain that they also claim are boundaries. In their Hi-C data (see for example Figure 1B), once can observed interactions between these small domains that resemble compartmental interactions in mammalian Hi-C data, suggesting that these small domains do not correspond to boundaries, they correspond to active domains. These features are not observed in the Hi-C data presented in this manuscript, which may be due to the fewer reads in the Hi-C data. Authors need to consider this and try to differentiate between effects of transcription and effects of insulator proteins.

The nature of interbands has long been controversial, and it has indeed been proposed that interbands correspond to active genes. A number of pieces of evidence argued against this model, including the identification of gene bodies that reside in polytene bands, dry weight measurements that indicated that interbands represented an insufficient fraction of genomic sequence cover actively transcribed loci, and recently, fine mapping of interbands to smaller genomic regions that largely separate promoters. The Zhimulev 2014 PLoS One paper provides a fairly detailed review of much of this evidence.

We agree with the reviewer that there are additional levels of complexity in Hi-C data. We observe strong compartment effects, as well as hierarchical interactions between domains of various sizes. These are interesting phenomena for future work, and we plan on pursuing many of them. We do not attempt to address all of these potential levels of organization in this manuscript. Rather, we feel that the strong correspondence between the domain boundaries that we have identified and both the molecular features of insulators and locations of mapped polytene interbands make a compelling case that these regions represent a significant contribution to the architecture of fly chromosomes.

4) Authors cannot directly compare features of nc12 and nc14 embryos unless they have the same number of Hi-C contacts in both samples. I think it is critical to obtain more reads for the nc12 libraries.

We agree that these comparisons require deeper complexity nc12 libraries. We included nc12 in the original manuscript because we had the data and thought it would be of interest, but nc12 was not central to any of our conclusions and we have decided not to include it here.

https://doi.org/10.7554/eLife.29550.041

Article and author information

Author details

  1. Michael R Stadler

    Department of Molecular and Cell Biology, University of California, Berkeley, CA, United States
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Investigation, Visualization, Methodology, Writing—original draft, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon 0000-0002-3333-4184
  2. Jenna E Haines

    Department of Molecular and Cell Biology, University of California, Berkeley, CA, United States
    Contribution
    Software, Formal analysis, Investigation, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
  3. Michael B Eisen

    1. Department of Molecular and Cell Biology, University of California, Berkeley, CA, United States
    2. Department of Integrative Biology, University of California, Berkeley, CA, United States
    3. Howard Hughes Medical Institute, Berkeley, CA, United States
    Contribution
    Conceptualization, Software, Formal analysis, Supervision, Funding acquisition, Methodology, Writing—review and editing
    For correspondence
    mbeisen@berkeley.edu
    Competing interests
    No competing interests declared
    ORCID icon 0000-0002-7528-738X

Funding

Howard Hughes Medical Institute

  • Michael Eisen

American Cancer Society (126730-PF-14-256-01-DDC)

  • Michael R Stadler

National Science Foundation

  • Jenna E Haines

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We are especially thankful to Emily Brown for her assistance in adapting Hi-C to fly embryos, to Xiao-Yong Li for help with embryo sorting and with optimizing the fixation and chromatin isolation protocols, and to Steven Kuntz for assistance with developing embryo sectioning protocols. We thank Mustafa Mir, Xavier Darzacq and members of the Eisen and Darzacq labs for critical discussions and advice supplied throughout the work. MS was supported by an American Cancer Society postdoctoral fellowship (126730-PF-14-256-01-DDC), JH was supported by the National Science Foundation Graduate Research Fellows Program, and the work was supported by an HHMI investigator award to ME.

Reviewing Editor

  1. Allan C Spradling, Reviewing Editor, Howard Hughes Medical Institute, Carnegie Institution for Science, United States

Publication history

  1. Received: June 13, 2017
  2. Accepted: November 13, 2017
  3. Accepted Manuscript published: November 17, 2017 (version 1)
  4. Version of Record published: December 21, 2017 (version 2)

Copyright

© 2017, Stadler et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 2,189
    Page views
  • 465
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading