screening for 8 bp sequences with a matrix score of at least 70. This search motif ought to detect all HSEs, which will include this sequence regardless of whether NTTCN or NGAAN occupies the most distal position. HSEs longer than the search motif were obvious as they generated multiple hits in the same region. All HSEs identified in the initial motif search were further screened prior to retention for analysis based on the following criteria derived from experimental data: (1) they must contain at least 3 contiguous alternating 5 bp subunits matching the canonical HSE, (2) they must have no more than 2 substitutions and (3) they must have no mismatch to G in position 2 (NGAAN type) and C in position 4 (NTTCN type), unless the 5 bp mismatched subunit was interior to an HSE that otherwise contained at least 4 subunits. 10 bp gaps were not tolerated under any circumstance. A final screen was performed to quantify the likely binding strength of each HSE subunit from its sequence similarity with the canonical motif, and to define boundaries of HSEs. We generated WebLogos [76] for each of the alternative 5 bp subunits recovered in the motif search. Bit-scores were summed for each base at the 5 positions in each subunit, and expressed as a percentage of the sum calculated if the preferred base were at every position. Subunits scoring ,50% and at an external position in a putative HSE were removed from further analysis. We specifically retained gapped HSEs in the dataset by including subunits scoring ,50% but flanked by higher-scoring subunits. Lastly, we used the 12-species genome alignment at UCSC supplemented with multiz alignments of local regions to assess the orthology of individual HSEs. Many HSEs detectable in this search were homologous to previously identified and, in some cases, experimentally confirmed HSEs in D. melanogaster, adding confidence to their assignment. A final screen was performed to test for the possible presence of HSEs consisting of only 2 5 bp subunits [25]. To exclude false positives, we applied a stringent threshold: only motifs conserved across multiple species with a maximum of one substitution (excluding positions G2 and T4) from the canonical HSE motif were retained. This search yielded the final collection of HSEs. We compared rates of divergence among species in HSE sequence with a neutral proxy, synonymous changes at 3rd positions of adjacent heat-shock coding regions. We performed this calculation among all orthologous HSEs and coding regions for all single copy genes except Hsp27 to ensure adequate sample References 1. Wittkopp PJ (2006) Evolution of cis-regulatory sequence and function in Diptera. Heredity 97: 139–147. 2. Wray GA, Hahn MW, Abouheif E, Balhoff JP, Pizer M, et al. (2003) The evolution of transcriptional regulation in eukaryotes. Mol Biol Evol 20: 1377–1419. 3. Crocker J, Tamori Y, Erives A (2008) Evolution acts on enhancer organization to fine-tune gradient threshold readouts. PLoS Biol 6: 2576–2587. 4. Zhang Z, Gerstein M (2003) Of mice and men: phylogenetic footprinting aids the discovery of regulatory elements. J Biol 2: 11. 5. Crocker J, Erives A (2008) A closer look at the eve stripe 2 enhancers of Drosophila and Themira. PLoS Genet 4: e1000276. 6. Dermitzakis ET, Bergman CM, Clark AG (2003) Tracing the evolutionary history of Drosophila regulatory regions with models that identify transcription factor binding sites. Mol Biol Evol 20: 703–714. 7. Kim J, He X, Sinha S (2009) Evolution of regulatory sequences in 12 Drosophila species. PLoS Genet 5: e1000330. 8. Lindquist S (1980) Varying patterns of protein synthesis in Drosophila during heat shock: implications for regulation. Dev Biol 77: 463–479. 9. Feder ME, Hofmann GE (1999) Heat-shock proteins, molecular chaperones, and the stress response: evolutionary and ecological physiology. Annu Rev Physiol 61: 243–282. 10. Lis JT, Xiao H, Perisic O (1990) Modular units of heat shock regulator regions: structure and function. In: Morimoto R, Tissieres A, Georgopoulos C, eds. Stress Proteins in Biology and Medicine. Cold Spring Harbor, NY: Cold Spring Harbor Press. pp 411–428. 11. Wu C (1995) Heat shock transcription factors: structure and regulation. Annu Rev Cell Dev Biol 11: 441–469. Phylogeny Orders HSE Evolution sizes. To minimize the issue of substitutional saturation, we based our calculations on transversions only. We performed 2-way x2 tests on all possible position pairs for the 5 positions within a subunit, using the entire pool of 1,611 identified motifs. The test was also performed on a corrected pool of motifs, to account for overrepresentation of HSEs from Hsp70, and other duplicates. The nGAAn-type motifs were tested separately from the nTTCn-type motifs, as these variants differ in binding ability [26]. The expect values for two given bases in their respective locations was the product of their rate of occurrences at those two positions. Results indicated where excesses and deficiencies arose in the data set. Three-way x2 tests were not performed due to insufficient sample size. Supporting Information Table S1 Extracted sequences of all HSEs and HSE motifs. All identified HSEs are indicated along with gene and species of origin. Individual motifs that constitute each HSE are also listed. Found at: doi:10.1371/journal.pone.0010669.s001 (0.14 MB XLS) Table S2 Annotated HSEs in heat-shock promoters. All individual HSEs are indicated within promoters from all species and all genes. HSEs are color coded by distance to the transcription start site. Found at: doi:10.1371/journal.pone.0010669.s002 (0.23 MB DOC) Acknowledgments The authors would like to thank the members of the laboratory groups of Martin Kreitman and Kevin White at the University of Chicago for comments and discussion on matters relevant to the current study, Marcos Antezana for assistance with the analysis, and Patricia Wittkopp for comments on the manuscript. Author Contributions Conceived and designed the experiments: ST MEF. Performed the experiments: ST MEF. Analyzed the data: ST RAH MEF. Contributed reagents/materials/analysis tools: MEF. Wrote the paper: ST RAH MEF. 12. Yao J, Munson KM, Webb WW, Lis JT (2006) Dynamics of heat shock factor association with native gene loci in living cells. Nature 442: 1050–1053. 13. Kremer SB, Gross DS (2009) SAGA and Rpd3 chromatin modification complexes dynamically regulate heat shock gene structure and expression. J Biol Chem 284: 32914–32931. 14. Anckar J, Sistonen L (2007) Heat shock factor 1 as a coordinator of stress and developmental pathways. Adv Exp Med Biol 594: 78–88. 15. Takemori Y, Enoki Y, Yamamoto N, Fukai Y, Adachi K, et al. (2009) Mutational analysis of human heat-shock transcription factor 1 reveals a regulatory role for oligomerization in DNA-binding specificity. Biochem J 424: 253–261. 16. Jedlicka P, Mortin MA, Wu C (1997) Multiple functions of Drosophila heat shock transcription factor in vivo. EMBO J 16: 2452–2462. 17. Vuister GW, Kim SJ, Orosz A, Marquardt J, Wu C, et al. (1994) Solution structure of the DNA-binding domain of Drosophila heat shock transcription factor. Nat Struct Biol 1: 605–614. 18. Gasch AP, Moses AM, Chiang DY, Fraser HB, Berardini M, et al. (2004) Conservation and evolution of cis-regulatory systems in ascomycete fungi. PLoS Biol 2: 2202–2219. 19. Harrison CJ, Bohm AA, Nelson HCM (1994) Crystal structure of the DNA binding domain of the heat shock transcription factor. Science 263: 224– 227. 20. Amin J, Ananthan J, Voellmy R (1988) Key features of heat shock regulatory elements. Mol Cell Biol 8: 3761–3769. 21. Perisic O, Xiao H, Lis JT (1989) Stable binding of Drosophila heat shock factor to head-to-head and tail-to-tail repeats of a conserved 5- bp recognition unit. Cell 59: 797–806. PLoS ONE | www.plosone.org 14 May 2010 | Volume 5 | Issue 5 | e10669