Ultra-long-range interactions between active regulatory elements

Screening for factors bound at sites of enriched CRE–CRE interactions

To quantify interactions between CREs, we merged DNase-accessible ENCODE (Luo et al. 2020) CREs from mESCs within 5 kb of each other and generated observed over expected (O/E) contact frequencies for each CRE pair from Hi-C and Micro-C data at both short range (0.1–1 Mb) and long range (1–10 Mb) (Fig. 1A). The 10-Mb cut-off was chosen to limit the number of combinations and to avoid sparse data at larger distances. We determined the overlap of CREs with binding sites (peaks) in the Cistrome DB (Mei et al. 2017) and ReMap2022 (Hammal et al. 2022) databases. For each peak data set representing one ChIP-seq experiment, we divided CREs into overlapping and nonoverlapping. We then compared the contact frequencies of CRE pairs overlapping the factor on both sides, or not at all, by calculating the Mann–Whitney U adjusted P-value and effect size (F = U/(n₁ × n₂)), where F > 0.5 means increased contact frequencies at bound compared with unbound (Fig. 1B). We annotated the factors into six classes: cohesin-associated, polycomb-associated, transcription cofactors, TFs, repressive factors, and other. We use transcription cofactors broadly here to mean chromatin and transcription-associated proteins that are not known TFs or repressors.

Figure 1.

Computational screening for factors associated with enriched interactions between CREs. (A) Schematic of the computational pipeline. CRE pairs are split into short range (0.1–1 Mb between pairs) and long range (1–10 Mb between pairs). These regions are overlapped with peaks from different data sets and contact frequencies at 10-kb resolution compared between pairs overlapped by the factor and between pairs not overlapped by the factor to calculate enrichment values. (B) Example of the analysis for a RING1B ChIP-seq data set. Observed over expected (O/E) contact frequencies for CREs at long range were split into overlapping RING1B peaks on neither or both sides, and Mann–Whitney F and adjusted P-values were calculated: (****) P = 3.8 × 10⁻¹⁸⁴. (C) Effect sizes for factors with significantly enriched chromatin interactions compared with unbound CREs in mESCs. The x- and y-axes show enrichment at short range and long range. Colors represent the group to which the factor belongs. (D) Pileup analysis for 5000 regions for each of the indicated hits from C for interaction pairs at different distances of genomic separation in Micro-C data from mESCs. For CFP1 and KLF4, peaks overlapping RING1B binding sites were excluded.

Open in new tabDownload PowerPointLink to figure

We tested our approach by measuring changes in enrichment in the six different classes in Micro-C and Hi-C data sets from cells in which cohesin-associated factors had been degraded (Hsieh et al. 2022) or Rnf2 (also known as Ring1B) (polycomb) knocked out (Boyle et al. 2020). As expected, RAD21 and CTCF degradation led to a decrease (Supplemental Fig. S1A,B) and WAPL degradation led to an increase (Supplemental Fig. S1C) in enrichment between cohesin-associated binding sites at short range. Rnf2 knockout (KO) led to an expected decrease in enrichment for polycomb-bound regions at both short and long range (Supplemental Fig. S1D).

After confirming that our approach was able to detect factors associated with enriched chromatin interactions, we performed our analysis on Micro-C data from WT mESCs (Fig. 1C; Supplemental Fig. S2A,B; Supplemental Table S1; Hsieh et al. 2020). Cohesin-associated factors were enriched for interactions at short but not long range, and polycomb-associated factors were similarly enriched at short and long range. Besides these already known mechanisms, binding sites for transcription cofactors and TFs tended to be enriched at both short and long range, but with higher enrichment at long range. We correlated the number of non-polycomb-bound TSSs per CRE bound by each factor and its enrichment. TFs and cofactors had a high correlation between enrichment and TSS overlap, showing that they are separate from polycomb and related to the presence of genes (Supplemental Fig. S2C). This is not a feature specific to Micro-C data or to mESCs as we also saw a similar pattern of enrichment for the different classes in Hi-C data from the GM12878 human lymphoblastoid cell line (Supplemental Fig. S2D–F; Supplemental Table S2; Rao et al. 2014). We used pileup analysis (Flyamer et al. 2020) to confirm some of the top hits. Sites bound by the CpG island (CGI) binding protein CFP1 (also known as CXXC1) or by KLF4, a pluripotency TF active in mESCs, showed enrichment at all distances up to 100 Mb (Fig. 1D). We call this noncohesin, non-polycomb-associated category of interactions between active regions at large distances “ultra-long-range interactions” (ULIs).

ULIs are seen as enriched stripes and central pixels in pileups, meaning that these active elements interact with each other and with the surrounding chromatin more than with other regions at similar distances. Note that although the level of enrichment above expected is similar across distances, this does not represent similar absolute contact frequencies. We considered potential technical artifacts that could lead to the appearance of ULIs. We divided accessible regions into quartiles based on DNase-seq signal, split them into those within 1 kb of a TSS or at least 5 kb from the nearest TSS, and generated pileups (Supplemental Fig. S3A,B). Although ULIs scaled with accessibility, TSSs in lowly accessible regions were much more enriched for ULIs than more accessible sites without TSSs, arguing against this signal simply being a result of increased digestion or cross-linking efficiency in accessible chromatin. We also considered that large protein complexes bound at active genes could lead to increased cross-linking efficiency. However, contact frequencies did not correlate with the number of Cistrome DB and ReMap2022 peaks overlapping CREs, suggesting this is not the case (Supplemental Fig. S3C,D). We also excluded normalization artifacts, as ULIs can be seen in both unbalanced and balanced (iterative convergence and eigenvector decomposition [ICE] normalized) (Imakaev et al. 2012) data, whether balanced on all, or only on cis, contacts (Supplemental Fig. S3E). ULIs are also seen regardless of normalizing by expected values, shifted controls, or not at all. We therefore conclude that ULIs represent bona fide enriched contacts between active regulatory elements.

