ChIP-seq and ATAC-seq data analysis

All short-read alignments were performed against our D. pseudoobscura genome using Bowtie 2 v2.2.7 (Langmead and Salzberg 2012). ChIP-seq peak calling was performed using MACS2 (version 2.0.10) (Zhang et al. 2008) with the default parameters. ChIP-seq normalization was performed using bamCompare from the deepTools suite (version 3.2.1) (Ramírez et al. 2016) with the following setting: “‐‐binSize 10 ‐‐operation log₂ ‐‐minMappingQuality 30 ‐‐skipNonCoveredRegions ‐‐ignoreDuplicates”. Read coverage of ATAC-seq was computed using deepTools bamCoverage for a bin size of 10 bp. To generate metaregion plots (Fig. 1D–F) of ChIP-seq/ATAC-seq signals or frequency of insulator binding sites surrounding TAD boundaries, a matrix A_ij was generated for each data set using deepTools computeMatrix and Perl scripts, in which each row represents a boundary and each column (j ∈[−40,40]) represents the signal value in a 1-kbp nonoverlapping bin within 40 kbp of the downstream and upstream flanking regions of that boundary. For CTCF and BEAF-32 binding sites, we summed values from columns of A_ij into a vector in which each element represents the signal value for the corresponding 1-kbp bin. For ChIP-seq and ATAC-seq data, we averaged values from columns of A_ij into a vector. To assess the significance of each signal at TAD boundaries, we generated 10,000 random samples of simulated TAD boundaries with the number and chromosome distribution confined by the observed data set using BEDTools shuffle (version 2.25.0) (Quinlan and Hall 2010). We then computed the sampling distribution of each signal value around TAD boundaries in the same way as described above for actual boundaries and determined the P-values.

The preprocessed ChIP data for D. melanogaster were obtained from the modENCODE Consortium (The modENCODE Consortium et al. 2010; http://www.modencode.org/).