Comparative Hi-C Reveals that CTCF Underlies Evolution of Chromosomal Domain Architecture

Sequence-Driven Evolution of CTCF Binding Profiles

CTCF binding is strongly correlated with the topological architecture of mammalian chromosomes and participates in long-range chromatin loops, thereby underlying global contact insulation. We analyzed mouse (Mus musculus [Mmus]), dog (Canis familiaris [Cfam]), and macaque (Macaca mulatta [Mmul]) CTCF chromatin immunoprecipitation sequencing (ChIP-seq) profiles from primary liver cells (Schmidt et al., 2012), aiming to define how conservation and divergence of the insulator binding landscape co-evolve with chromosomal topology. Pairwise CTCF ChIP-seq analysis identified conserved or divergent CTCF binding sites within syntenic chromosomal regions (Figures 1A, 1B, and S1). Sites with the strongest CTCF binding intensities were highly conserved (77% of the top 0.1 percentile), while lower-intensity CTCF binding sites were enriched for divergent binding (57% of mouse-divergent sites) (Figure 1B). We computed the sequence affinity of the different classes of binding sites to the canonical CTCF consensus motif in mouse and found that the levels of motif affinity for conserved sites were overall higher than the level for the mouse-divergent sites (Figure S1).

To understand the relationship between sequence affinity and CTCF binding at conserved or divergent sites, we correlated changes in CTCF binding with changes in CTCF sequence motif affinity among species. For this analysis, we used the canonical consensus motif from mouse that is the same in the other species (Schmidt et al., 2012). Remarkably, we found a direct association between sequence divergence and CTCF binding divergence. Conserved CTCF binding sites showed overall high motif affinities and a high degree of affinity conservation. Conversely, motifs underlying the divergent sites were evolutionarily dynamic and diverged in strong correlation to divergent binding intensity (Figures 1C and S1). The data show that when strong motifs in CTCF binding sites diverge, CTCF binding itself is concomitantly gained or lost. Interestingly, 65% of the sites that were conserved between mouse and dog were also conserved in macaque, while macaque-specific and dog-specific sites constituted another two populations of 775 and 891 sites, respectively, with weaker, more evolutionarily plastic motifs. Together, these data suggest that the CTCF insulator landscape is evolving under two regimes: the first involves a tight conservation of both sequence and binding landscape, and the second shows a dynamic interplay between divergence of specific cis elements and consequential evolution of the CTCF binding trait. The relatively direct influence of motif divergence on CTCF binding forms a potential link between sequence evolution and large-scale genome evolution.

CTCF Binding Site Evolution Is Correlated with the Mouse Hi-C Domain Structure

To investigate how the different classes of CTCF binding conservation correlate with chromosomal structure, we prepared Hi-C datasets on mouse liver cells (Figure S2). Filtering and normalization of the Hi-C ligation products was performed as before (Sofueva et al., 2013), revealing the characteristic chromosomal domain structure in these cells. Visualization of the CTCF occupancy groups with the Hi-C contact maps suggested that conserved CTCF binding sites were found at the borders of large-scale Hi-C domains while species-specific CTCF sites are located internal to domains (Figure 1D). This observation was supported with a genome-wide analysis whereby the relative position of conserved and divergent CTCF sites was determined with respect to all domains in the mouse genome (Figure 1E).

To further characterize chromosomal contacts around conserved and divergent CTCF sites, we analyzed the average contact distribution around these sites globally, measuring “contact insulation” by quantifying the decrease in contact probability between multiple elements separated by a CTCF site (Sofueva et al., 2013). Analysis of the composite contact insulation at multiple distance ranges, indicated strong insulation profiles for conserved CTCF sites, further supporting the idea that these conserved, high-intensity CTCF sites were co-occurring with the borders of large-scale domain (reminiscent of topological chromosomal domains) (Figure 1F, left panel). In comparison, the lower-intensity mouse-divergent sites showed a significantly weaker, more localized insulation profile (Figure 1F, right panel). Similar trends were also observed when classifying CTCF sites according to their conservation in macaque (Figure S1). In summary, we found strong correlation between the evolutionary dynamics of CTCF binding sites and mouse chromosome topology, indicating the possibility of a direct link between insulator site divergence and the evolution of topological domain structure.

Comparative Hi-C Reveals the Evolution of Chromosome Topologies