ULIs between active regions are independent of cohesin and polycomb

The strong enrichment for CFP1 binding sites in our screen (Fig. 1C; Supplemental Fig. S2B) prompted us to test the relationship between CGIs and ULIs. This showed a correlation between interactions and the density of CpG dinucleotides at CGIs devoid of polycomb (Fig. 2A). Interaction strength at TSSs also scaled strongly with the level of transcription (Fig. 2B). Notably, TSSs in quartile 4 (highly transcribed) showed enrichment with quartiles 3 and 2 but not with quartile 1, meaning that inactive genes do not interact at all (Supplemental Fig. S4A). Non-CGI promoters also interact but much less than CGI promoters, suggesting CGIs are not required but contribute to ULIs at TSSs (Supplemental Fig. S4B). Analysis of Micro-C data at 100-bp resolution within the 10-kb region surrounding TSSs showed that enriched interactions are centered at the TSS (Fig. 2C).

Figure 2.

Ultra-long-range interactions (ULIs) between active promoters and distal regulatory elements. Pileup analysis. (A) Non-RNF2 overlapping CGIs split by quartiles based on CpG density (2352–3015 peaks per group). (B) Non-RNF2 overlapping TSSs (4800–5128 peaks per group) divided into quartiles by expression level, based on 4sU-seq data, in Micro-C data from mESCs. (C) The 10 kb surrounding the Q4 TSSs (5128 peaks) at 100-bp resolution in Micro-C data from mESCs. (D,E) CGI Q4 (2651 peaks) regions in Micro-C data from RAD21-AID mESCs (D) and Hi-C data from WT or Rnf2 KO mESCs (E). (F) TSS-distal H3K27ac peaks not overlapping RNF2 or CGIs split by quartiles based on H3K27ac ChIP-seq signal (6161–6452 peaks per group) in Micro-C data from mESCs. (G) CGI Q4 regions (2651 peaks; +CGI) and Q4 promoter-distal H3K27ac (6268 peaks; −CGI) in Micro-C data from mESCs.

Open in new tabDownload PowerPointLink to figure

We confirmed that similar to polycomb-dependent contacts but unlike cohesin-driven interactions, interactions between active regions at short distances (<1 Mb) span TAD boundaries (Supplemental Fig. S4C). This indicates that the processes that drive ULIs also operate at shorter distances and are independent of loop extrusion. To formally show that ULIs are not dependent on polycomb and cohesin, we analyzed data from mESCs with Rnf2 deleted (KO) or RAD21 degraded (AID) for 3 h, which disrupts polycomb- and CTCF-mediated interactions, respectively (Boyle et al. 2020; Hsieh et al. 2022). These interventions did not affect ULIs (Fig. 2D,E) and neither did 3 h of depletion of CTCF or the cohesin unloader WAPL (Supplemental Fig. S4D). These results are in line with many shorter-range enhancer–promoter interactions being largely independent of cohesin and CTCF (Hsieh et al. 2022). However, we did note a strong reduction of interactions after 6 h of RAD21 depletion in an independent mESC data set (Supplemental Fig. S4E; Rhodes et al. 2020) but not in data from Rad21 KO noncycling mouse thymocytes (Supplemental Fig. S4F; Seitan et al. 2013). Furthermore, in human HAP1 cells (Haarhuis et al. 2017), we did not detect a loss of ULIs upon deletion of the cohesin subunit MAU2 (also known as SCC4), whereas WAPL deletion led to slightly decreased interactions (Supplemental Fig. S5A). It is possible that the stiffening of the chromatin fiber that occurs upon prolonged loss of WAPL (Tedeschi et al. 2013) leads to a less dynamic chromatin structure, preventing active regions from interacting. Overall, of the cohesin subunit perturbations examined, only 6 h of RAD21 degradation in mESCs led to a large reduction of ULIs, whereas shorter (3-h) degradation or complete KO in HAP1 and thymocytes did not. This may be because of the indirect effects of long-term RAD21 depletion in mESCs including, but not limited to, the cell cycle. We cannot exclude that cohesin has an impact in certain contexts, but the above results strongly suggest loop extrusion is not required for ULIs.

Analysis of a Hi-C data set from cell cycle–synchronized erythroblasts (Zhang et al. 2019) showed, as expected, that ULIs are lost in mitosis and reappear as early as anaphase/telophase (Supplemental Fig. S5B). This is before the formation of both CTCF loops and compartments (Zhang et al. 2019). Analysis of cell cycle–phased merged single-cell Hi-C data (Nagano et al. 2017) shows that ULIs persist throughout the G₁, S, and G₂ phases (Supplemental Fig. S5C).

Promoter-distal H3K27ac-marked CREs, namely, putative enhancers, which do not overlap with CGIs, also show enriched ULIs that scale with the level of H3K27ac enrichment (Fig. 2F). ULIs also occurred between these distal H3K27ac regions and CGIs, as well as expressed genes (Fig. 2G; Supplemental Fig. S5D). Taken together, these results show that ULIs scale with the activity of regulatory elements, are independent of polycomb and cohesin, are present in cycling and noncycling cells, are reformed quickly after mitosis, and persist throughout the cell cycle.

Rare interactions between distal CREs

