The role of insulators and transcription in 3D chromatin organisation of flies

The DNA in many organisms, including humans, is shown to be organised in topologically associating domains (TADs). In Drosophila, several architectural proteins are enriched at TAD borders, but it is still unclear whether these proteins play a functional role in the formation and maintenance of TADs. Here, we show that depletion of BEAF-32, Cp190, Chro and Dref leads to changes in TAD organisation and chromatin loops. Their depletion predominantly affects TAD borders located in heterochromatin, while TAD borders located in euchromatin are resilient to these mutants. Furthermore, transcriptomic data has revealed hundreds of genes displaying differential expression in these mutants and showed that the majority of differentially expressed genes are located within reorganised TADs. Our work identifies a novel and functional role for architectural proteins at TAD borders in Drosophila and a link between TAD reorganisation and subsequent changes in gene expression.


Introduction
Topologically associating domains (TADs) provide a fundamental unit for chromosome organisation 1,2 and are widely conserved across species 3 as well as during different developmental stages 4,5 , suggesting that they have a functional role. Furthermore, in Drosophila cells, changes in the 3D organisation of DNA after heat stress have been found to correlate with transcriptional changes 6 . Recent evidence points to defective 3D architecture as a major contributor for diseases, developmental defects and even ageing [7][8][9][10][11][12][13] . These results suggest that 3D organisation of the DNA is important in gene regulation.
There has been significant progress in generating empirical data on chromatin organisation in different organisms and tissues, but, despite this, the mechanisms that drive the formation of TAD borders remain unclear. Previous research has shown that TAD borders are enriched in housekeeping genes 6 , developmental enhancers 14 and highly conserved genomic regulatory blocks 15 . In addition, architectural proteins and insulators are enriched at TAD borders 16,17 . Two different mechanisms were proposed to be responsible for TAD formation: (i) compartment domains, which are formed by interactions among sequences that contain active or inactive histone modifications and (ii) loop domains that are flanked by CTCF binding sites and are formed by a cohesin driven loop extrusion mechanism [18][19][20] . The latter displays a strong loop localised at the top of the TAD, while the former lacks this chromatin loop. In mammalian systems, CTCF and cohesin are the main architectural components that are located at TAD borders and their depletion has been shown to disrupt TADs [21][22][23] . By contrast, in Drosophila, several insulator proteins occupy TAD borders, such as CTCF, BEAF-32, Chro and Cp190 16,[24][25][26][27] , but the majority of TADs lack the chromatin loop at the top of the TAD suggesting a prevalence of the compartment domains 27,28 . In particular, previous research has identified strong enrichment of BEAF-32 at TAD borders in Drosophila 16,25,26,29 , but this was more pronounced in cell lines derived from the embryo (Kc167 derived from dorsal closure stage and S2 derived from late embryonic stage) or whole embryos.
Interestingly, there are negligible changes in 3D chromatin organisation following BEAF-32 RNAi knockdown in Kc167 cells 25 despite the strong enrichment of BEAF-32 at TAD borders.
Kc167 cells display saturating levels of BEAF-32 at TAD borders, suggesting that upon RNAi knockdown, there is potentially still sufficient protein present in the cell to maintain TAD borders 30 . Furthermore, BEAF-32 displays the same binding motif as another architectural protein in Drosophila called Dref 31,32 . When BEAF-32 is depleted, one possibility is that Dref replaces it at TAD borders and this could explain the lack of changes in 3D organisation observed in Kc167 cells.
Two additional proteins, Cp190 and Chro, are enriched at TAD borders 14,24,29 . These proteins cannot bind independently to DNA, but are recruited mainly by BEAF-32 33 , with up to 91% of TAD borders in a Drosophila cell line (S2) displaying presence of BEAF-32 together with either Cp190 or Chro 29 . Like BEAF-32, the role of Cp190 and Chro at TAD borders is currently unclear.
Recently, the role of TADs in gene regulation has been challenged 34,35 . In one example, it was shown that changes in TAD borders and changes in transcription are not coupled when investigating a Drosophila balancer chromosome containing chromosome re-arrangements 34 .
However, the balancer chromosomes display a very small number of rearrangements that result in changes at only a few TAD borders. It is less likely that effects on gene expression will be observed when sampling only a few re-arrangements and one possibility is that more and stronger changes in TADs (e.g., more TAD borders are lost) would allow the observation of changes in gene expression that correlate with reorganisations of TADs.
We depleted BEAF-32 in BG3 cells (derived from the larval central nervous system) using RNAi knockdown and measured the changes in 3D chromatin organisation at sub-kilobase resolution together with changes in transcription to dissect the mechanism at TAD borders and evaluate the functional role of TADs. In BG3 cells, BEAF-32 has reduced levels at TAD borders 26 , which raises the question of whether a strong depletion combined with the low levels of BEAF-32 is sufficient to affect the borders of the TADs. We also performed double knockdowns of Cp190/Chro and BEAF-32/Dref using RNAi to disentangle the interactions between different architectural proteins at TAD borders.

