Enhancer–promoter interactions and transcription are largely maintained upon acute loss of CTCF, cohesin, WAPL or YY1

It remains unclear why acute depletion of CTCF (CCCTC-binding factor) and cohesin only marginally affects expression of most genes despite substantially perturbing three-dimensional (3D) genome folding at the level of domains and structural loops. To address this conundrum, we used high-resolution Micro-C and nascent transcript profiling in mouse embryonic stem cells. We find that enhancer–promoter (E–P) interactions are largely insensitive to acute (3-h) depletion of CTCF, cohesin or WAPL. YY1 has been proposed as a structural regulator of E–P loops, but acute YY1 depletion also had minimal effects on E–P loops, transcription and 3D genome folding. Strikingly, live-cell, single-molecule imaging revealed that cohesin depletion reduced transcription factor (TF) binding to chromatin. Thus, although CTCF, cohesin, WAPL or YY1 is not required for the short-term maintenance of most E–P interactions and gene expression, our results suggest that cohesin may facilitate TFs to search for and bind their targets more efficiently.


250)
using input DNA as a control. To create heatmaps we used deepTools 2.4.1 10 . We first ran bamCoverage function (--binSize 50 --normalizeTo1x 2150570000 --extendReads 300 --ignoreDuplicates -of bigwig) and normalized read numbers to 1x sequencing depth, obtaining read coverage per 50-bp bins across the whole genome (bigWig files). For spiked-in libraries, the number of deduplicated, unique reads aligning to the yeast genome in each sample was used to compute the bamCoverage scale factor to perform spike-in normalization. We then used the bigWig files to compute read numbers across 6 kb centered on either CTCF, RAD21, or YY1 peak summits as called by MACS2 in UT cells (computeMatrix reference-point --referencePoint=TSS --upstream 3000 --downstream 3000 --missingDataAsZero --sortRegions=no). We sorted the output matrices by decreasing UT enrichment, calculated as the total number of reads within a MACS2 called ChIPseq peak. Finally, heatmaps were created with the plotHeatmap function (--averageTypeSummaryPlot=mean --colorMap='Blues' --sortRegions=no). To identify the differential peaks between UT-and IAA-treated cells, we used MAnorm2 11 , which employs a hierarchical strategy for normalization of ChIP-seq data and assesses within-group variability of ChIP-seq signals based on an empirical Bayes framework. We note that the total ChIP-seq peak numbers in MAnorm2 are combined from UT-and IAA-treated cells and may differ from the number of MACS2 calling.