Enrichment in the pileups represents averages of many pairs of regions. High average enrichment could come from a few highly interacting pairs or from a tendency of many or all pairs to interact. To test this, we assessed the distribution of contact frequency values using two approaches, focusing on TSS quartile pairs at 10- to 25-Mb distance. First, we quantified the contact frequency in the central 5-kb bin containing the TSS. Higher expression quartiles had higher contact frequencies, although most values were zero in all quartiles, reflecting the sparsity of data at these large distances (Fig. 3A). We also looked at the signal of the individual “corner stripes” in the pileups (Fig. 3B). This revealed a higher number of nonzero values across the stripe with increased expression level (Q3 and Q4). We did not find evidence of a few particularly highly interacting pairs dominating the pileups, showing rather that the observed enrichments are a property of increased interaction frequencies between many pairs.

Figure 3.

Distribution of contact frequencies and DNA FISH. (A) Empirical cumulative distribution frequency of O/E contact frequency values in the central 5-kb pixel between TSSs in different expression quartiles separated by 10–25 Mb. (B, left) Depiction of the “corner stripe” from a pileup. (Right) Heatmaps of O/E contact frequency values in individual corner stripes for TSS–TSS quartile interaction pairs at 10- to 25-Mb distance. The line plot on top represents the mean signal. (C) Hi-C matrices for regions on Chromosomes 1, 9, and X used in DNA FISH experiments. CGI marks CGIs without RNF2. (D, top) Interprobe distances measured by DNA FISH between distal CREs on Chromosomes 1, 9, and X or a CRE and an equidistant control region in the same A compartment. N = 289 (Chr 1), 292 (Chr 9), and 208 (Chr X) measurements each for CRE–Control and CRE–CRE. Boxplot lines denote maximum (excluding outliers), interquartile range (upper and lower bound of box), median (center of box), and minimum values. (Bottom) Example images of DNA FISH. Blue indicates DAPI; green, CRE1 probe; magenta, CRE2 probe; and cyan, control probe (adjacent to CRE2 in the genome). Scale bar, 2 μm.

Open in new tabDownload PowerPointLink to figure

The analyses above show that ULIs represent higher average contact frequencies between active CREs compared to the surrounding chromatin. Contact frequencies are related to, but are not a direct measure of, spatial distances (Fudenberg and Imakaev 2017; Finn et al. 2019). Nevertheless, we expected that active CREs separated from each other by large genomic distances would be closer together in the nucleus on average than with regions at similar distances that are not active CREs. To test this, we performed DNA FISH on mESC nuclei using probes between pairs of CREs at 9.1, 8.9, and 6.8 Mb of genomic separation or between one CRE and an equidistant control region in the same A compartment (Fig. 3C). CRE–CRE pairs had significantly smaller distances separating them than CRE–Control pairs (linear mixed-effect model, P = 0.03); although medians were lower for all individual CRE–CRE pairs, only one of these pairs (Chr X) was statistically significant (Fig. 3D; Supplemental Fig. S5E). This indicates that although the average enrichment between distal CREs is high in pileups, this only translates to marginally smaller distances between individual distal CREs in nuclei (median distance in all cases was >1 micron). We only detected two instances of colocalization (<200 nm distance; both in CRE–CRE pairs) between 3156 total measured distances, showing that the enriched contacts seen between many regions in Hi-C/Micro-C are exceedingly rare between individual pairs (see Discussion).

ULIs occur primarily in cis and in the A compartment

Because ULIs are enriched at such large genomic distances along a chromosome, we tested if they are also enriched in trans between chromosomes. Indeed, high CpG-density regions devoid of polycomb are enriched in trans interactions (Supplemental Fig. S6A), but although the level of enrichment compared with that expected is comparable to that within chromosomes, the absolute trans contact frequencies are very low. We also considered that interactions might occur between homologs and tested this using Hi-C data from hybrid mESCs with phased SNPs (Han et al. 2020). This showed that interactions occur primarily within chromosomes and not between homologs (Supplemental Fig. S6B).

Active regions in the A compartment at distances of several tens of megabases, such as those we look at here, will be separated by multiple intervening A and B compartments. We wondered if the stripe seen in the pileups would span the whole intervening chromatin or be constrained to A compartments. We picked regions close to A/B compartment switches and saw that the interaction enrichment was confined to the A compartment; namely, the stripe does not span into the B compartment (Supplemental Fig. S6C).

Interaction dynamics and DNA methylation

Their relationship to transcriptional activity and H3K27ac suggests that ULIs are dynamic between tissues. To investigate this, we used data from mESCs and differentiated neural progenitor cells (NPCs) (Bonev et al. 2017) and selected H3K27ac peaks enriched in either of the cell states (Fig. 4A). These regions showed correspondingly higher enrichment of ULIs in the cell state with higher H3K27ac. To investigate more rapid dynamics, we examined Hi-C data from human macrophages stimulated with lipopolysaccharide (LPS) and interferon (IFNG) (Fig. 4B; Reed et al. 2022). Regions that gained H3K27ac upon stimulation also gained ULIs, and regions losing H3K27ac, although relatively lowly enriched to begin with, lost ULIs over the time course. Although the small number of regions affected and the time resolution of the experiment precludes delineating which changes come first, it appears that changes in H3K27ac and ULIs accompany one another.

Figure 4.

Interaction dynamics in mouse, human, zebrafish, and Drosophila melanogaster. (A, left) H3K27ac signal in mESCs and neural progenitor cells (NPCs) at mESC-enriched (1509 peaks) and NPC-enriched (2491 peaks) regions. (Right) Pileup analysis for the differentially enriched H3K27ac peaks in Hi-C from mESCs and NPCs. (B, top and bottom) Pileup analysis for regions with increasing (top; 1058 peaks) or decreasing (bottom; 553 peaks) H3K27ac across the time course of LPS and IFNG stimulation of THP-1-derived macrophages. (Middle) H3K27ac signal in THP1-derived macrophages at regions with increasing (top) or decreasing (bottom) H3K27ac. (C, left) H3K27ac signal in zebrafish brain and muscle at brain-enriched (2315 peaks) and muscle-enriched (1994 peaks) regions. (Right) Pileup analysis for the differentially enriched H3K27ac peaks in brain and muscle. (D) Average whole-genome bisulfite (WGBS) signal in the differentially enriched H3K27ac peaks in brain and muscle. (E) Pileup analysis for D. melanogaster TSSs (2943 regions per group), divided by expression level based on RNA-seq, in Hi-C data from D. melanogaster eye-antennal imaginal discs.

