Computation-Based Discovery of Related Transcriptional Regulatory Modules and Motifs Using an Experimentally Validated Combinatorial Model

The eve MHE

We recently described a comprehensive model for the transcriptional integration of multiple intercellular signals that act together to establish the identities of a subset of muscle and cardiac progenitor cells in the dorsal mesoderm of the Drosophila embryo (Fig. 1A–C; Halfon et al. 2000). In this model, expression of the progenitor identity gene even skipped (eve) is regulated via the action of at least five TFs binding to multiple sites in a 312 bp enhancer located approximately 6 kb downstream of the eve coding region, the Muscle and Heart Enhancer (MHE). Three of the TFs, that is, dTcf, Mothers Against Dpp (Mad), and Pointed (Pnt) function downstream of and are activated by the Wingless (Wg), Decapentaplegic (Dpp), and Ras/MAP kinase signaling pathways, respectively, while the other two TFs, Twist (Twi) and Tinman (Tin), are mesodermally restricted selector proteins. The transcriptional code comprised by these five transcription factors can account for all of the signaling events and genetic data known to affect eve regulation in the mesodermal progenitors (although it is likely that additional factors remain to be identified; see below).

The MHE Sequence is Well Conserved in a Distantly Related Drosophila Species

As part of our continuing analysis of the MHE, we cloned and sequenced the corresponding region downstream of the eve gene in the distantly related (∼40 Mya; Kwiatowski et al. 1994) species Drosophila virilis. We found that the D. virilis element (vMHE) maintains extensive sequence conservation with the MHE (Fig. 1D,E), and is functional when used to drive reporter gene expression in D. melanogaster, albeit more weakly than the MHE (data not shown). Consistent with our model of MHE function, at least one representative of each relevant TF binding site can be found in the vMHE, although not every functional MHE sequence motif is conserved (e.g., Ets4, Fig. 1E; Halfon et al. 2000). Importantly, several blocks of conserved sequence can be observed in which no TF binding sites have been characterized to date (Fig. 1E, gray shading), suggesting that additional but presently unidentified factors may be involved in MHE-derived gene regulation. Of note, sequences included in a proposed variant eve mesodermal enhancer (Knirr and Frasch 2001) are not conserved in D. virilis (Fig. 1D, blue arrow).

A Computational Search for Dorsal Mesodermal Enhancers

The combinatorial Wg+Dpp+Ras/MAP kinase signaling code necessary for eve transcription appears to be necessary not only for Eve expression but also for that of other dorsal mesodermal genes (Carmena et al. 1998, 2002; Halfon et al. 2000). We thus hypothesized that the same transcriptional code—Twi, Tin, dTcf, Mad, and Pnt—might also be responsible for the regulation of other genes. Moreover, we reasoned that, were this the case, this code might act via a cis-regulatory element resembling the MHE in that it would contain clustered binding sites for all five TFs. To test this idea comprehensively, we undertook a computational search of the entire Drosophila genome to identify elements that, based on their related composition, would be predicted to function as MHE-like dorsal mesodermal enhancers (DMEs).

The genome was searched to locate regions which contain predicted binding sites for all five TFs found within the eve MHE. The search was conducted essentially as follows (see Methods): the ScanACE program (Hughes et al. 2000) was used to identify each occurrence in the genome of a predicted binding site for each of the five TFs. An algorithm was then run to detect all instances in which all five TFs are found within a 500 bp window; these regions were termed “elements.” The 500 bp window size was chosen in recognition of the ∼300 bp size of the model (the MHE) and in an effort to limit the number of positive returns from the search so as to make feasible testing of our predictions.

