Introduction

The mammalian genome is organized at multiple length scales into compartments, topologically associating domains (TADs) and chromatin loops^1–4. This three-dimensional chromatin folding contributes to essential cellular processes such as transcription and DNA replication and is actively maintained by protein cofactors ^5–8. The architectural protein CTCF establishes the anchors of many DNA loops and TADs likely by blocking loop extrusion by cohesin^7,9–16. Consequently, perturbation of CTCF binding sites at domain boundaries can lead to aberrant cross-domain enhancer-promoter contacts, which have been linked to gene misexpression in cancer and developmental disorders ^17–21. Furthermore, CTCF binding and looping on or near genes can directly enforce enhancer-promoter loops and facilitate differential splicing and polyadenylation ^22–29. Although these examples demonstrate that CTCF can modulate transcription in specific contexts, acute depletion of CTCF only leads to moderate transcriptional defects, highlighting the complex and nuanced roles of CTCF in regulating gene expression ^6,7,30,31.

CTCF binding is sensitive to DNA methylation, and consequently global DNA demethylation ‘reactivates’ hundreds to thousands of previously masked CTCF peaks ^32–36. However, these reactivated CTCF sites comprise only a small subset of all potential methylated CTCF sites, suggesting that DNA methylation serves as a specialized mechanism to regulate specific CTCF binding sites ³⁵. While methylation-sensitive CTCF peaks have been linked to gene expression and mRNA processing in some instances,^{20,21,24–26,37} the functions of most methylation-sensitive CTCF peaks are unknown. Moreover, the broader impacts of reactivated CTCF binding upon global DNA demethylation on chromatin looping have been largely unexplored.

In this study, we employed a selective inhibitor of DNMT1, the maintenance DNA methyltransferase, to induce genome-wide DNA demethylation and investigated the chromatin looping patterns of the resulting reactivated CTCF peaks. We found that reactivated CTCF peaks within genes form loops to highly-looping partner peaks located in active chromatin regions. Remarkably, reactivated CTCF peaks and their highly-looping partners are located close to nuclear speckles, condensate bodies of protein and RNA that have been implicated in transcription and splicing ^38–44. Subnuclear localization relative to nuclear landmarks such as nuclear speckles is linked to genome structure and function, so we explored this relationship further by developing acute degron systems for both CTCF and speckles. Through these experiments, we clarify the directionality of the relationship between CTCF and speckles and glean novel insights into how methylation-sensitive CTCF contributes to chromatin structure and gene expression.

Results

DNMT1 inhibition reactivates CTCF peaks and loops on gene bodies

Loss of DNA methylation leads to increased CTCF occupancy at hundreds to thousands of CTCF sites depending on the cell type³⁵, but a global understanding of their looping patterns and functions is lacking. Recently, selective, noncovalent DNMT1 inhibitors (e.g., GSK3484862) have become available, enabling dynamic control of DNA methylation with fewer undesirable effects such as DNA damage and cytotoxicity ^45,46. Consequently, we employed these new DNMT1 inhibitors to probe the functions and looping patterns of methylation-sensitive CTCF sites.

We treated K562 and HCT116 cells with GSK3484862 (here after referred to as DNMT1i, 10 µM) for three days, and performed CTCF ChIP-seq and HiChIP to investigate CTCF-associated changes in chromatin looping (Fig 1A). Under these conditions, DNMT1i treatment led to minimal growth defects in both cell types (Fig S1A). Consistent with prior studies ³⁵, DNMT1i treatment increased CTCF binding at thousands of sites, previously termed “reactivated peaks” due to their overlap with CTCF peaks in other cell types^35,47 (Fig 1B,C, 2,079 up- and 2 down-regulated in K562, 6,686 up- and 328 down-regulated in HCT116). As expected, reactivated CTCF sites were largely methylated in wild-type cells (Fig S1B), and many CTCF binding motifs enriched within these peak regions were methylated at two key CpG sites implicated in CTCF binding (Fig S1C) ^36,47.

Figure 1

We next evaluated whether these reactivated CTCF peaks form chromatin loops. CTCF HiChIP revealed that a large subset of reactivated peaks resides within at least one loop anchor and these are henceforth called “looping” peaks (see Methods) ^48,49. An example of reactivated looping CTCF peaks is shown at the PCDHA locus (Fig 1D). On average, loops that are anchored by reactivated peaks increase in strength upon DNMT1 inhibition (Fig 1E). We next sought to evaluate these findings using our published Hi-C dataset following DNMT1 inhibition ⁵⁰. Specifically, we analyzed Hi-C interaction strength using standard aggregate pileup analysis over all reactivated CTCF sites identified in our ChIP-seq data. Consistent with our Hi-ChIP data, this analysis revealed an increase in Hi-C interaction signal at reactivated peaks following 3 days of DNTM1 inhibition (Fig 1F).

We next characterized the genomic annotations of looping and nonlooping reactivated CTCF peaks. Consistent with previous findings using azacitidine ³⁵, loss of DNA methylation through DNMT1i treatment led to the preferential reactivation of CTCF peaks on gene bodies, particularly in non-intron regions (exons, promoters, and transcription termination sites (TTS), Fig 1G, S1D, S1E). Notably, this enrichment on gene bodies was particularly pronounced looping reactivated CTCF peaks (Fig 1G, Fig S1D, S1E, S1F), highlighting an underappreciated propensity for CTCF peaks on gene bodies to engage in chromatin loops. Altogether, selective DNMT1 inhibition leads to the emergence of thousands of reactivated CTCF peaks on gene bodies that engage in new chromatin interactions.

Reactivated CTCF peaks interact with highly-looping partners near nuclear speckles

We next scrutinized the CTCF partners that loop to reactivated CTCF peaks. Notably, the CTCF partners of reactivated peaks tend to engage in a greater total number of different loops compared to all CTCF peaks (i.e., have more looping partners, Fig 2A), particularly for K562 cells. This trend is in part driven by the tendency of non-intronic gene body CTCF peaks to more frequently participate in loops with highly-looping CTCF partners in comparison to intronic or intergenic CTCF peaks (Fig S2A). To characterize these highly-looping partners, we ranked all CTCF peaks by the number of different partners they interact with and called those with the most partners ‘highly-looping’, analogously to how superenhancers or supersilencers were previously characterized (Fig 2B, 2C, S2B) ^51,52. HiChIP loops between reactivated peaks and highly-looping peaks increased in strength upon DNMT1 inhibition, consistent with increased interactions between these sites (Fig S2C). The propensity for a given CTCF site to engage in many contacts may reflect architectural stripes, a chromatin structural feature identified from Hi-C data that is postulated to result from persistent loop extrusion at a single region ^53–55. Stripe anchors are enriched for chromatin accessibility and histone modifications associated with active enhancers ⁵³. Similarly, regions near highly-looping CTCF peaks are accessible and enriched for active chromatin marks such as H3K9ac and H3K27ac (Fig 2D). Furthermore, based on previously published Hi-C data ⁵⁰, highly-looping CTCF sites have a stronger propensity to engage in chromatin interactions compared to ordinary CTCF sites (Fig 2E). In addition, stripe anchors identified from these Hi-C data also reside closer to highly-looping CTCF peaks than all peaks, consistent with the similarities between these features (Fig S2D). Thus, these results show that upon DNA demethylation, reactivated CTCF sites form loops to highly-looping partners associated with active regulatory marks and stripe anchors.