Open in new tabDownload PowerPointLink to figure

To examine dynamic changes in ULIs between tissues in vivo, we took advantage of data from Danio rerio (zebrafish) and selected brain- or muscle-enriched promoter-distal H3K27ac peaks (Yang et al. 2020). These were accompanied by enrichment of ULIs in the respective tissue (Fig. 4C). We also saw a corresponding higher level of DNA methylation in the tissue where the regulatory elements were inactive and no interactions were seen (Fig. 4D). Because of the relationship between CGIs and ULIs (Fig. 2A), we reasoned that the focal demethylation at these regulatory elements may be responsible for the interactions and that gain in DNA methylation would lead to their loss. To test if methylation is required for ULIs, we used Hi-C data from D. melanogaster, a species with little DNA methylation and lacking CGIs (Deaton and Bird 2011). We could detect ULIs between expressed TSSs and H3K27ac-positive, non-polycomb-bound (as well as polycomb bound) regions in D. melanogaster eye-antennal imaginal discs (Fig. 4E; Supplemental Fig. S6D; Loubiere et al. 2020). This shows that ULIs are not specific to vertebrates and that DNA methylation is not required for their formation.

ULIs are not directly dependent on transcription

In our initial screening, we had seen enrichment for the binding sites of many sequence-specific TFs at ULIs (Fig. 1C; Supplemental Fig. S2A,B). To determine whether homotypic interactions between TFs could drive ULIs, we compared enrichment values for sites where a particular TF is bound at both interacting sides to where that TF is bound at only one side but there is at least one other TF bound at the other side. This is to avoid comparison with regions devoid of binding altogether at one side, namely, inactive regions, which would skew the enrichment. Binding at both sides yielded a stronger enrichment than at one side, but the two were correlated with significant enrichment for every TF also when bound at only one side (Fig. 5A). As an example, we split H3K27ac peaks into those bound by MYC or not in the GM12878 lymphoblastoid cell line and saw that MYC-bound regions interacted with other active regions not bound by MYC (Supplemental Fig. S7A). Although we cannot exclude some contribution, this analysis indicates that homotypic TF–TF interactions are unlikely to drive ULIs.

Figure 5.

Independence of BRD4, Mediator, EP300 activity, and transcription. (A) Effect sizes of TFs for long-range (1- to 10-Mb) contact enrichment compared with TF-unbound CREs in mESCs. The x- and y-axes show enrichment for TF binding at both sides and only one side, but at least one other TF bound at the other side. (ρ) Spearman's correlation coefficient. (B) Pileup analysis between CGI Q4 regions in dTAG-BRD4 mESCs with or without BRD4. (C) Pileup analysis between CGI Q4 regions in WT, Med15^−/− Med16^−/− Med23^−/− Med24^−/− Med25^−/− (tailless), and MED14− SMASh degron CH12 cells. (D, left) Mean H3K27ac signal in CGI Q4 regions in mESCs treated with DMSO or A-485. (Right) Pileup analysis between CGI Q4 regions in mESCs treated with DMSO or A-485. (E) Pileup analysis between CGI Q4 regions in Micro-C data from WT mESCs or treated with flavopiridol (FLV) or triptolide (TRP) for 45 min. (F) Pileup analysis between CGI Q4 regions in RPB1-AID mESCs with or without RPB1. The top two rows represent asynchronous cells; the bottom two rows represent cells arrested in mitosis (Mit +) and released into G₁.

Open in new tabDownload PowerPointLink to figure

We next tested if transcription cofactors may be required for ULIs. BRD4 has been implicated in chromatin organization (Rosencrance et al. 2020; Linares-Saldana et al. 2021), and both Mediator and BRD4 have been shown to be unevenly distributed in the nucleus and to be enriched in so-called transcriptional condensates (Cho et al. 2018; Sabari et al. 2018). However, our analysis shows that neither BRD4 degradation (Linares-Saldana et al. 2021) nor Mediator disruption (El Khattabi et al. 2019) affects ULIs (Fig. 5B,C). YY1 was implicated as a regulator of enhancer–promoter contacts (Weintraub et al. 2017), but a later report showed that degradation of YY1 had no effect on these (Hsieh et al. 2022). Using data from the latter study, we also found no effect of YY1 on ULIs (Supplemental Fig. S7B). Treatment with the EP300 catalytic inhibitor A-485 (Pelham-Webb et al. 2021) leads to loss of H3K27ac but has no effect on ULIs (Fig. 5D).

We considered whether ULIs may be driven by association with nuclear speckles, as this is strongly correlated with gene activity (Chen et al. 2018). Knockdown of the splicing component Srrm2 leads to a partial disruption of speckles (Hu et al. 2019b), but this did not grossly affect ULIs (Supplemental Fig. S7C). Nuclear speckle association is primarily seen in the A1 subcompartment (Chen et al. 2018). We split active regions (H3K27ac peaks) in GM12878 cells by A1 and A2 (Rao et al. 2014) association and saw a higher level of ULI enrichment in A2 (Supplemental Fig. S7D). These data suggest that speckles cannot fully explain ULIs, although we cannot exclude that they may have some contribution.

