Analysis of variation at transcription factor binding sites in Drosophila and humans
© Spivakov et al.; licensee BioMed Central Ltd. 2012
Received: 28 March 2012
Accepted: 8 June 2012
Published: 5 September 2012
Advances in sequencing technology have boosted population genomics and made it possible to map the positions of transcription factor binding sites (TFBSs) with high precision. Here we investigate TFBS variability by combining transcription factor binding maps generated by ENCODE, modENCODE, our previously published data and other sources with genomic variation data for human individuals and Drosophila isogenic lines.
We introduce a metric of TFBS variability that takes into account changes in motif match associated with mutation and makes it possible to investigate TFBS functional constraints instance-by-instance as well as in sets that share common biological properties. We also take advantage of the emerging per-individual transcription factor binding data to show evidence that TFBS mutations, particularly at evolutionarily conserved sites, can be efficiently buffered to ensure coherent levels of transcription factor binding.
Our analyses provide insights into the relationship between individual and interspecies variation and show evidence for the functional buffering of TFBS mutations in both humans and flies. In a broad perspective, these results demonstrate the potential of combining functional genomics and population genetics approaches for understanding gene regulation.
Gene expression is tightly controlled by transcription factors (TFs) that are recruited to DNA cis-regulatory modules (CRMs). Many TFs have well-documented sequence preferences for their binding sites (transcription factor binding sites (TFBSs)) . However, in contrast to the startling simplicity of the amino acid code, the 'regulatory code' at CRMs has a more ambiguous relationship between sequence and function. Chromatin immunoprecipitation (ChIP) coupled with genome-wide analyses have made it possible to map TF binding positions globally in vivo, which in some cases can serve as good predictors of CRM transcriptional outputs [2–4]. At the same time, these analyses often cannot explain the exact rules underlying TF binding to a given sequence, and functional prediction based on sequence alone has had limited success, in particular in mammalian systems .
Evolutionary analyses across species have proven to be a powerful approach in elucidating the functional constraints of DNA elements, in particular protein-coding genes, but are less interpretable in the context of CRM architecture [6, 7]. In part, this is due to the fact that CRMs often have a 'modular', rather than 'base-by-base', conservation that may escape detection by conventional alignment-based approaches . Moreover, conservation in DNA binding profiles can be detected even without apparent DNA sequence constraint . Even at the level of individual TFBSs, differences in sequence may be hard to interpret - as such differences, for example, may reflect evolutionary 'fine-tuning' to species-specific factors to preserve uniform outputs rather than signifying a lack of functional constraint [6, 10–12].
A complementary way to analyze the relationship between sequence and function is to explore intra-species (that is, polymorphic) variation of functional elements. Variation at DNA regulatory elements is relatively common and at least a fraction of it falls directly at TFBSs [13, 14]. While some regulatory variants have been associated with major changes in transcription factor binding [15–17], gene expression [18, 19] and disease phenotypes , many others do not result in apparent aberrations in function. This difference in itself suggests that analyzing TFBS variability in the context of the same species may lead to insights into cis-regulatory logic. For example, high tolerance of a binding site to deleterious variation may indicate that such variation is effectively 'buffered', either at the level of the same regulatory module or elsewhere in the system.
Until recently, large-scale population genomics studies of metazoan TFBSs were unthinkable because of the limited number of available genotypes and global TF binding profiles. However, advances in sequencing technology have paved the way for high-throughput efforts, such as the human 1000 Genomes project  and Drosophila Genetic Reference Panel (DGRP) , that are making available an increasing number of individual genomes originating from the same population. Combining these data with the binding maps of dozens of TFs in both species generated by the Encyclopedia of DNA Elements (ENCODE) for human , and modENCODE and other published sources in Drosophila [2, 24–30] has provided an unprecedented resource for analyzing TFBS functional constraints.
Here we use three different approaches to take advantage of variation data in this context. First, we analyze TFBSs position-by-position to confirm that the levels of variation are generally consistent with TFBSs functional constraints predicted by their position weight matrix (PWM) models and highlight some intriguing exceptions. Next, we draw inspiration from Haldane's  and Muller's  genetic load model to devise a metric of TFBS variation that takes into account the loss of PWM match score associated with a mutation and makes it possible to investigate per-instance TFBS functional constraints. Finally, we take advantage of per-individual binding maps for a human transcription factor (CTCF) to highlight the 'buffering' of genetic variation at TFBSs at the level of binding, particularly in evolutionarily conserved regions.
We aim to analyze TFBS functional constraints using the binding data generated by the ENCODE, modENCODE and published sources. Prior to these global analyses, however, we first examined the relationship between binding sites' match to consensus, their conservation and variation using three well-characterized Drosophila TFs, Twist (Twi), Biniou (Bin) and Tinman (Tin), which have large numbers of TFBSs whose general occupancy is predictive of specific spatio-temporal activity . The discovered PWMs for these TFs from both in vitro and in vivo studies are good predictors for their binding  and their binding sites show an appreciable level of variation, presumably much of which is deleterious but not lethal. For these TFs, 24 to 28% of the bound sites overlapped with SNPs identified by the DGRP  in 162 isogenic lines of Drosophila melanogaster (hereafter we refer to variation across these lines as 'individual variation'). As expected, variation at the same sequences detected outside of TF-bound regions (that is, at potentially random motif matches) was even higher, with 35% of them containing known SNPs (Fisher test, P < 1e-50 compared to the bound sites).
PWMs are an established way of representing the sequence preferences of TFBSs, with PWM match scores reflecting the similarity of a given sequence to the hypothetical 'ideal' binding site for a given TF . To study the relationship between PWM scores and variation, we compared the variation properties of Twi, Bin and Tin motifs at three score ranges ('strong', 'medium' and 'weak' scoring). Weaker (that is, potentially 'less optimal') motifs generally showed higher levels of individual variation (Figure 1c), as further confirmed using only the strongest scoring sites from each bound region to reduce the contribution of non-functional motif matches (Figure S1B in Additional file 1). This result is consistent with the expectation that selection would predominantly work towards increasing TFBSs' match to consensus . We revisit this question more formally later in the study.
As well as looking across the entire PWM, we can consider each motif position in turn. Consistent with previous findings for other TFs in yeast  and Drosophila , cross-species variation at Twi, Bin and Tin motif positions strongly anti-correlated with their information content (Figure 1a; Figure S2 in Additional file 1). Variation across individuals also anti-correlated with positional information content (Figure 1d), confirming the general link between evolutionary conservation and population diversity . There are, however, some interesting exceptions. For example, positions 6, 7 and 12 of the Twi motif are less varied in the population than would be expected from their information content (Figure 1d, left panel). These positions correspond to the 'spacer' region of the CANNTG E-box consensus motif recruiting basic helix-loop-helix (bHLH) proteins, for which specific sequence preferences were documented depending on specific dimerization partners . Similarly, we found the first two positions of the Bin motif to be highly constrained despite their very low information content (Figure 1d, middle panel), suggesting that these positions may also be subject to specific restrictions depending on the cis-regulatory context of each motif instance. From this analysis we conclude that PWMs that have a strong correlation between information content and cross-species conservation are likely to be good descriptors of TF sequence binding preferences in a population context.
We do not assume that TFBS load at a given site reduces an individual's biological fitness. Rather, we argue that binding sites that tolerate a higher load are less functionally constrained. This approach, although undoubtedly a crude one, makes it possible to consistently estimate TFBS constraints for different TFs and even different organisms and ask why TFBS mutations are tolerated differently in different contexts. Conceptual and statistical considerations associated with TFBS load are discussed at length in Materials and methods; here we will only outline several major points. First, since binding events limited to minor alleles are likely to be overlooked by a single-genome ChIP analysis, we compute the decrease in PWM match score relative to the major and not the highest-scoring allele as in the 'classic' genetic load metric. In addition, since we focus on the deleterious effects of variation, we have assumed that mutations yielding increased PWM match scores have a load of zero. We avoided the use of negative load values for these 'gain-of-score' mutations, as it is possible that such mutations will often be near-neutral, while in some cases they may even be deleterious.
Most of the analyzed TFBSs have no detected variation, in particular in human, and therefore a zero load. This affects the statistical power, making it challenging to examine many TFs one-by-one. However, analyzing the data globally for all included TFs in each organism has allowed us to identify a number of significant trends, as presented below. Technically, the high proportion of sites with no detected variation also leads to a considerable zero-inflation of TFBS load distributions, which violates the assumptions of conventional significance tests. Therefore, instead we estimate significance by using permutation tests, as further described in Materials and methods. For the same reason, we also chose to present average (more precisely, trimmed mean) TFBS load values in many comparative analyses as a metric that reflects both the frequency of variation (that is, zero versus non-zero load) and the intensity of its effect (that is, the distribution of non-zero load).
We first asked whether motif load would be able to detect the expected link between evolutionary and individual variation. We used a published metric, branch length score (BLS) , to characterize the evolutionary conservation of a motif instance. This metric utilizes both a PWM-based model of the conservation of bases and allows for motif movement. Reassuringly, mutational load correlated with BLS in both species, with evolutionarily non-conserved motifs (BLS = 0) showing by far the highest degree of variation in the population (Figure 3b). At the same time, approximately 40% of human and fly TFBSs with an appreciable load (L > 5e-3) still mapped to reasonably conserved sites (BLS > 0.2, approximately 50% percentile in both organisms), demonstrating that score-reducing mutations at evolutionarily preserved sequences can be tolerated in these populations.
Earlier in the study we have shown evidence that 'weaker' motifs (that is, those with a poorer PWM match) are more prone to variation, suggesting that they are less functionally constrained. Weaker sites have many more possible variants with similar match scores, while mutations at stronger sites are less likely to preserve their match. Motif load is based on the decrease in PWM score associated with mutations and not sequence variation per se and is therefore more 'protected' from this bias. Using this metric, we confirmed our original findings, suggesting that TFBSs with higher PWM scores are generally more functionally constrained compared to 'weaker' sites (Figure 3c). The fraction of detected sites mapping to bound regions remained similar across the whole analyzed score range, suggesting that this relationship is unlikely to be an artifact of higher false-positive rates at 'weaker' sites (Figure S4A in Additional file 1). This global observation, however, does not rule out the possibility that a weaker match at some sites is specifically preserved to ensure dose-specific TF binding. This may be the case, for example, for Drosophila Bric-à-brac motifs, which exhibited no correlation between motif load and PWM score (Figure S4B in Additional file 1), consistent with the known dosage-dependent function of Bric-à-brac in embryo patterning .
We then used motif load to address whether TFBSs proximal to transcription start sites (TSSs) are more constrained compared to more distant regulatory regions. We found this to be the case in human, but not Drosophila (Figure 3d; see Discussion). CTCF binding sites in both species were a notable exception, tolerating the lowest mutational load at locations 500 bp to 1 kb from TSSs, but not closer to the TSS (Figure 3d, bottom panel), suggesting that the putative role of CTCF in establishing chromatin domains  is particularly important in proximity of gene promoters.
Taken together, CTCF binding data for multiple individuals show that mutations can be buffered to maintain the levels of binding signal, particularly at highly conserved sites, and this effect cannot be explained solely by the flexibility of CTCF's sequence consensus. We asked whether mechanisms potentially accountable for such buffering would also affect the relationship between sequence and binding in the absence of mutations. Training an interaction linear model across the whole set of mapped CTCF binding sites revealed that conservation consistently weakens the relationship between PWM score and the binding intensity (P = 1.9e-7; Figure 5e). Thus, CTCF binding to evolutionarily conserved sites may generally have a reduced dependence on sequence.
Deciphering the cis-regulatory 'logic' of gene regulation is one of the biggest challenges genomics faces today. Understanding the functional constraints of regulatory elements across species has been the focus of much 'evo-devo' research, leading to many exciting insights, such as the preservation of CRM function without a base-to-base preservation of sequence [9–11] and the impact of protein-protein interactions . Variation across individuals presents a snapshot of 'evolution in action', giving access to potentially suboptimal alleles without having to resort to artificial perturbation, and are a promising resource for population functional genomics studies as well as more formal association analyses. Such 'pop-fun' approaches will complement the insights obtained from 'evo-devo' studies.
Here we have used three different approaches to investigate TFBS functional constraints based on variation data. In the first one, using position-by-position comparisons, we have found that variability at TFBS positions generally correlates with information content, consistent with previous findings based on cross-species comparisons in Drosophila and human for other TFs [35, 36] and population studies in yeast . It should be noted that the majority of PWMs used in this study have been derived from comparing the sequences across all binding sites in one genome detected by genome-wide ChIP studies. Variation analyses look at sequence diversity in a different 'dimension': that is, across individuals at a particular point in the genome for each given binding site. That these two dimensions generally correlate with each other (and often also with in vitro biochemical data such as SELEX and protein binding microarrays [46, 47]) has been a reassuring confirmation of the general validity of PWM models to describe the sequence 'code' for the analyzed TFs. This, in turn, is an important prerequisite for using PWM scores to compute TFBS mutational load, a per-instance metric that combines the penetrance of a motif mutation with the loss of the PWM match it causes.
Cis-regulatory variation is accountable for serious deleterious effects, and yet it is common [14, 20]. Understanding TFBS functional constraints is therefore interesting for at least two reasons. First, it may shed light on the regulatory architecture of the genomes. For example, our finding that CTCF motifs tolerate the lowest load a short distance away from TSSs underlines the importance of chromatin architecture at the distal ends of promoter regions. In addition, TFBS constraints are indicators of how the system deals with noise in cis-regulatory networks, and the variation analyses presented here support such phenomena as homotypic redundancy . Interestingly, it was previously shown that homotypic clustering does not affect Drosophila TFBS turnover rate in the phylogenetic context , but the dynamics of selection inside a population need not correspond to that observed between species. For example, retaining multiple instances of neighboring homotypic sites in a given species may in itself bear the selective advantage to provide robust buffering to variation and other perturbations.
Genetic load, the concept that lies at the foundation of our constraint metric, was initially put forward by J Haldane  and HJ Muller , primarily in the context of the debate on hard versus soft selection. Here, however, we use this metric outside of such context and fully acknowledge that this is a crude, albeit computable parameter. We do not imply that a high TFBS load weakens the fitness of the individual bearing it, as would be the case in the 'classic' application of this concept. Rather, we take advantage of this concept to inquire why this probably does not occur - that is, why mutations at TFBSs are tolerated differently in different genomic contexts, likely without causing a significant reduction of an individual's fitness.
There is no doubt that mutational load is an imperfect metric. More sophisticated models linking fitness to the PWM score have been developed for cross-species phylogenetic analyses [49, 50] and their adaptation to population studies, although likely not straightforward, would be interesting to explore in the future. In addition, we know that the basic assumption of PWM models - that the frequency of nucleotide N at motif position K is proportionate to its positive impact on the binding affinity - does not always hold and even when it does, the amplitude of this effect may not be fully consistent across the TFs. Differences between motif sequences at different genomic locations may reflect TFBS optimization for a specific context rather than a lack of constraint. It was shown, for example, that differences at just two positions of the glucocorticoid receptor motif affect the choice of binding partners , while different k-mers of the apparently degenerate RACRYNNNNNACG motif in yeast are associated with the regulatory regions of genes with different functions . It is possible, therefore, that some mutations resulting in a loss of PWM match are, in fact, beneficial rather than deleterious and may be indicative of positive selection that was recently shown to occur at a fraction of Drosophila TFBSs by He et al. . However, in line with the assumption of He et al., we believe that the predominant direction of positive selection would be towards increasing PWM scores, and such mutations will have a zero load according to our definition.
These limitations, however, are universal for the problem of modeling functional constraints based on sequence alone. The predictive power of PWMs is probably comparable with our ability to predict the impact of mutations on RNA and protein structure. The rapidly increasing bulk of genotyping data will increase the statistical power of these analyses, but only experimental validation of the effects of TFBS mutations can give a definitive answer. This is why direct analyses of TF binding across individuals hold much promise. Using multi-individual CTCF binding maps [16, 44], it was reassuring to confirm that the loss of CTCF binding associated with a TFBS mutation is generally proportionate to its impact on motif PWM match. But perhaps more importantly, using these data has allowed us to observe that this relationship does not always hold, suggesting that variation at many sites, and in particular the most evolutionarily conserved ones, can be efficiently buffered at the binding level. We do not know the exact nature of these buffering mechanisms, and whether their prevalence at highly conserved sites is evolutionarily driven or is merely a side effect of the increasing complexity of regulatory networks [53, 54]. It can be expected that such buffering effects would be, at least in part, due to interactions with heterologous proteins. Given the multifaceted functions of CTCF, it is very likely that such interactions will involve different partners, depending on specific regulatory context. Studies of more 'specialized' TFs may therefore be more appropriate to address these questions. For example, analyses of individual variation at human NFκB  and yeast Ste12  pinpointed candidate interaction partners that affect the binding in the absence of mutations at the analyzed TF's own binding sites. We attempted to use the NFκB data to ask the reverse question, that is, look for factors that may help maintain the binding when mutations at conserved TFBSs are present; unfortunately, the number of such cases was extremely low, prohibiting this analysis. It is possible that mutations at conserved NFκB sites are poorly tolerated, implying that they are less efficiently 'buffered'. However, studies involving a larger number of individuals and/or using organisms with higher variation rates, such as Drosophila, will be required to adequately address this question.
Theoretically, TFBS mutations can be buffered at many different levels - starting from the motif itself that may 'absorb' a number of mutations due to a permissive consensus, to the level of CRMs (for example, homotypic motifs and protein interaction partners), cis-regulated genes (involving possible 'backup' by shadow enhancers ) as well as further along the regulatory network  - which may potentially explain the apparent redundancy that is often observed in the network architecture, both at the level of cooperative TF binding to enhancers and multiple 'cross-talking' pathways . Consistent with previous observations at individual CRMs , our observations suggest that much variation is buffered immediately in cis, via the redundancy of TFBS consensus sequences, neighboring homotypic motifs or other factors preserving regulator binding (or at least the overall CRM output). If true, this model may explain two of our preliminary observations that we initially found puzzling: that the levels of tolerated load did not significantly vary depending on the functional annotation of regulated genes (not shown) and that candidate Drosophila enhancers with seemingly very deleterious mutations at Bin, Tin and Twi binding sites were still able to drive reporter gene expression in vitro (Figure S7 in Additional file 1). It is clear, however, that this phenomenon requires further investigation, perhaps drawing more input from the biology of individual TFs. Finally, it is worth noting that a number of disease-causing mutations are located in regulatory regions, and presumably are either not buffered or inappropriately buffered. A well-studied example of this is the regulatory mutations in Pax6 regulatory regions associated with neurodevelopmental abnormalities . In addition, the majority of genome-wide association studies do not implicate a protein-coding variant . To fully understand these diseases we must gain a more complete knowledge of how variation impacts regulatory function.
Integrating genome-wide TF binding profiles with individual variation data in Drosophila and humans, we show that TFBSs are functionally constrained and yet mutations at them can be tolerated, providing evidence for possible 'buffering' effects. Beyond their direct biological implications, these results highlight the potential of integrating functional genomics and population genetics approaches for understanding cis-regulatory function.
Materials and methods
Data sources and basic analysis
Motif discovery data were from the modENCODE and ENCODE repositories [23, 24, 60, 61], with the exceptions of Bin, Tin and Twi that were from Zinzen et al. . Drosophila ChIP data were from Zinzen et al., modENCODE and other published sources [2, 24–30]; human ChIP data were from ENCODE  (see Tables S1 and S2 in Additional file 2 for details). CTCF multi-individual data were from [16, 44]. EPO alignments for 12 mammals were from Ensembl [62, 63]; phastcons scores  and multiz alignments for 12 Drosophila species were from Flybase [65, 66]. Drosophila variation data were from the DGRP , additionally filtered as described below. Human variation data were from the 1000 Genomes Pilot Project . Motif matches were detected using patser  (in case of overlapping matches, only the strongest-scoring motif was included) and overlaps with ChIP regions ('bound' motifs) were called using bedTools . Analysis was performed in R, Python and Perl with Ensembl API.
Filtering of DGRP data
DGRP SNPs were additionally filtered according to the following criteria: ε ≤ 0.02 (per SNP); p × ε ≤ 0.01 (per allele); coverage ≥ 3 (per allele); median coverage ≤20 (across strains); number of strains with detected homozygous alleles ≥100; number of strains with calls scored as 'heterozygous' ≤5%. The combination of these filters removed 31.3% low-confidence SNPs and increased the overlap with the SNPs detected by the Drosophila Population Genomics Project  based on a subset of the same Drosophila lines (not shown).
Motif selection for the analysis
For each modENCODE and ENCODE TF, a single combination of motif and cell type was chosen based on appreciable enrichments at TF-bound versus unbound regions, the total numbers of TF-bound motifs and a correlation between per-position evolutionary conservation and information content. Motif PWM score thresholds for human TFs were determined using TFM_PVALUE (P = 4e-8) , consistent with the thresholds used in ENCODE integrative analyses . For Drosophila TFs, thresholds were defined based on balancing the number of detected instances and motif enrichment at bound compared to unbound regions. Near-identical PWMs were removed based on Pearson correlation analyzed with STAMP [71, 72]. See Supplementary note on TF selection in Additional file 2 for more detail. The properties of selected motifs are listed in Tables S1 and S2 in Additional file 2. PWMs are listed in the data/motifs.txt files at  and , respectively. The positions, sequences, PWM scores and variation properties of all TFBSs included in this study are listed in Additional file 3 (Drosophila) and Additional file 4 (human).
Position-wise motif analysis
Reshuffled PWMs were generated by ten per-position permutations of the 'real' PWMs. Reshuffled motif matches were detected within the 200 bp proximity of real TF binding sites at the same PWM score thresholds as the real motifs. Position-wise variation data obtained for each permuted motif instance was then 'de-reshuffled' to match the positions of the real PWM to compute the total diversity per permuted motif position. For human motifs, the score thresholds used elsewhere in the study resulted in very low numbers of reshuffled motif instances detected near the corresponding TF binding sites. To overcome this, analyses in Figure 2 used slightly relaxed score thresholds for both real and reshuffled human motifs, adjusted such that the total number of motif instances detected with the 10 reshuffled PWMs was at least 1.5-times higher than the number of real instances for each TF.
Branch length score
BLS calculation was reimplemented in Perl for distributed computation on an LSF compute farm according to , allowing for a 50 bp motif movement either way along the alignment and a drop of motif score ≤1. Branch lengths are given relative to 12 eutherian mammals or 12 Drosophila species, respectively. Tree lengths were computed using Ensembl API.
TFBS mutational load
where w0 is the PWM score of the major allele, and wi and pi are the score and frequency of each allele, respectively. Classically, genetic load is expressed with respect to the maximum observed value (w0 = wmax). However, we have instead chosen to express it relative to the major allele (w0 = wmaj). The main reason for this is that, in the absence of ChIP data for each individual or isogenic line, TFBSs whose minor alleles have a higher PWM score than the major allele are subject to a significant ascertainment bias. Indeed, only TF-bound TFBS instances are included in the analysis, and we are much more likely to detect TFBSs as 'bound' when their weaker major alleles are also strong enough to ensure TF binding. Additionally, for reasons explained in the main text, we have postulated that TFBSs with stronger-scoring minor alleles have a zero load irrespective of frequency. Using the human data presented an additional challenge of interpreting heterozygous genotypes. Since the immediate phenotypic trait associated with TFBS's match to consensus (that is, TF binding) occurs in cis, we have taken the decision to consider each human allele separately. We did not focus exclusively on homozygous genotypes, as this approach would further reduce the statistical power of the analysis that was already limited by the low variation rates in the human genome.
Significance testing of TFBS load
Significance testing on TFBS load data was non-trivial, as their distributions are sparse (especially in the case of human data), with the majority of TFBSs having a load of zero. In statistical terms, these data present a case of zero-inflation, in which the observed zeros are a mixture of missing data (that is, mutations that are not observed due to a limited number of available genotypes) and 'real' zeroes (mutations that never occur because their deleterious effect is prohibitively strong). To overcome this problem, we have initially used generalized additive models (gam) based on zero-inflated distributions of the response variable (ZAGA for Drosophila and BEINF0 for human implemented in the R package gamlss ; not shown). However, gam P-values may be difficult to interpret, especially when the model includes random effects  (in our case, the TF identity). We have therefore eventually turned to permutation tests, permuting motif load values separately for each TF to avoid bias associated with specific properties of individual factors. To test the significance of trends, we used a permutation statistic based on : the dot product of the normalized data vector × and the index vector (1,..., N), where N is the length of X.
CTCF per-individual ChIP analysis
The analysis was based on lymphoblastoid lines, for which genotypes were available from the 1000 Genomes Pilot Project . We focused on the CTCF-binding data from McDaniell et al.  (Gm12892, Gm19239, Gm19238 and Gm19240) and confirmed the results using an independently generated dataset (Gm12872, Gm12873 and Gm12874)  processed through quantile normalization using the R/Bioconductor package preprocessCore. The remaining two datasets from  (Gm12878 and Gm12891) were excluded due to highly inconsistent overall binding score distributions. Global major allele data were from [75, 76]; assuming all reference alleles as major gave consistent results (not shown). Interaction models were plotted using the R package effects . The sequences, PWM scores and ChIP binding signals for all TFBSs included in these analyses are listed in Additional files 5 (individuals from ) and 6 (individuals from ).
branch length score
Drosophila Genetic Reference Panel
Encyclopedia of DNA Elements
position weight matrix
transcription factor binding site
transcription start site
The authors wish to thank David Garfield and Ian Dunham for valuable comments on the manuscript, and Eric Stone, John Marioni, David Thybert and all members of Birney and Furlong labs for helpful discussions. MS is supported by an EMBL interdisciplinary fellowship (EIPOD). This work was partially funded by a grant to EEF from the DFG (FU 750/1-1).
- Portales-Casamar E, Thongjuea S, Kwon AT, Arenillas D, Zhao X, Valen E, Yusuf D, Lenhard B, Wasserman WW, Sandelin A: JASPAR 2010: the greatly expanded open-access database of transcription factor binding profiles. Nucleic Acids Res. 2010, 38: D105-1010. 10.1093/nar/gkp950.PubMedPubMed CentralView ArticleGoogle Scholar
- Zinzen RP, Girardot C, Gagneur J, Braun M, Furlong EEM: Combinatorial binding predicts spatio-temporal cis-regulatory activity. Nature. 2009, 462: 65-70. 10.1038/nature08531.PubMedView ArticleGoogle Scholar
- Junion G, Spivakov M, Girardot C, Braun M, Gustafson E, Birney E, Furlong E: A transcription factor collective defines cardiac cell fate and reflects lineage history. Cell. 2012, 148: 473-486. 10.1016/j.cell.2012.01.030.PubMedView ArticleGoogle Scholar
- Visel A, Blow M, Li Z, Zhang T, Akiyama J, Holt A, Plajzer-Frick I, Shoukry M, Wright C, Chen F, Afzal V, Ren B, Rubin E, Pennacchio LA: ChIP-seq accurately predicts tissue-specific activity of enhancers. Nature. 2009, 457: 854-858. 10.1038/nature07730.PubMedPubMed CentralView ArticleGoogle Scholar
- Vavouri T, Elgar G: Prediction of cis-regulatory elements using binding site matrices-the successes, the failures and the reasons for both. Curr Opin Genet Dev. 2005, 15: 395-402. 10.1016/j.gde.2005.05.002.PubMedView ArticleGoogle Scholar
- Lusk RW, Eisen MB: Evolutionary mirages: selection on binding site composition creates the illusion of conserved grammars in Drosophila enhancers. PLoS Genet. 2010, 6: e1000829-10.1371/journal.pgen.1000829.PubMedPubMed CentralView ArticleGoogle Scholar
- Weirauch MT, Hughes TR: Conserved expression without conserved regulatory sequence: the more things change, the more they stay the same. Trends Genet. 2010, 26: 66-74. 10.1016/j.tig.2009.12.002.PubMedView ArticleGoogle Scholar
- Hare EE, Peterson BK, Iyer VN, Meier R, Eisen MB: Sepsid even-skipped enhancers are functionally conserved in Drosophila despite lack of sequence conservation. PLoS Genet. 2008, 4: e1000106-10.1371/journal.pgen.1000106.PubMedPubMed CentralView ArticleGoogle Scholar
- Schmidt D, Wilson MD, Ballester B, Schwalie PC, Brown GD, Marshall A, Kutter C, Watt S, Martinez-Jimenez CP, Mackay S, Talianidis I, Flicek P, Odom DT: Five-vertebrate ChIP-seq reveals the evolutionary dynamics of transcription factor binding. Science. 2010, 328: 1036-1040. 10.1126/science.1186176.PubMedPubMed CentralView ArticleGoogle Scholar
- Crocker J, Tamori Y, Erives A: Evolution acts on enhancer organization to fine-tune gradient threshold readouts. PLoS Biol. 2008, 6: e263-10.1371/journal.pbio.0060263.PubMedPubMed CentralView ArticleGoogle Scholar
- Crocker J, Potter N, Erives A: Dynamic evolution of precise regulatory encodings creates the clustered site signature of enhancers. Nat Commun. 2010, 1: 99-10.1038/ncomms1102.PubMedPubMed CentralView ArticleGoogle Scholar
- He BZ, Holloway AK, Maerkl SJ, Kreitman M: Does positive selection drive transcription factor binding site turnover? A test with Drosophila cis-regulatory modules. PLoS Genet. 2011, 7: e1002053-10.1371/journal.pgen.1002053.PubMedPubMed CentralView ArticleGoogle Scholar
- Garfield D, Haygood R, Nielsen W, Wray G: Population genetics of cis-regulatory sequences that operate during embryonic development in the sea urchin Strongylocentrotus purpuratus. Evol Dev. 2012, 14: 152-167. 10.1111/j.1525-142X.2012.00532.x.PubMedView ArticleGoogle Scholar
- Zheng W, Gianoulis TA, Karczewski KJ, Zhao H, Snyder M: Regulatory variation within and between species. Annu Rev Genomics Hum Genet. 2010, 12: 327-346.View ArticleGoogle Scholar
- Kasowski M, Grubert F, Heffelfinger C, Hariharan M, Asabere A, Waszak SM, Habegger L, Rozowsky J, Shi M, Urban AE, Hong M-Y, Karczewski KJ, Huber W, Weissman SM, Gerstein MB, Korbel JO, Snyder M: Variation in transcription factor binding among humans. Science. 2010, 328: 232-235. 10.1126/science.1183621.PubMedPubMed CentralView ArticleGoogle Scholar
- McDaniell R, Lee B-K, Song L, Liu Z, Boyle AP, Erdos MR, Scott LJ, Morken MA, Kucera KS, Battenhouse A, Keefe D, Collins FS, Willard HF, Lieb JD, Furey TS, Crawford GE, Iyer VR, Birney E: Heritable individual-specific and allele-specific chromatin signatures in humans. Science. 2010, 328: 235-239. 10.1126/science.1184655.PubMedPubMed CentralView ArticleGoogle Scholar
- Zheng W, Zhao H, Mancera E, Steinmetz LM, Snyder M: Genetic analysis of variation in transcription factor binding in yeast. Nature. 2010, 464: 1187-1189. 10.1038/nature08934.PubMedPubMed CentralView ArticleGoogle Scholar
- Chen K, van Nimwegen E, Rajewsky N, Siegal ML: Correlating gene expression variation with cis-regulatory polymorphism in Saccharomyces cerevisiae. Genome Biol Evol. 2010, 2: 697-707.PubMedPubMed CentralGoogle Scholar
- Majewski J, Pastinen T: The study of eQTL variations by RNA-seq: from SNPs to phenotypes. Trends Genet. 2011, 27: 72-79. 10.1016/j.tig.2010.10.006.PubMedView ArticleGoogle Scholar
- Manolio T: Genomewide association studies and assessment of the risk of disease. N Engl J Med. 2010, 363: 166-176. 10.1056/NEJMra0905980.PubMedView ArticleGoogle Scholar
- The 1000 genomes project consortium: A map of human genome variation from population-scale sequencing. Nature. 2010, 467: 1061-1067. 10.1038/nature09534.PubMed CentralView ArticleGoogle Scholar
- Mackay TFC, Richards S, Stone EA, Barbadilla A, Ayroles JF, Zhu D, Casillas S, Han Y, Magwire MM, Cridland JM, Richardson MF, Anholt RRH, Barran M, Bess C, Blankenburg KP, Carbone MA, Castellano D, Chaboub L, Duncan L, Harris Z, Javaid M, Jayaseelan JC, Jhangiani SN, Jordan KW, Lara F, Lawrence F, Lee SL, Librado P, Linheiro RS, Lyman RF, et al: The Drosophila melanogaster Genetic Reference Panel. Nature. 2012, 482: 173-178. 10.1038/nature10811.PubMedPubMed CentralView ArticleGoogle Scholar
- The ENCODE Consortium: An integrated Encyclopedia of DNA Elements in the human genome. Nature. 2012, 489: 57-74. 10.1038/nature11247.View ArticleGoogle Scholar
- Roy S, Ernst J, Kharchenko PV, Kheradpour P, Negre N, Eaton ML, Landolin JM, Bristow CA, Ma L, Lin MF, Washietl S, Arshinoff BI, Ay F, Meyer PE, Robine N, Washington NL, Di Stefano L, Berezikov E, Brown CD, Candeias R, Carlson JW, Carr A, Jungreis I, Marbach D, Sealfon R, Tolstorukov MY, Will S, Alekseyenko AA, Artieri C, Booth BW, et al: Identification of functional elements and regulatory circuits by Drosophila modENCODE. Science. 2010, 330: 1787-1797.PubMedPubMed CentralView ArticleGoogle Scholar
- Nègre N, Brown CD, Ma L, Bristow CA, Miller SW, Wagner U, Kheradpour P, Eaton ML, Loriaux P, Sealfon R, Li Z, Ishii H, Spokony RF, Chen J, Hwang L, Cheng C, Auburn RP, Davis MB, Domanus M, Shah PK, Morrison CA, Zieba J, Suchy S, Senderowicz L, Victorsen A, Bild NA, Grundstad AJ, Hanley D, MacAlpine DM, Mannervik M, et al: A cis-regulatory map of the Drosophila genome. Nature. 2011, 471: 527-531. 10.1038/nature09990.PubMedPubMed CentralView ArticleGoogle Scholar
- Li X-Y, MacArthur S, Bourgon R, Nix D, Pollard DA, Iyer VN, Hechmer A, Simirenko L, Stapleton M, Hendriks CLL, Chu HC, Ogawa N, Inwood W, Sementchenko V, Beaton A, Weiszmann R, Celniker SE, Knowles DW, Gingeras T, Speed TP, Eisen MB, Biggin MD: Transcription factors bind thousands of active and inactive regions in the Drosophila blastoderm. PLoS Biol. 2008, 6: 24-10.1371/journal.pbio.0060024.View ArticleGoogle Scholar
- MacArthur S, Li X-Y, Li J, Brown JB, Chu HC, Zeng L, Grondona BP, Hechmer A, Simirenko L, Kernen SVE, Knowles DW, Stapleton M, Bickel P, Biggin MD, Eisen MB: Developmental roles of 21 Drosophila transcription factors are determined by quantitative differences in binding to an overlapping set of thousands of genomic regions. Genome Biol. 2009, 10: R80-10.1186/gb-2009-10-7-r80.PubMedPubMed CentralView ArticleGoogle Scholar
- Bushey AM, Ramos E, Corces VG: Three subclasses of a Drosophila insulator show distinct and cell type-specific genomic distributions. Genes Dev. 2009, 23: 1338-5010. 10.1101/gad.1798209.PubMedPubMed CentralView ArticleGoogle Scholar
- Jakobsen JS, Braun M, Astorga J, Gustafson EH, Sandmann T, Karzynski M, Carlsson P, Furlong EEM: Temporal ChIP-on-chip reveals Biniou as a universal regulator of the visceral muscle transcriptional network. Genes Dev. 2007, 21: 2448-2460. 10.1101/gad.437607.PubMedPubMed CentralView ArticleGoogle Scholar
- Sandmann T, Girardot C, Brehme M, Tongprasit W, Stolc V, Furlong EEM: A core transcriptional network for early mesoderm development in Drosophila melanogaster. Genes Dev. 2007, 21: 436-449. 10.1101/gad.1509007.PubMedPubMed CentralView ArticleGoogle Scholar
- Haldane JBS: The cost of natural selection. J Genet. 1957, 55: 511-524. 10.1007/BF02984069.View ArticleGoogle Scholar
- Muller HJ: Our load of mutations. Am J Hum Genet. 1950, 2: 111-176.PubMedPubMed CentralGoogle Scholar
- Stormo GD, Zhao Y: Determining the specificity of protein-DNA interactions. Nat Rev Genet. 2010, 11: 751-760.PubMedGoogle Scholar
- Nuzhdin SV, Rychkova A, Hahn MW: The strength of transcription-factor binding modulates co-variation in transcriptional networks. Trends Genet. 2010, 26: 51-53. 10.1016/j.tig.2009.12.005.PubMedPubMed CentralView ArticleGoogle Scholar
- Moses AM, Chiang DY, Kellis M, Lander ES, Eisen MB: Position specific variation in the rate of evolution in transcription factor binding sites. BMC Evol Biol. 2003, 3: 19-10.1186/1471-2148-3-19.PubMedPubMed CentralView ArticleGoogle Scholar
- Kim J, He X, Sinha S: Evolution of regulatory sequences in 12 Drosophila species. PLoS Genet. 2009, 5: e1000330-10.1371/journal.pgen.1000330.PubMedPubMed CentralView ArticleGoogle Scholar
- Huston M: Biological Diversity: The Coexistence of Species on Changing Landscapes. 1994, Cambridge, UK: Cambridge University PressGoogle Scholar
- Kophengnavong T, Michnowicz JE, Blackwell TK: Establishment of distinct MyoD, E2A, and twist DNA binding specificities by different basic region-DNA conformations. Mol Cell Biol. 2000, 20: 261-272. 10.1128/MCB.20.1.261-272.2000.PubMedPubMed CentralView ArticleGoogle Scholar
- Charlesworth B: Fundamental concepts in genetics: effective population size and patterns of molecular evolution and variation. Nat Rev Genet. 2009, 10: 195-205.PubMedView ArticleGoogle Scholar
- Kheradpour P, Stark A, Roy S, Kellis M: Reliable prediction of regulator targets using 12 Drosophila genomes. Genome Res. 2007, 17: 1919-1931. 10.1101/gr.7090407.PubMedPubMed CentralView ArticleGoogle Scholar
- Godt D, Couderc JL, Cramton SE, Laski FA: Pattern formation in the limbs of Drosophila: bric à brac is expressed in both a gradient and a wave-like pattern and is required for specification and proper segmentation of the tarsus. Development. 1993, 119: 799-812.PubMedGoogle Scholar
- Ohlsson R, Bartkuhn M, Renkawitz R: CTCF shapes chromatin by multiple mechanisms: the impact of 20 years of CTCF research on understanding the workings of chromatin. Chromosoma. 2010, 119: 351-360. 10.1007/s00412-010-0262-0.PubMedPubMed CentralView ArticleGoogle Scholar
- Fiston-Lavier A-S, Singh ND, Lipatov M, Petrov DA: Drosophila melanogaster recombination rate calculator. Gene. 2010, 463: 18-20. 10.1016/j.gene.2010.04.015.PubMedView ArticleGoogle Scholar
- Maurano M, Wang H, Kutyavin T, Stamatoyannopoulos J: Widespread site-dependent buffering of human regulatory polymorphism. PLoS Genet. 2012, 8: e1002599-10.1371/journal.pgen.1002599.PubMedPubMed CentralView ArticleGoogle Scholar
- Bradley RK, Li X-Y, Trapnell C, Davidson S, Pachter L, Chu HC, Tonkin LA, Biggin MD, Eisen MB: Binding site turnover produces pervasive quantitative changes in transcription factor binding between closely related Drosophila species. PLoS Biol. 2010, 8: e1000343-10.1371/journal.pbio.1000343.PubMedPubMed CentralView ArticleGoogle Scholar
- Bulyk ML: Protein binding microarrays for the characterization of DNA-protein interactions. Adv Biochem Eng Biotechnol. 2007, 104: 65-85.PubMedPubMed CentralGoogle Scholar
- Hallikas O, Palin K, Sinjushina N, Rautiainen R, Partanen J, Ukkonen E, Taipale J: Genome-wide prediction of mammalian enhancers based on analysis of transcription-factor binding affinity. Cell. 2006, 124: 47-59. 10.1016/j.cell.2005.10.042.PubMedView ArticleGoogle Scholar
- Gotea V, Visel A, Westlund JM, Nobrega MA, Pennacchio LA, Ovcharenko I: Homotypic clusters of transcription factor binding sites are a key component of human promoters and enhancers. Genome Res. 2010, 20: 565-577. 10.1101/gr.104471.109.PubMedPubMed CentralView ArticleGoogle Scholar
- Halpern AL, Bruno WJ: Evolutionary distances for protein-coding sequences: modeling site-specific residue frequencies. Mol Biol Evol. 1998, 15: 910-917. 10.1093/oxfordjournals.molbev.a025995.PubMedView ArticleGoogle Scholar
- Moses AM: Statistical tests for natural selection on regulatory regions based on the strength of transcription factor binding sites. BMC Evol Biol. 2009, 9: 286-10.1186/1471-2148-9-286.PubMedPubMed CentralView ArticleGoogle Scholar
- Meijsing SH, Pufall MA, So AY, Bates DL, Chen L, Yamamoto KR: DNA binding site sequence directs glucocorticoid receptor structure and activity. Science. 2009, 324: 407-410. 10.1126/science.1164265.PubMedPubMed CentralView ArticleGoogle Scholar
- Swamy KBS, Cho C-Y, Chiang S, Tsai ZT-Y, Tsai H-K: Impact of DNA-binding position variants on yeast gene expression. Nucleic Acids Res. 2009, 37: 6991-7001. 10.1093/nar/gkp743.PubMedPubMed CentralView ArticleGoogle Scholar
- Gibson G: Epistasis and pleiotropy as natural properties of transcriptional regulation. Theor Popul Biol. 1996, 49: 58-89. 10.1006/tpbi.1996.0003.PubMedView ArticleGoogle Scholar
- Bolouri H, Davidson EH: Transcriptional regulatory cascades in development: initial rates, not steady state, determine network kinetics. Proc Natl Acad Sci USA. 2003, 100: 9371-9376. 10.1073/pnas.1533293100.PubMedPubMed CentralView ArticleGoogle Scholar
- Barolo S: Shadow enhancers: Frequently asked questions about distributed cis-regulatory information and enhancer redundancy. BioEssays. 2012, 34: 135-141. 10.1002/bies.201100121.PubMedPubMed CentralView ArticleGoogle Scholar
- Hartman JL, Iv JLH, Hartwell L: Principles for the buffering of genetic variation. Science. 2001, 291: 1001-1004. 10.1126/science.291.5506.1001.PubMedView ArticleGoogle Scholar
- Costanzo M, Baryshnikova A, Myers CL, Andrews B, Boone C: Charting the genetic interaction map of a cell. Curr Opin Biotechnol. 2011, 22: 66-74. 10.1016/j.copbio.2010.11.001.PubMedView ArticleGoogle Scholar
- Ludwig MZ, Bergman C, Patel NH, Kreitman M: Evidence for stabilizing selection in a eukaryotic enhancer element. Nature. 2000, 403: 564-547. 10.1038/35000615.PubMedView ArticleGoogle Scholar
- Sisodiya SM, Free SL, Williamson KA, Mitchell TN, Willis C, Stevens JM, Kendall BE, Shorvon SD, Hanson IM, Moore AT, Van Heyningen V: PAX6 haploinsufficiency causes cerebral malformation and olfactory dysfunction in humans. Nat Genet. 2001, 28: 214-216. 10.1038/90042.PubMedView ArticleGoogle Scholar
- ENCODE Motif Browser. [http://www.broadinstitute.org/~pouyak/motif-disc/human]
- modENCODE Motif Browser. [http://www.broadinstitute.org/~pouyak/motif-disc/fly]
- Paten B, Herrero J, Fitzgerald S, Beal K, Flicek P, Holmes I, Birney E: Genome-wide nucleotide-level mammalian ancestor reconstruction. Genome Res. 2008, 18: 1829-1843. 10.1101/gr.076521.108.PubMedPubMed CentralView ArticleGoogle Scholar
- Ensembl Genome Browser. [http://www.ensembl.org/index.html]
- Siepel A, Bejerano G, Pedersen JS, Hinrichs AS, Hou M, Rosenbloom K, Clawson H, Spieth J, Hillier LW, Richards S, Weinstock GM, Wilson RK, Gibbs RA, Kent WJ, Miller W, Haussler D: Evolutionarily conserved elements in vertebrate, insect, worm, and yeast genomes. Genome Res. 2005, 15: 1034-1050. 10.1101/gr.3715005.PubMedPubMed CentralView ArticleGoogle Scholar
- Tweedie S, Ashburner M, Falls K, Leyland P, McQuilton P, Marygold S, Millburn G, Osumi-Sutherland D, Schroeder A, Seal R, Zhang H: FlyBase: enhancing Drosophila Gene Ontology annotations. Nucleic Acids Res. 2009, 37: D555-559. 10.1093/nar/gkn788.PubMedPubMed CentralView ArticleGoogle Scholar
- Flybase. [http://www.flybase.org]
- Hertz GZ, Stormo GD: Identifying DNA and protein patterns with statistically significant alignments of multiple sequences. Bioinformatics. 1999, 15: 563-577. 10.1093/bioinformatics/15.7.563.PubMedView ArticleGoogle Scholar
- Quinlan AR, Hall IM: BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010, 26: 841-842. 10.1093/bioinformatics/btq033.PubMedPubMed CentralView ArticleGoogle Scholar
- Drosophila Population Genomics Project. [http://www.dpgp.org]
- Touzet H, Varré J-S: Efficient and accurate P-value computation for position weight matrices. Algorithms Mol Biol. 2007, 11 (2): 15-View ArticleGoogle Scholar
- Mahony S, Benos PV: STAMP: a web tool for exploring DNA-binding motif similarities. Nucleic Acids Res. 2007, 35: W253-W258. 10.1093/nar/gkm272.PubMedPubMed CentralView ArticleGoogle Scholar
- STAMP: A Tool-kit for DNA Motif Comparison. [http://www.benoslab.pitt.edu/stamp]
- Stasinopoulos DM, Rigby RA: Generalised additive models for Location Scale and Shape (GAMLSS) in R. J Stat Software. 2007, 23: 1-46.View ArticleGoogle Scholar
- Cox D, Hinkley D: Theoretical Statistics. 1974, London: Chapman & Hall, 188-View ArticleGoogle Scholar
- Dewey FE, Chen R, Cordero SP, Ormond KE, Caleshu C, Karczewski KJ, Whirl-Carrillo M, Wheeler MT, Dudley JT, Byrnes JK, Cornejo OE, Knowles JW, Woon M, Sangkuhl K, Gong L, Thorn CF, Hebert JM, Capriotti E, David SP, Pavlovic A, West A, Thakuria JV, Ball MP, Zaranek AW, Rehm HL, Church GM, West JS, Bustamante CD, Snyder M, Altman RB, et al: Phased whole-genome genetic risk in a family quartet using a major allele reference sequence. PLoS Genet. 2011, 7 (9): e100228010-View ArticleGoogle Scholar
- Human Synthetic Major Allele Data from Dewey. [http://datadryad.org/handle/10255/dryad.34659]
- Fox J: Effect displays in R for generalised linear models. J Stat Software. 2003, 8: 1-27.View ArticleGoogle Scholar
- Nei M: Analysis of Gene Diversity in Subdivided Populations. Proc Natl Acad Sci USA. 1973, 70: 3321-3323. 10.1073/pnas.70.12.3321.PubMedPubMed CentralView ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.