Figure 2

To understand the features that distinguish reactivated CTCF peaks and their highly-looping partners from other CTCF peaks, we next considered their subnuclear localization. To do so, we investigated their distribution across genomic subcompartments that exhibit variable transcriptional activity and association with nuclear landmarks, such as the nuclear lamina or nuclear speckles. In particular, we found a striking enrichment of reactivated CTCF peaks and highly-looping CTCF peaks in the A1 subcompartment of the genome, which is differentiated from the A2 subcompartment largely due to its localization near nuclear speckles (Fig S2E) ⁵⁶. Studies have shown that nuclear speckles are hubs of DNA-DNA contacts ^{42,44,57–61}. Additionally, CTCF can form stress-responsive clusters near nuclear speckles in some cell types, and genomic subcompartments associated with speckles have been correlated with increased CTCF binding and looping ^38,41,62,63. However, a direct connection between nuclear speckles and CTCF binding and looping has not been extensively explored.

To determine if reactivated peaks and highly-looping CTCF peaks are preferentially located near nuclear speckles, we performed Cut&Tag for SON, a core speckle protein, in K562 and HCT116 cells with and without DNMT1 inhibition (Fig 2F) ^38,40,58,64. Consistent with localization near speckles, reactivated peaks appeared in regions of higher SON signal compared to all CTCF peaks (Fig 2G, S2F). Notably, this trend is strong for looping reactivated CTCF peaks (Fig 2G). This enrichment of looping reactivated peaks near speckles is further corroborated by publicly available SON TSA-seq data in wildtype cells (Fig S2G) ⁴⁰. Thus, reactivated CTCF peaks reside closer to nuclear speckles than non-reactivated peaks.

Strikingly, relative to all CTCF peaks, highly-looping CTCF peaks are enriched at nuclear speckles (Fig 2H, S2H). More broadly, the number of different loops a CTCF peak engages in and their corresponding strength are positively correlated with SON Cut&Tag signal (Fig 2I, S2I). In addition, we observed that CTCF peaks with strong Cut&Tag SON signal preferentially engage in strong shorter-range interactions in published Hi-C interaction data (Fig 2J) ⁵⁰. Overall, our results reveal that reactivated CTCF peaks interact with highly-looping partners located near nuclear speckles.

DNMT1i-induced gene activation and speckle association depend on CTCF

Our data revealed that CTCF reactivation and looping are highly correlated with nuclear speckle association. As such, we next sought to study the nature of this relationship through degrading either CTCF or nuclear speckles. We first generated a HaloTAG-CTCF homozygous knock-in HCT116 cell line, which enabled the inducible recruitment of CTCF to the VHL E3 ubiquitin ligase upon addition of the small molecule HaloPROTAC3, causing subsequent proteolytic degradation within 8 hours (Fig 3A, 3B, S3A). We then performed SON Cut&Tag after pretreating our CTCF-HaloTag knockin cells with HaloPROTAC3 or vehicle for 8 hours before co-treating with DNMT1i or vehicle for 24 hours (32 hours total for HaloPROTAC3 treatment). This shorter timepoint was used to capture direct effects of CTCF degradation and DNA demethylation as well as avoid toxicity caused by CTCF depletion. We found that reactivated CTCF peaks gain SON Cut&Tag signal upon DNMT1i treatment (Fig 3C). This gain in SON signal is blocked when CTCF is degraded by HaloPROTAC3 treatment, suggesting that speckle association at these reactivated sites is dependent on CTCF (Fig 3C). Thus, SON associates with reactivated CTCF peaks in a CTCF-dependent manner.

Figure 3

Because reactivated peaks fail to engage with nuclear speckles if CTCF is not present, we hypothesized that gene activation induced upon loss of DNA methylation may be dampened in the absence of CTCF. Prior work has suggested that cell-type specific DNA methylation may protect a “labile” set of CTCF peaks from reactivation ³⁵ and that methylation-dependent CTCF occupancy on gene bodies has been linked to changes in gene processing, but a global understanding of the functions of methylation-sensitive CTCF peaks on transcription remains lacking ^24,26,28. Thus, we tested whether reactivated CTCF peaks promote gene activation upon DNMT1 inhibition.

To assess the role of CTCF in enabling DNMT1i-dependent gene activation, we pretreated our CTCF-HaloTag knockin cells with HaloPROTAC3 or vehicle for 8 hours and then co-treated the cells with DNMT1i or vehicle for 24 hours before performing RNA-seq (Fig 3A, 3D). Consistent with prior studies ⁷, acute CTCF degradation altered the expression of hundreds of genes with approximately equal numbers up- and down-regulated (488 up-regulated, 338 down-regulated, Padj <0.05, |log₂fold-change| > 0.5) (Fig S3B), with downregulated genes possessing higher baseline levels of CTCF at their promoters (Fig S3C). As expected, DNMT1i treatment in this cell line for 24 hours primarily resulted in gene upregulation, with 135 genes upregulated and 8 genes downregulated (Padj <0.05, |log₂fold-change| > 0.5) (Fig S3B).

We focused our analysis on genes that were differentially expressed in DNMT1i versus DMSO treatment in the presence of CTCF (Fig S3B, left). Using K-means clustering, we identified four gene clusters within this group defined by unique z-score normalized expression patterns across the conditions (Fig 3D). Clusters 1, 2, and 4 comprise genes that increase in expression upon DNMT1 inhibition, whereas Cluster 3 contains the few genes that decrease in expression upon DNMT1 inhibition. Moreover, genes in Cluster 1 are upregulated by CTCF depletion, whereas genes within Clusters 2, 3, and 4 are unaffected by CTCF depletion alone.

We focused on Clusters 2 and 4, which contain lowly-expressed genes that are upregulated upon DNMT1 inhibition and unaffected by CTCF depletion in the absence of DNMT1i (Fig 3D, S3D). Significantly, while genes in Cluster 4 are upregulated upon DNMT1 inhibition with or without CTCF present, genes in Cluster 2 depend on CTCF for DNMT1i-mediated upregulation. To determine how CTCF binding might be related to the observed changes in gene expression, we performed CTCF ChIP-seq in our HaloTag-CTCF HCT116 cells following 24 hours of DNMT1i or vehicle treatment. Notably, genes in Cluster 2 preferentially gained CTCF occupancy upon DNMT1 inhibition compared to genes in other clusters (Fig 3D, right, S3E). Thus, these data demonstrate that a subset of genes enriched for CTCF reactivation depend on CTCF for expression upon DNMT1i treatment.

We next assessed the correlation between CTCF-dependent genes and speckle association. Intriguingly, genes in Cluster 2 have higher SON TSA-seq levels than genes in Cluster 4 (Fig 3E), while genes in cluster 1, which are upregulated by CTCF depletion alone, are closest to nuclear speckles. Moreover, genes in cluster 2 gain a modest level of SON Cut&Tag signal upon DNMT1i treatment (Fig S3F). Taken together, these experiments show that CTCF is necessary for both gene activation and speckle protein association upon DNMT1 inhibition.

Acute disruption of nuclear speckles alters gene expression without disrupting CTCF