Active regulatory elements, including CGIs, tend to be preferentially located toward the nuclear interior (Boyle et al. 2001; Beck et al. 2018). To test whether ULIs are a reflection of a more central nuclear position of these regions, we used GPSeq data from HAP1 cells, which measures radial positioning genome-wide (Girelli et al. 2020). We divided GPSeq regions into three bins (1–3) from the periphery to the center of the nucleus. As expected, CGIs lacking H3K27me3 in HAP1 cells were enriched in the most central GPSeq bin (Supplemental Fig. S7E). Analysis of pileups between regions in the three bins revealed no trend toward more interactions between centrally occupying regions, in cis or in trans (Supplemental Fig. S7F). However, this analysis normalizes the signal to the corners of the pileup to see enrichment compared with surrounding chromatin. When excluding this step and looking simply at average O/E contact frequencies in the different bins, interactions are highest between regions in the most central bin at the largest genomic distances (25–100 Mb) and in trans (Supplemental Fig. S7G). This likely reflects the higher level of chromosome intermingling in the center of the nucleus (Kalhor et al. 2011). Our analyses indicate that (1) interactions between CGIs happen more often in the center of the nucleus because they tend to be localized there, (2) focal enrichment in contact frequencies between CGIs over the surrounding chromatin is similar regardless of nuclear radial positioning measured by GPSeq, and (3) trans or very distal cis interactions are higher between regions in more centrally occupying regions. We therefore conclude that radial positioning contributes to, but does not fully explain, ULIs.

As levels of transcription are correlated to ULIs, we investigated whether loss of transcription would disrupt ULIs. Data from mESCs treated with flavopiridol or triptolide (Hsieh et al. 2020), which block transcription elongation and initiation, respectively, for 45 min did not show a loss of ULIs (Fig. 5E). Neither did degradation of RPB1, the largest subunit of Pol II, for 6 h (Fig. 5F; Jiang et al. 2020). ULIs are lost in mitosis and regained before cells enter G₁ (Supplemental Fig. S5A). When analyzing data from cells that were synchronized in mitosis and then released into G₁ with or without RPB1, we saw ULIs in both conditions (Fig. 5F). These results for long-range interactions are similar to the conclusion that inhibition of transcription has a modest global impact on shorter-range promoter–promoter and enhancer–promoter contacts (Hsieh et al. 2020; Jiang et al. 2020; Goel et al. 2023), although prolonged Pol II depletion can lead to the reduction of some enhancer–promoter contacts owing to changes in cohesin binding (Zhang et al. 2023).

Simulations of multivalent binding are consistent with observed long-range interactions

In an attempt to understand potential microscopic mechanisms underlying ULIs, we turned to simulations of chromatin. We considered a model in which binding and unbinding of multivalent factors to chromatin lead to “bridging-induced clustering,” without any affinity between the factors themselves (Fig. 6A; Brackley et al. 2016). This occurs because multivalent binding causes looping of chromatin and increased concentration of chromatin and binding sites near the already bound sites (Fig. 6B). This causes a feedback loop whereby more factors bind, further increasing the concentration of chromatin and binding sites. To test whether such a model would lead to increased interactions also between very distal binding sites, we modeled a chromatin fiber consisting of 5000 beads (3 kb per bead, representing 15 Mb of chromatin in total) using a bead-and-spring polymer model. In the fitting-free model, 172 of the beads represent randomly distributed binding sites for five different multivalent binding factors. Multivalent binding factors (200 in total, 40 per type) diffuse and bind with low affinity to all beads (nonspecific binding) and with high affinity to their corresponding binding sites.

Figure 6.

Simulations of chromatin and multivalent binding proteins. (A) Schematic depiction of the model with different multivalent binding beads (top) that bind specifically to colored binding sites and nonspecifically to other chromatin beads (bottom). (B) Schematic illustration of bridging-induced clustering. Binding of a multivalent factor leads to bridging and increased concentration of chromatin and binding sites in the vicinity of the bound site. This leads to increased likelihood for more factors to bind, causing a feedback loop and clustering. (C) Snapshot of a configuration in the model. Left side represents chromatin (polymer colored from red to blue), with specific binding sites represented by differently colored beads. Right side represents only the specific binding sites. Instances of clustering are highlighted with arrows. (D) Cumulative contacts in virtual Hi-C between binding sites in the simulations at distances >1 Mb (left) and >5 Mb (right). (E) Proposed model for interactions between CREs. Inactive genes in A or B compartments do not interact. In the A compartment, CREs overlapping loop extrusion boundaries interact via cohesin at short range, and polycomb-repressed CREs interact at short and long range. Active regulatory elements interact with each other at both short and long range, possibly through bridging-induced clustering

Open in new tabDownload PowerPointLink to figure

We ran 180 simulations and detected clustering, as described previously (Fig. 6C). We then generated virtual Hi-C plots (Supplemental Fig. S8A) and analyzed total interactions between binding sites. We detected enrichment between binding sites at >1 Mb and >5 Mb distance, as well as a stripe emanating from the binding sites, similar to what is observed in Hi-C/Micro-C data (Fig. 6D). We considered what could lead to the stripe pattern seen in both the Hi-C/Micro-C data and our simulations. When we performed the simulations without the nonspecific affinity of binding factors to chromatin beads, retaining only specific binding, the central enrichment between sites was weaker, and the stripe pattern was almost completely gone (Supplemental Fig. S8B). This indicates that nonspecific binding is a strong contributor to the overall interactions and the stripe pattern.