Micro-C assay for mammalian cells
We briefly summarize the Micro-C experiment here. The detailed protocol and technical discussion are available in our previous study 12 . Mouse embryonic stem cells (JM8.N4) and the derivative genome-edited lines were cultured in the recommended conditions 13 . When cells grew to ~70% confluency, we resuspended them in 0.05% of trypsin, inactivated with cell culture media, and resuspended in 1% formaldehyde crosslinking media (without FBS). Cells were crosslinked for 10 min at room temperature (RT) while nutating. We then added 1 M Tris-HCl pH 7.5 (final concentration 375 mM) and incubated for 5 min at RT to quench the crosslinking reaction. Cells were spun down and washed twice with cold PBS. Cell pellets were crosslinked again with the freshly prepared 3 mM DSG crosslinking solution (in base media without FBS) for 45 min at RT while nutating. The crosslinking reaction was quenched and washed following the same steps as above. We routinely split cells into 1 million cells per vial after fixation and perform MNase titration and Micro-C with 1 million cells. Crosslinked cell pellets can be snap frozen in liquid nitrogen and stored at -80 °C for months or used immediately for the next step. We note that 1) cells directly crosslinked on the dish typically yield a similar result; 2) using freshly made formaldehyde and DSG solution is critical to obtain high-quality Micro-C data; and 3) to avoid loss of cells, low-retention tubes and tips are strongly recommended. Crosslinked cell pellets were resuspended in MB1 (50 mM NaCl, 10 mM Tris-HCl pH 7.5, 5 mM MgCl2, 1 mM CaCl2, 0.2% NP-40) at a concentration of 1x10 6 cells/100 µL and were incubated for 20 min on ice. Cells were spun down and washed once with MB1. We then added the appropriate amount of Micrococcal nuclease (MNase) and incubated the tube for 10 min at 37 °C while shaking in a thermomixer at ~850 rpm. The optimal digestion condition results in ~90% of mono-nucleosome and ~10% of di-nucleosome. The MNase reaction was inactivated by adding 4 mM EGTA and incubated for 10 min at 65 °C. Digested chromatin was spun down and washed twice with 1 mL of cold MB2 (50 mM NaCl, 10 mM Tris-HCl pH 7.5, 10 mM MgCl2). We note that and MNase titration that yields 90% monomer/10% dimers substantially reduces contamination with un-digested (un-ligated) dimers in Micro-C data. MNased-fragmented chromatin was then subjected to the three-step end-repair protocol to generate ligatable ends filled with biotin-dNTPs. First, the pellet was resuspended in the end-repair buffer (50 mM NaCl, 10 mM Tris-HCl pH 7.5, 10 mM MgCl2, 100 µg/mL BSA, 2 mM ATP, 5 mM DTT) and the 5'-ends were phosphorylated with 25 units of T4 Polynucleotide Kinase (NEB #M0201) for 15 min at 37 °C while shaking in a thermomixer at 1000 rpm for an interval of 15 sec every 3 min. Second, to convert the mixed types of nucleosomal ends to cohesive ends, we added 25 units of DNA Polymerase I, Large (Klenow) Fragment (NEB #M0210) directly to the reaction and incubated the tube for an additional 15 min at 37 °C while shaking in a thermomixer at 1000 rpm for an interval of 15 sec every 3 min. Third, to repair the nucleosomal DNA ends to the blunt and ligatable ends, we supplemented 66 µM of dNTPs (dTTP, dGTP (NEB #N0446), biotin-dATP (Jena Bioscience #NU-835-BIO14), biotin-dCTP (Jena Bioscience #NU-809-BIOX), and 1X T4 DNA ligase reaction buffer (NEB #B0202) directly into the reaction and incubated for 45 min at 25 °C while shaking in a thermomixer at 1000 rpm for an interval of 15 sec every 3 min. The end-repair reaction was then inactivated with 30 mM EDTA for 20 min at 65°C without shaking. Next, chromatin was pelleted by centrifugation for 5 min at ~10,000xg at 4 °C and washed once with cold MB3 (50 mM Tris-HCl pH 7.5, 10 mM MgCl2). End-repaired nucleosomes were then subjected to proximity ligation with ~5000 cohesive end units (CEU) of T4 DNA ligase (NEB #M0202) in 1X T4 DNA ligase reaction buffer (NEB #B0202) for at least 2 hrs at room temperature with slow rotation on an orbital shaker. Biotin-dNTPs at the unligated DNA termini were removed by ~500 units of Exonuclease III (NEB #M0206) in 1X NEBuffer 1 (NEB #B7001) for 15 min at 37°C while shaking in a thermomixer at 1000 rpm for an interval of 15 sec every 3 min. Samples were then reverse crosslinked with 1X proteinase K solution (500 ug/uL Proteinase K (ThermoFisher #AM2542), 1% SDS, 0.1 µg/µL RNaseA) at 65°C overnight. DNA was extracted by the standard phenol:chloroform:isoamyl alcohol (25:24:1) and ethanol precipitation procedure. DNA was then purified again with the ZymoClean DNA Clean & Concentrator-5 Kit (Zymo #D4013). Purified DNA was separated on a 3% NuSieve GTG agarose gel (Lonza #50081). The gel band corresponding to the size of dinucleosomal DNA (~250 to 350 bp) was cut and purified with the Zymoclean Gel DNA Recovery Kit (Zymo #D4008). We note that size selection for DNA larger than 200 bp greatly reduces the ratio of unligated monomers in Micro-C data. Micro-C sequencing libraries were generated by using the NEBNext Ultra II DNA Library Prep Kit for Illumina (NEB #E7645) with the following modifications. We first repaired the purified DNA again using the End-it DNA End-Repair Kit (Lucigen #ER0720) following the manufacturer's suggested conditions. The mix was incubated for 45 min at 25 °C and then inactivated the enzyme reaction for 10 min at 70 °C. This step is optional, but we find it increases the library yield. Biotinylated DNA was captured by Dynabeads MyOne Streptavidin C1 beads (ThermoFisher #65001) in 1X BW buffer (5 mM Tris-HCl pH 7.5, 0.5 mM EDTA, 1 M NaCl) on a nutator for 20 min at room temperature. Beads were washed twice with 1X TBW buffer (0.1% Tween20, 5 mM Tris-HCl pH 7.5, 0.5 mM EDTA, 1 M NaCl) for 5 min at 55 °C while shaking in a thermomixer at 1200 rpm, rinsed once with Tris buffer (10 mM Tris-HCl pH 7.5), and then resuspended in Tris buffer. We then performed 'on-bead' end-repair/A tailing and adapter ligation following the NEB protocol. After adapter ligation, beads were washed once with 1X TBW and rinsed once with Tris buffer. The Micro-C library was generated by using the KAPA HiFi HotStart ReadyMix (Roche #KK2601) or the Q5 High-Fidelity 2X Master Mix (NEB #E7645) with the manufacturer's suggested conditions. We recommend using a minimal PCR cycle to reduce PCR duplicates, typically 8 to 12 cycles, which can generate high-quality Micro-C data. Prior to sequencing, purifying the library twice with 0.85X AMPure XP beads (Beckman #A63880) can eliminate primer dimers and adapters. We used Illumina 100 bp paired-end sequencing (PE100) to obtain ~400M reads for each replicate in this study.

Micro-C data processing and analyses
Valid Micro-C contact read pairs were obtained from the HiC-Pro analysis pipeline (v2.11.3) 14 , and the detailed description and code can be found at https://github.com/nservant/HiC-Pro. In brief, paired fastq files were mapped to the mouse mm10 genome independently using Bowtie 2.3.0 with 'very-sensitive-local' mode 8 . Aligned reads were paired by the read names. Pairs with multiple hits, low MAPQ, singleton, dangling end, self-circle, and PCR duplicates were discarded using Samtools (v.1.9). Paired reads with distances shorter than 100 bp (e.g., unligated mononucleosome) were also removed. Output files containing all valid pairs were used in downstream analyses. We recommend running a pilot sequencing run (~10M reads) and checking the following quality control statistics before moving forward to a high-coverage sequencing: (1) bowtie mapping rate; (2) reads pairing percentage; (3) ratio of sequencing artifacts; (4) ratio of cis/trans contacts; (5) unligated monomer percentage. If any of the above statistics is not optimal, one might consider checking mapping and filtering parameters or further optimizing the Micro-C experiment. The summary of Micro-C interactions in this manuscript is available in Supplemental Table 2.
We evaluated the reproducibility and data quality for the Micro-C replicates using two published methods independently (https://github.com/kundajelab/3DChromatin_ReplicateQC) 15 . In brief, QuASAR-QC (v1.5.7) calculates the correlation of values in two distance-based transformed matrices. GenomeDISCO (v1.0.0) measures the difference in two graph diffusion-smoothed contact maps. We computed the matrix similarity scores between the biological replicates or between the untreated and IAA-treated cells for the 10-kb, 25-kb, 50-kb, and 250-kb Micro-C data.  20 . Briefly, we analyzed the insulation profile by using a 1-Mb sliding window that scans across Micro-C contact matrices at 20-kb resolution and assigns an insulation intensity to its corresponding bin. The insulation scores were obtained and normalized as the log2 ratio of the individual score to the mean of the genome-wide averaged insulation score. Chromatin boundaries can be identified by finding the local minima along with the normalized insulation score. Boundaries overlapping with low mappability regions were removed from the downstream analysis. The arrowhead analysis defines Ai,i+d = (M*i,i-d-M*i,i+d)/(M*i,i-d+ M*i,i+d), where M* is the normalized contact matrix. Ai,i+d can be thought of as the measurement of the directionality preference of locus i, restricted to contacts at a linear distance of d. Ai,i+d will be strongly positive/negative if either i,i-d or i,i+d is inside the domain and the other is not, but Ai,i+d will be close to zero if both loci are inside or outside the domain. By assigning this query across the genome, the edges of a domain will be sharpened and TADs can be detected. For aggregate domain analysis (ADA), each domain was rescaled to a pseudo-size by Ni,j=((Ci-Dstart)/(Dend-Dstart), (Cj-Dstart)/(Dend-Dstart)), where Ci,j is a pair of contact loci within domain D that is flanked by Dstart and Dend, and Ni,j is a pair of the rescaled coordinates. The rescaled domains can be aggregated at the center of the plot with ICE or distance normalization. Coolpup (v0.9.5) (https://github.com/open2c/coolpuppy) 21 has implemented a handy function to perform ADA with the COOL file. The lists of TAD called by the insulation score analysis or Arrowhead are available in Supplemental Table 3

or 4.
To identify loops/dots, we tested two novel algorithms, Mustache (v1.0.1) (https://github.com/aylab/mustache) 22 and Chromosight (v0.9.8) (https://github.com/koszullab/chromosight) 23 , for the high-resolution Micro-C data. We found that both approaches outperform the HICCUPS algorithm in the JUICER package 18 and the 'Call-dots' function in the Cooltools package in sensitivity and specificity to discover focal contact enrichment. In Mustache analysis, we called loops with balanced contact matrices at resolutions of 400 bp, 600 bp, 800 bp, 1 kb, 2 kb, 4 kb, 10 kb, and 20 kb using the calling options --pThreshold 0.1 --sparsityThreshold 0.88 -octaves 2. We then combined all loops at different resolutions. If an interaction was detected as a loop at different resolutions, we retained the precise coordinates in finer resolutions and discarded the coarser resolution. In Chromosight analysis, we used the 'detect' function to call loops with balanced contact matrices at resolutions of 400 bp, 600 bp, 800 bp, 1 kb, 2 kb, 4 kb, 10 kb, and 20 kb using calling options listed in Supplemental Table 5. We then combined all loops at different resolutions by the same approach as described above. We applied the lists of loop anchor to many downstream analyses by using Bedtools (v2.30.0) 24 , R (v4.0), Python (v3.8), or MATLAB (v2021a), including (1) comparison of loop anchors between Micro-C and Hi-C or between different chromatin states; (2) distribution of loop strength or length; (3) cross-correlation with ChIP-seq, RNA-seq, and mNET-seq data; (4) ratio of boundary crossing, etc. For analysis of paired genomic loci (e.g., paired ChIP-seq peaks, genetic features, etc.) within a distance ranging from 2 kb to 2 Mb, we used Chromosight's 'quantify' function to measure the probability of loop pattern for all intersections quantitatively. The loops were filtered by the following parameters: loop score >0.35 for 10-kb resolution, loop score >0.3 for 4-kb resolution, loop score >0.2 for 2kb resolution, and the q-value lower than 10 -5 for all resolutions (Supplemental Table 5). For aggregate peak analysis (APA) to assess genome-wide loop intensity, loops were centered and piled up on a 20-kb x 20-kb matrix with 400-bp resolution balanced data or 50-kb x 50-kb matrix with 1-kb resolution balanced data. Contacts close to the diagonal were excluded and normalized by a random shift matrix to avoid distance decay effects. The ratio of loop enrichment was calculated by dividing normalized center contacts in a searching window by the normalized corner submatrices. We used the same approach and normalization method to analyze the genome-wide target-centered loop intensity. Instead of aggregating at the intersection of loop anchors, the matrix is centered at the paired ChIP-seq peaks or genomic features. Coolpup (v0.9.5) (https://github.com/open2c/coolpuppy) 21 has implemented the APA function for the COOL file. The lists of loops called by Mustache or Chromosight are available in Supplemental Table 6 or 7. The lists of loop quantification for cohesin, E-P, and P-P loops are available in Supplemental Table 8 -10.

Single-particle imaging experiments
All single-molecule imaging experiments were performed with a similar setting as described in our previous studies 3, 26 and detailed below. In brief, the experiments were performed with a custom-built Nikon TI microscope equipped with a 100X/NA 1.49 oil-immersion TIRF objective (Nikon apochromat CFI Apo TIRF 100X Oil), an EMCCD camera (Andor iXon Ultra897), a perfect focus system to correct for axial drift and motorized laser illumination (Ti-TIRF, Nikon), and an incubation chamber maintaining a humidified 37 °C atmosphere with 5% CO2 for the sample and the objective. Excitation was achieved using the following laser lines: 561 nm (1 W, Genesis Coherent) for JF549/PA-JF549 and TMR dyes; 633 nm (1 W, Genesis Coherent) for JF646/PA-JF646 dyes; 405 nm (140 mW, OBIS, Coherent) for all photo-activation experiments. Laser intensities were controlled by an acousto-optic Tunable Filter (AA Opto-Electronic, AOTFnC-VIS-TN) and triggered with the camera TTL exposure output signal. Lasers were directed to the microscope by an optical fiber, reflected using a multi-band dichroic (405 nm/488 nm/561 nm/633 nm quad-band, Semrock) and focused in the back focal plane of the objective. Emission light was filtered using single band-pass filters placed in front of the camera (Semrock 593/40 nm for TMR and JF549/PA-JF549 and Semrock 676/37 nm for JF646/PA-JF646). The angle of incident laser was adjusted for highly inclined laminated optical sheet (HiLo) conditions 27 . The microscope, cameras, and hardware were controlled through the NIS-Elements software (Nikon). For 'fast-tracking' stroboscopic illumination (spaSPT) at ~133 Hz, the excitation laser (633 nm for PA-JF646 or 561 nm for PA-JF549) was pulsed for 1 ms at maximum (1 W) power at the beginning of the frame interval, while the photoactivation laser (405 nm) was pulsed during the ∼447 µs camera transition time. Each frame consisted of a 7-ms camera exposure time followed by a ~447-μs camera inactive time. The camera was set for frame transfer mode and vertical shift speed at 0.9 μs. With this setup, the pixel size after magnification is 160 nm and the photon-tograyscale gain is 109. Typically, 30000 frames with this sequence were collected per nucleus, during which the 405-nm intensity was manually tuned to maintain an average molecule density of ~0.5 localizations per frame, corresponding to ~15,000 localizations per cell per movie. Maintaining a very low density of molecules is necessary to avoid tracking errors. For 'slow-tracking' (slow-SPT) experiments, long exposure times (50 ms, 100 ms, and 250 ms) and low constant illumination laser intensities (0.5% -2% of 0.5 W power) were used to measure residence time. The camera was set for normal mode and vertical shift speed at 3.3 μs. We generally recorded each cell with 1200 frames for a 250 ms exposure time, 3000 frames for a 100 ms exposure time, or 6000 frames for a 50 ms exposure time. We included H2B-HaloTag cells for the photobleaching correction for each experiment.

spaSPT analysis
To estimate the likelihood of diffusion coefficients for each trajectory we used the SASPT package (v1.0) (https://github.com/alecheckert/saspt) 28 . The detailed discussion is available in Heckert et al 28 . In brief, we applied "State Array (SA)", a grid of state parameters that span a range of diffusion coefficients (0 to 100 μm 2 /s), to calculate the posterior occupations of each point in the grid. The SA method conceptually produces a similar result as the Dirichlet process mixture models (DPMM) and retains its ability to model complex diffusive mixtures, while mitigating the issue of expensive likelihood functions. Instead of allowing infinite number of states (K→∞), the method fixes the number of states at a large but finite value. For each state j = 1, …, K, the algorithm chooses a fixed set of state parameters θj. The model simplifies to: As a result, the SA method can compute more complex likelihood functions that incorporate localization error. Next, the algorithm infers the posterior p (Z, τ | X) with the Dirichlet prior p(τ) and corrects the systematic overestimation for the occupations of slow states due to defocalization. By marginalizing the posterior distribution on localization error, the method naturally incorporates uncertainty about the localization error of different states. Alternatively, we analyzed the spaSPT data with the kinetic modeling framework implemented in the Spot-On package (v.1.0.4) 26 . Briefly, the model infers the diffusion constant and relative fractions of two or three subpopulations from the distribution of displacements computed at increasing lag times (1∆τ, 2∆τ,…). This is performed by fitting a semi-analytical model to the empirical histogram of displacements using non-linear least squares fitting. Defocalization is explicitly accounted for by modeling the fraction of particles that remain in focus over time as a function of their diffusion constant. We used the following setting for Spot-On analysis in the manuscript:

Slow-SPT analysis
To extract residence times from slow-SPT data, we used long exposure times (50 ms, 100 ms, or 250 ms) to motion-blur freely diffusing molecules into the background 3,29-31 . We then recorded the trajectory length of each 'bound' molecule and used these to generate a survival curve (1-CDF), and performed double-exponential fitting to estimate the unbinding rates for non-specific binding (Kns) and specific binding (Ks). We note that localization errors can cause both falsepositive and false-negative detections. The Kns is likely to be contaminated by localization errors (e.g., from molecules close to being out-of-focus) and experimental noise. To filter out contributions from tracking errors and slow-diffusing molecules, we applied an objective threshold as previously described to consider only particles tracked for at least Nmin frames. To determine Nmin, we plotted the inferred residence time as a function of Nmin and observed convergence to a single value after ~2.5 s. We thus used this threshold to determine the value of Ks. To correct the biases from photobleaching, cell drifting, and background fluctuating, we assume that all these factors should affect H2B-HaloTag to the same extent as those affecting YY1-HaloTag. We can use an apparent unbinding rate for H2B-HaloTag as Kbias, consistent with our FRAP analysis. Thus, we performed the slow-SPT experiments for YY1-HaloTag and H2B-HaloTag with the same camera and laser settings on the same day. We then obtained the residence time as:

FRAP image analysis
To quantify FRAP movies, we wrote a pipeline in MATLAB (v2021a). Briefly, our algorithm automatically detects the bleached spot, the background spot, and the nucleus segments by Gaussian smoothing, hole-filling, and segmenting a nucleus in a FRAP movie. Cell drift is also automatically corrected by the optimal linear translation in x and y. Next, we quantify the bleach spot signal as the mean intensity of a slightly smaller region, which is more robust to lateral drift. The FRAP signal is corrected for photobleaching using the measured reduction in total nuclear fluorescence and internally normalized to its mean value during the 20 frames before bleaching. .

Inferring parameters related to YY1's target search mechanism
We used the parameters inferred from our spaSPT and the residence time measurements from our FRAP or slow-SPT analysis. The detailed discussion is available in Hansen et al 3 . Briefly, from the Spagl State Array analysis, we determined that the total bound fraction for YY1 is ~31%. However, the total bound fraction (0 -0.1 μm 2 /s) contains both YY1 molecules bound specifically to their target sites and non-specific interactions (e.g., sliding or jumping on DNA). We previously estimated the fraction that is non-specifically bound using a mutant CTCF with a His-to-Arg mutation in each of the 11 zinc-finger domains 3 . This CTCF mutant is virtually unable to interact specifically with its binding sites. The Spagl analysis estimated the bound fraction to be ~8.1% for this mutant in mESCs. Since we did not perform the spaSPT experiments for YY1's DNA binding domain mutants, we thus inferred the FBOUND,specifc ~= FBOUND,total in this manuscript. We next determined the average time for YY1 to find its cognate site after dissociating from the previous site. We will use 's' and 'ns', as abbreviations for specific and non-specific binding, respectively, in the following discussion. The pseudo-first-order rate constant for specific binding sites, k * ON,s, is related to the fraction bound by:

I. Technical controls of Micro-C assay
To address if the fine-scale chromatin loops/dots detected in Micro-C are due to higher accessibility of micrococcal nuclease (MNase) to promoters and enhancers, we here cite other studies and present results of additional analyses on our Micro-C data that demonstrate accessibility biases are largely negligible.

Micro-C data normalization
The 3C-based assays involve a series of biochemical reactions that may introduce noise to the data analysis and the final results. Such noise must be eliminated before interpreting the data. Early normalization methods for Hi-C data focused on explicit factors causing noise. Factors like sequence uniqueness, GC content, uneven distribution of restriction sites across the genome may introduce biases into the 3C-based analysis 1 . As a solution to these issues, the Mirny and Aiden labs proposed using matrix balancing methods to handle all noise sources "implicitly" 2,3 . Two fundamental assumptions underlie the methods: 1) visibility across all genomic regions should be equal, and 2) all Hi-C biases are one-dimensional and factorizable. A balanced matrix should allow equal visibility to any genomic locus regardless of system biases, that is, a normalized matrix whose rows and columns sum to the same amount. In Micro-C/Hi-C analysis, three primary matrix balancing algorithms are commonly used, including "the square root of vanilla coverage (VC)", "Knight-Ruiz (KR)", and "Iterative Correction and Eigenvector decomposition (ICE)". The idea of VC is similar to coverage normalization as it divides each matrix element by its row sum and column sum to eliminate different sequencing coverage of each locus. ICE and KR repeat the VC process until all rows and cols sum to the same value and converge. Results from KR and ICE are typically nearly identical and highly reproducible. Additionally, we often performed distance normalization (e.g., observed/expected) to remove the bias caused by 1D distance-dependent interactions in 3C-based data.

MNase digests chromatin uniformly
Recent studies from the Henikoff and Längst labs have experimentally and computationally addressed this issue 4,5 , indicating that MNase equally accesses DNaseI-hypersensitive and nonhypersensitive sites, sites of active and inactive transcription, and euchromatin as well as heterochromatin (see Fig.9 in Chereji et al, Fig.4 in Schwartz et al, and a snapshot from Chereji et al in Supplementary Fig. 1 below). They also state that "although heterochromatin and euchromatin appear different when observed cytologically at low resolution, at the molecular level, MNase and other proteins can access heterochromatin regions at rates similar to those of accessing euchromatin". We thus conclude that MNase can digest chromatin uniformly across the genome regardless of its accessibility level.

Micro-C is not biased toward chromatin accessible regions
Previously, we and others have independently demonstrated that Micro-C contacts are not biased toward highly accessible chromatin regions [6][7][8] . Notably, a recent Micro-C-based technique, Micro-Capture-C (MCC), that focuses on the interactions around promoter regions also did not find substantial bias toward DNaseI hypersensitive sites 8 . To further support this point, we provide key results below demonstrating that Micro-C data is not affected by chromatin accessibility levels.
We thus conclude that the E-P interactions captured by Micro-C represent bona fide proximal chromatin contacts.
• When we analyzed Micro-C reads as single-end data like ChIP-seq or MNase-seq, we do not observe any apparent correlation with chromatin accessibility. In Supplementary Fig.  2, panel i, we show that neither raw (r 2 = 0.12) nor balanced (r 2 << 0.1) Micro-C data show substantial correlation with ATAC-seq signal. While the raw Micro-C data shows an anticorrelation to ATAC-seq around ATAC-seq peaks, this trend disappears after matrix balance (Supplementary Fig. 2, panel ii). We note that all the Micro-C data in this manuscript were normalized either by matrix balancing (detailed in the Methods section). • If contacts at accessible promoters/enhancers were artifactually amplified by Micro-C, we should expect them to be far more prevalent than contacts at polycomb/heterochromatin regions. However, genome-wide P(s) analysis shows no evidence of significant differences in contact probability for different types of chromatin regions, including accessible promoters and enhancers, Polycomb-bound and heterochromatic domains, within a 500-kb range (Supplementary Fig. 2, panel iii). • As we showed in our previous work 6 , very deeply sequenced Hi-C data 9 can recapitulate some E-P loops detected by Micro-C (Supplementary Fig. 2, panel iv). While these dots in Hi-C are much fuzzier and their precise locations are often indiscernible, the example below clearly indicates that they are present in both data. It is widely accepted in the 3Dgenome field that Hi-C has minimal bias towards highly accessible chromatin, so contact enrichment around E-P interactions present in both methods cannot be an artifact of Micro-C. Critically, Micro-C can further resolve what appears as a large blurry dot in Hi-C into multiple finer-scale dots and pinpoint the precise location of their interactions (Supplementary Fig. 2, panel iv, zoomed-in box), indicating that Micro-C is a superior method to Hi-C in identifying fine-scale interactions.

1D coverage does not explicitly affect loop calling in the high-resolution data
To test whether loop calling is affected by coverage in our data, we analyzed "1D coverage vs. the number of significant loops per locus". Supplementary Fig. 3 below shows a similar trend across all resolutions, suggesting that 1D coverage does not explicitly affect dot calling in the highresolution data (Wilcoxon Rank Sum and Signed Rank Tests). We also included the balanced data that demonstrates the number of loops is not dependent on the coverage. Supplementary Fig. 2. Quality control of Micro-C analysis. i) Scatter plots with density heatmap compare the correlation between Micro-C reads coverage per kb (x-axis) and ATAC-seq signal per kb (y-axis). Neither raw Micro-C data (r 2 = 0.11) nor the balanced data (r 2 << 0.01) show substantial correlation with ATAC-seq signal. ii) Histogram shows the genome-wide averaged signal enrichment of ATAC-seq (left), Micro-C raw signals (middle), or Micro-C balanced signals across a ±3-kb region. Data are grouped by the intensity of ATAC-seq from high to low. Raw Micro-C signals are slightly higher in regions with low ATAC-seq signals than in regions with high ATAC-seq signals. The balanced data shows homogeneous coverage throughout the genome. iii) Genomewide contact decaying P(s) analysis across the regions enriched with different chromatin features, including promoter, enhancer, polycomb, heterochromatin, and random loci. iv) Snapshots of Micro-C vs. Hi-C around Igf1r region.