BEAF-32, Cp190 and Chro have a functional role at TAD borders in BG3 cells.
We performed single knockdown of BEAF-32 and combinatorial knockdown of Cp190 and Chro in BG3 cells followed by in situ Hi-C ( Figure S1 and Tables S1 and S2). The RNAi knockdowns lead to specific and strong reduction in both the mRNA levels and protein levels and do not affect the cell cycle ( Figure S1). High-resolution contact maps were generated for both knockdown mutants. The biological replicates displayed high similarities and were merged for the downstream analysis ( Figure S2). There was a noticeable re-organisation in the contact maps as a result of the knockdowns when compared to wild type BG3 cells ( Figure S2 and Figure 1A). BG3 BEAF-32 resulted in loss of long-range interactions and showed an increase in short-range interactions ( Figure 1A). Likewise, BG3 Cp190 -Chro also exhibited reduced long-range interactions and increased short-range interactions, but the loss of longrange interactions were less pronounced compared to BG3 BEAF-32 -( Figure 1A).
Several papers have proposed that BEAF-32, Cp190 and Chro control the borders of TADs 16,25,29 . We then investigated the TAD border classification as performed previously 26 . In particular, we used HiCExplorer 25 and identified between 2000 and 2600 TADs (see Table   S2 and Methods), which is consistent with other studies 14,25,26 . TAD borders were classified into weak and strong borders depending on whether they can be detected with increasing stringency of the TAD calling algorithm, with strong borders being detected even with the more stringent parameters (see Methods). To investigate the robustness of these TAD borders, we downsampled all Hi-C libraries by 20% and repeated the analysis (see Figure   S3). 706 of the 989 strong TAD borders in WT cells are robust, meaning they are recovered in both full and downsampled datasets ( Figure 1B).  Figure 1C).
Next, in order to distinguish between direct and indirect effects, we evaluated how many of the maintained and lost robust borders have BEAF-32, Chro or Cp190 ChIP peaks in their vicinity. Figure 1E shows that majority of maintained TAD borders (94%) are direct targets of the three proteins, but only half of the lost TAD borders (47%) are direct targets (also see Figure S4A). Furthermore, the majority of maintained TAD borders (70%) retain BEAF-32 or Cp190 upon knockdown, but most of the lost borders (70%) lose binding of these architectural proteins after knockdown ( Figure S4B-C). This further confirms that the direct maintained and direct lost TAD borders are indeed controlled by the three architectural proteins.
Some regions displayed high conservation of the TAD structure organisation ( Figure 1F), while others showed reorganisation ( Figure 1G-H). We observed that a loss of a TAD border could result in either movement of the TAD borders (see top panel in Figure 1G) or aggregation of two TADs (see bottom panel in Figure 1G).
We also found new border formation in both knockdowns, ranging between 400 to 600 weak borders and 200 to 300 strong borders. The majority of these new borders moved more than 2Kb in the mutants compared to WT ( Figure 1G-H and Figure S5). A small proportion of the new TAD borders result in splitting the original TAD in two separate TADs (see bottom panel of Figure 1H and Figure S5). Out of all the new borders, only 43 were common between both knockdowns ( Figure 1D). This may be explained by the fact that Chro and Cp190 are able to bind chromatin independent of BEAF-32 36 . Interestingly, the majority of these new borders have BEAF-32, Chro or Cp190 ChIP peaks in their vicinity and retain BEAF-32 or Cp190 upon knockdown ( Figure 1E and Figure S4A-C). To identify the roles of BEAF-32, Cp190 and Chro at TAD borders, we focused on two groups: (i) maintained borders (robust TAD borders that are strong in WT cells and are maintained strong in both mutants) and (ii) lost borders (robust TAD borders that are strong in WT cells and are lost in the two mutants).