Enriched interactions in the simulations occur between binding sites of both the same color (homotypic) and different colors (heterotypic) (Supplemental Fig. S8C). This is expected, as bridging-induced clustering exploits nonspecific interactions between all different binding factors and chromatin, so that mixed clusters form (see Fig. 6C). If we simulate the system with a single binding factor and, concomitantly, more separated binding sites, the enrichment and stripe pattern is weaker (Supplemental Fig. S8D). One prediction of the model is therefore that clustering depends on the distance to other binding sites and that isolated CREs would undergo less bridging-induced clustering than CREs with many nearby CREs. To test this, we took H3K27ac peaks in mESCs with different number of H3K27ac peaks within 100 kb, and sampled peaks to have similar enrichment for H3K27ac in the different groups (Supplemental Fig. S8E). We generated pileups for these regions with different numbers of “neighbor” sites at 1- to 100-Mb distances using Micro-C data. In agreement with the prediction, more isolated peaks had lower enrichment of interactions than did peaks with more neighboring sites (Supplemental Fig. S8F).

We next determined how often two distal regions in the model were in close proximity. We found that although interactions between distant binding sites for factors are enriched, occurrences of colocalization (<200 nm) between simulated FISH probes situated around such sites are rare as in our experimental data (1% and 4% for distances of 5.9 and 1.7 Mb, respectively) (Supplemental Fig. S8G,H).

In summary, our chromatin simulations with multivalent binding factors are consistent with several aspect of ULIs. These include enrichment for long-range (>1 Mb) interactions represented by rare colocalization events, a stripe pattern, no requirement for homotypic binding, and a correlation with the number of nearby active sites.

Hi-C data processing