Micro-C analysis is not biased by global changes in chromatin interactions
To test whether these normalization approaches can address global changes in genome organization and correct systemic biases, we computationally injected 3x more long-range interactions (> 2Mb) and inter-chromosomal interactions into wild-type Micro-C data (Supplementary Fig. 4, top panel). In this case, the artificially injected reads generated an extreme case of reads being redistributed to longer ranges (Supplementary Fig. 4, top panel). Importantly, the results of compartment strength, TADs, cohesin loops, and E-P/P-P loops are not affected by the global changes caused by artificial reads (Supplementary Fig. 4, bottom panels). We therefore conclude that matrix balancing and distance normalization can reveal the bona fide chromatin structures and differential chromatin interactions regardless of systemic biases and global changes. Supplementary Fig. 3. Quality control of loop calling with Micro-C data (raw and balanced signals) at different resolutions. The curves show the relationship between the number of loops being identified per locus (x-axis) and Micro-C reads coverage (y-axis; Top group=raw; Bottom group=balanced). A similar trend across all resolutions suggests that 1D coverage does not explicitly affect loop calling in the high-resolution data (Wilcoxon Rank Sum and Signed Rank Tests; n=37 biologically samples; Data are presented as mean values SD).
Supplementary Fig. 4. Micro-C analysis is not biased by global changes in chromatin interactions.

