Naturally Occurring Deletions of Hunchback Binding Sites in the Even-Skipped Stripe 3+7 Enhancer

Changes in regulatory DNA contribute to phenotypic differences within and between taxa. Comparative studies show that many transcription factor binding sites (TFBS) are conserved between species whereas functional studies reveal that some mutations segregating within species alter TFBS function. Consistently, in this analysis of 13 regulatory elements in Drosophila melanogaster populations, single base and insertion/deletion polymorphism are rare in characterized regulatory elements. Experimentally defined TFBS are nearly devoid of segregating mutations and, as has been shown before, are quite conserved. For instance 8 of 11 Hunchback binding sites in the stripe 3+7 enhancer of even-skipped are conserved between D. melanogaster and Drosophila virilis. Oddly, we found a 72 bp deletion that removes one of these binding sites (Hb8), segregating within D. melanogaster. Furthermore, a 45 bp deletion polymorphism in the spacer between the stripe 3+7 and stripe 2 enhancers, removes another predicted Hunchback site. These two deletions are separated by ∼250 bp, sit on distinct haplotypes, and segregate at appreciable frequency. The Hb8Δ is at 5 to 35% frequency in the new world, but also shows cosmopolitan distribution. There is depletion of sequence variation on the Hb8Δ-carrying haplotype. Quantitative genetic tests indicate that Hb8Δ affects developmental time, but not viability of offspring. The Eve expression pattern differs between inbred lines, but the stripe 3 and 7 boundaries seem unaffected by Hb8Δ. The data reveal segregating variation in regulatory elements, which may reflect evolutionary turnover of characterized TFBS due to drift or co-evolution.


Evolution of Transcriptional Regulatory Sequences
The molecular basis for phenotypic divergence and standing variation is often attributed to differences in the regulation of transcription [1][2][3]. The mechanistic principles of regulatory DNA and factor structure and function such as; multiple transcription factor binding sites (TFBS), TFBS motif degeneracy, cooperativity and number of trans factors [3,4] and interactions between transcription factors (TFs), enhancers and promoters [5,6] impose unique rules on their evolution. Regulatory DNA has no single ''active-site'', since most regions consist of multiple transcription factor binding sites. Evolutionary analyses of experimentally verified TFBS demonstrate examples of conservation, but also reveal evolutionary turnover of TFBS, were some sites are lost and others gained [7][8][9].
It has been postulated that selection mainly acts on the transcriptional output of a gene (timing, location and amount) and does not preserve individual TFBS [10,11]. That is, changes in TFBS and even losses are permitted, if the transcriptional output is preserved. Such models of stabilizing selection acting on transcriptional output can account for both loss of functional binding sites and evolutionary fine-tuning of regulatory elements [12]. They also suggest that positive selection may sometimes play a role, acting on compensatory mutations in cis or trans. Several studies [13][14][15][16] have investigated the evolutionary origin of TFBS, including co-evolution within regulatory sequences. From first principles one would predict both co-evolution in cis (promoters, regulatory modules, more distantly located signals like insulators) and co-evolution of sequence elements with the trans environment (abundance of transcription factor, mediator or holenzyme components). The model of trans co-evolution is corroborated by studies of between-species hybrids [17], which e.g. reveal misexpression of genes in hybrids of D. melanogaster and D. simulans, two closely related species, most likely due to species-specific cis-trans compensatory evolution. Also, genome-wide changes in cis elements of co-expressed genes in two distantly related yeast species document the co-evolution of the TF repertoire of an organism and the regulatory elements of coordinately expressed genes [18,19]. Numerical models show how mutation and drift can generate binding sites, and predictably that selection can speed up fixation of new TFBS [20]. Crucially, functional polymorphism (both single nucleotide polymorphism: SNPs or insertion/deletion polymorphism: indels) in human enhancers, are shaped by positive selection [21].

Insertion and Deletion Polymorphism in Regulatory DNA
Population genetics studies have largely neglected indels, perhaps because they represent a minority of segregating variation in most genomes [22]. Deletion polymorphism in the intergenic region of Adh in Drosophila pseudoobscura does not conform to neutral evolution, but exhibits signatures of purifying selection, i.e. deletion (but not insertion) polymorphism was removed from introns over time [23]. On a larger scale, Comeron and Kreitman [24] revealed a bias in the insertion and deletion frequency distribution in D. melanogaster populations. While deletion events were more common and on average longer, insertions were at significantly higher frequency. This may reflect both mutational bias (because the mechanisms causing deletions are different from those causing insertions) and a difference in selection pressures, with purifying selection keeping a large fraction of deletions at low frequency in the population [24,25]. Ometto et al. [25], on the other hand, also concluded that weak positive selection might increase the population frequency of some insertions, which is supported by a genome-wide study in D. melanogaster [26]. Population genetic analyses of Bicoid response genes in D. melanogaster revealed single nucleotide polymorphism (SNPs) in 13 of 85 predicted Bicoid binding sites [27]. Most notable was the high frequency of SNP and indel polymorphism in the Orthodenticle (Otd) early head enhancer. These polymorphisms clustered on two haplotypes, both at intermediate frequency. Transgene tests showed that the Otd haplotypes differ in transcriptional output [27]. Similarly, studies of the Endo16 promoter and other sea-urchin enhancers [28][29][30] show that many TFBS are affected by segregating indel variation. In particular, in Endo16 two rare insertions affect the same part of the promoter. One of these generated a functional repressor module [29].

