GGGT, nullo-2: GTCGGAGCAGGCAACGGGCAT, sema1: TCGTCGGTAGGTAAAAGTTGT, tat: CCAGCCGCAGGTA- TTTAGTTC, zen1m: CACTATTTGAATAGACACTGT, tatc- gat: TCACTACTATCGATGACGATG (TAGteam or TATC- GAT sites are underlined). Tiled genomic microarray design A Drosophila melanogaster tiled genomic microarray set was designed by J.R.M. in concert with the bioinformatics team at Roche NimbleGen using Genome Release 5. This array set, which comprises two HD2 (2.1 million feature) microarrays, utilizes 50- mer oligonucleotide probes with up to 100 close matches per sequence tolerated, and a median probe spacing of 33 bp. For our design, we chose to tolerate a large number of close matches in order to include on the array more heterochromatin and repeat regions, including transposon sequences. The Design Names of the arrays are as follows: 081229_Dm_JM_ChIP_1_HX1 and 081229_Dm_JM_ChIP_2_HX1. Transcriptome analyses using gene arrays Total RNA was extracted from three independent collections of 2–3 hr yw and zld2 Drosophila embryos by TRIzol (Invitrogen). cDNA was prepared using the GeneChip HT One-Cycle cDNA Synthesis Kit (Manufactured by Invitrogen for Affymetrix) and labeled with the BioArrayTM HighYieldTM RNA Transcript Labeling Kit (Enzo). Labeled probes were hybridized to Drosophila Genome 2 Affymetrix arrays and processed by a GeneChip Fluidics Station 400. Data were acquired by the GeneChip Scanner 3000 and processed/normalized by the Affymetrix GeneChip Operating Software (GCOS). Genes were identified as present when at least two of the three replicates had present (P) assignment (p,0.05). t-test analysis was performed on the data from the three biological replicates. The fold change of each gene was determined by the ratio of yw mean/zld2 mean. Transcriptome analyses using tiling arrays Double-stranded cDNA was generated from total RNA isolated (TRIzol, Invitrogen) from staged (1–2 hrs and 2–3 hrs at 25uC, verified by DAPI staining of a portion of the collected embryos) wild-type and zld2 embryos, then amplified/labeled using either Cy3- or Cy5-coupled random nonomers, respectively. 15 ug of each labeled cDNA were used for competitive hybridization, coupling the same staged wild-type and zld2 samples in order to facilitate direct comparison (see Roche NimbleGen Gene Expression Protocol [http://www.nimblegen.com/products/lit/ expression_userguide_v5p0.pdf] for specific labeling and process- ing details). After hybridization for 16–20 hrs, the arrays were washed, dried and then scanned on an Axon GenePix 4000B microarray scanner from Molecular Devices. Intensity readings of probes were corrected according to their GC content using a set of random probes on the arrays, and normalized by the Lowess normalization method [49]. After normalization, a median filter was applied (using a sliding window of three probes with the center probe given the median value). The expression level of each exon was calculated by taking the median of all the probes covering the exon, and the expression level of each RNA isoform was calculated by averaging the expression levels of all the exons of the isoform without weighting. Using the ‘‘fdrtool package’’ in R [50], we set the background threshold as 5% FDR, which was previously shown to be an appropriate cutoff for tiling array data [51]. Genes were considered as expressed (present) if more than 50% of the probe signals were higher than the threshold, otherwise they were considered as not expressed (absent). Using this approach, we concluded that 49.5% of the Zelda Coordinates Early Gene Networks genome was represented as RNA in 1–2 hr embryos, which is similar to findings from other studies [4,52,53]. Since, in most cases, each gene was represented by multiple probes, and thus was interrogated multiple times in one experiment, we applied the t- test to obtain p-values for each gene in the yw and zld2. samples. The fold change of each transcript was determined as the ratio of yw/zld2. ChIP–chip 1–2 hr yw embryos were fixed in 2% formaldehyde for 20 min. Nuclei were harvested and sonicated to release and shear genomic DNA. The chromatin was immunoprecipitated by purified anti- Zld antibody and protein A beads. The chromatin immunopre- cipitated (ChIPed) DNA samples were amplified by GenomePlex WGA1 (Sigma) followed by GenomePlex WGA3 (Sigma) twice. ChIPed DNA was labeled/amplified with Cy3-coupled random nonomers and the corresponding input DNA was labeled with Cy5-coupled random nonomers; 34 mg of each labeled DNA was combined and co-hybridized to the Drosophila tiling array set (as described above) according to the Roche NimbleGen ChIP-chip Protocol (http://www.nimblegen.com/products/lit/chip_userguide_ v6p1.pdf). After hybridization for 16–20 hrs, the arrays were washed, dried and then scanned on an Axon GenePix 4000B microarray scanner from Molecular Devices. For each dataset, the ratio of ChIP to input intensities was obtained, log2 transformed, and standardized to z scores. The datasets of the dye swap replicates were averaged to eliminate dye bias. ChIP enrichment scores were calculated using the R package ‘‘Ringo’’ [54] as the sum of probe levels minus the threshold with settings of 0.05 p-value and 10 minimum probes. A higher score indicates a higher ChIP/input ratio in the region. Data was visualized using the Integrated Genome Browser [55], and a median filter was applied to the ChIP/input ratio for the visualization. Validation was performed by qPCR to assess the enrichment of the enhancer region of sc/sisB, a known target of Zld, and the coding region of CG18125, a negative control, in ChIPed versus input DNA. The sc/sisB region showed 23.53 fold enrichment over CG18125 (data not shown). Enrichment test for Zld binding site The enrichment index for a specific DNA motif was calculated as the density of the motif in Zld-bound regions divided by the density of that motif in the Drosophila genome. Enrichment indices of sites were calculated in 100 bp non-overlapping windows across the 5 kb flanking regions from the center of Zld-bound regions. The background was estimated as the average of enrichment indices of 20 random heptamers. TAGteam site PWM for Zld recognition We took seven CAGGTAG related heptamers, including five that were previously identified by ten Bosch et al. [2]: CAGGTAG, TAGGTAG, CAGGTAA, CAGGCAG and CAGGTAT, and 2 newly discovered related heptamers from the Zld-bound region upstream of nullo: CAGGCAA and CAGGTAC, and then weighted every heptamer according to their enrichment indices to generate a primary PWM using the ‘‘Biostrings’’ package [56] and visualized by the ‘‘seqLogo’’ package [57] in R. The enrichment analysis with this PWM defined sequences with a p- value larger than 0.0005 as not enriched in Zld-bound regions. We then generated a new PWM by using the enriched sequences (p#0.0005) with their surrounding nucleotides. We repeated the enrichment test using this new PWM (Figure 2B) and eight TAGteam motifs (p#0.0003). PLoS Genetics | www.plosgenetics.org 12 October 2011 | Volume 7 | Issue 10 | e1002339