BACKGROUND
the ccctc-binding factor  is an 11-zinc-finger protein that is functionally conserved in vertebrates, insects and nematodes  <cit> .
it has been shown to be involved in various critical biological processes and has long been regarded as a versatile regulator.
ctcf knockout mice exhibit early embryonic lethality prior to implantation  <cit> .
in adult organisms, the ctcf protein is ubiquitously expressed across most tissues, and its expression levels vary in a cell type-specific manner, suggesting a role in maintaining phenotypic diversity  <cit> .
furthermore, ctcf depletion results in dysregulated transcription of hundreds of genes in oocytes  <cit>  and dramatically deregulates cell-cycle progression during t lymphocyte lineage commitment within the thymus  <cit> .
other studies have suggested that ctcf impacts cellular activity by playing diverse roles in gene regulation, including promoter activation/repression, genomic imprinting, enhancer blocking, and, most recently, long-range chromatin interactions  <cit> .
recent genome-wide mapping of ctcf occupancy in multiple cell lines has enabled the identification of ctcf binding sequences .
several hypotheses have been proposed to explain the uniquely versatile characteristics of ctcf based on the underlying binding motif sequence.
the most widespread is the zinc-finger model, in which the capacity of ctcf to confer vastly different functions has been attributed to the interplay between the zinc-finger engagement and the underlying sequence  <cit> .
early studies in which zinc-fingers were depleted in a stepwise manner reported involvement of multiple zinc-fingers within a broad ~50 bp sequence  <cit> , indicating that ctcf binding may be partially stabilized by interplay between peripheral zinc-fingers and other cofactors.
following this, a ~11–15 bp core consensus sequence was identified as bound by 4– <dig> central zinc-fingers  <cit> .
a variation at the 12th bp position of the consensus sequence was found to be tightly associated with the dna methylation level of the binding site which in turn was associated with remodeling of ctcf binding in immortalized cells  <cit> .
several studies have suggested that the functional diversity of ctcf is associated with its sequence variations  <cit> , but the mechanism by which sequence variation determines its function remains unexplored.
consequently, we sought to identify and characterize the variation in ctcf binding sequences and its relation to the ctcf functional spectrum.
using a newly developed motif discovery workflow, we were able to distinguish three different ctcf binding motifs.
integrative analyses of data on transcription factors, histone modifications, chromatin conformations and gene expression across multiple cell lines suggest distinct functionalities of these three ctcf binding sequence variations.
in particular, our analysis revealed that cpg coverage and methylation status at the 12th position of the ctcf binding motifs have a marked effect on the colocalization of cohesin, which in turn implies that the variations in the ctcf binding sequence mediate different effects on chromatin interactions.
to test this assumption, we examined the effects of the three ctcf binding sequence variations in relation to chromatin interactions, and found that chromatin interactions between promoters and distal regulatory elements tends to be mediated by ctcfs that bind to motifs with higher cpg content.
RESULTS
ctcf binding sequence variations detected by a motif discovery workflow
it has been shown that the functional diversity of ctcf may be associated with the occupancy of the protein at its binding sites  <cit> .
this led us to the hypothesis that the functional diversity of ctcf may result from its binding affinity, and thus be influenced by the variations in the ctcf binding motifs.
to obtain an extensive sample of the variations of ctcf binding motifs, we revisited the genome-wide data of chromatin immunoprecipitation of ctcf followed by high-throughput sequencing  from the encode project  <cit> .
while many of the ctcf binding sites were observed to be bound across all or most cell types, ~20–50 % of ctcf sites showed some level of cell type specific binding  <cit> .
however, cell type specific ctcf binding sites have recently been shown to be less occupied than constitutive sites  <cit> , indicating that cell type specific binding of ctcf is less stable and weaker than constitutive binding, which implies the possibility that the detection of cell type specific ctcf binding events may have a higher false positive rate.
in order to obtain high confidence ctcf binding sites, we extracted ctcf chip-seq peaks that appeared in all the  <dig> tested cell lines   and regarded these as “constitutive” binding sites.
a total of  <dig>  constitutive ctcf binding sites were detected by macs  <cit>  with a threshold fdr <  <dig> .
sequences of 200 bp in length centered at the summit of each ctcf binding peak were extracted and pooled for motif discovery.
in the rest of this paper, if not stated otherwise, the analysis was performed on the aforementioned constitutive ctcf binding sites.
different from previous motif generation methods that attempt to obtain a maximal description for all pooled sequences as a whole  <cit> , our workflow searches for multiple motifs, balancing the number of sequences each motif represents and the divergent nature of the recognition motifs of a given binding protein.
in other words, our workflow detects motifs that represent the consensus of mutually exclusive subsets of pooled sequences, and the motifs consequently represent sequence subsets that could be very different.
our workflow consists of five major steps:  motif generation;  motif evaluation;  sequence elimination;  motif updating; and  a stopping criterion .fig.
<dig> motif discovery workflow and three detected ctcf sequence variations.
a cartoon of the motif discovery workflow .
b motif logo for three core ctcf consensus sequences .
the stars  in the first row indicate the most informative sites.
c average gc content at the most informative sites of the three ctcf binding sequence variations.
d distribution of three ctcf binding sequence variations
we applied this workflow to the ctcf binding sequences described above, and identified three distinct ctcf core motifs of high confidence .
the details concerning iteration and convergence of the process are shown in additional file 3: figure s <dig>  we named the three ctcf binding sequence variations as ctcf-a, ctcf-b and ctcf-c, respectively .
ctcf-a motifs constituted the largest fraction of the motifs , while  <dig> and 14 % of the ctcf binding sites contained the ctcf-b and ctcf-c motifs, respectively .
the remaining 5 % of the binding sites did not fall into any of the three categories above, and were excluded from subsequent analyses.
at the dna sequence level, the three motifs did not show substantial divergence , although the gc content in the informative sites of ctcf-a was moderately larger than in the ctcf-b and ctcf-c motifs.
we then extracted the flanking region  centered at each ctcf binding motif, and calculated the binding intensity by counting the number of chip-seq reads that mapped within the flanking region.
the binding affinity differed significantly among the three motifs , suggesting three ctcf motifs may confer distinct functionalities.
we therefore sought to investigate whether there is functional divergence among ctcf binding sites containing the three different motifs.
actively expressed genes are predominantly associated with ctcf-a
to investigate the association between gene expression activity and the occupancy of ctcf binding sites, we used a set of histone modification data from the broad institute  <cit>  available from the encode project , in combination with counts of histone marks in the flanking regions of each ctcf binding site .
we observed that regions close to ctcf-a binding sites were highly enriched for chromatin features that have been associated with active regulatory genome elements.
for example, an enrichment of h3k4me <dig> and h3k9ac is considered a strong indication of an active promoter  <cit> , and an enrichment of h3k27ac separates active from poised enhancers  <cit> .
in contrast, the vicinity of ctcf-a sites were depleted of repressive chromatin marks, such as h3k27me <dig>  <cit> .
we did not observed any significant enrichment of particular chromatin features in the vicinity of ctcf-b and ctcf-c binding sites.
this pattern was even more pronounced for tissue specific ctcf binding sites , in that tissue specific ctcf-a binding sites showed a higher frequency of active chromatin marks than constitutive ctcf-a binding sites .fig.
<dig> expression of ctcf-linked genes.
the distribution of associated gene types for the different ctcf binding sites is shown for all three cell types.
active genes are more likely to be linked with ctcf-a, compared with ctcf-b and ctcf-c binding sites.
.
the expression levels are showed in logarithmic transformation
fig.
<dig> genome-wide distribution of three ctcf binding sequence variations.
a the heatmap shows fold-enrichment of  <dig> histone modifications in the flanking regions  of the three ctcf binding site variations  .
b distribution of three ctcf binding sequence variations in promoters, enhancers and insulators as annotated by the encode project.
c the distributions of three ctcf binding sequence variations in tss flanking regions .
the pie charts show the distributions of the three ctcf binding sequence variations in the whole genome and in promoter regions ; the star  indicates a significant difference between the two by the hypergeometric test
these results prompted us to investigate the association between gene expression and the occupancy of ctcf binding sequence variations.
we divided all the coding genes into four groups based on ctcf occupancy in the flanking regions of their transcriptional start sites , referred to as ctcf-a-linked , ctcf-b-linked , ctcf-c-linked  or ‘other’ genes.
excluding chip-seq peaks that did not correspond to any of these three ctcf binding motifs, we assigned each of the  <dig> genes annotated by ensemble into a group corresponding to its nearest ctcf binding motif within a  region of its tss, or, if none, as ‘other genes’.
we examined the transcriptional activities of the genes of each group in the gm <dig>  k <dig> and helas <dig> cell lines based on rna-seq data available from encode  <cit> .
most of the examined genes were transcriptionally silent; however, in all three cell lines, we found that the expression levels of al genes were, in general, higher than those of bl  and cl genes .
taken together, these results suggest that ctcf-a binding sites are more frequently involved in active gene transcriptional regulation than the two other types of sites.
functional diversity of ctcf sequence variations
as the next step, we studied the overlap between the ctcf sites and annotations of regulatory regions.
twelve percent of all ctcf binding sites were located within 3 kb regions ranging from 2500 bp upstream to 500 bp downstream of annotated tsss of coding genes  while the remaining 88 % of the sites were evenly distributed  between intergenic and intragenic locations , which is consistent with a previous mapping  <cit> .
we further considered the overlap with regulatory regions annotated by the encode chromatin state pattern  <cit> .
specifically, we studied the overlap of ctcf binding sites with  <dig> types of genomic annotations, including active promoters, weak promoters, poised promoters, strong enhancers, poised enhancers, and insulators  <cit> .
to further refine the annotation of the chromatin states, we pinpointed promoter, enhancer and insulator elements by adding additional restrictions of location or cofactors.
for example, a dna segment with a promoter-like chromatin state is regarded as a true promoter only if it also locates within  of an annotated tss.
likewise, enhancers and insulators needed to be colocalized with p <dig>  <cit>  and cohesin  <cit> , respectively, the data on p <dig> and cohesin binding sites were obtained from the encode chip-seq dataset  <cit> .
in all three cell lines examined, the three ctcf motifs exhibited pronounced differences in inferred functionality.
particularly, ctcf-a binding sites exhibited a strong tendency to overlap with promoters , while ctcf-c binding showed enrichment for overlap with insulators .
on the other hand, ctcf-b did not show enrichment for any of the chromatin states we tested , and we have for this reason, mainly compared ctcf-a and ctcf-c in the rest of the paper.
we can make several predictions based on the inferred functionality distributions.
first, considering the enrichment of ctcf-a binding sites in promoter regions  and the strong association with gene expression activity , we would expect a higher incidence of colocalization between ctcf and transcription factor  binding.
to test this, we integrated chip-seq data for  <dig> different transcription factors available from encode   <cit>  and calculated the fold-enrichment of colocalization of these transcription factors with the ctcf binding sites compared to the input signal .
as shown in fig.
<dig>  in all three cell lines we tested, we observed that binding of ctcf to ctcf-a sites colocalized with binding of most of the tfs we tested in all three cell lines.
second, the engagement of ctcfs in mediating transcriptional insulation has been found to be coupled with cohesin  <cit> ; therefore, we can expect a subset of ctcf binding sites to be colocalized with cohesin.
in fact, the results show a strong association between ctcf-c bindings and cohesins  .
we also observed a strong association between ctcf-c bindings and lamina associated domains , which also represent a repressive chromatin environment  <cit> .fig.
<dig> fold-enrichment of functional genome elements within binding regions of the three ctcf variations in a gm <dig>  b k <dig> and c helas <dig>  cell lines
in summary, ctcf-a binding regions were enriched for active histone modifications and tended to appear within promoter regions.
as such, they presumably consist of active regulatory elements in enhancers and promoters.
in the absence of active histone marks and enrichment for insulator-like genome segments, ctcf-cs binding regions possibly function as enhancer-blocking insulators.
the sequence and dna methylation variation at the 12th position of the ctcf binding motif
we next asked what characteristics of the motifs might influence the divergence of their functionality.
given the differences in gc content among the three motifs, it is possible that differential dna methylation levels at the cpg dinucleotide sites may result in different ctcf binding affinities.
two positions in the ctcf recognition motif have been reported to show an enrichment of cpg dinucleotides, with a strong association between dna methylation level and ctcf occupancy  <cit> .
to investigate the association between dna methylation levels at the two sites and the potential functions of the three sequence variations, we extracted the cpg methylation status in  regions centered at each ctcf binding motif  <cit> .
ctcf-a sequence regions have a relatively high overall cpg content and high dna methylation levels at the 12th position.
since the overall gc content varies among the ctcf-a, −b and -c sequences, it was not surprising to find that of all cpg dinucleotides detected in the three examined cell lines , ctcf-a sequences had much higher overall cpg levels than ctcf- b and -c sequences .
in agreement with wang et al.
<cit> , we also observed ultrahigh enrichment of cpg dinucleotides at the 12th position of ctcf recognition sequences particularly in the ctcf-a subgroup .
given the high cpg level at the 12th position of the ctcf-a sequences, we would also expect a correspondingly high dna methylation level.
indeed, the dna methylation level at the 12th position of ctcf-a binding sites was relatively high compared to other cpg sites .fig.
<dig> cpg coverage and methylation status at the 12th position of the ctcf motifs.
a cpg coverage  distribution within regions  of the center of ctcf-a and ctcf-c binding sites in gm <dig> .
see additional file 11: figure s <dig> for the corresponding dna methylation distribution.
b depletion of cohesin  at ctcf binding sites.
the black and gray bars represent cohesin coverage at all ctcf binding regions and at ctcf binding sites that are methylated  at the 12th position, respectively.
c correlation  between the chip-seq read counts for cohesin  and the dna methylation level at the 12th position of the ctcf binding sites across the three cell lines.
the stars  indicate the statistical significance in each test
one puzzling observation is that both gc content and dna methylation levels appear relatively high in and around ctcf-a binding sites located in transcriptionally active regions.
however, this may in part be explained by the fact that the overall dna methylation level at the constitutive ctcf binding sites is quite low compared to inactive or unbound ctcf sequences.
when comparing the dna methylation levels of ctcf-a binding sites with control regions, which were defined as the genome segments with high confident ctcf-a motifs but with little ctcf binding signal , we found that the dna methylation levels at ctcf-a binding sites were much lower than in the control regions.
if the regulation of functional diversity among ctcf binding sequence variations is influenced by the relative dna methylation level at the 12th position, we would expect different associations between the dna methylation levels at this position and key cofactors.
to test this hypothesis, we divided the ctcf binding sites into two groups labeled as dna “methylated” and “unmethylated” according to the dna methylation level at their 12th position.
thus, if the dna methylation level at the 12th position of a ctcf binding site was found to be greater than 20 %  <cit> , this binding site would be classified as “methylated” otherwise the binding site is classified as “unmethylated”.
as previously reported  <cit> , we observed that the binding intensity of ctcf is negatively correlated with the methylation level at the 12th position of the binding sites .
we also examined the colocalization of cohesin with methylated and unmethylated ctcfs.
cohesin is one of the most important cofactors of ctcf, and exhibits distinct functionalities in the presence or absence of ctcf  <cit> .
because the peak numbers and genome wide distributions of chip-seq reads from the two subunits of cohesin  are substantially different, we took only the overlapping binding regions of rad <dig> and smc <dig> to be representative of cohesin binding sites, investigating whether the colocalization of cohesin with ctcf was associated with dna methylation level at the 12th position, we found that cohesin was highly depleted from methylated ctcf binding sites , the finding that is consistent across the three cell lines we examined .
we also tested the effect of methylation at other sites within the ctcf core binding region , however, the 12th position was the only site that showed statistically significant depletion of cohesin occupancy in all three cell lines.
moreover, estimates of the cohesin binding affinity  in gm <dig> and helas <dig> cell lines were negatively correlated with the dna methylation level at the 12th position of the associated ctcf binding sequences .
the data thus indicate that the co-binding of cohesin with ctcf is, to some degree, negatively related to dna methylation at the 12th position of the ctcf binding site.
renda and colleagues demonstrated that high-affinity binding to a 12 bp variation of the cfct consensus sequence involved only 4– <dig> specific zinc-fingers of the ctcf  <cit> .
ren and colleagues reported that 17 % of all evolutionary nucleotide changes in the ctcf binding sites took place as c-to-t substitutions as a unique nucleotide change at the 12th position  <cit> .
our results, when combined with those of others  <cit> , imply that the methylation level at the 12th position of ctcf binding sequence may alter the binding environment, resulting in different zinc-finger binding and, in turn, recruiting different cofactors that ultimately leads to divergent functionalities.
functional diversity of ctcf binding sequence motifs is reflected in the genome 3d structure
because the ctcf protein is important for mediating chromatin-chromatin interactions  <cit> , we next asked what connection might exist between ctcf binding site variation and the dna loops in which they are involved.
to address this question, we revisited the chromatin interactions database in encode, which is based on data obtained with the chromosome conformation capture carbon copy  method  <cit> .
these data contain chromatin interactions between  <dig> tss-containing fragments and  <dig> distal restriction fragments covering the encode pilot project regions representing 1 %  of the genome for the three cell lines examined above .
we first used the sequence scanning tool fimo with default settings to search each distal fragment for ctcf-a, −b and -c motifs  <cit> .
each fragment with the ctcf binding type was tagged by which had the highest confidence for this fragment .
in all three cell lines, we observed that fragments tagged with ctcf-a were enriched for tss-distal chromatin interactions, while ctcf-c-linked fragments showed depletion of tss-distal interactions compare to ctcf-a-linked fragments .
ctcf-b, however, was enriched for tss-distal interactions in helas <dig> and gm <dig>  but not in k <dig>  in agreement with recent works showing that ctcf is involved in mediating long-range interactions  <cit> , distal fragments lacking ctcf occupation also showed a low tendency to form loops.
we further examined the correlation of  <dig> histone modifications ,  <dig> dna-binding factors  from encode and  <dig> chromatin domains  with loop formation , and observed that tss-distal looping is more strongly correlated with ctcf-a binding than with most of these other biological elements .fig.
<dig> association between the three ctcf binding sequence variations and functional elements.
a enrichment of tss-distal chromatin interactions in the three ctcf binding sequence variations.
the symbols, “*” and “**” indicate the significance of hypergeometric tests for p-value < 1e- <dig> and p-value < 1e- <dig>  respectively.
b distribution correlations between tss-distal chromatin interactions and  <dig> genomic features in gm <dig> cells.
.
the panel shows the distribution of pearson’s correlation coefficients
ctcf-a binding sites appeared more involved in the mediation of dna looping between tsss and enhancer-containing distal fragments than did ctcf-c binding sites.
we found that cohesin  was consistently and highly positively correlated with tss-distal looping in all three cell lines.
although cohesin was previously thought to act as insulator when colocalized with ctcf  <cit> , our results are in agreement with recent genome wide studies showing an import role for cohesin in establishing enhancer-promoter interactions  <cit> .
we also observed enrichment of the nfya and nfyb factors  around the ctcf-a binding regions.
these two proteins have been shown to interact with the ccaat/enhancer binding protein  which regulates gene expression  <cit> ; therefore we further examined the enrichment of regulatory genome elements in the distal fragments.
compared to all fragments, the distal fragments mainly occupied by ctcf-a sequences exhibited higher enrichment for enhancer-related marks , but not for h3k4me <dig>
CONCLUSIONS
to determine divergent dna recognition motifs of ctcf in various functional genomics contexts, we developed a novel motif discovery workflow focusing on the balance between the number of sequences a motif represents and the internal divergence within the sequence subgroups.
by applying the workflow to the encode chip-seq data set, we identified three ctcf core motif variations.
the expression activity patterns of the genes flanking the three ctcf motifs showed significant correlations with the gc content and the cpg enrichment of the motifs.
we also detected a strong association among the presence of the ctcf-a motif and enhancers and promoters defined by histone marks, and we further demonstrated that this association was supported by chromatin structural data from chromatin conformation capture based experiments.
the functional divergence of the motifs was further associated with possible genetic or epigenetic variations, in particular with the cpg dinucleotide coverage at the 12th position of the core binding motifs and the relative dna methylation level at this site, and the latter two features were also strongly and negatively correlated with the colocalization of ctcf with cohesin, the key cofactor of ctcf when it functions as an insulator protein.
these results suggest that the variation in dna methylation level at a single cpg site of the ctcf’s recognition motif has a determining influence on the divergence of its functions.
alternative preferences for critical cofactors, e.g.
cohesin, among the different motifs suggest potentially multiple molecular mechanisms for the ctcf functionalities.
the workflow we have introduced provides a new analytical tool for the studies of multifunctional dna binding proteins, particularly for those whose functional classification is not yet clearly defined.

 the ccctc-binding factor is an 11-zinc-finger protein that is functionally conserved in vertebrates, insects and nematodes. it has been shown to be involved in various critical biological processes and has long been regarded as a versatile regulator. ctcf knockout mice exhibit early embryonic lethality prior to implantation. in adult organisms, the ctcf protein is ubiquitously expressed across most tissues, and its expression levels vary in a cell type-specific manner, suggesting a role in maintaining phenotypic diversity. furthermore, ctcf depletion results in dysregulated transcription of hundreds of genes in oocytes and dramatically deregulates cell-cycle progression during t lymphocyte lineage commitment within the thymus. other studies have suggested that ctcf impacts cellular activity by playing diverse roles in gene regulation, including promoter activation/repression, genomic imprinting, enhancer blocking, and, most recently, long-range chromatin interactions. recent genome-wide mapping of ctcf occupancy in multiple cell lines has enabled the identification of ctcf binding sequences. several hypotheses have been proposed to explain the uniquely versatile characteristics of ctcf based on the underlying binding motif sequence. the most widespread is the zinc-finger model, in which the capacity of ctcf to confer vastly different functions has been attributed to the interplay between the zinc-finger engagement and the underlying sequence. early studies in which zinc-fingers were depleted in a stepwise manner reported involvement of multiple zinc-fingers within a broad ~50 bp sequence, indicating that ctcf binding may be partially stabilized by interplay between peripheral zinc-fingers and other cofactors. following this, a ~11–15 bp core consensus sequence was identified as bound by 4– central zinc-fingers. a variation at the 12th bp position of the consensus sequence was found to be tightly associated with the dna methylation level of the binding site which in turn was associated with remodeling of ctcf binding in immortalized cells. several studies have suggested that the functional diversity of ctcf is associated with its sequence variations, but the mechanism by which sequence variation determines its function remains unexplored. consequently, we sought to identify and characterize the variation in ctcf binding sequences and its relation to the ctcf functional spectrum. using a newly developed motif discovery workflow, we were able to distinguish three different ctcf binding motifs. integrative analyses of data on transcription factors, histone modifications, chromatin conformations and gene expression across multiple cell lines suggest distinct functionalities of these three ctcf binding sequence variations. in particular, our analysis revealed that cpg coverage and methylation status at the 12th position of the ctcf binding motifs have a marked effect on the colocalization of cohesin, which in turn implies that the variations in the ctcf binding sequence mediate different effects on chromatin interactions. to test this assumption, we examined the effects of the three ctcf binding sequence variations in relation to chromatin interactions, and found that chromatin interactions between promoters and distal regulatory elements tends to be mediated by ctcfs that bind to motifs with higher cpg content.  ctcf binding sequence variations detected by a motif discovery workflow it has been shown that the functional diversity of ctcf may be associated with the occupancy of the protein at its binding sites. this led us to the hypothesis that the functional diversity of ctcf may result from its binding affinity, and thus be influenced by the variations in the ctcf binding motifs. to obtain an extensive sample of the variations of ctcf binding motifs, we revisited the genome-wide data of chromatin immunoprecipitation of ctcf followed by high-throughput sequencing from the encode project. while many of the ctcf binding sites were observed to be bound across all or most cell types, ~20–50 % of ctcf sites showed some level of cell type specific binding. however, cell type specific ctcf binding sites have recently been shown to be less occupied than constitutive sites, indicating that cell type specific binding of ctcf is less stable and weaker than constitutive binding, which implies the possibility that the detection of cell type specific ctcf binding events may have a higher false positive rate. in order to obtain high confidence ctcf binding sites, we extracted ctcf chip-seq peaks that appeared in all the tested cell lines and regarded these as “constitutive” binding sites. a total of constitutive ctcf binding sites were detected by macs with a threshold fdr < . sequences of 200 bp in length centered at the summit of each ctcf binding peak were extracted and pooled for motif discovery. in the rest of this paper, if not stated otherwise, the analysis was performed on the aforementioned constitutive ctcf binding sites. different from previous motif generation methods that attempt to obtain a maximal description for all pooled sequences as a whole, our workflow searches for multiple motifs, balancing the number of sequences each motif represents and the divergent nature of the recognition motifs of a given binding protein. in other words, our workflow detects motifs that represent the consensus of mutually exclusive subsets of pooled sequences, and the motifs consequently represent sequence subsets that could be very different. our workflow consists of five major steps: motif generation; motif evaluation; sequence elimination; motif updating; and a stopping criterion.fig. motif discovery workflow and three detected ctcf sequence variations. a cartoon of the motif discovery workflow. b motif logo for three core ctcf consensus sequences. the stars in the first row indicate the most informative sites. c average gc content at the most informative sites of the three ctcf binding sequence variations. d distribution of three ctcf binding sequence variations we applied this workflow to the ctcf binding sequences described above, and identified three distinct ctcf core motifs of high confidence. the details concerning iteration and convergence of the process are shown in additional file 3: figure s we named the three ctcf binding sequence variations as ctcf-a, ctcf-b and ctcf-c, respectively. ctcf-a motifs constituted the largest fraction of the motifs, while and 14 % of the ctcf binding sites contained the ctcf-b and ctcf-c motifs, respectively. the remaining 5 % of the binding sites did not fall into any of the three categories above, and were excluded from subsequent analyses. at the dna sequence level, the three motifs did not show substantial divergence, although the gc content in the informative sites of ctcf-a was moderately larger than in the ctcf-b and ctcf-c motifs. we then extracted the flanking region centered at each ctcf binding motif, and calculated the binding intensity by counting the number of chip-seq reads that mapped within the flanking region. the binding affinity differed significantly among the three motifs, suggesting three ctcf motifs may confer distinct functionalities. we therefore sought to investigate whether there is functional divergence among ctcf binding sites containing the three different motifs. actively expressed genes are predominantly associated with ctcf-a to investigate the association between gene expression activity and the occupancy of ctcf binding sites, we used a set of histone modification data from the broad institute available from the encode project, in combination with counts of histone marks in the flanking regions of each ctcf binding site. we observed that regions close to ctcf-a binding sites were highly enriched for chromatin features that have been associated with active regulatory genome elements. for example, an enrichment of h3k4me and h3k9ac is considered a strong indication of an active promoter, and an enrichment of h3k27ac separates active from poised enhancers. in contrast, the vicinity of ctcf-a sites were depleted of repressive chromatin marks, such as h3k27me. we did not observed any significant enrichment of particular chromatin features in the vicinity of ctcf-b and ctcf-c binding sites. this pattern was even more pronounced for tissue specific ctcf binding sites, in that tissue specific ctcf-a binding sites showed a higher frequency of active chromatin marks than constitutive ctcf-a binding sites.fig. expression of ctcf-linked genes. the distribution of associated gene types for the different ctcf binding sites is shown for all three cell types. active genes are more likely to be linked with ctcf-a, compared with ctcf-b and ctcf-c binding sites.. the expression levels are showed in logarithmic transformation fig. genome-wide distribution of three ctcf binding sequence variations. a the heatmap shows fold-enrichment of histone modifications in the flanking regions of the three ctcf binding site variations. b distribution of three ctcf binding sequence variations in promoters, enhancers and insulators as annotated by the encode project. c the distributions of three ctcf binding sequence variations in tss flanking regions. the pie charts show the distributions of the three ctcf binding sequence variations in the whole genome and in promoter regions; the star indicates a significant difference between the two by the hypergeometric test these results prompted us to investigate the association between gene expression and the occupancy of ctcf binding sequence variations. we divided all the coding genes into four groups based on ctcf occupancy in the flanking regions of their transcriptional start sites, referred to as ctcf-a-linked, ctcf-b-linked, ctcf-c-linked or ‘other’ genes. excluding chip-seq peaks that did not correspond to any of these three ctcf binding motifs, we assigned each of the genes annotated by ensemble into a group corresponding to its nearest ctcf binding motif within a region of its tss, or, if none, as ‘other genes’. we examined the transcriptional activities of the genes of each group in the gm k and helas cell lines based on rna-seq data available from encode. most of the examined genes were transcriptionally silent; however, in all three cell lines, we found that the expression levels of al genes were, in general, higher than those of bl and cl genes. taken together, these results suggest that ctcf-a binding sites are more frequently involved in active gene transcriptional regulation than the two other types of sites. functional diversity of ctcf sequence variations as the next step, we studied the overlap between the ctcf sites and annotations of regulatory regions. twelve percent of all ctcf binding sites were located within 3 kb regions ranging from 2500 bp upstream to 500 bp downstream of annotated tsss of coding genes while the remaining 88 % of the sites were evenly distributed between intergenic and intragenic locations, which is consistent with a previous mapping. we further considered the overlap with regulatory regions annotated by the encode chromatin state pattern. specifically, we studied the overlap of ctcf binding sites with types of genomic annotations, including active promoters, weak promoters, poised promoters, strong enhancers, poised enhancers, and insulators. to further refine the annotation of the chromatin states, we pinpointed promoter, enhancer and insulator elements by adding additional restrictions of location or cofactors. for example, a dna segment with a promoter-like chromatin state is regarded as a true promoter only if it also locates within of an annotated tss. likewise, enhancers and insulators needed to be colocalized with p and cohesin, respectively, the data on p and cohesin binding sites were obtained from the encode chip-seq dataset. in all three cell lines examined, the three ctcf motifs exhibited pronounced differences in inferred functionality. particularly, ctcf-a binding sites exhibited a strong tendency to overlap with promoters, while ctcf-c binding showed enrichment for overlap with insulators. on the other hand, ctcf-b did not show enrichment for any of the chromatin states we tested, and we have for this reason, mainly compared ctcf-a and ctcf-c in the rest of the paper. we can make several predictions based on the inferred functionality distributions. first, considering the enrichment of ctcf-a binding sites in promoter regions and the strong association with gene expression activity, we would expect a higher incidence of colocalization between ctcf and transcription factor binding. to test this, we integrated chip-seq data for different transcription factors available from encode and calculated the fold-enrichment of colocalization of these transcription factors with the ctcf binding sites compared to the input signal. as shown in fig.  in all three cell lines we tested, we observed that binding of ctcf to ctcf-a sites colocalized with binding of most of the tfs we tested in all three cell lines. second, the engagement of ctcfs in mediating transcriptional insulation has been found to be coupled with cohesin; therefore, we can expect a subset of ctcf binding sites to be colocalized with cohesin. in fact, the results show a strong association between ctcf-c bindings and cohesins. we also observed a strong association between ctcf-c bindings and lamina associated domains, which also represent a repressive chromatin environment.fig. fold-enrichment of functional genome elements within binding regions of the three ctcf variations in a gm b k and c helas cell lines in summary, ctcf-a binding regions were enriched for active histone modifications and tended to appear within promoter regions. as such, they presumably consist of active regulatory elements in enhancers and promoters. in the absence of active histone marks and enrichment for insulator-like genome segments, ctcf-cs binding regions possibly function as enhancer-blocking insulators. the sequence and dna methylation variation at the 12th position of the ctcf binding motif we next asked what characteristics of the motifs might influence the divergence of their functionality. given the differences in gc content among the three motifs, it is possible that differential dna methylation levels at the cpg dinucleotide sites may result in different ctcf binding affinities. two positions in the ctcf recognition motif have been reported to show an enrichment of cpg dinucleotides, with a strong association between dna methylation level and ctcf occupancy. to investigate the association between dna methylation levels at the two sites and the potential functions of the three sequence variations, we extracted the cpg methylation status in regions centered at each ctcf binding motif. ctcf-a sequence regions have a relatively high overall cpg content and high dna methylation levels at the 12th position. since the overall gc content varies among the ctcf-a, −b and -c sequences, it was not surprising to find that of all cpg dinucleotides detected in the three examined cell lines, ctcf-a sequences had much higher overall cpg levels than ctcf- b and -c sequences. in agreement with wang et al., we also observed ultrahigh enrichment of cpg dinucleotides at the 12th position of ctcf recognition sequences particularly in the ctcf-a subgroup. given the high cpg level at the 12th position of the ctcf-a sequences, we would also expect a correspondingly high dna methylation level. indeed, the dna methylation level at the 12th position of ctcf-a binding sites was relatively high compared to other cpg sites.fig. cpg coverage and methylation status at the 12th position of the ctcf motifs. a cpg coverage distribution within regions of the center of ctcf-a and ctcf-c binding sites in gm. see additional file 11: figure s for the corresponding dna methylation distribution. b depletion of cohesin at ctcf binding sites. the black and gray bars represent cohesin coverage at all ctcf binding regions and at ctcf binding sites that are methylated at the 12th position, respectively. c correlation between the chip-seq read counts for cohesin and the dna methylation level at the 12th position of the ctcf binding sites across the three cell lines. the stars indicate the statistical significance in each test one puzzling observation is that both gc content and dna methylation levels appear relatively high in and around ctcf-a binding sites located in transcriptionally active regions. however, this may in part be explained by the fact that the overall dna methylation level at the constitutive ctcf binding sites is quite low compared to inactive or unbound ctcf sequences. when comparing the dna methylation levels of ctcf-a binding sites with control regions, which were defined as the genome segments with high confident ctcf-a motifs but with little ctcf binding signal, we found that the dna methylation levels at ctcf-a binding sites were much lower than in the control regions. if the regulation of functional diversity among ctcf binding sequence variations is influenced by the relative dna methylation level at the 12th position, we would expect different associations between the dna methylation levels at this position and key cofactors. to test this hypothesis, we divided the ctcf binding sites into two groups labeled as dna “methylated” and “unmethylated” according to the dna methylation level at their 12th position. thus, if the dna methylation level at the 12th position of a ctcf binding site was found to be greater than 20 %, this binding site would be classified as “methylated” otherwise the binding site is classified as “unmethylated”. as previously reported, we observed that the binding intensity of ctcf is negatively correlated with the methylation level at the 12th position of the binding sites. we also examined the colocalization of cohesin with methylated and unmethylated ctcfs. cohesin is one of the most important cofactors of ctcf, and exhibits distinct functionalities in the presence or absence of ctcf. because the peak numbers and genome wide distributions of chip-seq reads from the two subunits of cohesin are substantially different, we took only the overlapping binding regions of rad and smc to be representative of cohesin binding sites, investigating whether the colocalization of cohesin with ctcf was associated with dna methylation level at the 12th position, we found that cohesin was highly depleted from methylated ctcf binding sites, the finding that is consistent across the three cell lines we examined. we also tested the effect of methylation at other sites within the ctcf core binding region, however, the 12th position was the only site that showed statistically significant depletion of cohesin occupancy in all three cell lines. moreover, estimates of the cohesin binding affinity in gm and helas cell lines were negatively correlated with the dna methylation level at the 12th position of the associated ctcf binding sequences. the data thus indicate that the co-binding of cohesin with ctcf is, to some degree, negatively related to dna methylation at the 12th position of the ctcf binding site. renda and colleagues demonstrated that high-affinity binding to a 12 bp variation of the cfct consensus sequence involved only 4– specific zinc-fingers of the ctcf. ren and colleagues reported that 17 % of all evolutionary nucleotide changes in the ctcf binding sites took place as c-to-t substitutions as a unique nucleotide change at the 12th position. our results, when combined with those of others, imply that the methylation level at the 12th position of ctcf binding sequence may alter the binding environment, resulting in different zinc-finger binding and, in turn, recruiting different cofactors that ultimately leads to divergent functionalities. functional diversity of ctcf binding sequence motifs is reflected in the genome 3d structure because the ctcf protein is important for mediating chromatin-chromatin interactions, we next asked what connection might exist between ctcf binding site variation and the dna loops in which they are involved. to address this question, we revisited the chromatin interactions database in encode, which is based on data obtained with the chromosome conformation capture carbon copy method. these data contain chromatin interactions between tss-containing fragments and distal restriction fragments covering the encode pilot project regions representing 1 % of the genome for the three cell lines examined above. we first used the sequence scanning tool fimo with default settings to search each distal fragment for ctcf-a, −b and -c motifs. each fragment with the ctcf binding type was tagged by which had the highest confidence for this fragment. in all three cell lines, we observed that fragments tagged with ctcf-a were enriched for tss-distal chromatin interactions, while ctcf-c-linked fragments showed depletion of tss-distal interactions compare to ctcf-a-linked fragments. ctcf-b, however, was enriched for tss-distal interactions in helas and gm but not in k in agreement with recent works showing that ctcf is involved in mediating long-range interactions, distal fragments lacking ctcf occupation also showed a low tendency to form loops. we further examined the correlation of histone modifications, dna-binding factors from encode and chromatin domains with loop formation, and observed that tss-distal looping is more strongly correlated with ctcf-a binding than with most of these other biological elements.fig. association between the three ctcf binding sequence variations and functional elements. a enrichment of tss-distal chromatin interactions in the three ctcf binding sequence variations. the symbols, “*” and “**” indicate the significance of hypergeometric tests for p-value < 1e- and p-value < 1e- respectively. b distribution correlations between tss-distal chromatin interactions and genomic features in gm cells.. the panel shows the distribution of pearson’s correlation coefficients ctcf-a binding sites appeared more involved in the mediation of dna looping between tsss and enhancer-containing distal fragments than did ctcf-c binding sites. we found that cohesin was consistently and highly positively correlated with tss-distal looping in all three cell lines. although cohesin was previously thought to act as insulator when colocalized with ctcf, our results are in agreement with recent genome wide studies showing an import role for cohesin in establishing enhancer-promoter interactions. we also observed enrichment of the nfya and nfyb factors around the ctcf-a binding regions. these two proteins have been shown to interact with the ccaat/enhancer binding protein which regulates gene expression; therefore we further examined the enrichment of regulatory genome elements in the distal fragments. compared to all fragments, the distal fragments mainly occupied by ctcf-a sequences exhibited higher enrichment for enhancer-related marks, but not for h3k4me  to determine divergent dna recognition motifs of ctcf in various functional genomics contexts, we developed a novel motif discovery workflow focusing on the balance between the number of sequences a motif represents and the internal divergence within the sequence subgroups. by applying the workflow to the encode chip-seq data set, we identified three ctcf core motif variations. the expression activity patterns of the genes flanking the three ctcf motifs showed significant correlations with the gc content and the cpg enrichment of the motifs. we also detected a strong association among the presence of the ctcf-a motif and enhancers and promoters defined by histone marks, and we further demonstrated that this association was supported by chromatin structural data from chromatin conformation capture based experiments. the functional divergence of the motifs was further associated with possible genetic or epigenetic variations, in particular with the cpg dinucleotide coverage at the 12th position of the core binding motifs and the relative dna methylation level at this site, and the latter two features were also strongly and negatively correlated with the colocalization of ctcf with cohesin, the key cofactor of ctcf when it functions as an insulator protein. these results suggest that the variation in dna methylation level at a single cpg site of the ctcf’s recognition motif has a determining influence on the divergence of its functions. alternative preferences for critical cofactors, e.g. cohesin, among the different motifs suggest potentially multiple molecular mechanisms for the ctcf functionalities. the workflow we have introduced provides a new analytical tool for the studies of multifunctional dna binding proteins, particularly for those whose functional classification is not yet clearly defined. 
