Multiomic analysis of cohesin reveals that ZBTB transcription factors contribute to chromatin interactions

Abstract One bottleneck in understanding the principles of 3D chromatin structures is caused by the paucity of known regulators. Cohesin is essential for 3D chromatin organization, and its interacting partners are candidate regulators. Here, we performed proteomic profiling of the cohesin in chromatin and identified transcription factors, RNA-binding proteins and chromatin regulators associated with cohesin. Acute protein degradation followed by time-series genomic binding quantitation and BAT Hi-C analysis were conducted, and the results showed that the transcription factor ZBTB21 contributes to cohesin chromatin binding, 3D chromatin interactions and transcriptional repression. Strikingly, multiomic analyses revealed that the other four ZBTB factors interacted with cohesin, and double degradation of ZBTB21 and ZBTB7B led to a further decrease in cohesin chromatin occupancy. We propose that multiple ZBTB transcription factors orchestrate the chromatin binding of cohesin to regulate chromatin interactions, and we provide a catalog of many additional proteins associated with cohesin that warrant further investigation.

Structural maintenance of the chromosome (SMC) complexes (cohesin, condensin and SMC5 / SMC6) are crucial for the organization of interphase chromatin structures (28)(29)(30). Previous genomic distribution mapping efforts suggested that cohesin occupies promoters, enhancers and insulators (31)(32)(33). Cohesin is loaded by NIPBL and MAU2, stabilized by CTCF at insulators, and removed by WAPL (34)(35)(36). These direct regulators of cohesin have been shown to regulate global 3D genome organization (37)(38)(39)(40). Previous studies have identified multiple transcription factors involved in 3D genome organiza tion tha t interact directly or indirectly with cohesin. These factors, including CTCF, YY1, ZNF143, ADNP, MAZ, WIZ and BHLHE40 (41)(42)(43)(44)(45)(46)(47)(48)(49)(50)(51), usually occupy insulators or gene promoters, and their perturba tion modula tes the global 3D chroma tin organization or affects a subset of chromatin interactions. We hypothesized that there are additional transcription factors involved in cohesin-media ted chroma tin interactions. Ther efor e, the compr ehensi v e identification of cohesininteracting partners will provide a valuable r esour ce for the identifica tion of candida te regula tors of 3D genome organization and will support the idea that the redundancy between cohesin and various transcription factors modulates local chromatin interactions or specific interactions under different conditions.
The ZBTB transcription factors Mod (mdg4) and Cp190 are known to play critical roles in the regulation of insulator activities in Drosophila (52)(53)(54)(55)(56). The zinc finger domain interacts with specific DNA sequences, and the BTB domain mediates homo-and heteromeric interactions ( 57 , 58 ). These properties might be responsible for the ability of ZBTB transcription factors to mediate 3D chromatin structur es. Ther e ar e > 40 ZBTB transcription factors in mammals, and many of them have been reported to function in T-cell and B-cell lineage specification during de v elopment and diseases (59)(60)(61)(62). The major focus was previously on their roles in transcription regulation. Howe v er, whether and how they regula te 3D chroma tin structur es r emain unknown.
Her e, we r eport tha t ZBTB21 facilita tes the chroma tin occupancy of cohesin and modulates 3D chromatin interactions and transcriptional r epr ession. A number of ZBTB transcription factors (ZBTB7A, ZBTB7B, ZBTB11 and ZBTB35) have been shown to interact with cohesin. The double depletion of ZBTB21 and ZBTB7B further decreases the chromatin occupancy of cohesin. These results imply that various ZBTB transcription factors facilitate cohesin chromatin binding in addition to performing transcriptional regulatory activities that modulate chromatin interactions in mammalian cells.

Plasmid construction
The donors and single guide RNAs (sgRNAs) for the degron-tagged cells were generated as in the previous report ( 63 ). The homology arms were amplified from mouse genomic DNA with a length of 400 bp upstream or downstream of the targeted site. The homology arms were then ligated with mouse auxin-inducible degron (mAID)-green fluorescent protein (GFP) donor cassettes by Gibson assembly (NEB, E2611S). The ligated constructs were used to target the sequence near the TAG stop codon of the protein of interest by CRISPR / Cas9 gene editing to generate degron cell lines. The plasmid pMK243 (Addgene, 72835) was digested by restriction enzymes BglII and MluI (NEB, R0144S and R0198S). The coding sequences of the target protein and the GFP open reading frame (ORF) were ligated with digested pMK243 by Gibson assembly. The correct construct was used to generate Tet-on inducible expression of ZBTB stable HEK293T cell lines. For primers, see Supplementary Table S9.

Stable cell line generation
The degron donors and corresponding sgRNAs were transfected into the mES-OsTIR1 cell line generated in a previous study ( 64 ). The transfected cells were selected with neomycin (Gibco, 21810-031), and homogenous clones were obtained after the validation by genotyping and Sanger sequencing. The HEK293T cell lines expressing a doxy cy cline-inducib le GFP-tagged ZBTB proteins were generated by inserting the expression ORF into the AAVS1 locus by CRISPR / Cas9 gene editing. The clones were validated by genotyping and Sanger sequencing. In order to not express the ZBTB protein to very high le v els, ZBTB-GFP heterozygous clones were selected for downstream analyses.