Combined Dref and BEAF knockdown shows an enhanced effect on TAD border distribution
Dref is a DNA binding protein that shares a similar binding motif with BEAF-32, meaning that upon depletion of BEAF-32, Dref could potentially replace it at TAD borders. To investigate this, we performed a combinatorial knockdown of BEAF-32 and Dref ( Figure S1) followed by in situ Hi-C (Tables S1 and S2). Again, the combinatorial knockdown resulted in specific and efficient depletion at both mRNA and protein levels and does not affect the cell cycle ( Figure   S1). In the BEAF-32 Dref double knockdown (BG3 BEAF-32 -Dref -) we noticed a more pronounced effect in the reorganisation of the 3D interaction compared to the BG3 BEAF-32 or BG3 Cp190  Figure 2C and D). The majority of these new borders are movements of borders in the mutant compared to the closest WT border ( Figure S5 and Figure 2G). Overall, we found that there is a large overlap between TAD borders that are lost in the three mutants and also a large subset of TAD borders that disappear only in the BG3 BEAF-32 -Dref mutant, indicating that there is a subset of TAD borders that require Dref for maintenance ( Figure 2C).
To distinguish the direct targets from indirect, we aligned TAD borders with the protein occupancy (see methods and Figure S4D).  Figure S4G-H). This suggests that Dref displays redundancy to BEAF-32, by maintaining TAD borders when BEAF-32 is absent. When both architectural proteins are depleted then these TAD borders that were maintained after BEAF-32 single knockdown are also lost.

Reorganisation in TADs correlates with changes in gene expression
Several studies have shown that TAD reorganisation leads to changes in transcription that correspond to developmental defects or diseases [7][8][9][10][11][12] Figure 3 and Table S3). The majority of DEGs are upregulated in the mutants compared to WT. Interestingly, almost all of these are found inside robust TADs in both WT and mutants. Figure 3 shows  Figure 3) with no specific localisation near or away from TAD borders (see Figure S7).  Figure 4B and Figure S8B).
Noticeably, there is also strong divergent transcription at the maintained borders ( Figure 4B).
The lack of enrichment for Top2 at TAD borders that are maintained in the two mutants ( Figure 4B) indicates a potential role for supercoiling at these borders.
The RNA-seq signal around maintained and lost TAD borders, does not show noticeable changes in the two mutants ( Figure 4B and S8B). At maintained borders, given that there are negligible changes in gene expression, these results were expected (see Figure 3).
Nevertheless, given the large number of differentially expressed genes associated to reorganised TADs (see Figure 3), one could expect a change in the RNA-seq signal at lost TAD borders in the mutant. Since the differentially expressed genes are randomly located inside the TAD ( Figure S7), loss of TAD borders will often correlate with changes in gene expression at a larger distance from the TAD border and, this, cannot be captured in the analysis in the vicinity of TAD borders ( Figures 4B and S8B).
By contrast, at lost TAD borders (in BG3 BEAF-32 and BG3 Cp190 -Chro -), there is less DNA accessibility and transcription indicating that these borders are in a repressed chromatin state ( Figure 4B).