II. Technical controls of auxin-inducible degradation 1. The effect of basal degradation due to degron leakiness
To determine to what extent degron leakiness may affect our "untreated" cells, we quantified CTCF, RAD21, WAPL, and YY1 protein levels in each degron clone and compared them to wildtype cells by Western blot (Extended Data Fig. 2b-c). We also conducted correlation test, reproducibility tests and differential analyses of all untreated samples in ChIP-seq (Extended Data Fig. 3a-c), Micro-C (Extended Data Fig. 4a-d), and RNA-seq data (Extended Data Fig. 5a).
Although CTCF-AID, RAD21-AID, and YY1-AID cell lines show some basal degradation (~40%, ~30%, ~25% reduction of protein levels compared to wild type cells, respectively), we find no significant changes in chromatin association, 3D genome organization, or transcriptome in these untreated cells. The results of our quality controls are summarized below in Supplementary Fig.  5a.
Consistent with differential expression analysis identifying only a few significantly altered genes after 3-hour depletion, Principal component analysis (PCA) of RNA-seq and mNET-seq revealed that all the clonal cell lines derived from the C59 parental clone (see Methods for details; control (C59), ∆CTCF (C58), ∆RAD21 (F1), ∆WAPL (C40) clones) with or without IAA treatment are clustered together, except for cells depleted of RAD21 after 12 and 24 hours (DEGs > 1000). Additionally, we found that YY1-depleted cells are clustered separately from C59-derived cells and are more similar to wild-type JM8.N4 cells. The results of PCA are summarized below in Supplementary Fig. 5b.
deregulation (over thousands of DEGs), possibly because of the secondary effects such as cell cycle defects (Fig. 3h and Extended Data Fig. 5c-d). While our WAPL depletion clearly increases the fraction of cohesin retention on chromatin (~35%) (Extended Data Fig. 2f) and expands the grid of cohesin loops in Micro-C maps (Extended Data Fig. 4h), we do not detect a global effect on transcription across the time course (~10 DEGs, padj < 0.01 & > 2-fold change) (Extended Data Fig. 5c-d), contradicting the previous study that found ~779 DEGs (padj < 0.05) after 6-hour WAPL depletion using TT-seq 11 . We noticed that we have used a more stringent cut-off to call DEGs than other studies. By adjusting the threshold to padj < 0.05, we can detect more DEGs in the range comparable to previous studies (DEGs=~394).
Together, we conclude that 1) our mNET-seq method is sensitive to massive transcriptional alterations; 2) our mNET-seq results largely agree with previous reports; and 3) many of the DEGs detected using "padj < 0.05" fail to survive with "padj < 0.01 & 2-fold change", indicating that these changes to gene expression are mild (see Table below).