Enhancers of Eve as a Model of Regulatory Evolution
Early embryonic development in D. melanogaster is regulated by numerous genes through a complex network of activation and repression, resulting in segmental boundaries along the embryo length [31][32][33][34][35]. The accurate temporal and spatial expression of these genes is mainly achieved by integration of multiple TFs and their binding to regulatory sequences. Some regulatory functions (required for a given expression pattern) are aggregated in distinct modules like the eve stripe 2 enhancer (s2e) and the stripe 3 and 7 enhancer (s3+7e). These experimentally verified ''minimal'' enhancers [36,37] suffice to generate 4-7 cell-wide Eve stripes in early development. Not all regulatory sequences contain modular enhancers, and often spacer sequences (separating regulatory modules) have function, meaning that the length of these sequences matters for proper function of flanking cis-modules [10,38].
The cumulative effect of nucleotide changes in s2e between species is a turnover of functional motifs within enhancers [8,11,39]. Notably, the s2e from D. erecta is less effective than s2e from the more distantly related D. pseudoobscura at complementing a deletion of the s2e in D. melanogaster [11]. Quantitative analysis of the amount of Eve in stripe 2 illustrated the functional deficiency of the D. erecta s2e in the D. melanogaster genetic background. This means that, for a given enhancer, the spatial and temporal features of the expression pattern are highly conserved, but the quantity of gene product probably less so. The expression level of developmentally-specific gene products may exhibit changes over evolutionary time, possibly reflecting ''developmental system drift'' [40,41].
The aim of the current study was to gauge the level of polymorphism in the well-characterized regulatory regions in D. melanogaster, with particular focus on insertion and deletion polymorphism. Consistent with other studies and evolutionary theory, SNP and indel-polymorphism are rare in TFBS. However we find two peculiarly large and common deletions in and close to the eve stripe 3 and 7 enhancer. Both deletions remove binding sites for Hunchback, prompting analysis of the genetics and phylogeography of one of those polymorphisms and its potential phenotypic effects. The data provide insights into the nature of variation segregating in cis-regulatory elements.

Flies and Populations
Several populations of flies where studied. The population genetic surveys were done on collections of inbred lines derived from North Carolina, collected in 2000 and 2005 [27,42], and a Costa Rican sample from Peter Andolfatto, made isogenic for the second chromosome by three generations of crosses. Walter Eanes provided DNA from thirteen US East coast populations [43]; a total of 380 individuals used to test for clinal variation in the eve region. Jean-Claude Walser provided a sample of 46 cosmopolitan populations [44], in which DNA from 100 lines in each population was pooled.

PCR, DNA Sequencing and Genotyping
Primers were designed with primer 3 version 0.3 (frodo.wi.mit.edu [45]) for 13 well-characterized early developmental enhancers or promoters and several other non-coding regions (see Table  S1). The regions studied were several parts of the eve locus (the late element, s2e, s3+7e, and the promoter, along with two spacer sequences), Kruppel promoter and CD1, salm wing blade enhancer, ems abdominal enhancer, en regulatory region and promoters ofAntp, Ubx-bxd, tll, Act57B, RpL29/CG30390 and RpL30. The sequence variation in those regions was assessed by PCR followed directly by DNA sequencing. PCR was done as before [46] with Takara Taq and MJ Tetrad machines on 96 well plates. Products where purified by Qiagen purification columns or Exo-sap. DNA sequencing was done on purified PCR products, with the forward and reverse primer using Applied biosystems reagents. The ethanol purified reaction products where run in the University of Chicago sequencing facility or the ABI sequencing machine at the Institute of Biology, University of Iceland.
The deletion of the Hb8 site in s3+7e (see below) and the wild type allele were genotyped with PCR using allele-specific primers (Table S1). We ran separate reactions for both alleles on individuals from the East coast sample and on bulk DNA samples from the cosmopolitan sample. This was used to infer geographic distribution of specific variants, but does of course not yield information about frequency. All sequences were submitted as Popset data to NCBI (accession numbers: KJ465109-KJ465866), except two alignments that were shorter than 200 bp (provided in fasta format as Supporting information S1 and S2).