To further explore the observed connection between speckles and CTCF looping, we tested whether nuclear speckles play a causal role in promoting CTCF occupancy and/or looping. To accomplish this, we sought to acutely degrade speckles to evaluate whether speckles actively maintain CTCF loops while mitigating confounding factors associated with mitotic arrest caused by long-term SON depletion ⁶⁵. Consequently, we engineered a K562 knock-in cell line with the core speckle proteins, SON and SRRM2, both endogenously tagged with FKBP12^F36V for use with the dTAG degradation system ^58,66. Upon addition of dTAG-13, FKBP12^F36V-tagged proteins are recruited to the CRBN-Cul4-Ring E3 ligase and undergo acute proteasomal degradation (Fig 4A). Treatment of the K562 FKBP12^F36V double knock-in cell line with dTAG-13 resulted in near-complete degradation of both SON and SRRM2 within 6 hours (Fig 4B, S4A).

Figure 4

We next assessed if speckle disruption directly alters gene expression by performing RNA-seq after 6 hours of dTAG-13 treatment. Upon speckle degradation, 541 genes were down-regulated and 184 genes were up-regulated (padj < 0.05, |log₂fold-change| > 0.25) (Fig S4B). Consistent with a role of nuclear speckles in mediating active transcription, a large majority of down-regulated genes were located within the closest decile to nuclear speckles by TSA-seq (Fig S4C), and genes with the highest SON Cut&Tag signal exhibited the strongest decrease in expression (Fig S4D). These results show that disrupting nuclear speckles has immediate effects on RNA abundance.

To determine whether speckle disruption alters CTCF binding and/or looping, we performed CTCF ChIP-seq and HiChIP in the speckle dTAG cell line after treatment with dTAG-13 or DMSO for 6 hours. In contrast to the immediate transcriptional effects, CTCF occupancy remained unchanged upon speckle degradation, with only 2 significantly downregulated peaks (Fig 4C, S4E). CTCF HiChIP also revealed that loop strength was largely unchanged across speckle Cut&Tag and TSA-seq quantiles (Fig 4D, S4F) despite a positive correlation between speckle association and baseline loop strength (Fig 2I). Even after 12 hours of dTAG treatment, there was only a very subtle decrease in CTCF looping near speckles (Fig S4F). In addition, changes in gene expression were not associated with changes in loop strength or CTCF binding (Figure S4G, H). Conversely, depleting CTCF for 24 hours in our HaloTag-CTCF HCT116 knockin cell line did not noticeably change speckle morphology (Fig S4I). Altogether, acute disruption of nuclear speckles did not appreciably disrupt CTCF binding or looping, suggesting that maintenance of 3D CTCF loop structure is independent of speckle association. This complements our finding that conversely, speckle association depends upon reactivated CTCF.

Discussion

DNA methylation accounts for much of the variation in CTCF binding across cell types ³⁵, but the functions and features of methylation-sensitive CTCF peaks are incompletely understood. Here we leveraged a selective DNMT1 inhibitor ^45,46 to gain new insights into how reactivated CTCF peaks interface with chromatin looping and transcription. This new inhibitor is a valuable tool for reactivating CTCF without the confounding effects caused by the covalent nature and toxicity of earlier hypomethylating agents.

We found that reactivated CTCF peaks preferentially form loops on gene bodies and interact with promiscuous partner peaks at stripe anchors. Architectural stripes may contribute to transcription by promoting enhancer contacts with gene bodies ^53–55,67, and may also represent a mechanism to poise transcription ⁵⁵. Thus, we propose that reactivated peaks enable transcriptional activation by interacting with regulatory elements through highly-looping CTCF partner peaks (Fig 4E). However, DNA demethylation causes pleiotropic effects ranging from changes in gene expression to transcription factor binding ^45,68, and we cannot rule out that these other features could contribute to the effects we observe.

Strikingly, reactivated CTCF peaks and their highly-looping partners are located near nuclear speckles, consistent with prior studies showing that subcompartments near nuclear speckles exhibit increased CTCF binding and looping ^41,62. We first explored this relationship by degrading CTCF and intriguingly observed that reactivated CTCF peaks gain SON Cut&Tag signal upon DNMT1 inhibition in a CTCF-dependent manner, indicating a functional tie between reactivated CTCF and speckles. This is consistent with a recent preprint reporting that CTCF helps establish basal speckle-gene contacts⁶⁹. Additionally, our CTCF degron experiments demonstrate that a large fraction of genes that are differentially expressed upon DNMT1 inhibition depends on CTCF for up-regulation. These genes preferentially gain CTCF binding and are located near nuclear speckles. Consequently, CTCF peaks that are reactivated upon DNMT1 inhibition could facilitate the activation of genes near speckles, perhaps by mediating the engagement of genic CTCF sites with speckle-associated stripe anchors and/or splicing factors, although this model requires further study (Fig 4E). A role of CTCF in facilitating chromatin loops near speckles is synergistic with proposed models suggesting that intragenic looping contributes to gene splicing ^25,70. Thus, reactivated CTCF peaks may reorganize genes near nuclear speckles as a means of influencing transcription.

We next studied whether speckles impact CTCF binding or looping by achieving the first acute degradation of nuclear speckles. Previous work knocking down SRRM2 support a causative role of speckles in modulating chromatin structure but did not specifically study changes in CTCF binding or looping ⁷¹. Although acute disruption of speckles supports their direct involvement in transcription, we observed minimal impacts on CTCF binding or looping. This notion is consistent with recent work showing that CTCF clustering and looping is not dependent on phase separation or nuclear bodies ^72,73 and that transcriptional disruption does not immediately impact chromatin looping or CTCF binding ^53,74,75. Our findings also suggest that the genome interactions disrupted upon SRRM2 depletion ⁷¹ are not mediated by CTCF or only occur after a longer time. However, it is possible that nuclear speckles contribute to the establishment, but not maintenance, of CTCF loops, which may not be evident in our system due to the short time points employed for speckle degradation. Moreover, other features enriched near speckles, such as open chromatin, high GC content, and active gene expression, could instead contribute to increased CTCF binding and looping near speckles. Overall, our speckle degradation experiments help clarify the directionality of the relationship between speckles and CTCF.

Taken altogether, our findings highlight a potential mechanism through which DNA methylation modulates genome architecture — by protecting a subset of genic CTCF sites from engaging with promiscuous CTCF peaks located near nuclear speckles. In particular, genes near nuclear speckles may be susceptible to activation through CTCF-dependent mechanisms, and thus might require this layer of protective regulation. Furthermore, the GC-rich composition of non-intron gene bodies might increase the likelihood of CTCF motifs in comparison to other genomic regions, increasing the vulnerability of these sites to ectopic CTCF binding. More broadly, our work nominates new functions for the known antagonism between gene body methylation and CTCF occupancy in regulating transcription through remodeling genome architecture.

Methods

Cell culture

