nucleotide states at pairs of sites (Figure 12) suggests that certain combinations of nucleotides outside the highly conserved regions favor or disfavor HSF binding. This finding suggests that the influence of nucleotide composition at individual sites in a HSE on binding affinity may be non-additive. Interestingly, significant correlations among states at different sites occurred in only the nGAAn-type motif. In conjunction with the finding of a bias towards HSEs with nTTCn motifs at the most distal site (tail-to- tail), this suggests the intriguing possibility that HSF-binding and heat-shock gene expression could differ mechanistically depending on HSE orientation in Drosophila. While our sample of HSEs is both larger and phylogenetically more diverse than previously available, it reaffirms rather than alters the consensus sequence previously derived from more limited samples of HSEs. Importantly, our results overall suggest some caveats on the uses of this consensus sequence. For example, in silico searches for HSEs (e.g., [36] should accommodate gapped HSEs and recognize that sites outside the central 3 positions of each 5 bp unit may consequential. Such practices may reconcile the apparent contradiction between experimentally determined binding sites of HSF and predicted binding sites. While we have posited that pre-existing detailed understanding of a CRE’s function ought to facilitate the analysis of its evolution, we may have underestimated the level of detail that is necessary. Although the heat-shock promoter is among the best understood of eukaryotic promoters and piecemeal experimental investigations of numerous variants establish that most are consequential for heat-shock gene expression (see Introduction), a precise functional read-out of evolutionary variation in the heat-shock promoter [as, for example, [66] have attempted at the level of the cis-regulatory module] is currently difficult to impossible to obtain experimen- tally, even with recent advances in throughput [67]. As we have reviewed, the number, arrangement, sequence, and position of HSEs is consequential for heat-shock gene expression. Thus a single 365 bp-unit HSE represents 415 potential combinations of nucleotides, as even the highly-conserved positions sometimes vary when a canonical 5 bp-unit is adjacent to it. The HSE may be tail- to-tail or head-to-head. HSEs may be larger than 365 bp units, presenting both more combinations and increasing potential for cooperative binding of HSF. Separation between the HSE and the TSS can vary. Multiple HSEs can occur within a single heat-shock promoter, with identical HSEs having different impacts on heat- shock expression depending on their position and arrangement. Finally, mechanisms in addition to HSF binding to HSEs, both in cis and in trans, regulate heat-shock expression [68]. The specific patterns observed in the 145 promoters of the 13 species presumably are feasible for heat-shock gene expression. Which unobserved patterns are not feasible and which are feasible but have not evolved? Of the observed arrangements, which encode higher vs. lower levels of gene expression (and what levels), and which are functionally equivalent? Elucidation of the cis-regulatory code rests on the answers to such questions. The enormous complexity of even such a seemingly simple and well-understood promoter as the heat-shock promoter, however, suggests that this elucidation may be among the ‘‘grand challenges’’ of biology [69]. While ‘‘nothing in biology makes sense except in the light of evolution’’ [70], sometimes even this light is insufficient to make sense. Materials and Methods Sequences for 8 D. melanogaster heat-shock genes (Hsp22, Hsp23, Hsp26, Hsp27, DnaJ-1, Hsp68, Hsp70, and Hsp83) were obtained from FlyBase (http://flybase.org). We screened the highest scoring Phylogeny Orders HSE Evolution results of BLASTn and tBLASTn searches with D. melanogaster nucleotide coding and protein sequences using sequence homol- ogy, synteny, FlyBase orthology assignments, and reciprocal BLAST to D. melanogaster to identify orthologs in the other 11 species of Drosophila with sequenced genomes. Several putatively single-copy heat-shock genes had lineage- specific duplicate copies, which underwent further analysis. Six copies of Hsp70 are in the D. melanogaster genome. We used BLASTn with the coding sequence of D. melanogaster Hsp70Aa as a query to each of the other 11 sequenced genomes. All hits with an e-value of 0.0 were further examined to determine whether they constituted a putative copy of Hsp70, using a combination of Gbrowse orthology calls, synteny and relative position, and reciprocal BLAST to D. melanogaster. Hsp70 homologs were generally found in high-scoring clusters of .75% identity to D. melanogaster Hsp70, often in tail-to-tail configurations. Hsp68 and Hsp70 were distinguishable from diagnostic sites [71]. Analyses incorporated all verified copies found in each species [64], except where specifically noted. The identification of regulatory elements is dependent on the recognition of the promoter region upstream of the TSS. We chose a 1000 bp region that should encompass the majority of proximal promoter regulatory elements [72]. Transcriptional start sites were first determined in D. melanogaster using experimental data collected in the Eukaryotic Promoter Database [73] and extrapolated to other species by referencing the whole genome alignments via the UCSC Genome Browser (http://genome.ucsc.edu/). The 1000 bp region upstream of the conserved TSS was then extracted for each gene in each species. In several cases this process was complicated by gene duplication or similarity. First, assessment of TSS location was confounded for Hsp68, as several species appear to have copies of Hsp70 aligned to the D. melanogaster Hsp68 gene region. We located the UCSC aligned regions in the genomes of species where the alignment appeared incorrect, and confirmed that the predicted coding sequence of the adjacent gene was Hsp70 through comparisons to diagnostic sites [71]. Local multiz alignments obtained from the PromAn server [74] were exclusively used for the promoter region of Hsp68 to define the TSS in each species. Second, for several lineage-specific duplications or for the multi-copy gene Hsp70, only a single paralog is aligned at UCSC. We used local alignments to determine placement of the TSS in all paralogs. Third, Hsp22 is duplicate in D. grimshawi, but only one copy has identifiable HSEs in the region upstream of the putative TSS and was included in analysis. We obtained heat-shock promoter sequence from D. santomea strain CAR1600 [40] with primers designed from the closely related D. yakuba [63]. DNA was extracted from single individuals using a DNeasy kit (Qiagen, Rockville, MD), and PCR was performed using the following conditions: 94uC for 2 min., 35 cycles of 94uC for 30 sec., 50–55uC for 30 sec., and 72uC for 1 min, with a final extension at 72uC for 6 min. Sequence reads were trimmed and edited using Sequencher 4.5.In D. yakuba, Hsp70 occurs in two tail-to-tail clusters. One set of primers enabled sequencing of the entire intergenic region of one cluster (GE24511 and GE26203 in D. yakuba); the second cluster (GE26149 and GE24569) proved refractory to PCR amplification. As D. santomea is not included in the UCSC genome alignment, we aligned the TSS from this species based on the closely related D. yakuba for promoter region delineation. The defined promoter region upstream of the aligned TSS in all promoters was used as the search region for HSEs. First, we employed a motif finder, part of Mod Tools [75], to locate potential HSEs by searching for the 8-nucleotide sequence TTCNNGAA, allowing for a maximum of 4 mutations and PLoS ONE | www.plosone.org 13 May 2010 | Volume 5 | Issue 5 | e10669