Population Genetic Analysis
Metrics of population genetics (S, p, h, Haplotype number) were calculated for SNPs and indels with Tassel vs. 2.1 (www. maizegenetics.net [47]), either for individual regions or as a sliding window for the haplotype analysis. Tassel was also used to calculate LD, and R (www.r-project.org version 12.3 [48]) for testing of contingency tables. DNAsp vs. 4.1 (www.ub.edu/dnasp [49]) was also used to test for deviations of Tajima's D and Fu and Li's estimators. Furthermore Hudson's haplotype test [50] (utilizing the ms program and the psub option) was used to test for positive selection in four eve regions.

Phylogenetic Shadowing
A 2 kb region surrounding the stripe 3+7 minimum enhancer was blasted against the 12 finished genomes (insects.eugenes.org/ species), and the orthologous regions extracted (except D. willistoni which did not return a significant blast hit). The Drosophila species (abbreviated) and  [51]). Divergence in these sequences is considerable, requiring manual curating in Genedoc (www.psc.edu/biomed/genedoc [52]), with special devotion to characterized TFBS from redfly.ccr.buffalo.edu [53] and OR-egAnno [54]. In addition two additional Hb sites (Hb15 and Hb16) found by Stanojovic et al. [55] and two Stat binding sites discovered by Yan et al. [32] were included. We found that the D. melanogaster Stat binding sites differ from the genomic sequence, probably due to sequencing error (Stat-1 was reported to start with an A and stat-2 was reported as GTTCCCCGAAA, highlighted bases differ).
We also used (jaspar.genereg.net [56]) to predict Hb binding sites (score above 6) in the ,8000 bp upstream of eve, in D. melanogaster, D. sechellia, D. yakuba and D. pseudoobscura. Based on multiple alignments from Mavid, and Multiz alignments from the Santa Cruz genome browser (downloaded in December 2013), we mapped predicted Hb binding sites in orthologous and more rapidly evolving regions.

Testing the Effects of a Segregating Deletion on Adult Phenotypes
A set of 20-60 healthy inbred lines from NC [46] were used for the two experiments conducted to test the effects of a 72 bp deletion within s3+7e (called Hb8D, see below) on viability and developmental time. The first was a set of controlled crosses to lines deficient for eve, and the second was phenotyping of 60 genotyped inbred lines. All fly-rearing took place on cornmeal food at constant temperature, 25uC.
We first crossed the inbred lines to four stocks with characterized eve mutations. Ten inbred lines, homozygous for each allele (Hb8D or wt) were crossed to each eve mutant. . Three virgins of a mutant stock were crossed with 3 males from each of the 20 inbred lines, and allowed to lay eggs for 2-3 days. The offspring were counted and sexed, between 10 and 11 am, from day 10 to 18. The experiment was fully balanced and repeated three times, several weeks apart. The parents of all lines used in the crosses had been grown for 2 generations under controlled density (parents discarded between days 2 and 5 depending on visual assessment of egg number). We recorded both the total number of offspring (viability), and developmental time, summarized as the average time to eclosion for a given combination of, mutation, cross, genotype, sex and replicate.
For the association tests, 60 inbred lines where studied. The Hb8D/wt polymorphism was genotyped in three individuals of each line in the generation that was phenotyped. The rearing and measuring procedure was identical to the first experiment, except no crosses were required and only replicates were measured (two weeks apart).

Embryo Collections, Fixing and Staining
The embryos were collected, fixed and stained with standard protocols, as we have done before [8]. Four inbred lines with (NC25 and NC128) and without (NC006, NC017) the Hb8D laid eggs for 4-5 hrs at 22uC. Briefly, we collected embryos from each of the four lines, and they were fixed. Multiple embryo collections were pooled before staining with Eve primary antibody and a secondary antibody. The histochemical LacZ staining reaction was run for 12 minutes. The stained embryos were stored in 70% glycerol at 4uC, and photographed within a week.

Photography and Measurements
Each embryo in the appropriate developmental stage range was photographed three times at 20X magnification with water immersion on a Zeiss microscope. First a DIC sagittal section yielding maximum length of embryo and then two sections (DIC sagittal and bright-field) captured the stripes. Tiff photographs were saved and the X and Y coordinates of stripe boundaries assessed in ImageJ (rsb.info.nih.gov/ij/ [57]). First, a straight line was superimposed on the sagittal image, and the X-Y coordinates of anterior and posterior of the embryo recorded. Second, the same guideline was superimposed on the other two images and X-Y landmarks of the anterior and posterior boundary of each stripe were visually assessed and recorded. Third, the rotation of the embryo along the Dorsal/Ventral axis was scored. Finally, the stage of development was also visually assessed from eve pair-rule expression, in increments of 0.5 on the scale from 1 to 5, around cellularization [11]. The same investigator (AP) did all measurements.