III. Additional ChIP-seq analysis compared to previous studies 1. Proteins that enriched at Enhancer, promoter, Cohesin loop anchors
We plotted ChIP-seq data for various histone marks, TFs, and chromatin remodelers over the major types of loop anchors (Supplementary Fig. 7). Consistent with our previous characterization of the subtypes of chromatin structures below the scale of TADs 6 as well as a recent imaging study 13 , loop anchors enriched in CTCF and cohesin generally do not colocalize with sites of active transcription. In contrast, various TFs, coactivators, and Pol II are associated with E-P and P-P anchors.

Cohesin occupancy at active promoters
It has been reported that cohesin is also positioned by actively transcribing Pol II and by the transcription machinery 14,15 but the ChIP signal at these putative "CTCF-independent" loci is generally quite weak 16 (Supplementary Fig. 8a). Specifically, promoter-bound cohesin peaks were barely detectable in untreated cells, and though they slightly increased after CTCF depletion, we could still only detect a few thousand upregulated cohesin peaks upon CTCF degradation (Supplementary Fig. 8a-b) 14 . We thus suggest that in the absence of CTCF, the transcription machinery may halt cohesin extrusion but not as effectively as CTCF does. Whether the promoterbound cohesin peaks detected after CTCF depletion are functional remains to be determined.

