Searching for novel enriched sites in Zld-bound regions Enrichment indices of all possible heptamers were calculated for Zld-bound regions using a 500 bp window centered around the middle of Zld-bound peaks. Using an enrichment score of 3.5 as the cutoff, 11 heptamers were recovered (Figure S1B). The highest-ranking heptamer was CAGGTAG. Heptamers could be separated into two groups with different core sequences, CAGGTA and TATCGA (heptamers containing simple repeats were not further analyzed). The CAGGTA group was reminiscent of TAGteam sites. All of the Zld-bound regions that contain at least one of the new enriched sites that have the TATCGA core were analyzed by MEME4.4.0 [58] and a PWM was generated. Searching for genes associated with Zld-bound regions The distance from the center of every Zld-bound region to every annotated TSS was calculated (according to Drosophila melanogaster genome release 5.29). The closest TSS (for both strands) was assigned to each Zld-bound region. Thus each bound region typically has two assigned genes. Determination of overlap between Zld-bound regions and hotspots Hotspot data was obtained from Roy et al. [5]. Only the regions bound by eight or more transcription factors were analyzed in this study. Hotspots were ranked by their complexity scores [5]. In the 100 non-overlapping hotspot windows, Zld binding was calculated as the average of ChIP/input ratios of all probes located in the hotspot regions. To estimate the background of Zld binding in hotspots, the ChIP/input ratio for all probes on the tiling arrays were randomly shuffled 20 times, and the average of the ChIP/ input ratio from the shuffled probes corresponding to a hotspot was calculated. Determination of binding overlap between sets of transcription factors Target gene lists for AP factors (Bcd, Cad, Hb, Gt, Kr, Kni and Tll) were obtained from MacArthur et al. [14]. The target-gene list for Dl, Twi, and Sna (DTS) was obtained from Zeitlinger, et al. [35]. We limited the search to target genes that lie within 2 kb of the factor-bound regions. We identified the target genes bound by both Zld and each AP factor, Zld and DTS, and Zld and Mad, and distinguished whether or not they were expressed at 2–4 hrs (defined as bound by pol II, [14]). Pair-wise heat maps were generated to represent the overlapping percentages between factors. GO term enrichment analysis GO term analysis was performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.7 [59,60; http://david.abcc.ncifcrf.gov/] on the genes associated with Zld-bound regions. We limited this analysis to genes that are downstream of a Zld-bound region, except in cases where there was no downstream gene then we used the upstream gene. Supporting Information Figure S1 (A) Conservation of nucleotides around additional TAGteam and TATCGAT related sites. Conservation analysis (phastCons 15 ways) [62] of sequences encompassing the eight TAGteam sites, TATCGAT, TAGteam PWM (p#0.0003), TAT PWM (p#0.0006). The mean and median conservation scores of nucleotides were calculated for each individual motif within or outside the Zld-bound regions. The X-axis shows relative positions Zelda Coordinates Early Gene Networks to the site. Red and light red lines represent the median and mean conservation scores of sites within Zld-bound regions, respectively. Blue and light blue lines represent the median and mean conservation scores of sites outside Zld-bound regions, respectively. The TAGteam sites CAGGTAN that are within bound regions have higher conservation scores than those outside the bound regions, demonstrating evolutionary constraint. Interestingly, these correspond to the most over represented motif in hotspots, CAGGTA [5]. The TATCGAT site showed no significant difference. Sites that match the TAT PWM within Zld-bound regions are not more conserved than sites outside the bound regions. (B–D) Searching for novel enriched sites in Zld-bound regions. (B) Enrichment indices of all possible heptamers were calculated for Zld-bound regions using a 500 bp window centered over the middle of Zld-bound peaks. Eleven heptamers showed an enrichment index greater than 3.5. X-axis, -5 kb to +5 kb from the center of the bound regions. Y-axis, enrichment score. (C) Heptamer sequences are listed (left) with their enrichment indices in the 500 bp window at the center of the Zld-bound regions (middle), and the peak level of enrichment (right). (D) The heptamers could be separated into 4 groups, represented as different colors (see B) and aligned: CAGGTAG-related in blue, TATCGAT-related in red, GTCA- CAC-related in green, and CT repeats in purple (not shown in D). (E) Zld binding scores of the different enriched sites. The histogram shows the numbers of Zld-bound regions with at least one enriched motif counted in 10 windows according to the Ringo binding scores. X-axis, binding scores in 10 percentile non-overlapping windows with the highest scores in the 100-90 window. Y-axis, the number of Zld-bound regions. Bound regions containing TAGteam sites, especially CAGGTAG (black) and CAGGTAA (dark blue), tend to have higher binding scores. Bound regions with the newly discovered TATCGAT motif (red) were less enriched in the top 20% window. Bound regions without TAGteam or TAT sites (brown) tended to have lower binding scores. (F) Overlap of Zld- bound regions and hotspots. IGB views of the brk, slp1/slp2, odd, and Antp gene regions. Hotspots are shown as orange rectangles; height reflects the hotspot score, which depends on the number of factors bound and the strength of their binding [5]. Zld binding peaks are shown in blue with significance scores from the Ringo algorithm (see Materials and Methods) shown above as blue rectangles. Below the Zld peaks are: TAGteam sites (purple lines; limited to CAGGTAG, TAGGTAG and CAGGTAA), cis-regulatory modules (CRMs from REDfly [16]), and gene models using genome version BDGP R5/ dm3. Note the striking overlap of Zld-bound peaks and hotspots [5]. Also, several Zld-bound peaks that did not make the Ringo signficance cutoff (see Materials and Methods) overlapped with hotspots. (PDF) Figure S2 Expression of Zld-bound genes. The location of the center of all Zld-bound regions relative to the closest TSS (0 point on X-axis; numbers represent bp downstream (+) or upstream (2) of the TSS) was plotted against the log2-expression ratio wt/zld2 (A) or the expression level in wild-type embryos (B) of those genes. Dotted lines mark 22 kb and +2 kb. Zld often binds within 2 kb of the TSS of genes irrespective of how they are affected in the expression profiling assays (A). However, there was a significant difference (Fisher’s exact, p,0.0001) between the percentage of genes bound by Zld within 2 kb of their TSS for genes considered expressed (78.5% in blue) versus those considered not expressed (45.3% in red). (TIF) Figure S3 High Zld binding scores are associated with early developmental genes. (A) GO terms significantly enriched among PLoS Genetics | www.plosgenetics.org 13 October 2011 | Volume 7 | Issue 10 | e1002339