Summarizing the Expression Pattern
The raw landmark data indicating the length of the embryo and placement of stripes were processed in two ways. The relative positions of stripe boundaries were estimated by calculating distance of landmarks from the anterior and posterior end using standard geometric formula. First, the length of the embryos was estimated. Second, the relative distance from one embryo tip to the anterior and posterior boundary of each stripe was calculated.

Statistical Analysis of Adult and Embryonic Phenotypes
SAS version 8.2 [58] was used for analyses of phenotypes. The viability and developmental time analyses were conducted with mixed model ANOVA (proc MIXED). The model for the testcross was:

zL(CXG)zerror
Denoting the fixed effects of the mutation (M), that is the 4 different eve deficiencies or point mutations, the cross (C) designating the balancer (CyO) or the ''loss of function'' (LoF) eve mutation, the genotype (G) term which evaluates the effects of Hb8D, sex (S) and appropriate interaction terms. The effects of Line (L) and replicate vials (R) are considered random factors. Furthermore, the total number of offspring (O) was included as a covariate. As a large factorial model with 4 fixed terms runs the risk of being overly parametrized, higher order terms were evaluated and dropped if they were not significant at the 0.05 level. The association tests of the inbred lines data were simpler, with only terms denoting genotype, sex and total number of offspring, and not described here.
The relative location of histochemically detected Eve stripes was studied similarly. In order to remove the effects of orientation, a reduced model was fit, and the residuals were used in the subsequent analysis. The positioning of stripes was analyzed with a mixed model ANOVA. The dependent variables of interest are the relative positioning of stripe boundaries, with the anterior boundary of stripe 3 (S3A) and the posterior position of stripe 7 (S7P) being particular candidates given prior evidence on Hunchback distribution in the embryo [59]. The ANOVA model had the general form:

Y~GzTzGXTzL G ð Þzerror
Where G, indicating genotype (the presence or absence of Hb8), is a fixed main effect. The covariate T (for developmental time) captures the developmental progression and L is a random term for different inbred lines. The relative stripe position matrix (anterior/posterior boundary of all 7 stripes) was also summarized with Principle component (Proc PRINCOMP) on the correlation matrix. Only the first component, with eigenvalue 7.42, was analyzed for dependence on Hb8 genotype.