Maintained borders are associated with active promoters and enhancers whereas lost borders are located in heterochromatin.
Regulatory regions in the DNA can be defined by the presence of specific histone marks 37 .
Transcription has also been shown to be strongly implicated in the maintenance and formation of TADs 6,38,39 . The presence of Pol II and nascent transcription at maintained borders and their absence from lost borders indicate the existence of two classes of TAD borders in Drosophila, active and repressed borders, which display different mechanisms of maintenance. A similar classification into active and repressed domains in Drosophila has been previously proposed 40,41 . We investigated the presence of histone modifications to further dissect the potential factors and mechanisms that would be responsible for the maintenance of the TAD borders. We found that H3K4me3 (active promoter mark) and H3K4me1 (enhancer mark) together with H3K27ac were enriched at all maintained borders ( Figure 4C). Interestingly, depletion of BEAF-32 from promoters and enhancers is not sufficient to result in the loss of these TAD borders, which indicates the presence of a redundant mechanism with a different protein(s).
We observed strong enrichment of MOF (involved in maintenance of H4K16ac), JHDM1 (H3K36me3 demethylase), ISWI and NURF301 (nucleosome sliding) and WDS (involved in maintenance of H3K4me3) preferentially at maintained borders ( Figure S9). NURF301 was shown to co-localise together with Dref and Cp190 42 , which explains its enhanced level at the maintained TAD borders.
The lost borders were strongly enriched in H3K9me2, H3K9me3 and H3K27me2 (signatures for heterochromatin and Polycomb) suggesting a plausible association of these proteins with heterochromatin regions (Figures 4C and S8C). As we observed association of lost borders with heterochromatin and Polycomb, we further dissected and analysed the Polycomb complexes in detail at all borders. However, we did not observe enrichment of any Polycomb subcomplexes (Pc or dRING) at lost borders in the two mutants ( Figure S9). Nevertheless, we did find enrichment of Su(var)3-9 and HP2, which explains the enrichment of heterochromatin at lost TAD borders in the two mutants ( Figure S9). Note that, in Drosophila, Su(var)3-9 was previously reported to have a role in maintenance of TADs located in heterochromatin 43 .
While we observed heterochromatic signatures at the lost borders ( Figure 4C Figure S10D). This means that while the majority of borders are enriched in enhancers or active TSSs, maintained borders are located in euchromatin and lost borders in heterochromatin.
One possibility is that lost borders, while euchromatic, display higher levels of Pol-II pausing.
Using the Pol II pausing index definition from 25 (see Methods), we found only negligible differences in Pol II pausing for genes located within 5K windows around of maintained, lost and new borders ( Figure S10E). This indicates that Pol II pausing does not differentially affect maintained or lost borders.
A large proportion of maintained TAD borders in the knockdowns are also present in Kc167 cells and harbour housekeeping genes.
Previously, we showed that Kc167 cells display more short-range interactions and fewer long-range contacts when compared to BG3 cells, which was true also after down-sampling to control for library size differences 26 . Given that the three mutants we analysed here display  (Table S4 and Materials and Methods).  Figure 5B). We classified 140 loops that are maintained in both BG3 BEAF-32 and BG3 Cp190 -Chro mutants as maintained loops and 122 that are lost in both BG3 BEAF-32 and BG3 Cp190 -Chro mutants as lost loops ( Figure 5C). Figure 5D confirms that the strong level of interactions is maintained in the two mutants at maintained chromatin loops, but this is not the case at lost loops. We also found that there is no statistically significant difference in the size of the lost and maintained chromatin loops ( Figure 5E).

Majority of chromatin loops in
76% to 68% of these loops connect genes to each other or other genes ( Figure 5F), indicating that they are involved in the formation of gene domains 28 . Only 9% of the maintained and lost loops are promoter-enhancer loops ( Figure 5F), which indicates that this mechanism is less prevalent in Drosophila than previously proposed in mammalian systems 49 . When we select all genes that have their promoter located at one of the anchors of the loops, we found that only a small subset of genes (<10%) located at the lost or maintained loops display differential expression in the two mutants and this is true even when using a less stringent threshold to call differentially expressed genes (log 2 Figure 5H). Conversely, a lost chromatin loop can lead to no changes in gene expression of the target genes (bottom panel in Figure 5H). Thus, our results support a model where the presence or absence of a chromatin loop does not necessarily lead to regulation of the target gene.
Chro and Cp190 are known to be involved in long-range interactions 33 , but previous research identified the enrichment of Polycomb at Drosophila loops 48 . We found that both maintained and lost chromatin loops display high levels of BEAF-32 together with Chro and/or Cp190 at both anchors ( Figure 5I and S13A); i.e., 92% of maintained and 84% of lost loops have binding of BEAF-32, Cp190 and/or Chro ( Figure S13B). The maintained loops display higher levels of Chro at the anchors compared to lost loops, suggesting that the depletion of Chro is not sufficient to affect the maintained loops. In addition, 60% of lost loops lose binding of BEAF-32 and/or Cp190 upon their knockdown ( Figure S13C-D), thus, providing support that these loops are lost as a direct consequence of the depletion of the architectural proteins in our mutants.
Nevertheless, approximately half of the maintained loops lose BEAF-32 and/or Cp190 upon knockdown ( Figure S13C-D), which suggests that Chro is recruited by additional proteins at maintained loops or that other factors could help maintain these loops ( Figure S13A). We observed an enrichment of MED1 at the anchors of maintained and lost loops, but also 1 1 enrichment of CTCF and cohesin subunit Rad21. The majority of chromatin loops in our dataset are located near a MED1 ChIP peak ( Figure S13A), indicating that Mediator complex would be more important for chromatin loops in Drosophila. We also observed a small number of loops with enrichment of Polycomb peaks near their anchors ( Figure 5 and S13), but this is less pronounced than in the case of Mediator complex.