K562 was obtained from ATCC. HEK 293 T was obtained from Life Technologies. HCT116 was obtained as a gift from M. Shair. Cells were authenticated by Short Tandem Repeat profiling (Genetica) and tested for mycoplasma (Sigma-Aldrich, MP0035). All cell lines were cultured in a humidified 5% CO₂ incubator at 37 °C with media supplemented with 100 U/mL penicillin and 100 µg/mL streptomycin (Life Technologies, 15070063). K562 was cultured in RPMI-1640 (Life Technologies, 11875119) with 10% fetal bovine serum (FBS) (Peak Serum, PS-FB2). HEK 293 T was cultured in DMEM (Life Technologies, 1195073) with 10% FBS (Peak Serum, PS-FB2). HCT116 was cultured in McCoy’s 5A media (Thermo Fisher Scientific, 16600082) with 10% FBS (Peak Serum, PS-FB2).

Generating SON, SRRM2, and CTCF knock-in cell lines

To generate an N-terminal Puro-P2A-dtag SON knock-in donor plasmid, SON homology arms of 1 kb were amplified from K562 genomic DNA and cloned into the ph1 backbone Puro-P2A-dTAG, sub-cloned from pCRIS-PITChv2-Puro-dTAG, a gift from J. Bradner and B. Nabet (Addgene # 91793). The SON sgRNA, 5’-CACCGCTGCTCGATGTTGGTCGCCA-3’ that was previously used for N-terminal SON knock-in was cloned into PX459 by restriction enzyme digestion (Allen Institute,https://www.allencell.org/cell-catalog.html). To generate a C-terminal dTAG-P2A-blasticidin SRRM2 knock-in donor plasmid, SRRM2 homology arms of 1 kb were amplified from K562 genomic DNA and cloned into the ph1 backbone dTAG-P2A-BSD, sub-cloned from pCRIS-PITChv2-dTAG-BSD, a gift from J. Bradner and B. Nabet (Addgene # 91795). The sgRNA, 5’CACCGCCATGAGACACCGCTCCTCC-3’ that was previously validated for C-terminal SRRM2 knock-in^\was cloned into PX459 by restriction enzyme digestion ⁵⁸. Notably, this strategy eliminates the last four amino acids of SRRM2, which is necessary due to a sparsity of guides in the region and has been used to study SRRM2 in the context of speckles ⁵⁸. Cells were transfected with both the donor and guide plasmids by nucleofection (Lonza, program T016). Cells were selected 3 days after transfection with puromycin (2 µg/mL) (Thermo Fisher Scientific, A1113803) or blasticidin (8 µg/mL) (Thermo Fisher Scientific, A1113903) for a duration of 5-10 days. Single-cell clones were isolated from the bulk cell line through fluorescence-activated cell sorting. The SON knock-in was generated from wildtype K562 cells, and then the SRRM2 knock-in was introduced on top of the validated SON knock-in clone. To generate the HT-CTCF HCT116 knock-in cell line, the sgRNA 5’-CACCGTGCAGTCGAAGCCATTGTGG-3’ was subcloned into PX459 and HaloTAG-CTCF was subcloned into PDONR223. The guide and donor plasmids were transfected into cells using Fugene HD (Promega). After seven days, single cell clones were isolated from the bulk population through fluorescence-activated cell sorting using HaloTag TAMRA labeling to sort potential positive clones.

Immunoblotting

Pellets were harvested for western blot analysis and washed with PBS (Corning, 21-040-CV). For immunoblots of CTCF, cells were lysed on ice using RIPA buffer (Boston BioProducts) supplemented with fresh HALT Protease Inhibitor (Thermo Fisher Scientific) and the lysates were cleared through centrifugation. The protein concentration of the lysate was determined using the BCA Protein Assay Kit (Thermo Fisher Scientific) and immunoblots for CTCF were performed according to standard procedures with a 6% SureCast acrylamide page gel. For immunoblots of SON and SRRM2 and corresponding loading controls, 1 million cells per condition were lysed on ice for 30 minutes in 50uL RIPA buffer (Boston Bioproducts, BP-115) supplemented with 1mM PMSF, 2mM EDTA, and 1X Halt Protease inhibitors (Thermo Fisher Scientific, 78429). Subsequently, 50uL of wash buffer was added (50mM Tris pH 8.0, 2mM MgCl, 50mM NaCl, 1mM phenylmethylsulfonyl fluoride (PMSF), 1X Halt Protease inhibitors, 1:500 Benzonase (Millipore-Sigma, E1014-25KU) and the samples were rotated at room temperature for 20 minutes and then clarified by full speed centrifugation at 4 °C. The protein concentration of the lysates was determined using the BCA Protein Assay Kit (Thermo Fisher Scientific, 23225). Immunoblotting was performed according to standard procedures with the following modifications: a 3-8% tris acetate gel (Thermo Fisher Scientific, EA03752) was used and samples were run for 2 hours 15min at 80V at 4 °C in 1X tris-acetate running buffer (Thermo Fisher Scientific, LA0041) except for the loading control which was run for 90 minutes. Transfer was conducted in 1X transfer buffer (Research Product International, 501977679) containing 10% methanol and 0.05% sodium dodecyl sulfate (SDS) onto an activated 0.45um PVDF membrane (Millipore sigma IPVH09120) for 2 hours at 15V except for loading control which was transferred in 1X transfer buffer with 20% methanol at 12V for 90 minutes. The primary antibodies used for immunoblotting are as follows: CTCF (1:2000, Cell Signaling Technology D31H2); SON (1:1000, Abcam ab121759); SRRM2 (1:1000, ThermoFisher PA5-66827); GAPDH (1:2000, Santa Cruz Biotechnology sc-47724); Vinculin (1:2000, Cell Signaling Technology 13901S).

Immunofluorescence

Immunofluoresence was performed according to standard protocol. Briefly, K562 speckle knock-in cells treated with 500 nM dTAG-13 or DMSO for the indicated amount of time were cytospun onto glass coverslips for 5 minutes at 800 x rcf and fixed in 4% PFA for 15 minutes. HCT116 cells treated with 330 nM HaloPROTAC3 or DMSO were instead grown directly on glass coverslips and fixed in 4% PFA for 15 minutes. Cells were then washed 3X in PBS, permeabilized with 0.25% Triton-X for 10 minutes, and subsequently washed 3X with PBS. Cells were then blocked in 10% goat serum (Millipore Sigma, G9023) in PBS-T for 1 hour and incubated overnight in 5% goat serum in PBS-T with primary antibody at 4°C (SON (1:300, abcam ab121759) or SRRM2 (1:300, ThermoFisher PA5-66827)). Following 3X PBS-T washes, cells were incubated with secondary antibody for 1 hr (1:2000, Invitrgoen alexafluor488 A-11008), washed 3X in PBS-T, and mounted on a coverslip with DAPI (Prolong Diamond Anti-Fade, P36961). Cells were imaged on an LSM180 Confocal microscope. Image quantification was performed using CellProfiler on raw Tiff images to identify nuclei and then quantify average fluorescence signal within each nucleus⁷⁶.

Cell growth assays

K562 or HCT116 cells were plated in triplicate in a 96-well plate. Drug or vehicle (DMSO, Sigma Aldrich, D8418) were dosed at the specified concentration for 3 days. After 3 days of treatment, cell viability was measured using CellTiter-Glo (Promega, G7572) with the luminescence detector on the SpectraMax i3x (Molecular Devices) plate reader with the SoftMax Pro (version 6.5.1) software, and data were processed using Prism.

ChIP-seq

CTCF ChIP-seq was performed in duplicate. For the CTCF reactivation experiments, K562 or HCT116 wildtype or HT-CTCF HCT116 knock-in cells were treated with 10 µM GSK3482364 or 0.1% DMSO control for the indicated time. For the speckle depletion experiments, K562 speckle dTAG knock-in cells were treated with 500 nM dTAG-13 or 0.1% DMSO control for 6 or 12 hours. After drug treatment, cells were resuspended and fixed at a density of 10⁶ cells/mL in media containing 10% FBS (Peak Serum, PS-FB2) and 1% formaldehyde (Sigma Aldrich, F8775) for 10 min with rotation and quenched with 0.2 M glycine and washed in PBS. 5 million cells per condition were lysed on ice for 10 min in 600 uL of ChIP lysis buffer (50 mM Tris-HCl pH 7.5, 1% SDS, 0.25% NaDOC) with HALT protease inhibitors (Thermo Fisher Scientific, 78429). Cells were then diluted to 2mL with ChIP dilution buffer (50 mM Tris-HCL pH 7.5, 0.01% SDS, 150 mM NCl, 0.25% triton X-100, 1X HALT protease inhibitor) and sonicated at 4 °C with a Branson Sonicator (45% amplitude, 0.7s on, 1.3s off, total time 5 min). The supernatant was clarified by centrifugation at 16000 x rcf for 10 min at 4 °C. The sample was then diluted into 6 mL total with ChIP Dilution Buffer and the final concentration of triton was brought to 1%. 50 µL was removed as an input sample and kept at 4 °C overnight. Immunoprecipitations were carried out on the remaining sample overnight at 4 °C with rotation with 5 uL of antibody (Cell Signaling Technology D31H2 anti-CTCF, 3418S). 50uL of Protein G Dynabeads (Thermo Fisher Scientific, 10004D) were added to each sample and incubated at 4 °C with rotation for 3 hours. Supernatant was discarded on a magnet and beads were then washed twice with ice-cold RIPA wash buffer (10 mM Tris-HCl pH 8.1, 0.1% SDS, 150 mM NaCl, 0.1 NaDOC, 1% Triton X-100, 1 mM EDTA), once with ice-cold LiCl wash buffer (10 mM Tris-HCl, pH 8.1, 250 mM LiCl, 0.5% Triton X-100, 0.5% NaDOC), and once with ice cold TE buffer (10mM Tris-HCL pH 8.0, 50mM NaCl, 1mM EDTA). Chromatin was eluted in 100 µL of Elution buffer (10 mM Tris-HCl pH 8.0, 0.1% SDS, 150 mM NaCl, 1 mM EDTA) at 65 °C for 1 hour at 1000rpm. Input samples were brought to 100µL with Elution buffer and incubated alongside ChIP samples. Supernatant was removed and 50 µg of RNAse A (Sigma Aldrich, 10109142001) was added and incubated for 30 min at 37 °C. Next, 25 µg of Proteinase K (Invitrogen, 25530015) was added, and the samples were incubated 63 °C for 3 hours with 1000rpm agitation. DNA was subsequently purified by a 2X SPRI (Omega Bio-Tek, M1378-01). Sequencing libraries were then prepared from 5 ng of DNA by performing end-repair with the End-it DNA End-Repair Kit (Lucigen, ER81050), followed by A tailing with NEB Klenow Fragment (3’−5’ exo-) (New England Biolabs, M0212), adapter ligation with NEB DNA Quick Ligase (New England Biolabs, M2200), and PCR using NEB ultra mastermix (M0544S) and 10X KAPA Library Amp Primer Mix (Kapa Biosystems, KK2623). Samples were sequenced using a NovaSeq SP kit (Illumina) with 50 bp paired-end reads at a sequencing depth of approximately 25 million reads per sample.

HiChIP

Cells were treated in duplicate. For the DNMT1i experiment, K562 or HCT116 cells were treated with 10 µM GSK3482364 (ChemieTek) or 0.1% DMSO control for 3 days. For the speckle degron experiment, K562 speckle dTAG knock-in cells were treated with 500 nM dTAG-13 (Sigma Aldrich, SML2601) (or 0.1% DMSO control for 6 or 12 hours. CTCF HiChIP was performed on each sample according to a published protocol ⁷⁷. Specifically, cells were resuspended and fixed at a density of 10⁶ cell/mL in media containing 10% FBS (Peak Serum, PS-FB2) and 1% formaldehyde (Sigma Aldrich, F8775) for 10 min with rotation and quenched with 0.2 M glycine and washed in PBS. 10 million cells per condition were resuspended in ice-cold Hi-C Lysis Buffer (10mM Tris HCl pH 8.0, 10mM NaCl, 0.2% NP-40, 1X protease inhibitor) and lysed at 4 °C for 30 minutes. Tubes were centrifuged and the pellets were washed with 500 µL ice-cold Hi-C lysis buffer. Pellets were then resuspended in 100µL of 0.5% SDS and incubated at 62 °C for 10 minutes, after which they were quenched with Triton X-100 and incubated at 37 °C for 15 minutes. The samples were then digested with 375 units of DpnII (New England Biolabs, R0543L) for 2 hours at 37 °C with shaking at 900 rpm and then enzyme was heat inactivated at 65 °C for 20 minutes. Digested ends were filled in with biotinylated dATP (Life Technologies, 19524-016) for 1 hour at 37 °C at 900rpm, and in situ ligation was performed with T4 DNA ligase (New England Biolabs, M0202L) for 4 hours at room temperature with rotation. Nuclei were pelleted by centrifugation and resuspended in Nuclear Lysis Buffer (50mM Tris-HCL pH 7.5, 10mM EDTA, 1% SDS, 1X protease inhibitor) and sheared in Covaris Millitubes (cat # 520135) on a Covaris S220 sonicator with Fill Level 10, Duty Cycle 5, PIP 140, Cycles/Burst 200, for 4 minutes each. Samples were clarified by centrifugation and supernatant was split into two tubes and diluted with 2X the volume of ChIP Dilution Buffer (0.01% SDS, 1.1% triton X 100, 1.2mM EDTA, 16.7mM tris pH 7.5, 167mM NaCl). 30ul Pre-washed Protein G Dynabeads (Thermo Fisher Scientific, 10004D) were then added to each tube and samples were rotated for 1 hour at 4 °C. Supernatant was transferred to new tubes and antibody was added and allowed to bind overnight at 4 °C with rotation (15uL antibody per tube, Cell Signaling Technology D31H2 anti-CTCF, 3418S). 30uL of pre-washed Protein G Dynabeads were added to each tube and rotated at 4 °C for 2 hours. Beads were washed 3X each with Low Salt Wash Buffer (0.1% SDS, 1% triton X 100, 2 mM EDTA, 20 mM Tris HCL pH 7.5, 150 mM NaCl), High Salt Wash Buffer (0.1% SDS, 1% triton X 100, 2 mM EDTA, 20mM Tris HCl pH 7.5, 500 mM NaCl), and LiCl Wash Buffer (10mM tris pH 7.5, 250 mM LiCl, 1% NP-40, 1%Na-DOC, 1mM EDTA) at room temperature. Tubes for each sample were recombined and chromatin was eluted into 100 µL of Elution Buffer (50mM NaHCO₃, 1% SDS) per sample by incubating at room temperature for 10 minutes with rotation followed by 3 minutes at 37 °C shaking. Elutions were repeated for 200 µL total volume. 10 µL of Proteinase K were added per sample (Invitrogen, 25530015 10mg/mL) and the samples were incubated at 55 °C for 45 minutes at 900rpm followed by 67 °C for 2 hours at 900rpm. Samples were purified (Zymo #D4014) and eluted in 10 µL water and quantified by dsDNA qubit quantification. Biotinylated DNA was then enriched with 5 µL pre-washed Streptavidin C1 beads per sample (Life Technologies, 65001) in biotin binding buffer (5mM Tris-HCL pH 7.5, 1mM EDTA, 1M NaCl) for 15 minutes at room temperature with rotation. Beads were washed twice with Tween Wash Buffer (5mM Tris-HCL pH 7.5, 0.5mM EDTA, 1M NaCl, 0.05% tween 100) at 55 °C for 2 minutes at 900rpm and were then washed with TD buffer (10mM Tris-HCL pH 7.5, 5mM MgCl2, 10% Dimethylformamide). Beads were resuspended in TD buffer and Tn5 (Illumina, 20034198) was added at the suggested dilutions and tagmentation was allowed to occur at 55 °C for 10 minutes with 900 rpm shaking. Supernatant was removed and beads were washed in 50 mM EDTA for 30 minutes at 50 °C and twice more for 3 minutes each. Beads were then washed twice in Tween Wash Buffer at 55 °C for 2 minutes and once in 10mM Tris. DNA was amplified by PCR with Phusion HF 2X MM and Nextera universal barcodes described previously and was purified by two-sided SPRI purification (Omega Bio-Tek, M1378-01). Samples were sequenced using a NovaSeq SP kit (Illumina) with 100bp or 50 bp paired-end reads at a sequencing depth of approximately 100 million reads per sample.

RNA-seq

In triplicate separate cultures, cells were treated with the indicated compounds. For the CTCF degron experiment, HaloTag-CTCF HCT116 knock-in cells were pretreated with 330 nM HaloPROTAC3 (AOBIOUS AOB36136) or 0.1% DMSO control for 8 hours, and then were treated with 10µM GSK3482364 or 0.1% DMSO in addition to +/-330 nM HaloPROTAC3 for 24 more hours. For the speckle knockdown experiment, K562 speckle dTAG knock-in cells were treated with 500 nM dTAG-13 or 0.1% DMSO control for 6 hours. Total RNA was isolated from 1 million cells per condition using the RNeasy Plus Mini Kit (Qiagen, 74104). mRNA was isolated using the Poly(A) RNA Selection Kit V1.5 (Lexogen, 157.96) and total RNA seq libraries were prepared using the CORALL Total RNA-Seq Library Prep Kit (Lexogen, 095.24). Samples were sequenced using a NovaSeq SP kit (Illumina) 50 bp paired-end reads at a sequencing depth of approximately 30 million reads per sample for the speckle degron experiments. For the HaloTAG CTCF degron experiments, samples were instead sequenced to a depth of approximately 15 million reads per sample with 100bp paired end reads.

SON Cut&Tag

For the wildtype cell line experiments, K562 or HCT116 wildtype cells were treated in duplicate with either 10 uM GSK3482364 or 0.1% DMSO control for 3 days. For the CTCF degron experiment, HaloTag-CTCF HCT116 knock-in cells were pretreated with 330 nM HaloPROTAC3 (AOBIOUS AOB36136) or 0.1% DMSO control for 8 hours, and then were treated with 10 µM GSK3482364 or 0.1% DMSO in addition to +/− 330nM HaloPROTAC3 for 24 more hours. Cut&Tag was performed according to published protocol ⁶⁴. 1 million cells were harvested and washed in PBS. Cells were resuspended in ½ volume ice-cold NE1 buffer (20 mM HEPES-KOH pH 7.9, 10 mM KCl, 0.5mM spermidine, 0.1% Triton-X-100, 20% glycerol with Roche Complete Protease Inhibitor EDTA-free) and lysed on ice for 10 minutes. Nuclei were resuspended in PBS and fixed with 0.1% formaldehyde (Sigma Aldrich, F8775) for 2 minutes at room temperature and quenched with 75 mM glycine. Nuclei were centrifuged, drained, and resuspended in wash buffer (20 mM HEPES-NaOH pH 7.5, 150 mM NaCl, 0.5 mM spermidine, Roche Complete Protease Inhibitor EDTA-free) to ∼1 million cells per mL, and 250,000 nuclei per condition (as determined by hematocytometer) were used per condition moving forward. Concavalin A beads (Cell Signaling Technologies, 93569S) were then activated with binding buffer (20 mM HEPES-KOH pH 7.9, 10 mM KCl, 1 mM CaCl₂, 1 mM MnCl₂) and 10 µL were added per sample and incubated for 10 minutes at room temperature with rotation. Liquid was removed on a magnet and beads were resuspended in 50uL ice-cold antibody buffer (wash buffer with 2 mM EDTA and 0.1% BSA) and 1 µL of antibody (anti-SON Abcam #ab121759) was added and incubated upright with gentle nutating overnight at 4 °C. Beads were resuspended in secondary antibody in wash buffer (1:100 guinea pig anti-rabbit ABIN101961) and tubes were nutated at room temperature for 1 hour. Liquid was removed and beads were washed three times with wash buffer. pA-Tn5 (purified in-house according to published protocol ⁷⁸) was mixed 1:50 in Dig-300 buffer (without digitonin, 20 mM HEPES-NaOH pH 7.5, 300 mM NaCl, 0.5mM spermidine, Roche Complete Protease Inhibitor EDTA-free) and beads were resuspended in 100 uL and nutated at room temperature for 1 hour. Beads were washed three times in Dig-300 buffer. Beads were resuspended in 300 uL Tagmentation Buffer (Dig-300 buffer with 10 mM MgCl₂) and incubated at 37 °C for one hour. 16 mM EDTA, 0.1% SDS, and 50 µg of Proteinase K (Invitrogen, 25530015) was added per sample and incubated at 55 °C for one hour and then DNA was purified by PCI-chloroform extraction. DNA was amplified by PCR with NEBNext HiFi 2x PCR Master Mix (NEB M0541S) and library was purified by 1.3X SPRI purification (Omega Bio-Tek, M1378-01). Samples were sequenced using a NextSeq1000 (Illumina) with 50bp paired end reads at a sequencing depth of approximately 10-20 million reads per sample.

CTCF ChIP-seq data processing and peak calling

Reads were aligned to hg38 using Burrow-Wheeler Aligner (BWA-mem, v 0.7.15)⁷⁹. Sam files were converted to bam files and sorted using Samtools view and sort (v 1.5)⁸⁰. Identical ChIP-seq sequence reads were collapsed to avoid PCR duplicates using Picard Tools MarkDuplicates(v 2.18.5) with Java (v 1.8.0) with parameters VALIDATION_STRINGENCY=LENIENTASSUME_SORTED=true REMOVE_DUPLICATES=true. Reads were converted to tdf files for viewing on IGV using Igvtools (v 2.3.95) count with -w 25⁸¹. Tag directories were made from deduplicated bam files with Homer (v. 4.9) makeTagDirectory with Perl (v 5.26.1) ⁸². Peaks were called using Homer (v. 4.11) getDifferentialPeaksReplicates where -t was replicate ChIP tag directories, -i was the corresponding replicate input directories, with parameters -F 10 and -L 4. Reactivated peaks were called using Homer (v. 4.11) getDifferentialPeaksReplicates where -t was replicate ChIP tag directories for DNMT1i treatment, -b was the replicate ChIP directories for DMSO treatment, and -p was the peaks called in the corresponding DNMT1i treatment condition. These peaks were further restricted to those only called in DNMT1i and not DMSO by using bedtools (v 2.26.0) intersect with the -v and -wa options ⁸³. For aggregate analyses, DeepTools was used to convert deduplicated bam files to RPKM normalized bigwigs and to compare bigwigs to one another.

HiChIP data preprocessing and loop calling

HiChIP reads for each sample were aligned to hg38 and processed with HiC Pro (v 3.0.0) using the default settings with a DpnII digested genome generated with the digest_genome script ⁸⁴. For most analyses, HiCPro output .allvalidpairs files for replicate samples were then pooled. Loops were called using Hichipper (v 0.7.7) in call mode with default parameters and the flag --input-vi, and loops for each condition were kept if they had FDR< 0.01 and > 4 PETs⁴⁹. For downstream foldchange and CPM analysis using Diffloop, Hichipper was instead used to call loops from HiCPro outputs without pooling replicates and the unfiltered loop list was used as input to Diffloop. For all Hichipper analyses, merged peak files containing all peaks from DMSO and DNMT1i for a given cell line were used as the input peak file.

Hi-C analysis

Published LIMe-Hi-C data for three replicates of DMSO and DNMT1 inhibitor treatment conditions in K562 were employed for the analysis ⁵⁰. The same treatment conditions were used as employed in this study. LIMe-Hi-C samples were generated as previously described⁵⁰ using JuiceMe (v 1.0.0) to generate merged_nodups.txt.gz pairwise contact files. Contacts from all three replicates were subsequently pooled together using the JuiceMe CPU “mega” script (v 1.0.0). These files were then converted to the cool format at 10 kb resolution using cooler (v 0.8.11). These cool files were normalized using the cooler (v 0.8.11) “balance” function and were directly utilized for aggregate pileup analysis. To perform pileup analysis, CTCF peaks, identified as described above, were filtered based upon the presence of a CTCF motif using the motif file MA0139.1.tsv.gz downloaded from the JASPAR database and assigned to the motif with the strongest signal using the bioframe (v 0.3.3) “overlap” function. In addition, CTCF sites that were located within 1000 bp of a previously identified blacklisted region (ENCFF356LFX.bed.gz) were removed using the bioframe “subtract” function. Pileup analysis was then performed using the cooltools (v 0.5.1) “pileup” function on peaks deduplicated through the bioframe “cluster” function with a “min_dist” threshold set to a resolution of 10 kb to remove duplicate peaks located within the same genomic bin⁸⁵. Peaks for which motifs mapped to the – strand were reflected across the diagonal and the pileup matrix was then plotted in python (v 3.7) using standard matplotlib (v 3.5.2) functions. Highly-looping and reactivated peaks were defined as specified below. Speckle associated peaks were defined to be those with an average speckle SON Cut&Tag signal greater than or equal to the 90^th percentile genome-wide 20 kb bin average speckle signal.

The stripenn (v 1.1.65.7)⁵⁴ “compute” function was used to identify stripes from the merged DNMT1 inhibitor 10 kb cool file with the following parameters: --norm KR -m 0.95,0.96,0.97,0.98,0.99. The distances of highly-looping peaks and normal-looping peaks to stripe anchors tabulated in the stripenn “result_filtered.tsv” output were calculated using pybedtools (v 0.8.1) “closest” function and visualized in python (v 3.8) using standard plotting functions.

CTCF binding scatterplots and aggregate profile plots

To make scatterplots of CTCF normalized tags, RPKM normalized bigwigs were generated for ChIP samples and matched input controls from deduplicated bam files (see above) using deepTools (v 3.4.3) bamCoverage. ChIP files were normalized to their matched input controls using deepTools (v 3.4.3) bigwigCompare with the log2 option. A merged peak list was made by combining peak files from DNMT1i and DMSO and collapsing the list with bedtools (v 2.26.0) merge, and normalized ChIP-seq signal per condition was annotated onto these peaks using deepTools multiBigwigSummary. To make aggregate profile plots of CTCF ChIP-seq signal over genes or peak subsets, bigwig files of CTCF ChIP-seq signal were generated for each condition from deduplicated bam files (see above) with deepTools (v 3.4.3) bamCoverage with flag --normalizeUsing RPKM. To make aggregate plots over peak subsets, deepTools (v 3.4.3) computeMatrix and plotProfile were used with the following computeMatrix parameters: reference-point --referencePoint center -bs 50 -a 1000 -b 1000 --skipZeros. To make aggregate plots over genes, deepTools (v 3.4.3) computeMatrix and plotprofile were used with the following parameters: --scale-regions --regionBodyLength 1000-bs 50 -a 1000 -b 1000 --skipZeros. Downstream processing was completed in python (v 3.7.0) (www.python.org) with pandas (v 1.0.3), matplotlib (v 3.2.1) and seaborn (v 0.11.2).

Annotating CTCF peaks

CTCF peaks were annotated using Homer (v 4.11) annotatePeaks using the -gtf option with an Ensembl gtf file (v 104) with the standard priority order assigned ⁸². Motifs were assigned to CTCF peaks using gimmemotifs gimme scan function with -b and -f 0.1 and the JASPAR motif pfm as input. Strand-aware motif methylation heatmaps were made using Deeptools. To assess overlap of categories of CTCF peaks with a bed file of predefined subcompartments from SNIPER ⁵⁶ bedtools (v 2.26.0) map was used with the – o distinct flag, and data were further processed in python (v 3.7.0) (www.python.org) and any peaks overlapping multiple subcompartments were discarded. To assess overlap of categories of CTCF peaks with normalized TSA-seq counts, UCSC’s function bigWigToBedGraph was first used to convert a bigwig file of normalized SON TSA-seq counts to a bedgraph in 20kb bins, and then bedtools (v 2.26.0) map was used with – o mean flag to map TSA-seq counts onto CTCF peak list and the data were further processed in python (v 3.7.0). The following published files were used from 4D Nucleome for TSA-seq figures: 4DNFIVZSO9RI.bw 4DNFIBY8G6RZ.bw.

Assigning number of loops (partners) per CTCF peak

CTCF peaks were annotated with the number of different loops (partners) per peak by using pgltools (v 2.2.0) intersect1D to list each loop-peak interception followed by bedtools (v 2.26.0) merge on the list of peaks to sum the number of different loops in which each peak falls within an anchor ^83,86. To annotate peaks with the number of different loops their maximally-looping partner makes, pgltools (v 2.2.0) intersect1D was again used to list each loop-peak interception. Bedtools (v 2.26.0) map was used to annotate the number of different loops the peak in each partner anchor makes from the previous file using the -o max flag. This resulted in one entry per loop that lists the “original” peak and the number of loops its partner anchor makes. The maximally-looping partner per original peak was mapped back onto the original peak list with bedtools (v 2.26.0) map using the -o max flag so that each peak was listed once.

Calling and characterizing highly-looping peaks

Highly-looping peaks were called by adapting code used previously to define methylation rich regions (MRRs) and superenhancers ⁸⁷. Specifically, peaks were ordered by rank with respect to the number of different loops they make. The tangent line was drawn where the slope becomes 1 and any peaks above this cutoff were considered highly-looping for a given condition. To generate aggregate profile plots of public K562 histone modification bigwig files (see source data) over CTCF peaks, deepTools (v 3.4.3) computeMatrix was used with the following parameters: reference-point --referencePoint center -bs 50 -a 5000 -b 5000 --skipZeros and results were visualized with deepTools (v 3.4.3) plotProfile. Overlap of peaks with subcompartments and TSA-seq was performed as described above for all CTCF peaks. The following published files were used from the ENCODE consortium for the peak characterization pileups: ENCFF239EBH.bigWig, ENCFF469JMR.bigWig, ENCFF600FDO.bigWig.

Calculating looping changes across conditions

To calculate foldchange and mean CPM loop strength between conditions, mango Hichipper output files for individual replicates were input into Diffloop (v 1.16.0) and processed with the loopsmake mango function ⁴⁸. The full loop list was then filtered for loops with FDR <0.01 and with at least 8 PETs across samples for a given cell type and loops that appeared strongly (>5 PETs) in only one replicate but 0 PETS in the other were removed. The quickAssoc function was then used to get loop strength foldchange and CPM values for each loop. Diffloop output files obtained from the summary() function were then intersected with genomic features (such as reactivated peak overlap, subcompartment, Cut&Tag signal) for further analyses using pgltools (v 2.2.0) intersect1D and redundant loops were merged using pgltools (v 2.2.0) merge and files were further processed in python (v 3.7.0) (www.python.org) with pandas (v1.0.3), matplotlib (v 3.2.1) seaborn (v 0.11.2) and statannot (v 0.2.3) ⁸⁶.

Cut&Tag data processing

Paired end reads were aligned to hg38, deduplicated and converted to bam and tdf files as described above for ChIP-seq. To generate aggregate profile plots of SON Cut&Tag signal over CTCF peaks, normalized bigwig files averaging both replicates were generated from deduplicated bam files with deepTools (v 3.4.3) bamCompare with --operation mean --normalizeUsing RPKM --scaleFactorsMethod None. Then, deepTools (v 3.4.3) computeMatrix and plotProfile were used to plot SON Cut&Tag signal over categories of CTCF peaks with the following computeMatrix parameters: reference-point --referencePoint center -bs 50 -a 1000 -b 1000 –skipZeros. To compare Cut&Tag binned signal to peak and gene sets, replicate-averaged bedgraph files were generated with deepTools (v 3.4.3) bamCompare with --normalizeUsing RPKM --operation mean -of bedgraph --binSize 20000 --scaleFactorsMethod None. Binned Cut&Tag signal was mapped onto genes/peaks with bedtools map (v 2.26.0) with -o mean. To compare loop strength and loop foldchange from Diffloops summary file (described above) to binned Cut&Tag signal, the replicate-averaged Cut&Tag bedgraphs in 20 kb bins were intersected with loop files using pgltools (v 2.2.0) intersect1D with -wa -wb -allA followed by pgltools merge with -o mean ⁸⁶. For all plots with deciles, equal deciles for each plot were determined in python (v 3.7.0) with pandas (v 1.0.3) qcut with duplicates =’drop’.⁸⁸ Subsequent processing for all plots was completed in python (v 3.7.0) (www.python.org) with pandas (v1.0.3), matplotlib (v 3.2.1), seaborn (v 0.11.2) and statannot (v 0.2.3).

RNA-seq processing and differential expression analysis

RNA-seq data were analyzed following Lexogen’s recommended protocol. First, UMIs were extracted with umi-tools (v 1.1.1) extract ⁸⁹. Adapters were then trimmed with cutadapt (v 2.7) ⁹⁰. Trimmed reads were aligned to hg38 using Star aligner (v 2.6.0) with a genome generated using the Ensembl hg38 primary assembly and Ensembl GTF file v 104 with the following settings: --outFilterMultimapNmax 20 --alignSJoverhangMin 8 --alignSJDBoverhangMin 1 --outFilterMismatchNmax 999 --outFilterMismatchNoverLmax 0.6 --alignIntronMin 20 --alignIntronMa× 1000000 --alignMatesGapMa× 1000000 --peOverlapNbasesMin 40 --peOverlapMMp 0.8 --outSAMattributes NH HI NM MD⁹¹. Reads were then deduplicated with umi-tools dedup (v 1.1.1) with flags multimapping-detection-method=NH and unpaired-reads=discard.

For gene expression analysis, HTseq – count (v 0.11.2) was used to count reads in genes ⁹². Count files were then processed in DEseq2 (v 1.28) with R (v 4.0.2) using the standard workflow to get foldchange values and using alpha=0.05 ^93,94. Hierarchical clustering was performed on the regularized log₂-transformed counts using Euclidean distance as a metric and the complete clustering method. All genes with both an adjusted P < 0.05, calculated using the Benjamini-Hochberg correction, and |log₂(fold-change)| > 0.5 between wildtype DNMT1i and DMSO were included in the clustered heatmap. The regularized log₂-transformed counts for the variable genes were grouped by k-means clustering (k = 4) and the Z-scores for the variable genes, ordered according to their cluster identity, are depicted. For differential exon analysis, Dexseq (v 1.34.1) was used. Files were then imported into R and foldchange values were computed according to the standard Dexseq workflow. BioMart was used to map the genomic ranges for each gene onto RNA-seq differential expression output files. Gene expression changes were compared with Cut&Tag signal and CTCF peak overlap using bedtools (v 2.26.0) map with -o mean. Gene Ontology analysis was completed on the online DAVID portal by filtering for goterm_direct and goterm_BP_5 results with FDRj<0.2 ^95,96. Further analysis, including volcano plots, was completed in python (v 3.7.0) (www.python.org) with pandas (v 1.0.3), matplotlib (v 3.2.1) seaborn (v 0.11.2) and statannot (v 0.2.3).

Data Access

All sequencing data has been deposited on NCBI Gene Expression Omnibus under the super series accession code GSE217760. Custom code is available at https://github.com/liaulab.

Competing Interests

B.B.L. is a shareholder and member of the scientific advisory board of Light Horse Therapeutics. The remaining authors declare no competing interests.

Acknowledgements

We acknowledge members of the Liau lab for useful discussions. We acknowledge J. Nelson, C. Maesner and Z. Niziolek for assistance with cell sorting. We acknowledge the FAS Division of Science Research Computing Group at Harvard University for assistance with cluster computing. We acknowledge S. Johnstone for support with the HiChIP protocol. We acknowledge T. Lu for helpful discussions about speckle Cut&Tag. We acknowledge H.Roh for providing Protein A Tn5. A.P.S and C.L. were supported by the Herchel Smith Graduate fellowships. N.Z.L. was supported by a NSF Graduate Research Fellowship (grant no. DGE1745303). A.L.W was supported by a Simmons award from Harvard Center for Biological Imaging. This research was supported by startup funds from Harvard University.