Polymorphism in Regulatory DNA Includes Large Deletions of TFBS
First we surveyed the molecular variation, i.e. nature, frequency and distribution of polymorphisms, in 13 well studied Drosophila regulatory elements and several less well defined elements and spacer sequences. Few indel polymorphisms are found in the regulatory regions, 8 of the regions have no indels (Table 1). Purifying selection seems to affect both SNP and indel polymorphism, as there is a significant correlation between h for SNPs and indels (r = 0.48, p = 0.03, Figure S1A). The size and frequency of indels in characterized cis-elements was contrasted to those in noncoding regions surrounding two developmental genes, hairy and EGFR [46,60]. As was previously observed [46] most indels are short, and rarely do large indels (more than 10 bp) reach appreciable frequency ( Figure S1B). The notable exception is a 72 bp deletion in the stripe 3 and 7 enhancer (s3+7e) of eve ( Figure 1A and B). Interestingly this deletion removes a DNase I characterized Hunchback (Hb) binding site [55], and is henceforth called Hb8D. Bioinformatic analyses in Jaspar show that this site has a PWM score of 8.5, suggesting the notion that this a   transcription factor binding site presence/absence polymorphism. Oddly enough, less than 250 bp away (in the spacer separating s3+7e and s2e), another segregating large deletion also removes a putative Hunchback binding site ( Figure 1C). This site (here called Hbs1) is predicted with high PWM score, 11.2. That is the fourth highest score of 60 predicted Hb sites in the 8 kb region upstream of eve in D. melanogaster ( Figure S2A and Table S2). Most of the 21 DNaseI characterized Hb sites in s3+7e and s2e have lower scores than Hbs1. This 45 bp deletion in the spacer is referred to as Hbs1D. This putative Hb binding site has probably been unnoticed for two reasons. It sits outside the fragments tested for enhancer function, presumably because of restriction site locations [10,37]. Also, the D. melanogaster reference genome sequence contains the deletion. To iterate, the 45 deleted bases do not appear in the standard versions of the D. melanogaster genome and are only visible in genomic alignments with close Drosophila relatives or population genetic sequence data. The two deletions sit on distinct haplotypes, and are never found in the same inbred lines. They are both at appreciable frequency, in a sample of 55 Costa Rican chromosomes the Hb8D and Hbs1D are at 9% and 17% frequency respectively ( Figure 1B and C). This leads to the question, are these deletions harmful, neutral or beneficial?
Phylogenetic Footprinting of s3+7e shows the Hb Sites are Conserved Comparative genomic alignments of the s3+7e and the adjacent regions with 12 publicly-available Drosophila genomes [61] were used to assess the functional importance of these two predicted Hb binding sites, and other characterized Hb, Kni and Stat sites [32,55,59]. Similarly to the eve s2e, TFBS in s3+7e are highly conserved (Table S3); 3 of 13 Hb sites are identical from D. melanogaster to D. mojavensis and 9 have none or only one mutation between D. melanogaster and D. persimilis. The Hb8 site is found in all of the 12 species, except D. ananassae (most probably due to a gap in the genomic sequence), but has experienced several substitutions ( Figure 1D). The PWM score for Hb8 is 8.2 in D. melanogaster and D. simulans, but 9.9 in D. yakuba and D. pseudoobscura ( Figure S2 and table S2). On the other hand, the predicted Hbs1 site (with a PWM score of 11.2) is completely conserved between D. melanogaster and D. yakuba, but was not found in distantly related species. Those data suggest considerable evolutionary constraints on those sequences, arguing that they could indeed be functional Hb binding sites. But in the absence of functional tests they must regarded as putative Hb binding sites.
Additionally, the Hb8D also removes half of a putative Sloppy Paired 1 (Slp1) binding site. The putative Slp1 site is less conserved then the characterized Slp1 site in s2e [62] (Table S4), but no SNPs within either of these two (characterized and putative) Slp1 binding sites in eve, in 104 sequenced alleles, suggests selective constraint within D. melanogaster at least. The genome comparisons confirm that both Hb binding sites in eve affected by these two deletions have been protected by purifying selection. This prompts the question, why do these deletions of conserved TFBS occur at such high frequencies in populations? Here we focus mainly on studying the population genetics of Hb8D and assess its potential impact on development and fitness.

Polymorphism on the Hb8D and wt Haplotypes
How can a deletion removing a conserved binding site be at such high frequency in the population? One possibility is that the deletion of Hb8 is buffered by compensatory mutations (sitting on the same haplotype). To assess this, and to evaluate the polymorphism in the region, two strategies were deployed. One was deeper sequencing of four eve regions (the promoter, s2e, s37e and the late element) in inbred strains from North Carolina, and the other, a contrast of sequence diversity in alleles with or without the Hb8D in ,8 kb around s3+7e.
The Hb8D is at 32% frequency in the NC population (N = 63), and there is less variation on the Hb8D haplotypes compared to the wt haplotypes (Table 2). For instance p (which captures the number of substitutions and their frequency) is 25% to 100% lower on the Hb8D haplotypes. This is most extreme in the s3+7e, and notably weaker in flanking regions. This tendency was captured by other population genetic summary statistics (S, Haplotypes, haplotype diversity and Dxy -a measure of differences in nucleotide substitution rate between samples).  We next compared more extensively the sequence variation on the Hb8D and wt chromosomes and screened for variants that might possibly compensate for the loss of this Hb binding site. We estimated the polymorphism on two distinct haplotypes carrying either the wt or deletion polymorphism, by sequencing 16 (Hb8D) and 18 (wt) chromosomes of each type. The 8200 bp region we selected spans the eve neighborhood, from the 39UTR of CG12134 to the end of the transcript. There is reduced polymorphism (p and h) on the Hb8D haplotypes compared to wt haplotypes (Figure 2A and B), which is consistent with selection for the Hb8D bearing haplotype. Another indicator of long haplotypes is high LD between Hb8D and polymorphic sites in the region ( Figure 2C). Several sites more than 3 kb away from Hb8 are in high LD (r 2 . 0.7) with the deletion. Additionally, most polymorphism in the region shows perfect coupling or repulsion LD to Hb8D (data not shown). (The Hbs1D was only found in 3 (wt) lines. Omission of those 3 lines did not affect the outcome of the polymorphism analyses -data not shown). Furthermore, no variants are unique to the Hb8D haplotype. Finally, no potential compensatory mutations that strengthened or generated other Hb sites were observed. The data do reveal less diversity on the Hb8D haplotype, compared to the wt haplotype. Note however, standard tests of natural selection can not be deployed on these data because the sampling was not random from a population; lines were picked for sequencing to get similar representation of wt and Hb8D chromosomes.