We used comparative Hi-C to examine conservation and divergence of chromosome topology and to test how evolution of CTCF binding sites might underpin this. We collected liver cells from macaque, rabbit (Oryctolagus cuniculus [Ocun]), and dog and processed them using the same approach applied to mouse, yielding chromosomal contact maps for each of these species (Figures S2 and S3). Evaluation of the overall topological structure within the three newly profiled species first indicated the integrity of their reference genome structures and generated a resource for future refinement of such assemblies. More importantly, the data showed that the chromosome topologies in macaque, dog, and rabbit are characterized by a chromosomal domain structure that is similar to the one inferred before for human and mouse (Dixon et al., 2012). For example, comparison of a 9-Mb syntenic region highlighted the extensive conservation of chromosomal structure across all species (Figure 2A). The maps also revealed evidence of intra-domain differences between species (Figure 2B). We quantified the extent of structural conservation genome-wide using a computational approach that allowed us to comprehensively describe domain structure at multiple scales. This pairwise approach revealed extensive genome-wide interspecies conservation of chromosome structure (Figures 2C and S3). A systematic analysis of paired domains in mouse and dog revealed that conserved domains are smaller in size compared to other domains and are classified as both active and passive clusters (Figure S4). Together, these data facilitated extensive analysis of the evolution of chromosomal topologies within regions that did not go through substantial genome rearrangement, allowing examination of the evolution of both large-scale domain borders and the insulation structure within domains.

Divergent CTCF Binding Drives Local Structural Change within Domains

Hi-C maps from liver cells of different species allowed us to ask how the evolutionary dynamics of CTCF correlate with conservation or divergence of domain structure. Analysis of specific loci showed that conserved CTCF sites were typically located at the borders of large-scale chromosomal domains that were themselves conserved between mouse and dog (Figure 3A). To test these observations globally, we computed the contact insulation profiles from either the mouse or dog Hi-C maps around conserved CTCF sites, showing that these sites indeed globally served as conserved insulation points (Figure 3B). Similar results were derived using a comparison of mouse and macaque (Figure S5). We also observed that conserved CTCF sites were strongly enriched for Rad21 in mouse (79% of conserved sites compared to 51% mouse-divergent sites co-localize with Rad21) and that CTCF/cohesin co-occupied sites exhibited strong contact insulation in all three species (data not shown).

In contrast to these highly stable sites, our data showed that divergent CTCF sites were located primarily within domains and exhibited local contact insulation. Comparative analysis of contact insulation at divergent CTCF sites revealed that indeed these sites correlated with divergent contact insulation profiles. For example, dog-divergent CTCF sites (Mmus−/Cfam+) exhibited local contact insulation specifically in the dog genome, whereas these same sites exhibited background levels of contact insulation when examined in the mouse Hi-C data (Figures 3C and S5). Importantly, the change in insulation following CTCF binding site evolution was stronger at the local (20-kb) scale, but was not significant at the higher (80-kb) scale (Figure 3D), suggesting that large-scale domain changes either are not affected by CTCF evolution or are under strong negative selection and are therefore not observed. These observations were further strengthened when we examined CTCF binding sites that were “partially” conserved. CTCF sites that were bound in mouse and dog, but lost in macaque, were associated with reduced contact insulation in the macaque genome. Thus, the data demonstrate a relationship between CTCF binding divergence and divergence of local insulation structure and therefore point to a role for CTCF in driving structural change in the genome.

The continuous evolutionary dynamics of intra-domain looping can play a key role in tuning promoter-enhancer contacts within domains. Consistent with this, we observe long-range contacts between divergent CTCF sites and enhancers or transcription start sites (TSSs) (Figure S6). Furthermore, analysis of transcription data from mouse, dog, and macaque liver reveals that divergent CTCF sites are contacting differentially expressed genes with a greater frequency than non-differentially expressed genes (Kolmogorov-Smirnov test, p < 0.05) (Figure S6). Together, these data support the hypothesis that emergence of divergent CTCF binding sites can contribute to changes in gene expression.

Conserved CTCF Sites Are Directional and Interact with Other Conserved Sites