BEAF-32, Cp190, Chro and Dref knockdown does not affect A/B compartments
The checkerboard pattern seen on Hi-C maps led to the identification of A and B compartments which mark active and inactive regions of chromatin 50 . A/B compartments were also identified in Drosophila 39 using 10 Kb bins, and we showed that compartmentalisation changes between cell lines 26 Figure S14A). Nevertheless, we identified some switching between the A and B compartments (less than 5%) ( Figure S14B). When we zoomed in, we observed that the majority of these compartments are robust and consistent in the WT and mutants ( Figure   S14C). One interesting observation is that there is some rare local spreading of the B compartment (heterochromatin) into the A compartment (euchromatin) (e.g., yellow stripe in Figure S14C).
Saddle plots confirm that regions belonging to the same compartments (lower right corner A- Altogether, our results indicate that organisation of compartments in Drosophila is independent of the organisation of TADs. When investigating in which compartments are TAD borders localised, we found that most of the maintained borders are localised in the A compartment, while most of the lost or new borders are localised in the B compartment ( Figure S15A). This is not surprising since most of the lost borders are located in repressed chromatin, while the maintained ones are in active chromatin.
The majority of compartments that switch do not harbour any DEGs, even when using a lower threshold to call differential gene expression (log 2 FC threshold of 1) ( Figure S15B).

2
Furthermore, the fact that a compartment contains DEGs does not mean that all genes in that compartment change expression in the same direction (either upregulated or downregulated). For example, spreading of B compartment in Figure S14C corresponds to three genes displaying different behaviours: ine gene is downregulated, Dp is upregulated and FIG4 maintains expression in all three mutants (all three genes are located within the yellow stripe in Figure S14C). The relationship between changes in gene expression and compartment switching is complex and often compartment switching cannot be explained by a majority of genes changing expression in the same direction. Note that RNA-seq libraries capture only polyA transcripts and do not include other transcripts such as eRNAs or lncRNAs which could potentially contribute to compartment switching.

Discussion
The enrichment of architectural proteins at TAD borders raises the question of whether they have a functional role in TAD organisation or whether their co-localisation with borders is correlative in nature. In mammalian systems, depletion of CTCF or Cohesin disrupts TADs [21][22][23] . In Drosophila, several architecture proteins (including BEAF-32, Chro and Cp190) are enriched at TAD borders, but their functional role at TAD borders has not previously been investigated 16,[24][25][26][27]45 . Our results confirm that the architectural proteins are essential for TAD borders and their depletion results in reorganisation of TADs. In particular, we found that

3
To investigate that the effects we observe in 3D chromatin organisation are not a reflection of cell cycle arrest 46,51 , but are due to the knockdown of architectural proteins, we have performed a FACS analysis. This showed that none of our knockdowns lead to changes in the cell cycle progression ( Figure Figure S15A).
The enrichment of divergent transcription at BEAF-32 enriched TAD borders that are maintained in the two mutants, when coupled with the lack of enrichment for Top2 at these borders, possibly indicates that negative supercoiling accumulates at these TAD borders, which may be due to active transcription. This negative supercoiling is not relaxed due to lack of Top2. When negative supercoiling accumulates at these borders, positive supercoiling may accumulate inside TADs, which indicates a role for supercoiling in TAD borders 52,53 .

TAD reorganisation and transcription
We identified between approximately 600 and 800 differentially expressed genes in the three mutants and the majority of those are located within TADs that lose one or both borders or shifted the position of the borders (more than 89%). We also found that there are more statistically significant DEGs than expected by chance in reorganised TADs, however this is mainly the case when TAD borders move more than 2 Kb away from their WT position. This with a regulatory sequence, less than 10% of genes display differential expression when the contact is lost, but the same is true at maintained loops. This suggests that the presence of chromatin loops would not be essential for controlling gene transcription in the majority of cases 58,63-65 .

6
Cells were pelleted, washed in PBS and resuspended in 50% ethanol in PBS and stored until analysis at 4 0 C. On the day of the analysis cells were pelleted, washed in PBS and resuspended in FACS PI buffer (PBS, 01% triton, 100µg/ml RNase and 50µg/ml propidium iodide) at a concentration of 10 6 cells/ml. The cell cycle profile was analysed with the Guava easycyte HT flow cytometer using the Incyte software and FlowJo. For each sample 15000 cells were analysed.

In situ Hi-C protocol
Hi-C libraries were generated from 10 million cells by following the in situ Hi-C protocol as mentioned in 26  replicates were multiplexed and further sequenced at Oxford Genomics Centre and Edinburgh Genomics (Genepool) using HiSeq4000.

Hi-C analysis
Each pair of the PE reads was aligned separately to Drosophila melanogaster (dm6) genome 69,70 using BWA-mem 71 (with options -t 20 -A1 -B4 -E50 -L0). HiCExplorer was used to build and correct the contact matrices and detect TADs and enriched contacts 25 Table S1). Note that the number of reads and valid pairs used in this study are within values successfully used for previous work in Drosophila cells to detect TADs, chromatin loops and compartments (e.g., 14,26,46 ). In addition, we also showed that these libraries are robust to downsampling ( Figure S3) 26 . The matrices were corrected using the thresholds in Table S2, where values were selected from the diagnostic plots ( Figure S16). By using the corrected contact matrices, we detected TADs of at least 5 Kb width using a P-value threshold of 0.01, a minimum threshold of the difference between the TAD-separation score of 0.04, and FDR correction for multiple testing (--step 2000, --minBoundaryDistance 5000 --pvalue 0.01 --delta 0.04 --correctForMultipleTesting fdr). We selected these parameters to ensure that we recover a similar number of TADs as previously reported 26 . Finally, we called strong TAD borders using a stringent value of the threshold of the difference between the TAD separation score of 1 7 0.08. This value ensured that we retrieved the strongest half of TADs. The enriched contacts were extracted with HiCExplorer using the observed/expected ratio method.