WAPL and cohesin occupancy near SOX2/OCT4 binding sites
WAPL depletion leads to cohesin retention on chromatin 17 , appreciable both by ChIP-seq ( Fig.  3c-e, Extended Data Fig. 3e-g) and by salt fractionation, (Extended Data Fig. 2f, orange box). A recent study reported that thousands of cohesin peaks are lost near SOX2 and OCT4 binding sites after WAPL depletion 11 . Nonetheless, we did not observe such change (Supplementary Fig.  9), perhaps due to differences in WAPL degradation levels or the duration of depletion. Taken together, our ChIP-seq data confirmed effective degradation of loop extrusion factors within 3 hours and largely recapitulated previous observations 10, 16,18 . Supplementary Fig. 7. Heatmaps of mNET-seq and ChIP-seq signal enrichments in a 3-kb window around the four primary types of loop anchors. Supplementary Fig. 8. a) Snapshot of an example showing the position and level of cohesin at promoters (red shade) after CTCF degradation. b) MA plots showing the differential ChIP-seq peaks of RAD21 and SMC3 in the CTCF-depleted cells. The differential peaks (adjusted p-value < 0.05) are colored in red. X-axis: mean observations of UT and IAA cells. Y-axis: log2 fold-change comparing the UT and IAA-treated cells.

Loop extrusion model
Elegant experiments combining acute protein depletion of CTCF, cohesin, and cohesin regulatory proteins with Hi-C or imaging approaches have revealed the role of CTCF and cohesin in regulating the first two levels, TADs and compartments 10,[18][19][20][21] . These studies have shown that while CTCF plays only a minor role in compartmentalization, CTCF and cohesin removal largely eliminates TADs and chromatin loops anchored by these proteins across the genome. CTCF and cohesin are thought to form TADs and loops through DNA loop extrusion 22,23 . This model posits that cohesin extrudes bidirectional loops until it encounters convergent and occupied CTCF binding sites. When averaged across cell populations, the extruded chromatin appears to be spatially organized into a self-interacting domain (TAD or loop domain), and CTCF binding sites constitute the domain boundaries that restrict inter-domain interactions. Halting of cohesin extrusion at CTCF sites sometimes gives rise to sharp corner peaks in contact maps, known as loops or corner dots. WAPL, a cohesin unloader, releases cohesin from chromatin and WAPL depletion therefore increases cohesin residence times as well as the amount of cohesin on chromatin 19,21 .