Geographic Distribution of the Hb8D
What is the geographic distribution of Hb8D and does it correlate with geographic attributes? To study the geographic distribution, bulk DNA samples from 51 cosmopolitan samples, from Europe, Africa, Asia and South America [44] were genotyped with allele specific primers. There was evidence of Hb8D in 43 of the 51 populations (Table S5), consistent with an evolutionarily old and broadly distributed polymorphism. The cosmopolitan distribution of the Hb8D is unlikely if it was strongly deleterious.
Does this binding site deletion show any relationship with geographic attributes? To assess this we genotyped Hb8D in 13 east coast samples, from Maine to Florida [43]. The frequency ranged from 5% to 35% (Table 3) but there was not a significant relation between latitude and frequency of Hb8D (b = 20.006, p = 0.1). For comparison the s2e was also sequenced in the same individuals. Again, no unique SNPs are found on the Hb8D haplotype. Thus, nothing in in this broader N-American sample suggests complementary mutations in s2e. Curiously however, there is a significant reduction in s2e polymorphism with latitude (p = 0.02 for p and h). This does not explain the prevalence of Hb8D, but suggests geography (or history) affects variation in the regulatory regions of some developmental genes.

Testing for Effects of Hb8D on Viability and Developmental Time
Test crosses and analysis of inbred lines were used to gauge the putative impact of Hb8D on the number of offspring hatching and developmental time. Here developmental time is assessed as the time to eclosion (see methods).
Consistently with earlier studies [63,64] hemizygosity at eve reduces viability (Table 4) by about 20% in all crosses except to eve 5 (DF vs. Cy in Figure S3). However offspring number was not affected by the deletion of Hb8 binding site (Genotype term in Table 4). Number of hatching offspring differs between the four eve mutant stocks (Table 4) most likely due to varying genetic backgrounds. We also asked about factors influencing developmental time. The ANOVA's indicate difference among eve alleles, and potential effects of hemizygosity at the locus (Table 4). Most notably, Hb8D seems to reduce developmental time (Table 4)while hemizygosity at eve increases it. In three of the four crosses did Hb8D individuals develop significantly faster than the wt flies ( Figure 3). The Hb8D flies eclose on average 3.5 hours earlier, but again no effects are seen in eve 5 . This effect was also seen if the effect is estimated for sexes separately. In 13 of the 16 Mutation-Cross-Sex combinations Hb8D developed faster than flies with wt s3+7e, which is significant in a sign-test (binomial, p = 0.02). Note the Hb8D is tested in heterozygous form, thus in these crosses it appears to have dominant effects on developmental time.
We also examined the effects of Hb8D with association tests in 60 inbred lines. As before, Hb8D had no effect on offspring number. Peculiarly, the data do not confirm the association between Hb8D and developmental time (lower part of Table 4). The estimated developmental time is in the same range for both experiments suggesting they are not systematically different.
Together these data suggest an effect of Hb8D on developmental time, but further tests are needed to confirm or refute this.

Histochemical Staining of Eve Expression
Proximal phenotypes, like protein level at a specific time and location in the embryo, might be associated with functional variation in regulatory elements. To test this we stained for Eve in stage 14A embryos of four inbred lines, two Hb8D and two wt. Mixed model ANOVA shows that the relative positioning of the Eve stripe boundaries differs between the four inbred lines studied (Table S6). Both developmental stage and embryo orientation affect the anterior and posterior boundaries of stripes. Those sources of error were accounted for by i) working with the residuals after fitting the embryo orientation and ii) using developmental stage as a covariate. The average developmental stage does not differ between lines (p = 0.8), suggesting that rate of early development does not contribute to the line differences.
Hb repression establishes the anterior boundary of stripe 3 and posterior boundary of stripe 7 [62]. Thus, a priori, those features are most likely to be affected by Hb8D. However, the mixed model ANOVA does not indicate effects of the Hb8D on these stripe 3 and 7 boundaries (Figure 4). It is possible that this Hb site has broader function. The only putative signal in the data was with stripe 5; according to least square means stripe 5 is found more anteriorly in Hb8D. But this is not formally significant after Bonferroni correction for all 14 tests. A complementary analysis of principle components (PC) of the relative stripe positions does not implicate Hb8D in stripe positioning. The two largest principle components capture variation in (PC1) the central stripes and (PC2) the anterior -posterior axis of the embryo. The contribution of Hb8D to principle component 1 is not formally significant (F 1,10 = 4.25, p = 0.07). These results do not suggest that Hb8D affects Eve pattern in the early development.