Chromatin loops
Chromatin loops were called with the HICCUPS tool from the Juicer software suite 72 on all mutants as done previously 26

Compartments
Compartments were called as described in 26,28,50 . More specifically, we used Juicer 72 to compute the eigenvectors in 10 Kb bins for all conditions 26

Saddle Plot and Compartmentalisation Strength.
We use the procedure similar to 74 . We rank each genomic region by their eigenvector value

RNA extraction and sequencing
RNA extraction was carried out using Trizol according to manufacturer's instructions. RNA was further DNase treated and purified using RNeasy Mini kit (Qiagen) following the manufacturer's protocol. RNA was assessed qualitatively and quantitatively using Quibit and 1 8 Bioanalyzer 2100(Agilent). PolyA RNA selection, library preparation and sequencing were carried out by Novogene.
Finally, we used Picard tools (http://broadinstitute.github.io/picard/) to deduplicate reads, HTseq 81 to count reads and then DESeq2 82 to detect differential expressed genes. For DESeq2 we selected transcripts with at least 10 reads and used a p-value threshold of 0.05 and a log 2 FC threshold of 2.0 (for compartments and loops we reduced the log 2 FC threshold to 1.0). A previous work used Affymetrix GeneChip expression analysis to quantify changes in transcription upon BEAF-32 knockdown in BG3 cells and they observed negligible changes in gene expression 56 . Using RNA-seq, we found a larger number of genes displaying differential expression, but this can be explained by the increase sensitivity of RNA-seq.

Analysis of Differentially and Non-Differentially Expressed Genes
We removed all genes that were not expressed in WT or any of the mutants and then we

Pol-II pausing index
We followed the method from 46 and computed the pausing index as the ratio of the mean Pol-II ChIP signal over the promoter and over the gene body. Promoter region was selected from 200 bp downstream to 50 bp upstream of TSS and gene body from 50 bp upstream to gene end. Values of 0 and below were discarded.

Data
The full list of datasets used can be found in Supplementary Tables S6-S11.

DNase-seq:
We used pre-processed DNase-seq profiles from the modENCODE Consortium 37 .

Supplementary materials and methods
Information on the comparison of borders between WT and mutants, analysis of chromatin signals at TAD borders and clustering analysis can be found in Supplementary materials and methods.

Data access
All