BA T -Hi-C
The BAT Hi-C was perf ormed by f ollowing the protocol pub lished pre viously ( 66 ). The major difference between BAT Hi-C and in situ Hi-C is that BAT Hi-C employs a biotin-labeled bridge linker to mediate proximity ligation, while in situ Hi-C uses blunt end ligation. Briefly, the 1% formaldehyde-fixed cells were digested by AluI (NEB, #R0137S) for at least 12 h. Then, the restriction enzyme was inactiv ated b y heat. The DNA overhangs were added with a base dA by Kleno w (3 -5 ex o-) enzyme [40 l of NEB buffer 2, 8 l of 10 mM dATP, 40 l of 10% Triton X-100, 304 l of H 2 O and 8 l of Klenow (3 -5 exo-)]. The dA overhang DNA products were ligated with a bridge linker by T4 DNA ligase (NEB, #M0202L) at room tempera ture for a t least 6 h. Then the chroma tin fractions were isolated, and the unligated DNA was digested by lambda exonuclease. The ligated DNA was de-cross-linked at 65 • C overnight. The ligated DNAs were purified, and then fragmented by sonication. The fragmented DNA was examined by 1% DNA agarose gel electrophoresis. The size of the DNA product was ∼0.2-2 kb. The ligated DNA was purified by biotin-streptavidin pulldown. The DNA library was constructed with the NEBNext ®Ultra ™II DNA Library Prep Kit for Illumina ®. The 300-500 bp DNA products were purified by the Magen gel purification kit. The DNA was sequenced at Hiseq 4000, and 150 bp paired-end reads were obtained for downstream analyses.

4C-seq
The 4C-seq was described in a previous publication ( 67 ). Briefly, the cells were first treated a ppropriatel y and then cross-linked with 1% formaldehyde for 10 min at room temperatur e. The nuclei wer e isolated and treated with 25 l of 10 × cutsmart buffer and 100 U of NlaIII (NEB, #R0125S), and incubated overnight at 37 • C, 900 rpm to digest the chromatin. The next day, the restriction enzyme was inactiv ated b y heat. The pr oximal, digested chr omatin fragments were ligated by T4 DNA ligase (NEB, M0202) by incubating the reaction for at least 6 h at room temperature. After ligation, the DNA was re v erse cross-linked and extracted by an equal volume of phenol:chloroform:isoamyl alcohol (Solarbio, P1013). The purified DNA was fragmented by sonication with Biorupter at high energy 30 s on, 60 s off, for two cycles. The fragmented DNA was examined by 1% DNA agarose gel electrophoresis. The size of the DNA product was ∼0.2-2 kb. To construct a 4C-seq library, a biotin-labeled primer targeting the region of interest was used for primer extension. The amplified DNA was purified by biotin-streptavidin pulldown. The biotin-labeled singlestranded DNA was ligated to an adapter by bridge linkermedia ted liga tion. The liga ted DNA products were amplified by nested polymerase chain reaction (PCR), and 500-900 bp DNA products were purified by the Magen gel purification kit. The purified DNA was sequenced by Hiseq 4000 with 150 bp paired reads.

Immunoprecipitation and western blotting
A total of 20 million nati v e cells were used for immunoprecipitation experiments. A 1 ml aliquot of IP lysis buffer [20 mM Tris-HCl pH 7.5, 150 mM NaCl, 1% Triton X-100, sodium pyr ophosphate, ␤-glycer ophosphate, EDTA, leupeptin, phen ylmethylsuflon yl fluoride (PMSF) and protein inhibitor cocktail] was used to resuspend cell pellets by rotating the cell lysates at 4 • C for 30 min. The supernatant (5% cell lysis supernatant was saved as an input sample) was collected. Then, different antibodies were added to the pr e-clear ed w hole-cell l ysa tes and incuba ted overnight with rota tion a t 4 • C . The protein-antibod y complex es wer e incubated with Protein G for ∼2 h a t 4 • C . The beads were washed by sonification buffer four times and by high-salt buffer once. The supernatant was discarded, and a suitable volume of 1 × SDS loading buffer was added. The eluted IP samples were boiled at 100 • C for 5 min. The protein elution and input samples were used for western blotting analyses. Then the samples were run on dif ferent concentra tion SDSpolyacrylamide gel electrophoresis (PAGE) gels. They were transformed onto nitrocellulose membranes and the membranes were blocked with skim milk and then incubated with diluted primary antibody (1:1000-1:5000) overnight at 4 • C. On the next day, the membranes were washed and incubated with secondary antibodies (1:5000). The membranes were washed and analyzed on a G.E AI 600 RGB imaging system.

Live cell imaging
The cell lines that express proteins of interest with GFP tags were imaged on 35 mm glass-bottom dishes (MatTek, USA) by Nikon A1R-si. Before imaging, the cells were digested by 0.25% trypsin EDTA (Gibco, 25200-056) and transferred to the 35 mm glass-bottom dishes for ∼6 h. Then the cells were imaged by a confocal microscope (Nikon A1R-si). Hoechst 33342 (Beyotime, C1029) was used to stain the nucleus. Laser: 405 nm, 488 nm.
Nucleic Acids Research, 2023, Vol. 51, No. 13 6787 ChIP-qPCR and ChIP-mass spectrometry ChIP-quantitati v e PCR (qPCR) primers were designed and v alidated b y NCBI primer-blast. The Ct v alue of input DNA was used for normaliza tion. Calcula tion formula: %Input = 5%*2 ∧ (Ct Input -Ct Sample ). The results of ChIP-qPCR were displayed by prism 8. At least two replicates were run for each qPCR. The statistical significance was ev aluated b y Student's t -test. The primers for the ChIP-qPCR can be found in Supplementary Table S9.
The suitable concentrations of SDS-PAGE protein gels to run ChIP protein preparations were chosen according to the molecular weight of target proteins. Cut gels of the full lane were used for MS anal yses w hich were performed at the National Center for Protein Sciences at Peking Uni v ersity in Beijing, School of Life Science, Peking Uni v ersity) (Note: the part of the gels containing antibody was separated in to a single tube for MS analyses).

DNA library construction
ChIP DNA libraries were constructed according to the pr eviously r eported TELP library construction protocol ( 68 ) or the NEBNext ®Ultra ™II DNA Library Prep Kit (NEB, E765S). Briefly, DN A was end-r epair ed, and a poly(C) tail was added to the purified end-r epair ed DNA products. Poly(C) tails were extended with biotin-labeled adapters. Biotin-labeled products were captured by streptavidin C1 beads (Invitrogen, 11206D). Quick ligase (NEB, M2200S) was used to ligate another adapter to the biotinlabeled products. The first-round PCR amplification was performed on beads. The second-round PCR amplification was carried out with the first-round PCR products. The final PCR products were run on a 2% DNA agarose gel. DNA products of 200-400 bp were purified by a Magen gel extraction kit. The NEB library construction was according to the NEBNext ®Ultr a ™II DNA Libr ary Prep Kit for Illumina ® by following the manufacturer's instructions. The purified DNA was sequenced by Hiseq 4000 with 150 pair ed-end r eads.

RNA-sequencing
For RNA-seq, the ZBTB21-degron mESCs were induced by doxy cy cline f or 24 h, and treated with IAA f or 24 h. The untreated cells and IAA-treated degron cells were collected. RNA extraction was performed with Trizol by following the manufacturer's instructions. Our RNA-seq library was prepared by using the NEBNext ® Ultr a ™ RNA Libr ary Prep Kit for Illumina ® (NEB, USA), which is regularly used in our lab ( 64 , 69 , 70 ). The libraries were prepared by following the manufacturer's instructions and sent to Novogene and sequenced on Illumina Hiseq 4000; 150 bp paired-end reads were obtained for downstream analyses.

Chromatin fraction preparation
The chromatin was pr epar ed by following the protocol published previously ( 65 ). A total of 10 million cells were incubated with NP-40 cell lysis buffer (10 mM Tris-HCl pH 8.0, 10 mM NaCl, 1.5 mM MgCl 2 , 0.01% NP-40 and protease inhibitors) to lyse cell membranes. The cell suspension was transferred to a 1.25 ml sucrose cushion [24% sucrose (w / v) in NP-40 cell lysis buffer] to extract the nuclei. The supernatant (cytoplasmic fraction) was collected. The nuclei wer e gently r esuspended with 0.5 ml of glycerol buffer [20 mM Tris-HCl pH 8.0, 75 mM NaCl, 0.5 mM EDTA, 0.85 mM DTT, 50% (v / v) glycerol]. A 0.5 ml aliquot of nuclei lysis buffer (10 mM HEPES pH 7.6, 1 mM DTT, 7.5 mM MgCl 2 , 0.2 mM EDTA, 0.3 M NaCl, 1 M urea, 1% NP-40) was added to the resuspended nuclei. After centrifugation, the supernatant was saved as a nucleoplasmic fraction. The pellets were chromatin fractions used in this study; 1 × SDS loading buffer was added to the chromatin pellets and boiled at 100 • C for 5 min. Western blot was used to examine the target proteins.

Chromatin RNA-seq (ChAR-seq)
For ChAR-seq, the ZBTB21-degron mESCs were induced by doxy cy cline f or 24 h, and treated with IAA f or 6 h. The chromatin was pr epar ed by following the chromatin fraction preparation protocol as described above. The chromatin RNA was extracted from chromatin fractions with Trizol by following the manufacturer's instructions. The total RNA samples were firstly depleted of rRNAs and were constructed into next-generation sequencing libraries by following the experimental procedure for RNA-seq library prepar ation. The ChAR-seq libr aries were sequenced in the same way as for RNA-seq libraries.

ChIP-MS analysis
MS raw data were mapped to mouse or human proteins in the Uniprot database by the MS analysis platform in the National Center for Protein Sciences at Peking Uni v ersity, Beijing. Then, the peptide counts of all isoforms of a protein were combined to generate a single peptide-spectrum match (PSM) value for each protein. For mESCs, the candidate proteins were selected if the total peptide counts of the proteins in ChIP samples were at least twice more than those in input samples. For B16 and CH12F3 cells, the candidate proteins were identified by the R package DEP [ Pvalue < 0.05, log2 fold change (FC) > 1] ( 71 ). The missing data were imputed using the 'MinProb' method. The candidate proteins were annotated according to previous publications ( 72 , 73 ), Uniport, Genecards and g:Profiler databases. In addition, the Gene Ontology (GO) determination of candidate proteins was carried out by using R package Cluster Profiler ( 74 ).

ChIP-seq data analysis
Raw ChIP-seq reads were firstly trimmed using Cutadapt software ( 75 ) and then mapped to the genome using Bowtie2 software ( 76 ) in '-v ery-sensiti v e-local' mode. ChIP-seqs that were performed in mESCs were mapped to the mm10 genome. ChIP-seqs that were performed in HEK293T cells were mapped to the hg19 genome. ChIPseqs with spike-in control were mapped to both the mm10 and hg19 genome. Only uniquely and concordantly mapped r eads wer e kept for further anal ysis. The ma pped reads were converted to bigwig format using bamCoverage from deepTools2 software ( 77 ) in 10 bp bins with parameters '-normalizeUsing' set to RPGC. The mapped reads from ChIP-seq with spike-in control in mESCs were converted to bigwig format normalized by the spike-in reads mapped to the hg19 genome. Peaks were called by MACS2 software ( 78 ) using a q threshold of 0.05. The two replicates were filter ed by IDR softwar e ( https://github.com/nboley/idr ) with an idr threshold of 0.05 to obtain repeatable peaks. Peaks that overlapped with the mm10 or hg19 blacklist were removed from the peak lists.
The genomic distribution of peaks was annotated according to the priority of promoters , enhancers , CTCF sites , gene bodies and others. The cis-regulated elements were defined as below: Promoters: ±1 kb from transcription starting site (TSS) of Refseq genes.
Enhancers: H3K27ac peaks that did not overlap with promoters. Super enhancers wer e downloaded from a previous publication ( 79 ).
CT CF sites: CT CF peaks that did not overlap with enhancers and promoters.
Heatmaps and meta-gene profiles were calculated using computeMatrix of deepTools2 software centered around ±3 or ±5 kb of the peaks, SMC1 ChIA-PET anchors ( n = 24 036, in Figure 2 H) and genes. ChIP-seq signals of cohesin, NIPBL and ZBTB21 were counted at acti v e gene promoters. The correlation curves between cohesin, NIPBL and ZBTB21 were fitted and smoothed by LOESS r egr ession. Motifs surrounding ±250 bp of peak summits wer e pr edicted using Homer ( http://homer.ucsd.edu/homer/ ngs/index.html ) software. The genes regulated by each peak were assigned by the R package ChIPseeker ( 80 ). GO analyses were performed using the R package ClusterProfiler ( 74 ).
To define ZBTB21-bound and unbound genes, genes wer e first filter ed by r equiring the existence of ZBTB21 peaks for the bound ones and no peaks for the unbound ones. Then the top 1000 genes were selected as bound genes based on the ChIP-seq signal of ZBTB21 at gene promoters. The bottom 1000 no peak genes were selected as the unbound genes. The mm10 gene set was obtained from the GENCODE database.
The differential cohesin peaks were first identified by the R package Diffbind using a false discovery rate (FDR) threshold of 0.05 to obtain more robust changes in 0-1, 0-6 and 1-6 h. These differential peaks were then merged and clustered using the R package hclust. For clustering analysis, a heatmap was used to visualize the changes of cohesin binding. Due to the high variability of cohesin binding signals at different regions, the z-score, which was calculated by (x -mean) / standard deviation (SD), was used to normalize the data in order to make the change patterns clearer. To determine the optimal cluster numbers, we used the 'Elbow method' to test a range of cluster numbers and finally selected three clusters to reflect the biolo gicall y meaningful changes.
Peaks of ZBTB21, cohesin and CTCF ChIP-seq datasets in Figur e 2 C wer e first merged and then cluster ed by their ChIP-seq signal using the R package k-means. The control peaks were randomly selected from the mm10 genome (excluding the merged peak regions) to the same number of ZBTB21 peaks. ChIP-seq densities of cohesin, RAD21, NIPBL, WAPL and CT CF wer e then plotted according to the clustered and control categories. The heatmap in Figure 6 C was generated using the same method utilizing the merged peaks of the fiv e ZBTB transcription factors, cohesin and CTCF. The 1 / 2 / 3 / 4+ ZBTB peaks in Figure 6 F were defined as the regions co-occupied by one, two, three, four and more than four ZBTB transcription factors.
Statistical tests used in this study were all determined by two-side Wilco x on test or Student's test as mentioned in each figure legend.

Hi-C data analysis
The bridge linkers of raw Hi-C r eads wer e first trimmed using trimLinker of the ChIA-PET2 software ( 81 ). The processed reads were then mapped to the mm10 genome and filtered to obtain valid contact pairs using the HiC-Pro pipeline ( 82 ). RCP duplicated, dangling ended, re-ligated and self-cir cled r ead pairs wer e discarded for further analysis. Long-range ( > 20 kb) intrachromosomal valid pairs of all the samples were randomly sampled to the same depth ( n = 34 M) and merged together between the two replicates ( n = 68 M) to obtain higher resolution. The resolution of our Hi-C library is ∼25 kb, so the interactions < 20 kb are below our Hi-C resolution, which are potential false-positi v e interactions. Ther efor e, we r emoved the interactions < 20 kb in our analyses. Many previous publications also removed the short-distance interactions for subsequent analysis (83)(84)(85). The merged read pairs were converted to '.hic' format using hicpro2juicebox in HiC-Pro software and then converted to '.cool' format using hic2cool software ( https://github.com/4dn-dcic/hic2cool ). The '.cool' format Hi-C maps were balanced by cooler software ( 86 ). The snapshots of contact maps were visualized using Juicebox software ( 87 ) after Knight-Ruiz (KR, also known as balanced) normalization. The Hi-C maps were iterati v ely corrected using the 'balance' command from cooler software.
A / B compartments were called using call-compartments of Cooltools software ( https://github.com/open2c/ cooltools ) and adjusted by H3K27ac ChIP-seq signals. The saddle plots were computed by compute-saddle of Cooltools software at 50 kb resolution and normalized by the expected matrix in '-quantiles' mode. The enrichment score was calculated as (AA + BB) / 2AB, where AA refers to the mean normalized density of bins in the top 20% eigenvector scor es, BB r efers to the mean normalized density of bins in the bottom 20% eigenvector scores, and AB refers to the mean normalized density of bins between AA and BB.
To plot meta-TAD, we firstly calculated the average contact probability for all loci at a certain distance as the expected Hi-C contact value for each chromosome. Then we transformed the KR-normalized Hi-C matrix into an observ ed / e xpected matrix by dividing each normalized observed v alue b y its corresponding expected value at that distance. Second, we carried out a 2D meta feature analysis by piling up individual submatrices into an average matrix using Coolpup software (v0.9.5) with '-rescale' and '-local' par ameters. In gener al, this is similar to meta-gene analysis, commonly perf ormed f or ChIP-seq data. The positions of  (8) Transcription factors (46) Chromatin regulators (29) RNA binding proteins (29) DNA replication (11) Cofactors (15) Cell cycle ( Supplementary Table S9. ( E ) Heatmap showing the ChIP-seq densities of ZBTB21, SMC1, RAD21, CTCF, NIBPL, WAPL and GFP among merged peaks of ZBTB21, cohesin and CTCF. The peaks of ZBTB21, cohesin and CT CF wer e first merged together and then clustered using their ChIP-seq signals by the k-means method. The control peaks were randomly selected from the mm10 genome (excluding the merged peak regions) to the same number of ZBTB21 peaks ( n = 1827). GFP ChIP-seq was performed with an anti-GFP antibody using wild-type mESCs, which served as a negati v e control for ZBTB21. ChIP-seq data of cohesin (SMC1), RAD21 and WAPL were downloaded from pre vious pub lications ( 33 , 90 , 127 ).        Figure 4 E indica tes tha t IAA-induced loops (the loops that occurred after ZBTB21 degradation) distribute at the top longest loops among all the merged loops. Loops were ranked by their loop length and marked in red if they were hits in IAA-gained loop sets. The less responsi v e loops (see below) serve as negative controls. Statistics were computed by permuting the loop labels and recomputing the enrichment score of the IAA-induced or less responsi v e loop sets to generate a null hypothesis ( 89 ). Local Hi-C interactions of promoters , enhancers , CT CF sites and super enhancers were calculated by Coolpup software at 25 kb resolution and normalized by the expected matrix. The interaction frequency enrichment curves were fitted using all the bins that interacted with the bin containing the peaks. ZBTB21-responsi v e loops were defined as the top 500 loops that contained both ZBTB21 peaks and downregulated cohesin peaks; less ZBTB21-responsi v e loops were defined as the bottom 500 loops that did not contain ZBTB21 peaks and down-regulated cohesin peaks. Loops were ranked by the sum of the ZBTB21 ChIP-seq signal and the absolute value of log2FC of cohesin ChIP-seq signal upon ZBTB21 degrada tion a t the ATAC-Seq peaks inside each loop. Because the loop length distribution analyses indicated that loops > 500 kb showed more changes, loops > 500 kb were classified into ZBTB21-responsi v e and less-responsi v e loops. Down-regulated cohesin peaks were identified using cohesin ChIP-seqs before and after ZBTB21 degradation for 6 h by Diffbind software using a P -threshold of 0.05. ZBTB21, differential cohesin and GRO-seq signals were plotted by heatmaps centered around the ATAC-seq peaks in ZBTB21-responsi v e ( n = 7973) and less-responsi v e ( n = 1112) loops. GRO-seq signals were not classified into + or -strands due to the inability to distinguish ATAC-seq peaks between the + or -strands. APA plots of loops were calculated by Coolpup software at 25 kb resolution and normalized by the expected matrix. The enrichment score was computed as the mean normalized density of the central 3*3 bins.
The genes were assigned to ZBTB21-responsi v e and less-responsi v e loops and plotted using the absolute value of log2 ChAR-seq fold change upon ZBTB21 degradation for 6 h (Figure 5 D). The enrichment of downregulated, unchanged and up-regulated genes was shown as an observ ed / e xpected ratio in the same way as in a pre vious publication ( 90 ). The percentages of genes in each category of loops (observed) were normalized by the genome-wide percentage of the genes in each category (expected). The red dots indicate the observ ed / e xpected ratio of ZBTB21responsi v e and less-responsi v e loops. ZBTB21-responsi v e and less-responsi v e loops were randoml y perm uted 10 000 times across the genome to generate distributions of the ratio (the violin plot), and significance was determined by the distributions. The insulation scores of the ice-normalized Hi-C matrix of ZBTB21 were computed by cworld software ( https://github.com/dekk erlab/cw orld-dekk er ) at 25 kb resolution. The insulation scores of the CTCF Hi-C matrix were computed by cworld software at 20 kb resolution. CT CF Hi-C data wer e downloaded from a previous publication ( 25 ).

4C data analysis
Raw 4C reads were first filtered and trimmed by the enzyme sites using trimLinker from ChIA-PET2 ( 81 ) softwar e. The pair ed-end r eads wer e selected for downstr eam analysis by Cutadapt ( 75 ) software if the 5 end of reads1 contains bait sites. The processed reads2 were then mapped to the mm10 genome by bowtie2 ( 76 ) software using '-verysensiti v e-local' parameters. PCR duplica tes, self-religa tions and reads longer than the library size were discarded. The bigwig files were generated by bamCoverage of deepTools Aggrega ted insula tion score changes around ZBTB21-bound and ZBTB21-unbound promoters are shown at the bottom. The insulation scores of the ZBTB21 Hi-C maps were generated at 25 kb resolution. The insulation scores of the CTCF Hi-C maps were generated at 20 kb r esolution. CT CF Hi-C data were downloaded from previous publica tions ( 25 , 41 ). Quantifica tion results of the insulation scor es ar e shown on the right. Significance was determined by the Wilco x on test (**** P < 0.0001). ( D ) GO analysis of up-and down-regulated genes identified by ChAR-seq upon ZBTB21 depletion. The DEGs were identified by DEseq2. ( E ) Boxplot showing the expression changes re v ealed by ChAR-seq for ZBTB21-bound and ZBTB21-unbound genes before and after ZBTB21 degradation. Significance was determined by the Wilco x on test (**** P < 0.0001, ns, P > 0.05). ( F ) The violin plots (left) show the expression changes of the genes in ZBTB21-responsi v e and less-responsi v e loops. Significance was determined by the Wilco x on test (**** P < 0.0001). The violin plots (right) show that genes in ZBTB21-responsi v e loops ar e mor e enriched in up-regulated genes. The observ ed / e xpected ra tio was calcula ted by using the percentages of down-regulated, unchanged and up-regulated genes in loops (observed) normalized by the genome-wide percentage of the genes in each category (expected). The red dots indicate the observ ed / e xpected ratios of ZBTB21-responsi v e and less-responsi v e loops. ZBTB21-responsi v e and less-responsi v e loops were subjected to random permutation 10 000 times across the genome to illustrate the distributions of the ratios (violin plot), and significance was determined by the distributions. Up-or down-regulated genes were identified by DEseq2 (FDR < 0.05)    ( 91 ) software and normalized to the 1 × coverage of the mm10 genome using '-binSiz e 10 -normaliz eUsing RPGC' parameters.
For displaying 4C track, we smoothed the signals in the same way as in a previous publication ( 41 ), i.e. the genome to be displayed was binned using a 50 bp sliding window. At each bin, the smoothed signal was calculated as the mean signal of a 5 kb region around the bin. For 4C quantification, the signal at a specific region (chr8:88940000-89060000) was normalized by the signal of a 4 Mb downstream background region (chr8:84940000-85060000) using the non-smoothed normalized signals.

RNA-seq data analysis
Raw RNA-seq reads were first trimmed using Cutadapt software and then mapped to the mm10 genome using STAR software ( 92 ). Primaril y ma pped reads without PCR duplication were kept for further analysis. The signals of genes downloaded from GENCODE vM25 were obtained using featureCounts software. Differentially expressed genes (DEGs) were identified by the R package DE-seq2 ( 93 ) with an FDR threshold of 0.05.

ChAR-seq data analysis
Raw ChAR-seq reads were first trimmed using Cutadapt software and then mapped to the mm10 and dm6 genome using bowtie2 software (v2.3.5.1). Uniquely mapped reads without PCR duplication were kept for further analysis. The expression of the mm10 genes was counted by feature-Counts software at gene body regions [+300 from the TSS to the transcription end site (TES)] and normalized by the spike-in reads mapped to the dm6 genome in each sample. DEGs were identified by the R package DEseq2 ( 93 ) with an FDR threshold of 0.05.

Distribution per centage, CRISPR scor e and expr ession analysis
ChIP-seq peaks of transcription factors in Supplementary Figure S1 were downloaded from the ENCODE and Cistrome databases (Supplementary Table S8). ChIP-seqs of H3K27ac and CT CF in differ ent cell lines, as well as genes in the mm9 / mm10 and hg19 / hg38 genome, were downloaded to generate the genomic regions of promoters, enhancers and CTCF sites in corresponding cell lines. Peaks were annotated as these cis-elements and clustered based on the percentage of occupancy using the R package hclust. Gene expression data of the transcription factors in 53 different tissues were obtained from the GTEx database. All the samples were grouped and averaged by tissues. The expr ession per centage of each gene was defined as the per centage of tissues whose averaged transcripts per million (TPM) value was > 10. CRISPR scores of human genes were obtained from a pre vious pub lication ( 94 ), which were computed by log2FC in abundance of each sgRNA between the initial and final populations. Genes were ranked by their CRISPR scores and defined as essential if their CRISPR scor es wer e < 0 and the P -value was < 0.05.

Cohesin ChIP-MS identifies known and potential regulators of 3D chromatin organization
To comprehensi v ely identify protein factors that contribute to 3D genome organization, we performed cohesin ChIP-MS in mESCs (Figure 1 A) because cohesin plays a broad r ole in chr omatin ar chitectur e ( 84 , 95 ). Cohesin ChIP-MS identified almost all known components and direct regulators of cohesin complexes (Figure 1 B). The genomic distribution re v ealed by cohesin ChIP-seq among enhancers, promoters and CTCF sites was consistent with that obtained in previous studies (31)(32)(33) (Figure 1A; Supplementary Table S1). These results indicated the success of our cohesin ChIP experiment. The identified proteins were then grouped into different categories according to the UniProt, GeneCards and PubMed databases ( Figure 1B; Supplementary Table S2; see also the Materials and Methods). Many kno wn chromatin or ganization regulators were identified in our cohesin ChIP-MS data; for example, the known chromatin ar chitectur e proteins CT CF, YY1, ZFP143, ADNP and WIZ were among the transcription factors identified by cohesin ChIP-MS (Figure 1 C) ( 25 , 41-44 , 47 ). The exosome complex was predicted to regulate superenhancer chromosomal interactions ( 96 , 97 ). The exosome component DIS3 was among the RNA-binding proteins identified by cohesin ChIP-MS (Figure 1 C). Together, these analyses suggest that our cohesin ChIP-MS data in mESCs are reliable.
To provide further information on cohesin-interacting proteins, we performed genomic localization, expression specificity and essential score analyses. Genomic localization is expected to reveal the types of chromatin interactions regulated by these proteins, but not all these proteins have reliable ChIP-seq datasets in mESCs. We therefor e sear ched the ChIP-seq datasets for cohesin-associated pr oteins fr om different cell lines included in the ENCODE and Cistrome databases, and calculated the occupancy percentages of promoter -, enhancer -and CT CF-binding r egions in the corresponding cells (Supplementary Figure  S1A). We performed gene expression analyses across different tissues by using gene expression datasets from the GTEx database (Supplementary Figure S1B) and conducted CRISPR score analyses to evaluate the essentiality of the cohesin-associated proteins ( Supplementary Figure S1C) ( 94 ). We hypothesized that constituti v ely e xpressed, essential chromatin structural proteins, such as CTCF and YY1, play a general role in 3D genome organization ( 25 , 41 , 42 ), while less essential specifically expressed proteins, such as ADNP and AP-1, play a cell type-specific role or regulate a subset of chromatin interactions ( 44 , 48 ). Although these results were generated from different cells, they integrate useful inf ormation f or in vestigating these cohesin interactors in the future.
It would be informati v e to compare our cohesin ChIP-MS data with the cohesin interactome in different cell types. Pre vious studies hav e identified the cohesin interactome in human HCT116 and MOLT-4 cells ( 98 , 99 ). Howe v er, these MS data are not comparable with our data, as they were genera ted by dif ferent experimental protocols and in different species. We therefore performed cohesin ChIP-MS Nucleic Acids Research, 2023, Vol. 51, No. 13 6797 analysis in mouse B16 cells from skin and CH12F3 B lymphocytes by using a similar approach to the one we used for mESCs. We used the R package DEP (Differential Enrichment analysis of Proteomics data) to perform statistical analyses of these MS data ( P -value < 0.05, log2FC > 1) and identified 260 and 1332 cohesin-interacting proteins in B16 and CH12F3 cells, respecti v ely (the number of interacting proteins was quite different, possibly due to the different efficacy / accuracy of cohesin ChIP-MS in different cells). Our cohesin ChIP-MS data were quite reproducible and were strongly enriched in known cohesin components (Figures 1 D, E), indicating that the quality of these MS data was good. We further compared the cohesin ChIP-MS data of B16 cells and CH12F3 cells with the cohesin ChIP-MS data of mESCs. GO analysis showed that the factors identified by cohesin ChIP-MS were mainly enriched in chromosome organiza tion, histone modifica tion and non-coding RNA processing terms in different cell types (Figure 1 F). The factors that have previously been demonstrated to regulate cohesin, namely CTCF, WIZ and ADNP, were identified in three different cohesin ChIP-MS preparations ( 25 , 41-44 , 47 ). In addition to those known cohesin-interacting proteins , additional transcription factors , RNA-binding proteins and chroma tin regula tors were also detected in the cohesin ChIP-MS preparations (Figure 1 G). Moreover, many transcription factors wer e pr efer entially detected in different cell lines (Figure 1 H), which may be candidate regulators of cell type-specific 3D chromatin structures. These analyses provide a r esour ce listing constituti v e and cell typespecific cohesin-interacting proteins in different cell types.
Many of the proteins identified by cohesin ChIP-MS have not previously been reported to play roles in chromatin organization or to function together with cohesin. We detected four ZBTB transcription factors, ZBTB19, ZBTB21, ZBTB30 and ZBTB35, in mESCs (Figure 1 C, G), which were of considerable interest to us because ZBTB transcription factors are known to regulate insulator activities, as suggested by the functions of ZBTB factors in Drosophila (52)(53)(54)(55)(56). Insulator function has been proposed to involve determining 3D chromatin organization ( 4 , 100 , 101 ). Since ZBTB21 was detected in cohesin ChIP-MS preparations in thr ee differ ent cell types (Figur e 1 G), we primarily investigated the molecular functions of ZBTB21 in this study.

ZBTB21 interacts with cohesin at active promoters
ZBTB21 ChIP-MS and ChIP-seq were performed to valida te the associa tions between ZBTB21 and cohesin. We first introduced a degron-GFP tag into the C-terminus of endogenous ZBTB21 in mESCs ( Supplementary Figure S2A). Genotyping analyses confirmed the success of homozygous knockin (Supplementary Figure S2B). Then, widely used GFP antibodies were chosen to perform ChIP-MS and ChIP-seq in ZBTB21-degron-GFP cells. Consistent with the cohesin ChIP-MS results, cohesin components wer e pr efer entially identified in ZBTB21 ChIP-MS pr eparations ( Figure 2 A; Supplementary Table S2). Interestingly, condensin and SMC5 / SMC6 components were also detected in ZBTB21 ChIP-MS preparations (Figure 2 A). The pr otein-pr otein interactions between ZBTB21 and cohesin were further validated through IP-western blotting analy-ses of nati v e cell lysates but not IgG IP preparations (Figure 2 B), which served as a negati v e control. To gain deeper insight into ZBTB21-mediated cohesin regulation, we performed ZBTB21 immunoprecipitation after ZBTB21 depletion or treatment with different nucleases. The results showed that ZBTB21 forms specific pr otein −pr otein interactions with cohesin ( Figures S2C, D). The genomic localization of ZBTB21 was carefully characterized to obtain functional insights into ZBTB21-mediated regulation. Meta-gene analysis and a pie chart displaying the genomic distribution of the ZBTB21 peaks indica ted tha t ZBTB21 pr efer entially bound gene promoters (Figure 2 C). We identified ZBTB21 chromatin-binding sites in mESCs (Supplementary Tables S1 and S3). These binding sites could be independently validated by ChIP-qPCR ( Figure 2 D; Supplementary Figure S2E). Furthermore, genomic localization analyses indica ted tha t ZBTB21-occupied regions showed weak binding of cohesin and the cohesin loading factor NIPBL, which was significant compared with the signals in randomly selected regions and GFP ChIP-seq signals (Figure 2 E). ZBTB21 peaks belonged to Cluster 1 and were associated with weak cohesin ChIP-seq signals. Cluster 2 did not have ZBTB21 peaks and was associated with strong cohesin binding signals (Figure 2 F).
Furthermore, ZBTB21 and NIPBL preferentially occupy the gene promoters (Figure 2 G). As the ZBTB21 ChIPseq density increased, the chromatin binding strength of cohesin and NIPBL at gene promoters first increased and was then maintained at a fairly consistent le v el (Figure 2 H). ZBTB21-binding sites were enriched at the cohesin loop anchors identified by cohesin ChIA-PET (Figure 2 I). Motif analyses of the ZBTB21 ChIP-seq peaks re v ealed the enrichment of transcription factor [e.g. JUND (a component of the AP-1 family), ZNF143, CTCF and YY1] functions in 3D chromatin organization (Figure 2 J). GO analyses of the ZBTB21-bound gene indicated that these genes were involved in RNA metabolism pathways (Figure 2 K). These results suggested that ZBTB21 was associated with cohesin chromatin binding at promoters, indicating potential involvement in cohesin promoter occupancy.

ZBTB21 depletion leads to decreased chromatin binding of cohesin
To obtain functional insights into the interactions between ZBTB21 and cohesin, we performed cohesin (SMC1) and NIPBL ChIP-seq after ZBTB21 depletion (Suppelementary Tables S1 and S4). The western blotting results confirmed the acute depletion of ZBTB21 proteins in our ZBTB21degron mESCs (Figure 3 A). We also performed cohesin ChIP-seq with equal amounts of Drosophila cells spiked-in. The tornado plot displaying the spike-in and non-spike-in cohesin ChIP-seq signal changes at the ZBTB21-bound and unbound promoters consistently showed that ZBTB21 depletion decreased cohesin binding at ZBTB21-bound promoters but not for the unbound promoters (Figure 3 B). The chromatin binding of NIBPL at ZBTB21-bound gene promoters also displayed a weak decrease (Supplementary Figure S2F). Screenshots of specific genes showed that ZBTB21 depletion led to decreased cohesin and NIPBL binding at the Man1b1, Mmab and Znrd1 gene promoters (Supplementary Figure S2G). We also performed cohesin ChIP after ZBTB21 depletion, and ZBTB21 resumed after withdrawing the IAA. The results showed that cohesin chromatin binding decreased after ZBTB21 depletion and then resumed after the withdrawal of IAA ( Supplementary Figure S2H), suggesting that the roles of ZBTB21 in cohesin chromatin binding are specific. Collecti v ely, these results support that ZBTB21 specifically regulates cohesin chromatin occupancy, e v en though direct pr otein −pr otein interactions with recombinant proteins are currently lacking.
We next carried out time-course cohesin ChIP-seq after ZBTB21 depletion to determine the dynamic effects on cohesin. The clustering analyses of the differential peaks indicated the existence of three clusters (Figure 3 D, E): in Cluster 1, the cohesin signals ra pidl y decreased, and the peaks displayed high-density ZBTB21 signals; in Cluster 2, the cohesin signals first increased and then decreased, and the peaks displayed less ZBTB21 binding; in Cluster 3, the cohesin signals gradually decreased, and the peaks displayed an intermediate le v el of ZBTB21 binding. Moreover, there was an obvious decrease in NIPBL binding in Cluster 1, a small but significant decrease in Cluster 3 and no significant change in Cluster 2 after the rapid depletion of ZBTB21 (Figure 3 E). We then performed Pol II and CTCF ChIP-seq and ATAC-seq to obtain functional insights into ZBTB21-media ted cohesin regula tion (Figur e 3 F). The r esults re v ealed that ZBTB21 depletion leads to global transcriptional activation, as shown by the increased signals of Pol II ChIP-seq and ChAR-seq, suggesting that ZBTB21 functions in transcription r epr ession, which is consistent with the findings of a previous study ( 102 ), but could not simply explain its effects on cohesin chromatin binding. The ATAC-seq signals and CTCF chromatin binding did not show obvious changes after ZBTB21 depletion, indicating that CTCF and chromatin accessibility also could not explain ZBTB21-mediated cohesin binding.
To gain functional insights into ZBTB21 chromatin binding, we performed total RNA-seq and ChAR-seq to enrich nascent RNAs after ZBTB21 depletion (Supplementary Table S5). The differential analyses showed that 473 DEGs in ChAR-seq have ZBTB21 binding, but those are a subset of ZBTB21-bound genes. The number of those genes was higher than the number of genes with similar RNA-seq results (Supplementary Figure S2I). These results are consistent with a recent finding that there is a poor correlation between transcription factor binding and acutely affected gene expression ( 103 ). We also examined the expression changes of subunits of cohesin, condensin and SMC5 / 6, and found that the expression of most of the subunits did not change. Some of the subunits showed up-regulation, but less than our cut-off (Supplementary Figure S2J; Supplementary Tab le S5). Howe v er, we could not completely rule out other secondary effects of ZBTB21-mediated transcriptional r epr ession that might contribute to the reduction in cohesin binding.

ZBTB21 contributes to 3D chromatin interactions
To explore the direct roles of ZBTB21 in 3D chromatin organization, we performed BAT Hi-C after ZBTB21 depletion for 6 h. The BAT Hi-C method, which was pre-viously de v eloped by our labor atory ( 64 , 66 ), can gener ate high-resolution chromatin interaction maps in a more economical way than the classical in situ Hi-C method. The Hi-C data were processed via pre viously estab lished pipelines ( 64 , 66 ). The mapping rate was high and was quite reproducible between the two replicates (Supplementary Figure  S3A; Supplementary Tables S1 and S6). The saddle plot and meta-TAD analyses indicated that ZBTB21 depletion did not obviously affect large-scale chromatin interactions in mESCs, such as A / B compartments and TADs (Figures  4 A, B), and the decay curve of Hi-C interaction frequencies across different genomic distances likewise indicated no change (Figure 4 C), suggesting that ZBTB21 was not r equir ed for global, large-scale chromatin interactions in mESCs.
Chromatin loops were then compared before and after ZBTB21 depletion with IAA. Chromatin loops were identified with Mustache software ( 88 ), and Hiccups produced results with a similar tr end. Ther e wer e 5606 loops identified in the untreated condition and 6423 loops identified under IAA treatment. Overla p anal yses indica ted tha t 2867 loops were induced after ZBTB21 depletion. Density distribution analyses and loop strength changes across loops of different lengths indicated that loop str ength incr eased slightly after ZBTB21 depletion (Figure 4 D; Supplementary Figure  S3B). GSEA showed that IAA-induced loops were more enriched in larger loops (Figure 4 E). Accordingly, the average loop size increased from 417 kb to 455 kb after IAA treatment of ZBTB21-degron cells (Supplementary Figure S3C). We then investigated the features of the new Hi-C loops. The results showed that the new Hi-C loops were more enriched with regulatory elements (such as CTCF / cohesin peaks , enhancers , promoters or polycomb regions) than the less ZBTB21-responsi v e loops (Figure 4 F).
We performed chromatin interaction analyses focusing on cis-regulatory elements (enhancer, promoter and CTCFbinding regions) using the ZBTB21 depletion Hi-C dataset. The results showed that ZBTB21 depletion increased the Hi-C interaction signals at promoters , enhancers , CTCF sites and superenhancers (Supplementary Figure S3D). The aggregated analysis of Hi-C interaction frequencies centered on ZBTB21-bound promoters compared with promoters without ZBTB21 binding showed that ZBTB21 depletion significantly increased 3D chromatin interactions at ZBTB21-bound promoters but not at unbound promoters (Figure 4 G). For example, at highly transcribed gene loci (including Lrrc75a, Ncor1 and Pigl), the Hi-C interaction fr equencies incr eased, and these r egions also showed ZBTB21 binding and decreased cohesin binding after the depletion of ZBTB21 (Supplementary Figure S3E). At the Sall1 and Chd9 loci, we observ ed acti v e transcription, ZBTB21 binding and decreased binding of cohesin after ZBTB21 depletion. The interactions between Sall1 and Chd9 increased after IAA treatment of ZBTB21-degron cells, which could be independently validated by 4C-seq after ZBTB21 depletion (Figure 4 H). Although the effects of ZBTB21 on chromatin interactions are modest, they are statisticall y significant, w hich is reminiscent of a previous report that subtle changes in 3D chromatin structures can have a large functional impact, for example, on transcription ( 104 ). It seems that all cohesin peaks in these regions, Nucleic Acids Research, 2023, Vol. 51, No. 13 6799 independent of ZBTB binding, are smaller in IAA-treated samples than in untreated samples, suggesting that ZBTB21 may regulate cohesin loading or non-specifically modulate cohesin chromatin binding.
We next defined ZBTB21-responsive and less-responsive loops based on the existence of ZBTB21 ChIP-seq peaks and down-regulated cohesin peaks ( Figure 5 A; Supplementary Table S7). The ZBTB21-responsive loops showed higher transcriptional activity (Figure 5 A). Then, the loop interaction frequencies displayed in APA plots indicated that ZBTB21 depletion caused greater increases in the loop strength of the ZBTB21-responsi v e loops than in that of the less-responsi v e loops (Figure 5 B). These results suggested that ZBTB21 binding is associated with decreased cohesin occupancy and increased loop strength after the acute depletion of ZBTB21. Genomic regions with more ZBTB21 binding but less CTCF binding exhibited a greater increase in insulation scores after ZBTB21 depletion. Regions with higher densities of CTCF binding but lower densities of ZBTB21 binding displayed a greater increase in insulation scores after CT CF depletion (Figur e 5 C). The insulation scores at ZBTB21-bound promoters also showed a more significant increase than those at promoters not bound by ZBTB21 after ZBTB21 depletion (Figure 5 C). The insulation scores of ZBTB21-unbound promoters also appeared to be sensiti v e to ZBTB21 depletion, indica ting tha t the ef fects of ZBTB21 on 3D chroma tin interactions were not specific. The GO enrichment analysis of the DEGs showed that the up-regulated genes were enriched in terms such as ribonucleoprotein complex biogenesis, mRNA processing and histone modification, and the down-regulated genes were enriched in regulation of cellular component size, dendrite de v elopment and actin polymerization or depolymerization (Figure 5 D). We then explored the correlation between gene expression and the chromatin occupancy of ZBTB21. We found that the expression of ZBTB21-bound genes tended to increase more upon ZBTB21 depletion (Figure 5 E). Genes associated with ZBTB21-responsi v e loops were also more enriched for upregulated genes than those associated with less-responsi v e loops (Figure 5 F). These results collecti v ely suggested that ZBTB21 is associated with low le v els of CTCF and acti v e transcription, which involv es the genomic distribution of cohesin and changes in local chromatin interactions between promoters and their regulatory elements within the range of se v eral hundred kilobases. These changes in cohesin distribution and chromatin ar chitectur e have been extensi v ely associated with cohesin effects on cell type-specific programs in many different cell systems, including mESCs and human cancer cells, as reported previously (105)(106)(107). The biological relevance of these regulatory mechanisms is worthy of further investigation in various cell systems in the future.
Further analyses indicated that 3D chromatin interactions increased more at genes with increased Pol II occupancy, and the ZBTB21-responsi v e loops also ov erlapped more with the loops of genes with increased Pol II occupancy. These results indicated that the increased long-range interactions wer e corr ela ted with transcriptional activa tion after ZBTB21 activation (Figure 5 G). Howe v er, it is still unclear whether the increased long-r ange inter actions were the cause or consequence of transcription activation after ZBTB21 depletion.

Mor e ZBTB tr anscription factors ar e associated with the chromatin binding of cohesin
Ther e ar e > 40 ZBTB transcription factors in mammalian cells ( 57 , 58 ). Inspired by the results obtained for ZBTB21, we hypothesized that ZBTB transcription factors might generally function together with cohesin. We ther efor e chose more essential ZBTB transcription factors based on their CRISPR scores ( 94 ) or cell type-specific functions ( 59 , 108 ), including ZBTB7A, ZBTB7B, ZBTB11 and ZBTB35, and generated GFP-tagged stable HEK293T cell lines with inducible expression of these four factors as well as ZBTB21 (Supplementary Figure S4A). For side-by-side comparison of the protein interacting partners among different ZBTB factors, clones with relati v ely similar, reasonable (not too high or too low) expression levels were used for downstream analyses (Supplementary Figure S4B). We then performed GFP ChIP-MS analysis of these cell lines. The results showed that these ZBTB transcription factorinteracting partners shared many transcription factors and chroma tin regula tors (Supplementary Figure S4C). Surprisingly, ZBTB ChIP-MS also identified many components of cohesin, condensin and SMC5 / SMC6 complexes that were not found by GFP ChIP-MS performed in a cell line that expressed only GFP as a negati v e control (Figure 6A; Supplementary Tab le S2). IP-western b lotting e xperiments further validated the interactions of these components with cohesin (Figure 6 B). E2F8 IP did not detect SMC1 [E2F8 was shown not to connect to cohesin components in a previous study ( 109 )]. Furthermore, we expressed these ZBTB factors with a smaller tag (HA tag) in HEK293T cells, and HA antibody immunoprecipitation followed by western blotting showed that these ZBTB factors pulled down cohesin. Cohesin IP in HEK293T cells reciprocally verified these interactions ( Supplementary Figure S5A). Together, these results indicate that ZBTB transcription factors specifically interact with cohesin.
ChIP-seq experiments were then performed for ZBTB7A, ZBTB7B, ZBTB11, ZBTB21 and ZBTB35. Their genomic binding sites were distributed at promoters, enhancers and CTCF sites in the genome, and this binding could be validated by ChIP-qPCR with two different antibodies (Supplementary Figure S5B), suggesting that we captured specific binding sites of these ZBTB factors. The analyses of genomic distributions, motifs and GO functions indicated tha t dif ferent ZBTB factors recognized specific DNA sequences and were associated with di v erse biological functions (Supplementary Figure S6A-E). Interestingly, the ZBTB21 genomic distribution and GO functions seemed to differ between 293T cells and mESCs. Furthermore, the chromatin binding strength of these ZBTB factors was positi v ely correlated with the expression le v els of the target genes, which was consistent with previously known roles of ZBTB in recruiting transcriptional regulators (Supplementary Figure S6A-E) ( 108 , 110 , 111 ).
Low le v els of cohesin were observ ed at acti v e promoters where ZBTB factors co-localized, but the highest cohesin binding was observed at CT CF sites wher e ZBTB factor binding was very weak (Figure 6 C-E), and the chromatin binding strength was positi v ely associated with the occupancy of acti v e promoters by cohesin (Figure 6 D). Our k-means clustering analyses indicated that the fiv e ZBTB transcription factors appeared to co-localize in the genome but also re v ealed preferentially occupied sites (Figure 6 C). ZBTB11 and ZBTB35 showed pr efer ential binding clusters and stronger co-localization with Pol II, consistent with the recent identification of their functions in transcriptional r egulation ( 112 ). We pr edicted that if other ZBTB transcription factors affected cohesin similarly to ZBTB21, then sites showing higher le v els of ZBTB factor binding would also show more cohesin binding. Indeed, cohesin binding gradually increased as more ZBTB factors were bound (Figure 6 F).

Acute ZBTB21 and ZBTB7B depletion leads to a further decrease in cohesin occupancy
We chose to degrade ZBTB21 and ZBTB7B because they showed stronger ChIP-MS signals for the cohesin and condensin subunits than the other ZBTB proteins that we investigated (Figure 6 A). We first degraded ZBTB7B in mESCs (Figure 7 A) and then performed cohesin ChIP-seq after ZBTB7B depletion. The illustration of meta-and singlegene examples showed that ZBTB depletion led to a slight but significant decrease in cohesin binding at ZBTB7Bbinding sites (Figure 7 B, C). To directly investigate the functional relationships among ZBTB transcription factors, we generated a double degradation system to sim ultaneousl y deplete ZBTB21 and ZBTB7B by knocking degron tags into the C-termini of the two genes in mESCs (Figure 7 D). Western blotting confirmed the simultaneous degradation of ZBTB21 and ZBTB7B (Figure 7 E).
We first performed chromatin fractionation analyses after ZBTB21 and ZBTB7B depletion. Double degradation caused changes in the protein le v els of chromatin-bound SMC1 and SMC4 that were not statistically significant (Figure 7 F; Supplementary Figure S7A). Consistently, the total protein le v els of the cohesin and condensin components remained similar to those in the whole-cell extract (Supplementary Figure S7B), suggesting that the protein le v els of the cohesin and condensin components SMC1 and SMC4, at least, did not obviously change after ZBTB21 depletion. We then performed cohesin ChIP-seq after ZBTB21 and ZBTB7B depletion. The specific gene locus and meta-gene analyses indica ted tha t double depletion caused a greater decrease in cohesin binding than the depletion of ZBTB21 alone for the ZBTB21-bound promoters (Figure 7 G, H; Supplementary Figure S7C). Howe v er, cohesin binding did not show obvious changes among genes that were not bound by ZBTB21 (Supplementary Figure S7C), which is consistent with the local changes in cohesin binding after ZBTB21 depletion.

DISCUSSION
Cohesin plays critical roles in diseases and cell identity ( 113 , 114 ), and its direct regulators have been functionally linked to 3D chromatin structures ( 25 , 45 , 47 ). The proteomic profiling of cohesin-interacting partners can improve our understanding of the molecular mechanisms underlying 3D chromatin organization. Here, we report dozens of cohesin-interacting proteins associated with chromatin, and we re v eal that ZBTB transcription factors interact with cohesin. The depletion of ZBTB21 leads to a decrease in the occupancy of cohesin and to increased local 3D chromatin interactions for ZBTB21-bound promoters and transcriptional acti vation. Moreov er, the doub le depletion of ZBTB21 and ZBTB7B causes a further decrease in cohesin chromatin binding. Our results re v eal that cohesin interacts and cooperates with ZBTB transcription factors to contribute to 3D chromatin ar chitectur e and gene expression regulation.
The well-established CTCF-cohesin complex is responsible for the formation of large-scale chromatin structures, such as TADs. Howe v er, the molecular basis of smallscale chromatin interactions (such as chromatin loops associated with acti v e genes and sub-TAD structures) has remained unclear. Many cofactors and chromatin regulators have been implicated in the formation of biomolecular condensates to mediate local chromatin structures ( 20 , 115 , 116 ). Howe v er, recent e xperiments involving the chemical inhibition and rapid protein degradation of these factors did not re v eal obvious effects on small-scale chromatin structures ( 23 , 24 , 117 , 118 ). In this study, we provide evidence indica ting tha t ZBTB21 is both physically and functionally linked to cohesin: (i) cohesin ChIP-MS analysis identified ZBTB21 (Figure 1  Our findings were also consistent with the previous observation of strengthened interactions between superenhancers after cohesin depletion ( 26 ).
We obtained 25 kb resolution Hi-C maps in the current study, which is not high enough. This could be why we did not observe robust effects on 3D chromatin loop changes. Higher resolution chromatin interaction mapping (such as Micro-C) and elucidation of the detailed pr otein −pr otein interactions between cohesin and ZBTB21 would be valuable to clarify the molecular mechanism of ZBTB21 function in the futur e. A pr evious stud y showed tha t transcription factors are enriched at local cohesin-binding sites, and the pioneer transcription factors OCT4 and SOX2 create an open chromatin conf ormation f or cohesin binding ( 90 ). We identified many transcription factors in cohesin ChIP-MS preparations and showed decreased cohesin binding at ZBTB21-binding sites after acute depletion of ZBTB21, which is consistent with the results described by Liu et al ., indica ting tha t ZBTB21 may also function as a pioneer and ZBTB7B double-degron mESCs. ZBTB21 was detected with a GFP antibody because it was endo genousl y GFP tagged, ZBTB7B was identified with endogenous antibodies, and GAPDH served as a loading control. ( F ) Western blotting analysis of the chromatin fractions of targeted proteins in the ZBTB21-degron and ZBTB21 and ZBTB7B-degron cell lines under untreated and (6 h) IAA-treated conditions. SMC1 is a cohesin component. SMC4 is a condensin complex sub unit. ␣-Tub ulin is a cytoplasmic marker. Histone H3 is a chromatin marker. ( G ) Snapshots of cohesin ChIP-seq signals at the Rap2a locus in untreated, ZBTB21-depleted and ZBTB21-and ZBTB7B-double-depleted conditions. ( H ) Meta-gene analysis of cohesin ChIP-seq density at ZBTB21-bound genes upon ZBTB21 or ZBTB21 and ZBTB7B degradation. The lines indicate the mean cohesin ChIP-seq signal under each condition.
The boxplot (inset) shows the changes in the ChIP-seq signal at ±100 bp around TSS regions upon ZBTB21 or ZBTB21 and ZBTB7B degradation. Significance was determined by the Wilco x on test (**** P < 0.0001).
factor (similar to OCT4 and SOX2) to facilitate cohesin binding. ZBTB transcription factors have been reported to function in se v eral critical de v elopmental processes and diseases. For example, many ZBTB transcription factors, including ZBTB7A, ZBTB7B and ZBTB35, are required for the prolifera tion and dif ferentia tion of T cells ( 59 , 121 ); ZBTB7A is involved in tumorigenesis and pluripotency control through signaling pathways or chromatin remodelers ( 108 , 122 , 123 ); ZBTB7B dri v es brown fat de v elopment via long non-coding RNAs ( 111 ); ZBTB11 regulates neutrophil de v elopment, and its mutation causes intellectual disability ( 124 , 125 ); and ZBTB21 is associated with the pathogenesis of congenital heart defects in Down syndrome ( 126 ). Their effects on 3D chromatin interactions are probably cell type specific because the distributions of ZBTB ChIP-seq peaks in mESCs and HEK293T cells are different. It would be interesting to deplete these ZBTB factors during de v elopment or disease progression in future research. The molecular mechanisms of ZBTB functions have mostly been investigated in the context of transcriptional regulation. Even though the chromatin interaction analyses were performed within 6 h, these secondary effects on mature RNAs could be limited. We could not eliminate the possibility that misregulated gene expression may also contribute to their effects on cohesin binding and 3D chromatin interactions. ZBTB21 has zinc finger domains that bind DNA, and its BTB domains mediate pr otein −pr otein interactions. Moreover, Drosophila insulator proteins such as Mod (mdg4) and Cp190 are also ZBTB transcription factors and have been shown to function as regulators of insulator activities (52)(53)(54). Previous studies have shown that ZBTB factors are critical for the cell fate determination of specific lineages. Although the effects of ZBTB21 depletion on cohesin chromatin binding are minor, we envisaged that ZBTB21 proteins may function as insulator proteins, similar to previously reported ZBTB factors in Drosophila , which might pre v ent cohesin e xtrusion at critical chromatin regulatory elements. This potential mechanism, together with transcriptional r epr essor activity, would be important for the specific gene expression program during different de v elopmental processes.

DA T A A V AILABILITY
All the data that support the findings of this study are available from the corresponding authors upon reasonable request. Raw sequencing data can be found in the GEO database: GSE184272. Access to the mass spectrometry dataset is available in ProteomeXchange: PXD028860, PXD036004, PXD035893. Raw blot, gel and microscopy image data can be found in Mendeley data: DOI: 10.17632 / cwpb6kmxv7.2