Discussion
Sequence comparisons of close and more distantly related species show how TFBS emerge, change and get lost [8,65]. Is this turnover of functional sequences due to relaxed purifying selection, or does positive selection play a role [66][67][68]? There is substantial variation in gene expression among individuals and the bulk of expression QTLs map in cis [69][70][71]. The exact nature of those cis variants is rarely known, but a systematic review by Rockman and Wray [72] shows that SNPs, indels and length polymorphism in repeats can abolish TF binding and affect expression of neighboring genes.

Hunchback Site Polymorphisms are not Deleterious
Here we report that two large deletions segregating at moderate frequency remove predicted Hunchback binding sites in, and next to, the stripe 3 and 7 enhancer of eve. Both sites have high PWM scores and are evolutionarily conserved. One of them (Hb8) was characterized molecularly [55]. Three observations suggest that Hbs1, removed by a 45 bp deletion, is a true Hb binding site. It has among the highest PWM score of Hb sites in the eve region. It is evolutionarily conserved between D. melanogaster and D. erecta and resides less than 250 bp away from the Hb8 site. Stanojevic et al. [55] footprinted 4 Hb sites in the spacer between s2e and s3+7e, and recent thermodynamic models and quantitative measurements of TF abundances indicate that the spacer between s2e and s3+7e contains functional Hb motifs [73]. However functional assays are required to confirm that Hb binds to these two sites in vivo and modulates eve expression.
Our initial hypothesis was that these deletions of Hb binding sites are deleterious, as the loss or modulation of a single TFBS can have measurable effects [72,74,75]. This is refuted by several facts: 1) both mutations are at appreciable frequency, 2) individuals homozygous for each of those deletions survive as inbred stocks, 3) Hb8D has cosmopolitan distribution and 4) Hb8D does not seem to reduce viability and, if anything, it speeds up developmental time.
The genetic assays had sufficient statistical power to detect the effects of eve hemizygosity on offspring number (consistent with reported partial haplo-insufficiency at the locus [63,64]) and less so developmental time. Thus we conclude that the Hb8D is not strongly deleterious. The alternate scenarios are that the two deletions are either (nearly) neutral or favored by positive selection.
The most parsimonious explanation is that Hb8D is neutral and drifts in the population. This scenario is supported by haplotype tests, which do not point to the involvement of positive selection. However, the fact that the two deletions destroy binding sites for the same TF in the same enhancer is rather puzzling. Thus, it is tempting to hypothesize that the two Hb binding site deletions are favored by selection. Curiously, no other Hb sites in the s3+7e or s2e are affected, no substitutions are seen in more than 100 sequenced lines.

Variation in Early Development
Several studies have documented substantial variation in early Drosophila gene expression, with expression arrays [76], RNA seq [77] and in-situs [78]. As the deletions are found in s3+7e, it is most probable that they could affect Eve stripes 3 or 7. Hb is abundant in the anterior of the embryo, and drops adjacent to the anterior boundary of eve stripe 3. Hb is also produced in a narrower domain in the posterior, close to the posterior boundary of eve stripe 7 [62]. Hb demarcates the boundaries of those stripes (and stripes 4 and 6). Thus deletions of Hb sites would be expected to lead to an anterior shift of stripe 3 and posterior shift of stripe 7, because this regulatory module would be less sensitive to Hb repression (the absence of its full complement of binding sites). Our analysis of Eve expression in four inbred lines does not reveal effects of Hb8D on Eve stripe placement. Genetic and maternal factors affect the placement of expression boundaries; physical or environmental attributes like egg size do as well [78][79][80]. Note, lack of evidence does not prove the alternative. These results do not prove that the Hb8D does not affect Eve expression. The ideal test of the functionality of Hb8D and Hbs1D requires transgenic constructs in a common genetic background or homologous recombination into the eve locus of a particular line. It is unclear how such alterations would affect proximal or distal features of development. The quantitative tests suggest Hb8D acts dominantly, and speeds up development by , 3 hours. This seems unrealistic as the Eve pair rule pattern only takes ,50 minutes to mature [81], thus it is impossible that these effects (if real) are due to Hb and eve interaction during early development. But curiously both eve and hb also play a role in the developing neuronal system [82,83] but the functional interaction of Hb and eve in those tissues is largely unexplored. In the absence of functional or genetic confirmation we argue for cautious interpretation of the observed association of Hb8D and developmental time in the test-crosses. Finally, it is also possible that these deletions affect proximal developmental events, but that those effects are a minute or acceptable noise in the system.