Example of the consequences of perturbing CTCF/Cohesin
Paradigmatic experiments focusing on mouse development suggested that TAD disruption through inversion or deletion of CTCF sites around developmental genes such as Epha4, Kcnj, or Ihh can cause severe limb malformation 24 . Similar studies at the locus encoding the morphogen Sonic hedgehog have led to somewhat inconsistent effects on transcription and developmental phenotypes, perhaps due to manipulation of different CTCF sites 25,26 . Furthermore, CTCF and cohesin appeared to be crucial to some biological processes such as neuronal maturation 27 and lipopolysaccharide-induced inflammatory response 28 , but their presence seemed dispensable in other cases such as neuronal activity-dependent transcription 27 and immune cell transdifferentiation 28 .

E-P/P-P interacts across TAD boundaries
Our data reporting that ~30% of loops cross a TAD boundary are largely consistent with previous studies. First, a recent CRISPRi-based enhancer screening reported that ~30% of E-P interactions Supplementary Fig. 9. Profiles of genome-wide averaged ChIP-seq for CTCF, RAD21, SMC1A, and input around the 3kb region of SOX2 or OCT4 peaks in WAPL-depleted cells.
do not fall into the same TAD 29 . Second, promoter Capture Hi-C experiments 30 also observed that ~24% of promoter-centered interactions are not constrained within the same TAD, which aligns with our findings. Third, emerging evidence by single-cell Hi-C and super-resolution imaging indicates that TADs are dynamic chromatin structures that constantly form and dissolve [31][32][33][34][35] . Associations between regions in the same TAD are not necessarily pervasive and interactions are common between neighboring TADs 36,37 . Cohesin can traverse across boundaries frequently and mediate domain intermingling 38,39 . Fourth, the latest evidence indicates that the perturbations of TAD boundaries around the Sox2 and HoxD regions are mostly harmless to embryo development 40,41 , suggesting that high affinity enhancer-promoter interactions can overcome structural barriers. Finally, since the eQTL analysis in Delaneau et al., is based on the correlation of a set of ChIP-seq data for histone modifications, the cis-regulatory domains (CRDs) identified in this study most likely represent local chromatin domains segregated by chromatin states (e.g., active or inactive regions) rather than long-range E-P loops. It is uncertain whether this computational structure fully reflects 3D genome structures in situ. A recent study indicates that analysis of chromatin loops outperforms eQTL in explaining neurological GWAS results, suggesting there may be much ampler information in high-resolution 3D genome maps than in eQTL analyses 42 . Together, our findings of cross-TAD interactions are largely consistent with other studies, arguing that communication between enhancers and their cognate genes in another domain is not uncommon.

What are the 10-20% CTCF/cohesin-sensitive E-P/P-P loops?
Previous studies reported that CTCF and cohesin appear to mediate the long-range E-P loops but have a lesser impact on short-range loops 43 . Our genome-wide analysis reveals that long-range chromatin loops are more susceptible to disruption (Extended Data Fig. 6g), with CTCF/cohesinsensitive loops being much larger (median = ~85 kb for E-P and ~102 kb for P-P) than the insensitive loops (median = ~41 kb for E-P and ~71 kb for P-P). According to the loop extrusion model, cohesin might either directly bridge E-P interactions at TAD boundaries, or the process of extrusion might inherently increase the frequency of long-range E-P interactions 44 . To test this possibility, we compared the genomic features of the CTCF/cohesin-sensitive loops relative to TAD boundaries. CTCF and cohesin occupancy at the anchors of these loops is much higher than at the unaffected loop anchors (Extended Data Fig. 6h), suggesting that these loops may overlap with TAD boundaries, where CTCF and cohesin proteins are highly enriched. However, only ~5% of the CTCF/cohesin-sensitive loop anchors are located within a 1-kb window surrounding TAD boundaries (Fig. 4f-g). In fact, most of them are located even farther away than the unaffected loops ( Fig. 4f-g and Extended Data Fig. 6i) and tend to interact with another DNA locus within the same TAD without crossing the boundary (Fig. 4g).

Loop at specific CTCF sites
Interestingly, in comparison to a recent finding 12 , our data show only moderate CTCF persistence at architectural sites (e.g., loop anchors or TAD borders) upon CTCF degradation ( Supplementary  Fig. 10a). Interactions anchored by CTCF binding sites associated with mitotic bookmarking 45 and frequently occupied sites in single-molecule foot-printing (SMF) data 46 also did not persist upon CTCF degradation (Supplementary Fig.10b-c). Given the near-complete ablation of cohesin loops after CTCF degradation (Fig. 3f-g), we suggest that the persistent residual CTCF proteins, if any, are insufficient to halt or position cohesin.

Probing YY1 as a candidate regulator of E-P and P-P links
Our finding that E-P and P-P interactions remain largely intact in the absence CTCF, cohesin, and WAPL suggests that other proteins are likely responsible for E-P and P-P interactions. To address this, we searched for factors specifically enriched at E-P and P-P loop anchors. BRD2, BRD4, P300, ESRRB, SP1, Mediator, YY1, pluripotency TFs, and chromatin remodelers are all broadly enriched at enhancer and promoter loop anchors (Supplementary Fig. 7). Among them, Mediator has been proposed to function as a structural complex to bridge E-P interactions 47 and promote the folding of sub-TAD structures 48 , but two recent studies showed that loss of Mediator does not strongly affect E-P interactions 49,50 .
We therefore instead focused on YY1, a multifunctional zinc finger-containing TF (Extended Data Fig. 7a) that is ubiquitously expressed, highly conserved, and essential for embryonic development in mammals 51 . Heterozygous YY1 mutations cause Gabriele-de Vries syndrome, which is characterized by developmental delay and intellectual disability 52 . YY1 has been implicated in chromatin looping, particularly during early neural lineage commitment 53 , and a previous study suggested it to be a master structural regulator of E-P interactions 54 . However, YY1 redistribution on chromatin triggers little or no change in H3K27ac and H3K27ac-anchored HiChIP interactions upon rapid induction of erythroid differentiation 55 . These confounding results led us to investigate the role, if any, of YY1 in mediating E-P interactions using Micro-C and YY1 depletion. Supplementary Fig. 10 a) Profile of differential ChIP-seq signal for CTCF in a 3-kb window around loop anchors (gray) or non-loop anchors (blue). b) APA for loops where CTCF is bound (left) or lost (right) at the anchors during mitosis. c) APA for loops with either high (left) or low (right) CTCF occupancy characterized by single-molecular footprinting (SMF) assay at their anchors7.