supported by transcriptomic, epigenomic and chromatin-interactomic data, we show that the functional diversity of the ctcf binding motifs is strongly associated with their gc content, cpg dinucleotide coverage and relative dna methylation level at the 12th position of the motifs. further analysis suggested that the co-localization of cohesin, the key factor in cohesion of sister chromatids, is negatively correlated with the cpg coverage and the relative dna methylation level at the 12th position. however, the definitive characteristics of the ctcf binding motif required for its functional diversity still remains elusive. these results demonstrate the existence of definitive ctcf binding motifs corresponding to ctcf’s diverse functions, and that the functional diversity of the motifs is strongly associated with genetic and epigenetic features at the 12th position of the motifs. finally, we present evidences for a hypothetical model in which chromatin interactions between promoters and distal regulatory regions are likely mediated by ctcfs binding to sequences with high cpg. here, we describe a new motif discovery workflow by which we have identified three ctcf binding motif variations with highly divergent functionalities.
 the ccctc-binding factor has diverse regulatory functions. however, the definitive characteristics of the ctcf binding motif required for its functional diversity still remains elusive.  here, we describe a new motif discovery workflow by which we have identified three ctcf binding motif variations with highly divergent functionalities. supported by transcriptomic, epigenomic and chromatin-interactomic data, we show that the functional diversity of the ctcf binding motifs is strongly associated with their gc content, cpg dinucleotide coverage and relative dna methylation level at the 12th position of the motifs. further analysis suggested that the co-localization of cohesin, the key factor in cohesion of sister chromatids, is negatively correlated with the cpg coverage and the relative dna methylation level at the 12th position. finally, we present evidences for a hypothetical model in which chromatin interactions between promoters and distal regulatory regions are likely mediated by ctcfs binding to sequences with high cpg.  these results demonstrate the existence of definitive ctcf binding motifs corresponding to ctcf’s diverse functions, and that the functional diversity of the motifs is strongly associated with genetic and epigenetic features at the 12th position of the motifs. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords ctcfbinding motifdna methylationchromatin interactionissue-copyright-statement© the author 2015 