Already processed Hi-C and Micro-C data were, when required, converted to cooler format (Abdennur and Mirny 2020) using hic2cool (Dekker et al. 2017) or cooler cload pairs (see Supplemental Table S3). Hi-C pairs from allele-phased SNPs (Han et al. 2020) were downloaded from the NCBI Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) under accession number GSE132898 and split into separate files by genotype annotation, followed by conversion to cooler using cooler cload pairs. Other data sets (see Supplemental Table S4) were processed using the distiller-nf 0.3.3 (https://github.com/open2c/distiller-nf) pipeline with default settings, and coolers filtered for mapq ≥ 30 were used. ICE balancing (Imakaev et al. 2012) was performed using cooler balance, using “‐‐cis_only” or “‐‐trans-only” to generate cis/trans balanced weights. Replicate data sets were combined when available using cooler merge. Expected contact frequencies were generated using the cooltools 0.5.1 (Open 2C et al. 2022a) functions expected-cis or expected-trans with “‐‐clr_weight_name” set to the appropriate balancing (cis/trans/total). Hi-C data in Figure 3C were visualized using HiGlass in resgen.io (Kerpedjiev et al. 2018).

Contact screen

ReMap2022 peaks were downloaded from https://remap2022.univ-amu.fr/ (Hammal et al. 2022). Cistrome DB “factor” data for mouse and human were downloaded in batch from http://cistrome.org/db (Mei et al. 2017). Cistrome DB data sets corresponding to perturbation data sets (e.g., siRNA, KO, treatments) or misannotated data were discarded (see Supplemental Table S5). ReMap2022 and Cistrome DB data for mouse or human were combined and overlapped with ENCODE CREs (peaks within 5 kb merged using BEDTools merge) (Quinlan and Hall 2010) overlapping DNase-seq peaks in mESCs or GM12878, respectively (see Supplemental Table S6). O/E contact frequencies at 10-kb resolution (balanced using all contacts) and their respective coordinates for DNase-accessible CREs were extracted using coolpup.py (Flyamer et al. 2020) with the settings “‐‐store_stripes ‐‐flank 0 ‐‐expected expected_file ‐‐mindist 100000 ‐‐maxdist 10000000 ‐‐by_distance 100000 1000000 10000000.” For each data set in the Cistrome DB and ReMap2022 data, the contact frequency scores for short range and long range were split between pairs overlapping the data set at both ends or at neither end. Data sets with fewer than 500 regions or with fewer than 50 overlapping or nonoverlapping regions were discarded. A Mann–Whitney U test (using scipy.stats.mannwhitneyu with default settings) was performed to compare the two distributions and the effect size (F) calculated using F = U/n₁ × n₂, where n₁ and n₂ are the number of observations for overlapping and nonoverlapping regions. Adjusted P-values were calculated by multiplying the Mann–Whitney U P-value by the number of tests performed. TFs and transcription cofactors were derived from the AnimalTFDB 3.0 (Hu et al. 2019a) database followed by additional manual annotation (see Supplemental Tables S1, S2). For Supplemental Figure S1, the F-values were normalized to the mean F in each data set (e.g., Rnf2 KO), and this normalized value was used to divide treated over untreated or KO over WT. For Figure 5A, the analysis was performed in the same way as described above, except peaks were combined for all ChIP-seq data sets of the same factor. Comparison was made between regions unbound by any TF to those bound on both sides by the factor or only one side but at least one other TF bound on the other side.

Pileups

Pileups were generated with coolpup.py (https://github.com/open2c/coolpuppy, using v1.0.0 unless otherwise stated) (Flyamer et al. 2020). Most pileups are generated with the command “coolpup.py cool_file peak_file ‐‐flank 100000 ‐‐mindist 100000 ‐‐maxdist 100000000 ‐‐by_distance 100000 1000000 10000000 25000000 100000000 ‐‐expected expected_file.” Balancing (ICE normalization) was always performed on cis contacts, except for Supplemental Figure S6A, in which trans contacts were used, and Supplemental Figure S6B, in which all contacts were used. For all pileups, 5-kb resolution data were used, except for Figure 2C, in which 100-bp resolution was used. For Supplemental Figure S3E, mindist was set to 1,000,000, and clr_weight_name was altered based on if using total, cis, or no balancing, and “‐‐nshifts 0” “‐‐nshifts 5” or “‐‐expected expected_file” was used. Central pixel interactions and stripes in Figure 3, A and B, were generated using “‐‐mindist 10000000 ‐‐maxdist 25000000 ‐‐store_stripes.” For Supplemental Figure S4C, regions of interest were overlapped with TADs using bioframe.overlap. coolpup.py was run using “groupby = [“TAD”],” where TAD represents a unique value for each TAD in order on the genome. The number of TADs separating the two regions were calculated using a custom modify_2Dintervals_func in coolpup.py. For Supplemental Figure S6C, regions near compartment edges were annotated as + or − stranded depending on if the nearest boundary was to the left or right, and pileups were generated with “‐‐by_strand ‐‐mindist 1000000 ‐‐maxdist 100000000 ‐‐flip_negative_strand” and combining −− with ++ using a custom modify_2Dintervals_func in coolpup.py. For figures comparing interactions between different sets of regions, for example, Figure 2G, the different sets were annotated as + or −, and pileups were generated with “‐‐by_strand ‐‐flip_negative_strand” and combining +− and −+ using a custom modify_2Dintervals_func in coolpup.py. For Supplemental Figure S8C, pileups were generated on peaks annotated by the number of neighbors and binned (see Peaks and Coverage Files section) using the settings “‐‐expected expected_file ‐‐mindist 1000000 ‐‐maxdist 100000000 ‐‐groupby neighbor_bin1 neighbor_bin2 ‐‐ignore_group_order” and using coolpup.py v1.1.0. Pileup heatmaps were generated using the coolpuppy function plotpup.py with the settings “‐‐plot_ticks” and appropriate values for “‐‐cols” and “‐‐rows.” “‐‐norm_corners 10” was used in all cases except for Supplemental Figures S3E, S7G, and S8F. Stripe plots in Figure 3B were plotted using plotpup.py with the settings “‐‐stripe corner_stripe ‐‐plot_ticks ‐‐lineplot.”

Peaks and coverage files

Peaks (bed or narrowPeak) and coverage files (bigWig) used are listed in Supplemental Tables S6 and S7. Overlaps between peaks were generated using BEDTools intersect or bioframe overlap (Open2C et al. 2022b). Peaks were randomly sampled in some cases to get the same or similar number of peaks for comparisons. Distances to the closest peak were generated using BEDTools closest, and peaks were merged using BEDTools merge. For Supplemental Figure S3, A and B, and Figure 2F, ENCODE DNase-seq peaks or H3K27ac peaks not overlapping TSSs, CGIs, or RING1B were split into quartiles based on the signalValue column. CGIs and dinucleotide frequencies were derived from CAP-CGI data (Illingworth et al. 2010). TSSs were defined as the first TSS for all genes in refGene (O'Leary et al. 2016). For Supplemental Figure S4C, TAD boundaries were called at 250-kb resolution using cooltools insulation. Processed H3K27ac values from stimulated THP-1-derived macrophages (Reed et al. 2022) were split into quartiles based on the m0000_VST and m1440_VST columns, and regions with adjusted P-values < 0.05 and going from quartiles 1 or 2 to 4 (up-regulated) or 4 to 1 or 2 (down-regulated) were used. For Supplemental Figure S6C, compartments were called using cooltools eigs-cis at 50-kb resolution, and eigenvector E1 values above and below zero were called as A and B and adjacent bins in the same compartment merged. TSSs not overlapping RNF2 binding sites within 50 kb of compartment edges were compared to Q4 TSSs not overlapping RNF2 at least 100 kb from compartment edges. Zebrafish H3K27ac enrichment values were taken from Figure 2D source data from Yang et al. (2020) and brain/muscle enrichment defined as those with values greater than four in one and less than one in the other tissue. Processed D. melanogaster RNA-seq data from Loubiere et al. (2020) were split into quartiles based on the baseMean value. For GM12878, the top 20000 H3K27ac peaks (based on signalValue) not overlapping H3K27me3 from ENCODE were used and overlapped with ENCODE MYC peaks (Supplemental Fig. S7A) or subcompartments from Rao et al. (2014) (Supplemental Fig. S7D). GPSeq data at 100-kb resolution were downloaded from GitHub (https://github.com/ggirelli/GPSeq-source-data) (source data Fig. 2E), and values from two experiments were averaged and split into three bins, which were overlapped with CAP-CGI regions that do not overlap H3K27me3 in HAP1 (from ENCODE). For Supplemental Figure S8, E and F, the 50 nearest neighbors within ENCODE H3K27ac peaks were generated using BEDTools closest with the settings “-d -k 50.” Neighbors within 2.5 kb were discarded, and the number of neighbors within 100 kb were counted for each peak and divided into groups (zero, one to four, five to 10, 10 or more neighbors). The groups were sampled to generate equal numbers and similar distributions of H3K27ac based on the signalValue column using nullranges (Mu et al. 2023) in R with the settings “covar = ∼signalValue, method = ‘stratified.’” Peaks were converted between assemblies using UCSC liftOver (Hinrichs et al. 2006). Coverage heatmaps and lineplots were generated using deepTools computeMatrix with the settings “reference-point ‐‐referencePoint center -a 5000 -b 5000” and plotted using deepTools plotHeatmap or plotProfile (Ramírez et al. 2016).

Expression quartiles and H3K27ac differential abundance analysis

4sU-seq data from mESCs and differentiating toward NPCs (Boyle et al. 2020) were aligned to the mm10 genome in paired-end mode using STAR 2.7.1a (Dobin et al. 2013) with the settings “‐‐outFilterMultiMapNmax 1.” H3K27ac ChIP-seq data from mESCs (Luo et al. 2020) and NPCs (Bonev et al. 2017) were aligned to the mm10 genome in single-end mode using STAR 2.7.1a with the settings “‐‐outFilterMultiMapNmax 1 ‐‐alignMatesGapMax 2000 ‐‐alignIntronMax 1 ‐‐alignEndsType EndToEnd.” Duplicate reads were discarded using Picard (http://broadinstitute.github.io/picard/), and reads not aligning to autosomes, X, or Y were removed using SAMtools 1.10 (Li et al. 2009). Read counts were generated using HOMER 4.10 (Heinz et al. 2010) annotatePeaks.pl with the settings “-noadj -len 0 -size given” in the first exon of every refGene gene (for 4sU-seq) or merged H3K27ac peaks within 5 kb from mESCs and NPCs (for H3K27ac ChIP-seq). TMM normalization was performed using edgeR (Robinson et al. 2010) in R 4.1.3 (R Core Team 2022). For 4sU-seq, the normalized values from both replicates were averaged, divided by the length of the exon, and split into quartiles. For H3K27ac, DE analysis using limma (Ritchie et al. 2015) was performed with contrast “∼0 + Sample.” Peaks with adjusted P-values < 0.05 were split into those higher in mESCs or NPCs.

DNA FISH

mESCs were grown on 0.1% gelatin-coated slides in GMEM BHK-21 (Gibco 21710-025) supplemented with 15% fetal calf serum (Sigma-Aldrich F-7524), leukemia inhibitory factor (in-house), 2 mM L-glutamine (in-house), 1 mM sodium pyruvate (Sigma-Aldrich 58636), 1× nonessential amino acids (Sigma-Aldrich M7145), and 50 mM 2-β-mercaptoethanol (Gibco 31350-010). Cells were fixed with 4% paraformaldehyde for 10 min at room temperature, washed, permeabilized in PBS/0.5% Triton X-100, dried, and then stored at −80°C before hybridization. Slides were incubated in 100 μg/mL RNase A in 2× SSC for 1 h at 37°C, washed briefly in 2× SSC, passed through an alcohol series, and air-dried. Slides were incubated for 5 min at 70°C, denatured in 70% formamide/2× SSC (pH 7.5) for 40 min at 80°C, cooled in 70% ethanol for 2 min on ice, and dehydrated by immersion in 90% ethanol for 2 min and 100% ethanol for 2 min before air drying. One microgram of fosmid DNA was labeled by nick translation to incorporate green-dUTP (Enzo Lifesciences), Alexa fluor 594-dUTP (Invitrogen), or digoxigenin-11-dUTP (Roche). One hundred nanograms of each fosmid, 6 μL of Cot1 DNA per fosmid, and 5 μg of sonicated salmon sperm DNA were dried in a spin-vac and then reconstituted in 30 μL of hybridization mix. Probes were then denatured for 5 min at 80°C and reannealed for 15 min at 37°C. Fosmid probes were hybridized to slides under a sealed coverslip overnight at 37°C. Slides were washed the next day four times for 3 min in 2× SSC at 45°C and four times for 3 min in 0.1× SSC at 60°C, and the digoxigenin-labeled probe was detected with antidigoxigenin antibody (Roche 11333089001) and Alexa-Fluor 647 donkey antisheep antibody (Invitrogen A21448). Slides were stained with 4,6-diaminidino-2-phenylidole (DAPI) at 50 ng/mL, mounted in VectaShield (Vector Laboratories), and sealed with nail varnish. Slides were imaged on the SoRa spinning disk confocal microscope (Nikon CSU-W1 SoRa), and images were denoised and deconvolved using NIS deconvolution software (blind preset; Nikon). 3D images are shown in the figures as maximum intensity projections prepared using ImageJ (Schneider et al. 2012). The distances between the relevant spots were calculated using the Imaris v9.4 (Oxford Instruments) spot function. Statistical comparison was performed using a linear mixed-effect model using nlme::lme in R with “fixed = distance∼category + chromosome, random = ∼1|replicate,” where category refers to CRE–CRE or CRE–Control (P = 0.0325 for all data). For individual chromosome comparisons, only the values from those DNA FISH experiments were included and lme run with “fixed = distance∼category, random = ∼1|replicate” (P = 0.7777, P = 0.0686, P = 0.0431 for Chr 1, Chr 9, and Chr X, respectively). Probe coordinates are found in Supplemental Table S8.

Chromatin simulations

The general outline of the model is similar to that of Brackley et al. (2021). For a detailed description of the model, refer to the Supplemental Methods. Example code is provided at GitHub (https://github.com/efriman/Friman_etal_ULI). Briefly, LAMMPS (Thompson et al. 2022) was used to perform molecular dynamics simulations of a beads-and-spring polymer representing chromatin with certain beads in the polymer corresponding to high-affinity sites for binding factors. Each bead in the chromatin polymer represents ∼3 kb or a physical size of ∼30 nm. Additional beads representing multivalent binding factors were modeled with different affinities to specific chromatin beads and all chromatin beads. Binding beads undergo both binding and unbinding events. Simulations were initialized with the chromatin fiber as a random walk and binding proteins randomly distributed. Virtual Hi-C, pileups, and FISH were generated from configurations in the second half of the simulated trajectories. For virtual Hi-C, a contact was considered if two regions of the chromatin fiber were within 90 nm and binned at 21-kb resolution (seven chromatin beads). For virtual pileups, the contact distance was 54 nm, and binning resolution was 3 kb (one chromatin bead). Contacts in the pileups are cumulative, and the values represent counts (i.e., number of contacts). Virtual FISH distributions were generated based on the center of mass distances between 51-kb regions surrounding the two sites.

Software availability

The code used in the analysis is available at GitHub (https://github.com/efriman/Friman_etal_ULI) and as Supplemental Code. Custom analysis was performed in Python 3.8.12 (van Rossum and Drake 1995) in Jupyter notebooks (Kluyver et al. 2016) using pandas (McKinney 2010), numpy (Harris et al. 2020), and SciPy (Virtanen et al. 2020) for data analysis and Matplotlib (Hunter 2007) and seaborn (Waskom 2021) for plotting.