Spatial cluster analysis of YY1
To further link the spatial relationship of YY1 proteins with their different diffusion rates, we reconstructed the spatial distributions of trajectories by their likelihood of diffusion coefficients at 0.02, 0.1, 1, 5, and 20 µm²s -¹. Consistent with the live-cell imaging and PALM results, the immobile fraction showed apparent cluster-like structures with ~12 -30 trajectories per spot (from a total of ~15,000 trajectories per cell), while faster moving YY1 subpopulations showed weaker propensity to cluster within a constrained area (Supplementary Fig. 11a). Overlays of the immobile trajectories with the other fractions showed that the molecules within the bound regime had nearly complete overlap, but the faster-diffusing molecules were less likely to co-localize with the immobile clusters (Supplementary Fig. 11b). The averaged (Supplementary Fig. 11b, insets) and differential cluster signals (Supplementary Fig. 11b, bottom panel) between the immobile and the other fractions (Supplementary Fig. 11a) further confirmed that the faster-moving molecules more frequently travel to areas in the vicinity of the immobile clusters. These results are consistent with chromatin bound YY1 proteins forming clusters, while the diffusing YY1 molecules appear to search the 3D nuclear space outside of the clusters. We previously made similar observations for CTCF [56][57][58] . Supplementary Fig. 11. a) Spatial reconstruction of spaSPT data for YY1's trajectory densities. YY1 trajectories are binned by diffusion coefficients, as indicated (left). Segmentation of YY1 clusters with the reconstructed images from the spaSPT data. Individual YY1 clusters are colored (right). b) Spatial reconstruction of spaSPT data for YY1's trajectory densities. YY1 trajectories are binned by diffusion coefficients as indicated. The bound trajectories are colored in red. The diffusing trajectories are colored in green and overlaid with the bound trajectories. The insets show the averaged signals for all identified clusters aggregating at the center of the plot with the same imaging overlays. Differential signals comparing each diffusing fraction to the bound fraction are plotted at the bottom. The divergent colormap shows a higher signal than the bound in red and lower than the bound in blue.

Dynamics of YY1 vs CTCF
While YY1's average residence time of ~13-60 seconds is similar to that of many TFs, it is much shorter than the residence times of known structural factors such as CTCF (~1-4 min) and cohesin in G1 (~20-25 min) in mESCs 56 . These results may explain why YY1 depletion has a marginal effect on chromatin looping compared to CTCF and cohesin depletion. CTCF/cohesin loops are generally stronger and almost completely lost upon CTCF/cohesin depletion (Fig. 3f-g), whereas YY1 loops tend to be weaker and less affected by YY1 depletion (Fig. 5e&g). Furthermore, the chromatin bound fraction of YY1 (~30%) is considerably smaller than that for CTCF (~50-70%) 56 .

Consideration of phase-separation mechanism for gene regulation
Phase-separated transcriptional condensates have been proposed to mediate gene activation and E-P interactions 59 . Pol II, Mediator, BRD4, and many TFs containing intrinsically disordered regions (IDRs) that tend to aggregate into local-high concentration hubs in the nucleus [60][61][62] . Hub formation is thought to be a general property of TFs used to engage with regulatory elements and keep them in the spatial vicinity for gene regulation [59][60][61][62] . While not specifically addressed in this study, recent studies employing acute depletion/inhibition of Mediator and BRD4, which are supposed to be key players in condensate formation, also found no drastic effects on E-P interactions 49,50,63 . Thus, condensates also do not seem to be generally required for the maintenance of E-P interactions.

Evidence of CTCF and cohesin regulates gene expression
The evidence that CTCF and cohesin can directly or indirectly regulate E-P interactions and affect gene expression in many cases is overwhelming. This evidence includes the following: 1) Insertion of CTCF sites between an enhancer and a promoter can both reduce E-P interactions and strongly reduce gene expression [64][65][66] ; 2) CTCF binding site silencing 67,68 or genetic CTCF binding site loss 26,69,70 can cause aberrant E-P interactions and gene expression and drive disease; 3) Inversion or repositioning of CTCF sites can redirect E-P interactions that cause gene misexpression and diseases 71,72 . Two recent studies have also proposed variants of a time-buffering model based on mathematical modeling of E-P interactions and gene expression 73,74 . In both models, individual E-P interactions are memorizedeither as long-lived promoter states 66 or as long-lived "promoter tags" 73such that gene expression can be temporally uncoupled from E-P interactions, yet still be causally linked.

Consideration of other mechanisms by which cohesin regulates E-P/P-P loops
We probed the nuclear dynamics of two additional transcription factors, KLF4 and SOX2, by spaSPT and found that their bound fractions were reduced by ~20% 3 hours after cohesin degradation (Extended Data Fig. 9g). Furthermore, while our paper was under review an independent study from Gordon Hager's group 75 found precisely the same thing: for the TF glucocorticoid receptor, cohesin depletion leads to decreased chromatin binding. Thus, while we agree with the reviewer that our model was too strong a leap based on a single TF (YY1) in the original submission, since we now have evidence for 4 TFs (YY1, KLF4, SOX2, GR) in two different cell types (mESCs, human cancer cells) and from two different labs, we now believe that the evidence is strong enough to invoke a role for cohesin in facilitating TF binding.
Finally, the question that "How is an alternative model in which the act of loop extrusion acts to stir the nucleoplasm, resulting in higher diffusion and faster search times?". However, loop extrusion has the opposite effectit makes chromatin diffusion slower. Indeed, two recent studies ( 76,77 have both shown that cohesin depletion leads to a ~2-4-fold increase in chromatin diffusion dynamics. However, faster chromatin motion is extremely unlikely to noticeable accelerate the TF target search. In other words, the TF search time would be expected to increase by 1% due to cohesin depletion based on chemical kinetics. Instead, we observe the oppositethe 2.2-fold slower TF search for YY1. Thus, both the magnitude (1% increase vs. 54% decrease) and the direction (chemical kinetic predicts faster search, we observe much slower search) are inconsistent with the model proposed by the reviewer. In conclusion, we feel confident in ruling out changes in chromatin diffusion due to cohesin depletion as a possible explanation for the reduction in TF search kinetics.