This initial search resulted in 647 elements (Table 1 and Supplemental Table 1 available online at http://www.genome.org). To determine how likely it was that these would occur merely by chance, given the number of each predicted binding site in the genome, Monte Carlo simulations were performed in which the predicted sites were randomly distributed, and the cooccurrence search was run (see Methods). The results from these simulations indicate that cooccurrence of the five motifs happens significantly less than expected at random (P < 0.0001; Supplemental Fig. 1 available online at http://www.genome.org). This finding suggests that such cooccurrences have functional consequences that have acted to keep the individual motifs separated in the course of evolution except when needed for the expression of specific genes.

cis-Regulatory elements directing dorsal mesodermal expression should map in proximity to mesodermally expressed genes. We used the gene-expression annotations present in Flybase (Flybase 1999; http://flybase.bio.indiana.edu) to assess how many of the predicted elements were located either within the introns of or adjacent to genes with known mesodermal expression. The Flybase annotations indicate a strong bias for mesodermally expressed genes mapping near the elements identified by our search (P = 0.0001; Supplemental Table 2 available online at http://www.genome.org). This enrichment for mesodermal genes is consistent with our hypothesis that cooccurrence of the five TFs might be predictive for dorsal mesodermal regulatory elements.

cis-Regulatory modules frequently contain multiple occurrences of a given TF binding site (Arnone and Davidson 1997). As a first step in the analysis of our search results, therefore, we filtered the elements to select only those that occur in noncoding sequence and in which each TF binding motif is present at least twice, with the exception of that for dTcf, which is found only once in the model sequence, the MHE. We refer to this set as the “selected subset.” The selected subset contains 33 putative DMEs, in addition to the Eve MHE, for a total of 34 elements (Table 1).

DME-2 Functions as a Dorsal Mesodermal Enhancer for the hbr Gene

Element DME-2 maps to an intron of the heartbroken (hbr) gene (also known as stumps and dof; Michelson et al. 1998; Vincent et al. 1998; Imam et al. 1999). As this gene is known to be expressed in the embryonic dorsal mesoderm under control of the Wg, Dpp, and Ras signaling pathways (data not shown; Halfon et al. 2000; Carmena et al. 2002), we focused first on this element. Hbr is initially expressed broadly in the embryonic mesoderm, but at early stage 11 its mesodermal expression is dorsally restricted and undergoes dynamic modulation in subsets of dorsal cells, including the Eve-expressing muscle and heart progenitors (Fig. 2A,B; Halfon et al. 2000; Carmena et al. 2002). So as not to inadvertently delete potentially important sequences, we selected approximately 1.5 kb of intronic sequence, centered on the region identified by our search, and used it to drive reporter gene expression in transgenic embryos. This sequence element, hereafter referred to as the Hbr DME, drove expression of the reporter gene in a number of embryonic dorsal mesodermal cells (Fig. 2C; data not shown). No expression was detected in nonmesodermal Hbr-positive cells (i.e., the ectodermally derived tracheal precursors; Fig. 2C), although some expression was detected in epidermal cells overlying the reporter-expressing mesodermal cells (data not shown). Not all Hbr-positive mesodermal cells expressed the reporter gene, suggesting that mesodermal regulatory elements in addition to the DME are required to govern other aspects of mesodermal Hbr expression. Double-labeling with antibodies to Eve and to the DME-driven β-galactosidase revealed that at least one of the two Eve-positive mesodermal progenitors expressed the DME reporter (Fig. 2D). Expression in the other Eve progenitor was occasionally, although not consistently, detected (data not shown).

Ectopic mesodermal activation of the Ras/MAP kinase pathway causes an increased number of cells to express Hbr (Carmena et al. 2002). We therefore examined the activity of the Hbr DME in embryos that expressed a constitutively activated form of the Ras nuclear effector Pnt in the mesoderm. As predicted, activated Pnt induced Hbr expression in an expanded number of cells, and all of these cells expressed the DME reporter gene (Fig. 2E). A similar result was observed when the DME reporter construct was crossed into embryos mutant for the Ras-inactivated repressor, Yan (Fig. 2F). The DME is thus sufficient to drive reporter gene expression in a manner that recapitulates endogenous Hbr expression not only in wild-type conditions but also when the Hbr expression pattern has been altered by experimental manipulation of the pathways known to regulate this gene.

The Hbr DME was identified by virtue of its having predicted binding sites for a number of TFs necessary for inducing dorsal mesodermal expression, including both signal-activated factors such as Pnt and tissue-specific factors such as Twi. To determine whether these predicted sequences are functionally required for DME activity, the sites were mutated in the context of the entire DME. Mutation of the Twi binding sites led to a severe decrease in DME activity in the mesoderm, with the concomitant appearance of reporter gene expression in the Hbr-expressing cells of the developing trachea (Fig. 2G). Mutation of the Ets domain (i.e., Pnt and Yan) binding sites, which are predicted to mediate responses to Ras signaling, also caused an apparent derepression of reporter gene expression in tracheal tissue (Fig. 2H). DME-dependent mesodermal expression was not substantially affected by elimination of the Ets sites, although an occasional reduction/loss of reporter gene activity was observed (Fig. 2H). Both the Twi and Ets binding sites thus play a functional role in the DME, although unlike in the Eve MHE where the Ets binding sites are essential for enhancer activation, those in the DME appear to be required mainly for repression in ectodermally derived tissues. These data validate both our model of transcriptional regulation of dorsal mesodermal gene expression and our computational prediction of binding sites.

To further confirm that the Hbr DME is a genuine cis-regulatory element, the corresponding region from D. virilis was cloned and sequenced. A high degree of sequence conservation was observed spanning an approximately 800 bp portion of the 1.5 kb DME sequence, with complete conservation of many of the binding site motifs used in our search (Fig. 3 and Supplemental Fig. 2 available online at http://www.genome.org). The DME thus meets all of the requirements of a dorsal mesodermal enhancer element that functions in a manner similar to the Eve MHE: It drives reporter gene expression in the embryonic dorsal mesoderm in a pattern that overlaps that of the MHE; the reporter gene expression responds predictably to ectopic activation of signaling pathways; binding sites for TFs predicted to play a role in mediating DME-driven expression are functionally required; and the DME sequence has been conserved through evolution.

Additional Elements Map Close to Known and Novel Mesodermally Expressed Genes

Several additional elements in the selected subset map in proximity to genes with mesodermal expression. These include the known mesodermal genes dmef2 (DME-3), meso18e (DME-4), and rst (DME-7), as well as at least one gene, KP78b (DME-30), for which mesodermal expression has not been reported (Fig. 4; data not shown; Lilly et al. 1994; Taylor 2000; Strunkelnberg et al. 2001). We tested the ability of these elements, as well as elements DME-25 and DME-31, which map near to uncharacterized genes, to function as dorsal mesodermal enhancers in transgenic embryos. Surprisingly, none of these elements appeared able to drive reporter gene expression (data not shown). We cannot rule out the possibility that some or all of these elements are true DMEs and that our assay was insufficient to detect their activities—for instance, we may have failed to incorporate all of the necessary sequences into the reporter constructs, or the elements may be highly sensitive to specific promoter-enhancer interactions and thus fail in our construct. However, the simplest explanation is that despite their lying close to genes with the expected mesodermal expression pattern, these are, as has been seen in other computational regulatory element prediction approaches, false-positive results.

Identification of Additional Sequence Motifs

We next sought to determine whether we could identify additional sequence motifs that might represent binding sites for previously unrecognized TFs necessary for the generation of dorsal mesodermal gene expression. Even given the inclusion of false-positive elements, the presence of true DMEs among the search results should be sufficient to enable the identification of additional motifs common to a large number of the putative elements. We used the AlignACE program (Roth et al. 1998; Hughes et al. 2000) to look for such additional motifs within the selected subset of elements (see Methods).

AlignACE successfully identified four of the five motifs used to conduct the search; only the Mad motif, which is a highly degenerate sequence of low information content (see Supplemental Table 3, available online at http://www.genome.org) was not found. In addition, a large number of new motifs were identified (Supplemental Fig. 3, available online at http://www.genome.org; data not shown). In keeping with the expectation that important motifs are evolutionarily conserved, we screened the set of new motifs for their presence in phylogenetically conserved regions of either the Eve MHE or Hbr DME (see Methods). Twenty-five motifs, which cluster into 14 related groups and include the known dTcf, ETS, Tin, and Twi motifs, passed this filter. A search of the TRANSFAC database (Wingender et al. 2001; http://www.gene-regulation.com) revealed that one of these motifs, motif A, matched the binding site for the POU/homeodomain TF Oct-1 (Fig. 5A–C).

To test whether this motif represents a functional TF binding site, it was mutated in the context of an otherwise wild-type MHE. Analysis of transgenic embryos bearing this mutated MHE revealed an increased number of cells expressing the reporter gene, suggesting that this site is required for the binding of a transcriptional repressor (Fig. 5E). Although the actual repressor which binds to this site remains to be identified, preliminary results from a yeast one-hybrid screen suggest that the sequence is able to bind homeobox-containing TFs, consistent with the TRANSFAC data (S. Gisselbrecht and A. Michelson, unpubl.). In addition, available genetic data are compatible with the possibility that the Ladybird homeodomain proteins could act through this site (Jagla et al. 1997). These data show that analysis of the initial search output, even when confounded by the presence of false-positive results, can lead to the identification of previously unidentified binding sites important for dorsal mesodermal transcriptional regulation.

The MHE Model is Generalizable

The fact that we were able to identify a functionally related enhancer based on the MHE model shows that similar regulatory strategies can be used by coexpressed genes. As such, it may be possible to use variations on our model to discover enhancers for mesodermal genes subject to Wg and Ras but not Dpp control, or for nonmesodermal genes that also respond to combined Wg, Dpp, and Ras signaling. The MHE model is predicated on the presence of binding sites for the signal-responsive TFs dTcf, Mad, and Pnt, with tissue specificity provided by the mesodermal selector proteins, Twi and Tin (Halfon et al. 2000). Searching for the presence of the signal-activated factors in combination with selector proteins for other tissues may therefore be sufficient to detect signal-responsive genes in these other cell types. Furthermore, as additional DMEs are confirmed from analysis of our search results, shared structural features of the enhancers, such as a specific number, order, orientation, or spacing of binding sites may become apparent.

Interestingly, enhancers associated with a number of genes known to be coexpressed with eve and hbr in the dorsal mesoderm, such as heartless (Shishido et al. 1993) and Krüppel (Hoch et al. 1990) were not identified in the search, suggesting that at least some of these regulatory elements may depend on different mechanisms than those used by eve and hbr for their dorsal mesodermal expression. Thus, different regulatory mechanisms may lead to similar outcomes. A more complete understanding of these alternate strategies must await the isolation of regulatory elements for some of these additional genes.

Evaluating and Reducing the False-Positive Rate

Like those obtained with other approaches to the computational prediction of regulatory elements (e.g., Wasserman and Fickett 1998; Gailus-Durner et al. 2001; Krivan and Wasserman 2001; Berman et al. 2002; Markstein et al. 2002), our results include a number of false positives. These results are instructive, however; despite our analysis of the search results leading to the identification of a gene (KP78b) previously not known to have dorsal mesodermal expression, the correct cis-regulatory sequences for this gene were not correctly predicted. Indeed, four of the six demonstrated false positives mapped near genes with expression in the dorsal mesoderm. The presence of a computationally identified element near a gene with the expected pattern of expression thus does not guarantee that the element represents a true cis-regulatory module. Our data underscore the importance of empirical testing of computational predictions, which has not been performed comprehensively in other studies, and suggest that the efficacy of cis-regulatory prediction algorithms in general may be lower than initially estimated.

An important future refinement of our method will be to incorporate ways to better and more easily eliminate false-positive results. All of the elements we tested weredrawn from the selected subset of elements containing at least two of each TF binding site motif (except for the dTcf motif). Other filtering criteria may lead to a lower false-positive rate. A number of recent studies used a more general requirement that a minimum number of binding sites be present, without mandating a specific number for any particular site (Frith et al. 2001; Berman et al. 2002). However, it is not clear that adapting such an approach would significantly impact our success rate, as the site density of our apparent false positives is not significantly different from that of the true positives. Moreover, the fact that both of the elements we analyzed, the MHE and the hbr DME, contain only a single binding site for dTcf suggests that a simple minimum site number requirement such as that employed by Markstein et al. (2002) may not be widely applicable. We note as well that methods that are highly dependent on site number are highly sensitive to the information content of the given site motifs. Short or somewhat degenerate motifs will occur more frequently in the genome, meaning that elements composed of such binding sites will tend to be predicted at higher rates than those containing only longer or more invariant motifs. Use of a combinatorial strategy such as that described here that described here helps to reduce the leverage that low-information-content motifs have on the search results.

A considerable advantage will be gained by the ability to include additional types of information in the search algorithm. Both the Eve MHE and Hbr DME show significant degrees of sequence conservation in a distantly related species, and identification of conserved noncoding sequences should help in recognizing true-positive results (see also Loots et al. 2000; Wasserman et al. 2000). Extensive interspecific sequence comparisons are already possible for mammalian and nematode genomes, and a second Drosophila species will be sequenced in the near future. A second way of decreasing the number of false-positive results in future searches will be to incorporate genomic expression data from microarray, SAGE, high-throughput in situ hybridization, or similar expression profiling methods. A number of studies have already demonstrated that expression data, even when gleaned from separate experiments, provide a useful way of identifying potentially coregulated genes (Jensen and Knudsen 2000; Bussemaker et al. 2001; Pilpel et al. 2001).

Toward the Comprehensive Identification of cis-Regulatory Sequences

We were able to rapidly characterize the Hbr DME by starting with a single model and evaluating the search results with a combination of gene expression data, conserved sequence data, and alignment-mediated motif extraction. We propose that use of these several components together in an iterative searching strategy (Fig. 6), taking advantage of the rapidly growing availability of sequence, binding, and expression data in public databases, will provide an optimal approach to genomic-scale regulatory module discovery.

Fly Stocks and Reporter Gene Construction

Fly stocks and methods for reporter gene construction and mutagenesis are as described by Halfon et al. (2000). A minimum of three independent transgenic lines were analyzed for each construct. The following sequences were cloned by PCR and used for the reporter constructs: hbr/DME-2, GenBank accession no. AE003705, 113573-115127; DME-3, AE003831, 74207-74696; DME-4, AE003513, 21993-22696; DME-7, AE003426, 103046-103969; DME-25, AE003584, 51023-51921; DME-30, AE003690, 206275-207950; DME-31, AE003684, 115123-116875. For the motif A mutation, the sequence TATGCTAAT was changed to TATTATCAC.

D. virilis Cloning

Sequences from D. virilis were obtained by screening a genomic library (Blackman and Meselson 1986) with probes from the D. melanogaster Eve or Hbr coding regions (details available on request). Isolated clones were then digested and screened by Southern analysis for cross-hybridization with the MHE or DME sequences, respectively; cross-hybridizing bands were subcloned and sequenced. These sequences have been submitted to GenBank with accession nos. AF513981 (Eve MHE) and AF513982 (Hbr DME).

Sequence Alignments

Dot-plot alignments were created using MacVector (Accelrys Bioinformatics) with hash value = 1 and window = 10. Sequence alignments were performed using VISTA (Mayor et al. 2000) followed by manual adjustment for optimal alignment.

Computational Search

Binding sites for each of the five TFs were culled from the literature and used to construct position weight matrices (PWMs; Supplemental Table 3, available online at http://www.genome.org). Only sites with both in vitro binding and in vivo function were included. The program ScanACE (Hughes et al. 2000) was used to scan release 2 of the Drosophila melanogaster genome (Adams et al. 2000; http://www.fruitfly.org/sequence/download.html). Those sequences with matches that scored at a level equal to or greater than the lowest scoring functional sites in the Eve MHE were selected. We then ran a program (cooccur_scan.pl; available at http://arep.med.harvard.edu/Halfon_Grad_etal/supplemental.html) to determine regions of cooccurrence of predicted binding sites for all five TFs within at least 500 bp. Two special cases were addressed: (1) Both Twi and MAD sites can form palindromes; when this was observed, only one of the palindromes was included in further analysis. (2) Mad sites may match the initial bases of Ets sites; in these cases, the Ets sites were included to the exclusion of the overlapping Mad sites, as the Mad PWM has lower information content than the Ets PWM. The algorithm extends the windows for as long as the condition that binding sites for all five TFs appear within 500 bp holds. Using Drosophila genome annotations (http:// www.fruitfly.org/sequence/download.html), the program also assigns locations of cooccurrence regions with respect to known and predicted gene coding sequence. The results of this analysis were termed “elements.” A subset of the elements was then sorted to identify those in noncoding sequence with at least two predicted sites for each of the TFs except dTcf (termed the “selected subset”).

Monte Carlo Methods

To estimate the number of cooccurrences expected by chance given the number of motifs identified by ScanACE, a variant of the cooccur_scan.pl program was used. This program maintained the number of motifs found in the genome scan, randomized the locations of each motif, and then searched for cooccurrences as in the coccur_scan.pl program. As the D. melanogaster genome sequence is available in several FASTA format files, each representing a chromosomal arm or chromosome, the scans were performed per individual FASTA file and then summed to represent the whole genome. Gaps in the sequences denoted by strings of Ns greater than 50 letters were taken into account in location randomizations to ensure that the overall space for motif locations was nearly identical.

Motif Discovery and Evaluation

Sequence corresponding to the 34 elements in the selected subset was extracted from Release 2 of the Drosophila genome and searched for overrepresented subsequences using AlignACE (Roth et al. 1998; Hughes et al. 2000). AlignACE identified a total of 755 motifs from several runs using combinations of parameters intended to sample extensively the sequence space (all variations of –gcback 0.42, 0.45, 0.48 and –numcols 7, 10, 13). The ACE package clustering program revealed that the motifs cluster into 375 motif groups at the 0.8 correlation coefficient level.

The top 100 sites in the 34 elements for each of the 755 motifs were determined by running ScanACE with default parameters and –s 100 on the set of predicted element sequences. Using a Perl script, these sites were evaluated for phylogenetic conservation as determined by D. melanogaster versus D. virilis sequence alignments of the Eve MHE and Hbr DME. Twenty-five motifs mapped to perfectly conserved sites. Clustering using the ACE package clustering software at the 0.8 correlation coefficient level organized the 25 motifs into 14 motif groups. All of the 755 initial motifs were compared against TRANSFAC (Wingender et al. 2001) using CompareACE (Hughes et al. 2000). The ACE programs are available for download at http://arep.med.harvard.edu.