Can Co-evolution Explain the High Frequency Hb TFBS Deletions?
Co-evolution can occur via neutral changes (e.g. in the network neighborhood [84]) or via positive selection favoring compensatory changes in the genome. Here two co-evolution models that may account for these two Hb binding site deletions in eve are entertained. Those are i) cis-changes within eve or, ii) trans-changes in the function or abundance of activators and/or repressors.
First, the relatively high frequency of those two deletions could reflect co-evolution within eve. Hunchback acts both as a transcriptional activator and repressor during development [85][86][87]. Hb positively influences expression via the eve stripe 2 enhancer, but is part of two-tier repressor system that demarcates the boundaries of stripes 3, 4, 6 and 7 [62]. Stripes 3 and 7 are known to be activated by D-stat [32], an ubiquitously available activator (other agents may also play a role). The high frequency of Hb binding site deletions could be a co-evolutionary response to increased activation of stripes 3 and 7 expression, for instance via altered Dstat binding. This is unlikely as the two D-stat sites in s3+7e have not diverged between D. melanogaster and D. erecta (Table  S2) and no polymorphism is found in those sites within D. melanogaster. Binding sites for other agents activating eve stripes 3 and 7 may have changed; TFBS that could reside elsewhere in regulatory regions around eve. The eve regulatory region is 85-95% identical between D. melanogaster and D. simulans. We scanned the eve region of both species with Jaspar [56], and found hundreds of TFBS differing between the species (data not shown). Nonetheless, no changes in Hb or Dstat sites were found. It is also possible that miRNA docking sites or other regulatory elements in eve have changed, thus leading to selection for higher frequency of those two Hb site deletions.
Alternatively, changes in structure or function of trans-factors, like Hb itself, may have led to the increased frequency of those two Hb binding site deletions. It is improbable that a protein change is responsible, as the differences between the D. melanogaster and D. simulans Hb proteins are all on the D. simulans branch (unpublished results, Dagmar Yr Arnardottir and Arnar Palsson). We find it more plausible that the spatial or temporal amount of trans-factors has changed, for instance a lower amount of Dstat. The most intuitive scenario is, quantitative, temporal or even spatial changes in Hb expression in the embryo -which may have prompted co-evolution in regulatory elements sensitive to quantitative changes in Hb amount in development. The eve s3+7e might be such a critical Hb-target element. This is of course speculation, but in this scenario, one would expect that other Hb target enhancers, which produce expression overlapping the spatial and temporal patterns of eve s3+7e might also have experienced altered selection pressure. Thus, other Hb such target genes could also exhibit point mutations or deletions of conserved and presumably functional Hb binding sites. Note, we are not arguing positive selection is necessarily responsible; changes in Hb dose could lead to relaxation of selection for a subset of Hb target genes, and thus previously detrimental mutations in these genes could drift to higher frequency.

Conclusions
The genetic network governing early Drosophila development has been used to discover many of the basic principles of developmental genetics, regulatory DNA function and regulatory evolution [6,10,88,89]. Recent technical and analytical improvements have enabled quantitative analyses of enhancer function and logic [87,88,[90][91][92] and dosage compensation [77,93]. Developmental networks must cope with variation due to chance, the internal and external environment, and in the relevant genetic components. Studies point to the involvement of positive selection in the gain and loss of TFBS in Drosophila [66,94] and co-evolution within enhancers [39,95]. Furthermore, non-clocklike evolution of the s2e from four Drosophila species [11], indicates co-evolution of TF abundance and functional elements in cis-regulatory modules. The fact that two large deletions removing TFBS for Hb are found in close proximity in a regulatory element, might be an example of such co-evolution. However we favor the cautious explanation that these high frequency deletions reflect developmental system drift [40,41], i.e. permitted deviations in parameters of the Drosophila developmental regulatory network.     in these species. The order reflects approximately phylogenetic relationship available on http://insects. eugenes.org/species. There is length variation in T stretch between Kni5 and Hb11c; extra 1 and 2 bases in D. sim and D. gri respectively. As these are monomorphic stretches the core binding sites are presumably not affected. (DOC) Table S4 Little evolutionary conservation of a putative sloppypaired site in the eve stripe 3+7 enhancer. Full species names and accession numbers are listed in material and methods. Orthology of the sloppy-paired binding site region was determined by colinearity of binding sites in the stripe 3+7 region, were Hb8 and Hb9 flank the sloppy-paired binding site. Fewer than 50 bp separated Hb8 and Hb9 in all species. The exception is D. ananassae, were Hb8 was not detected.  Supporting information S1 Alignment of population sequencing of a part of evenskipped stripes 3+7 enhancer from North Carolina, in fasta format.

(TXT)
Supporting information S2 Alignment of population sequencing of a part of the hunchback regulatory region from North Carolina, in fasta format. (TXT)