While it is known that CTCF binding specificity greatly depends on its specific DNA consensus sequence, the specificity and directionality of CTCF/cohesin long-range contacts (Sofueva et al., 2013) and the way by which specific sites are assembled to define topological domains are not fully understood. As our data indicated that conserved CTCF binding sites have conserved motif affinities (Figure 1C), and because it is known that the CTCF consensus motif is non-symmetric, we asked whether conserved sites could also be conserved for the orientation of the CTCF motif. Indeed, 94% (3,265/3,483) of CTCF binding sites that are conserved between mouse and dog are also conserved in their orientation. To explore this further, we profiled contact insulation around conserved CTCF binding sites grouped according to the strand that the consensus motif was found on. We observed an asymmetric insulation behavior that was mirrored when the orientation of the motif was reversed (Figure 4A). This analysis uncoupled “insulation” (blue) from “preferential contacts” (red) and revealed that preferential contacts are made on one side of the oriented CTCF binding site, indicating that the orientation of the motif likely contributes to directionality of CTCFs long-range interactions. Consistent with this, we profiled the genome-wide relative position within chromosomal domains of Mmus+/Cfam+ conserved CTCF sites grouped according to the orientation of their binding motif as above. We observed that the conserved CTCF binding sites that are enriched at the edges of conserved domains (Figure 1E) have a specific orientation of their motif relative to chromosomal domains (Figure 4B). These observations were replicated when we compared mouse and macaque.

To characterize the contact relationship between evolutionarily stable or flexible CTCF sites and to further understand how they contribute to the evolution of chromosome domain structure, we performed a high-resolution high-throughput circular chromosome conformation capture (4C-seq) study. We designed four 4C-seq viewpoints to a series of neighboring conserved CTCF binding sites bordering conserved domains in the mouse and dog as well as to a mouse-specific site. The results showed that each conserved CTCF site engages in very strong and directional interactions with neighboring conserved CTCF sites (Figure 4C). Remarkably, the specific interactions mediated by conserved sites in the mouse genome were themselves precisely conserved in the dog genome (Figure 4D) and define the underlying domain structure. In each case, the long-range interaction was anchored by a pair of conserved CTCF sites whereby one CTCF site had an orientation on the “+” strand and the other on the “−” strand and could provide the basis for the observed directionality of CTCF-mediated interactions. Moreover, a viewpoint designed to a mouse-divergent site exhibited weak interactions within the mouse domain, analogous to the local insulation behavior observed in Figure 3B (Figure 4C). Importantly, the mouse-divergent viewpoint had no prominent interactions in the dog genome, confirming the specificity of its interaction network.

Global analysis of Hi-C contacts between pairs of CTCF binding sites stratified according to genomic distances in cis (Sofueva et al., 2013) confirmed the 4C-seq observation systematically (Figure 4E). Consistent with the high-resolution 4C-seq profiles, Hi-C trends showed that conserved CTCF sites strongly contacted one another within the same domain. Divergent CTCF sites engaged in significantly weaker contacts with other divergent sites, even when stratifying thoroughly for genomic distances. Importantly, little to no contact was observed in the mouse genome between dog-divergent sites. These results show that evolutionarily stable CTCF sites are engaged in strong contacts with one another and suggest that in so doing, they create an interaction network that may support the conservation of domain structure. On the other hand, divergent CTCF sites are involved in weaker interactions, perhaps reflecting the evolutionary flexibility of the binding sites themselves.

Domains Maintain Their Integrity during Chromosomal Rearrangements

Our data suggested that large-scale domain re-organization does not typically occur following insulator divergence. How then can it still be observed? Our interspecies comparative Hi-C data allowed us to ask what happens to the integrity of conserved chromosomal domains when genomes are challenged by structural rearrangements. If chromosomal domains act as modular units (e.g., to regulate gene expression), then large-scale rearrangements would be expected to occur at domain borders, so as to maintain the integrity of these structures. We scanned the mouse and dog genomes for differences in the distance between contiguous orthologous genes in the two species. Our analysis uncovered a number of complex rearrangements between the mouse and dog genomes involving insertions, inversions, and duplications. In each case, we discovered that the rearrangement occurred at the border between two chromosomal domains. This is exemplified in the Hi-C map from chromosome 15 in dog (Figure 5). Here, we found two domains, one containing the Slc5a9 gene and the other containing the Trabd2b gene (highlighted by red dots). Comparison of this region to the mouse genome revealed that a 2-Mb insertion occurred in the mouse genome that contains the Skint gene cluster, which is rapidly evolving and unique to the mouse lineage (Boyden et al., 2008). Remarkably, the insertion occurred directly between two neighboring dog domains in such a way as to perfectly maintain their integrity. A similar rearrangement event occurred at the Mrgpr gene cluster in the mouse genome (Dong et al., 2001), again preserving the structure of the neighboring domains (Figure S7). In another example, we observed a large-scale 5.5-Mb insertion in the dog genome containing multiple domains, and again, the domains on either side of the insertion have been maintained intact (Figure S8). These examples suggest that domains function as modular units and are selected against breakage during genome rearrangements.

Liver Homogenization and Fixation

Fresh or frozen liver from mouse, rabbit, macaque, and dog were processed for Hi-C or 4C-seq libraries. With the exception of mouse, the samples used for the Hi-C libraries were the same as the material used for CTCF ChIP-seq (Schmidt et al., 2012). Livers were fixed in 10% formalin for 20 min, and ∼1 g was cut and processed with a Dounce homogenizer (ten strokes with a loose pestle followed by ten strokes with the tight pestle). After filtration through a 70-μm nylon cell strainer, the sample was washed twice with PBS, spinning down at 852 rcf for 5 min at 4°C to collect the cells between washes. 1–5 × 10⁷ liver cells were then fixed for a second time in fixation buffer (1% formaldehyde, 750 μg/ml BSA in DMEM/Ham’s F12 [Invitrogen]) for 10–30 min at room temperature. The fixation reaction was quenched using 0.125 M glycine for 5 min at room temperature. Samples were washed twice with 10 ml PBS, pelleted into 1 × 10⁷ cells aliquots, and stored at −80°C. Mouse Hi-C libraries were prepared from fresh liver samples of biological replicates (9-week-old C57/BL6 mouse and the pooled livers from 2- to 4-week-old outbred mice. The libraries for the other three organisms were technical replicates.

Propidium Iodide Staining of Hepatocytes

Formaldehyde-fixed liver cells were lysed on ice in a hypotonic buffer (10 mM Tris-HCl [pH 8], 10 mM NaCl, 0.2% Igepal CA-640, EDTA-free protease inhibitors) for 30 min. Nuclei were stained with a propidium iodide (PI) staining buffer (100 μg/ml PI, 50 μg/ml RNase A, 0.05% Triton X-100) for 60 min on ice. Samples were analyzed on a MoFlo cell sorter (Beckman Coulter).

High-Throughput Mapping of Chromatin Interactions via Hi-C

The Hi-C method previously used (Sofueva et al., 2013) was modified to accommodate primary liver samples. Hepatocytes were lysed in Hi-C lysis buffer (10 mM Tris-HCl [pH 8], 10 mM NaCl, 0.2% Igepal CA-640, EDTA-free protease inhibitors) for 30 min on ice. The sample was transferred to Protein LoBind tubes (Eppendorf) and the nuclei were permeabilized by incubation with 0.1%–0.6% SDS for 1 hr at 37°C with 800 rpm shaking. The reaction was quenched with 0.67%–4% Triton X-100, 1 hr at 37°C, 800 rpm shaking. Nuclei were digested in 500 μl 1X NEBuffer 2 with 1500 U HindIII (New England Biolabs) and monitored for maximal digestion of the chromatin template, thus digestion times ranged from 24-72 hr. All other parts of the Hi-C protocol, including library preparation were performed as previously described. 75 bp paired-end sequencing was performed for each library according to manufacturers conditions using the Illumina Hi-seq platform.

Hi-C Interaction Matrix Generation and Domain Calling

Sequencing reads were aligned to the mouse (mm10), rabbit (oryCun2), macaque (rheMac2), and dog (canFam3) genome assemblies using Bowtie 0.12.8 (Langmead et al., 2009). The parameters used for the alignment allowed a maximum of three mismatches and strictly one alignment per read. Processing of the aligned reads and normalization of the interaction matrices were performed as previously described (Yaffe and Tanay, 2011; Sofueva et al., 2013). The pipeline produced normalized matrices of interactions binning the genome at different resolutions. Interaction matrices for each library were generated displaying seven different resolutions simultaneously (12,500, 25,000, 50,000, 100,000, 250,000, 500,000, and 1,000,000 bp). Domains were identified and clustered as described (Sexton et al., 2012) with the modification that scaling factors were inferred using fends 100–400 kb apart, to account for the lower resolution of the mouse map compared to the Drosophila map. Domain borders were called using the 95% percentile of the scaling track. A domain-level map was partitioned into two clusters, and clusters were assigned as passive/active according to Lamin B mouse embryonic fibroblast (MEF) data, as before. For the rabbit, macaque, and dog genome, the Lamin B MEF track for mouse was lifted over to the corresponding genome to label domain clusters. Domain calls in mouse and dog are available in Table S1.

ChIP-Seq Analysis

We used previously published ChIP-seq data for CTCF from mouse, macaque, and dog livers (Schmidt et al., 2012) and for Rad21 for mouse liver (Faure et al., 2012). Rad21 ChIP-seq data for macaque and dog was prepared as for CTCF. Mouse, macaque, and dog ChIP-seq reads were mapped using bowtie. Alignment was followed by extension of sequenced tags to 300-bp fragments and pileup into 50-bp bins. We normalized ChIP-seq coverage by computing the distribution of pile-up coverage on 50-bp bins and transforming each coverage value v into −log10 (1-quantile(v)). To define binding sites, we used a simple threshold on the sum of values from two biological replicates for each CTCF dataset and for the macaque Rad21 data. Rad21 ChIP data from mouse and dog were done in single, and the data were thresholded. Thresholds used were as follows: mouse CTCF = 2.2, macaque CTCF = 2.4, dog CTCF = 2.2, mouse Rad21 = 2.3, macaque Rad21 = 2.5, dog Rad21 = 3. Different thresholds did not change the results. Binding site width was standardized at 200 bp, and the ChIP-seq intensity for each site was calculated as the maximum value across the 200 bp. The relative distribution of CTCF within topological domains (Figures 1E and S5) was calculated as the distance of each CTCF site from the center of its domain. Half the size of the domain was added to convert it to a measure of distance from the edge of the domain, and this number was then divided by the size of the domain.

Interspecies Comparison of CTCF Sites

Macaque and dog CTCF ChIP-seq libraries were converted to mouse genome coordinates using the liftOver tool from UCSC. To reduce the chance of inaccurate liftOver, a number of filters were implemented: sites within low-mappability regions, repeats, or windows of 100 kb with insufficient synteny were excluded. To estimate mappability, each genome was broken into 50-bp reads and the whole-genome sequence was split into artificial reads and then mapped back to the genome. For each 50-bp bin, the mappability score was then defined to be the portion of artificial reads mapped uniquely to that bin. To estimate the level of synteny in the 100 kb around a CTCF site, the mappability tracks for macaque and dog were converted to the mouse genome using liftOver and all bins for which liftOver was not possible were converted to zeroes. The converted tracks were subsequently smoothed over 100 kb, and CTCF sites falling in regions below the top quartile of such smoothed tracks were excluded from all subsequent analysis. Divergent CTCF sites in mouse and dog are available in Table S1.

CTCF Binding Energy Function

A CTCF DNA-binding energy function from the Cortex CTCF binding sites (ENCODE Cortex CTCF mouse, GSM769019; Shen et al., 2012) was used to profile all genomes for their similarity to the CTCF consensus motif. The consensus motif is very highly conserved across all species (Schmidt et al., 2012). Given a set of genomic sites, we compute for each site the maximal energy value within a 200-bp window centered on the point.

Motif Orientation Analysis

Orientation of the motifs underneath conserved CTCF peaks was obtained using MEME (http://meme.nbcr.net/meme/), (Bailey and Elkan, 1994) with the parameters -revcomp -dna -nmotifs 1 -w 20 -mod zoops -maxsize 100,000.

Crossover Analysis

Crossover analysis was performed as described previously (Sofueva et al., 2013). The bands used were 5–7.5, 7.5–11.25, 10–15, 15–22.5, 20–30, 30–45, 40–60, 60–90, and 80–120 kb.

Distal Contact Analysis

To calculate the average interaction profiles for a group of genomic landmarks, HindIII fragment ends were grouped into classes by associating each end with a genomic element located within 5 kb and then grouping all fragment ends associated with an element of the same class. For the mouse, macaque, and dog genomes, three classes of CTCF sites (conserved, divergent present, divergent absent) and TSS sites were defined. These classes were further divided to sites within active or passive Hi-C domains. The remaining fragment end (not classified given other landmarks) was defined as the background.

4C-Seq

Preparation of 4C-seq samples, libraries, sequencing analysis, and normalization were all performed as previously described (Sofueva et al., 2013). Primer sequences were chosen to viewpoint sites that were as close as possible to CTCF ChIP-seq peaks (Table 1). Mouse primers were designed according to the genome-wide 4C-seq primer database from (van de Werken et al., 2012). For dog primers, a similar database was generated for the regions of interest.