The transformation of a developing epithelium into an adult structure is a complex process, which often involves coordinated changes in cell proliferation, metabolism, adhesion, and shape. To identify genetic mechanisms that control epithelial differentiation, we analyzed the temporal patterns of gene expression during metamorphosis of the Drosophila wing.
We found that a striking number of genes, approximately 50% of the Drosophila transcriptome, exhibited changes in expression during a time course of wing development. While cis-acting enhancer sequences clearly correlated with these changes, a stronger correlation was discovered between core-promoter types and the dynamic patterns of gene expression within this differentiating tissue. In support of the hypothesis that core-promoter type influences the dynamics of expression, expression levels of several TATA-box binding protein associated factors (TAFs) and other core promoter-associated components changed during this developmental time course, and a testes-specific TAF (tTAF) played a critical role in timing cellular differentiation within the wing.
Our results suggest that the combinatorial control of gene expression via cis-acting enhancer sequences and core-promoter types, determine the complex changes in gene expression that drive morphogenesis and terminal differentiation of the Drosophila wing epithelium.
Within many developmental contexts, cells assemble into epithelial sheets and coordinately differentiate to form organs and tissues. This differentiation process often involves changes in cell shape, cell-cycle arrest, and the emergence of cell-type-specific features. Whereas many factors (i.e., transcription factors, signaling pathways, adhesion molecules, and cytoskeletal elements) have been identified that control individual aspects of this differentiative process, it is less clear how these divergent changes in epithelial and cell biology are temporally coordinated to generate the adult structure.
Drosophila wing tissue is uniquely amenable to studies that concern epithelial differentiation. During larval stages of development, presumptive wing cells proliferate within a relatively uniform, mono-layered epithelium. When the animal pupariates and metamorphosis is initiated, however, this simple epithelium rapidly transforms into a complex wing structure. The most dramatic changes associated with wing differentiation involve cell shape, as evagination and elongation of the wing pouch gives rise to the much larger wing blade. This process requires wing-cell flattening, creation of a bi-layered epithelium, and the formation of tubes (i.e., wing veins)
[1, 2]. Coincident with this epithelial morphogenesis, each cell executes its final cell cycle
[3, 4], adopts a cell-type-specific shape
[5, 6], and forms a polarized actin-rich wing hair
, among other things. To identify regulatory mechanisms that initiate and coordinate the terminal differentiation of an epithelium, therefore, we have characterized the changes in gene expression that occur during a time course of Drosophila wing metamorphosis.
Our time-course study in the wing reveals a surprisingly complex pattern of gene regulation during metamorphosis, affecting nearly half of the genes in the Drosophila genome at one or more time points. We suggest that this complex regulatory pattern results from the combined influence of sequence-specific binding complexes within genetic enhancers, and temporal changes in core-promoter preference. As such, our focused analysis of a single epithelial tissue as it undergoes a dynamic developmental transition, suggests an important role for core-promoter sequences (and associated proteins) in coordinating epithelial differentiation.
Within a very short period of time (roughly 36 h at 25°C) the wing imaginal disc is transformed from a relatively simple epithelial sheet into a complex structure that resembles the adult wing (Figure
1A). To illustrate this fact, and to confirm the developmental timing of our genotype under study (w1118), wing tissue was dissected at several developmental time points between the late third larval instar (L3) and 36 h after puparium formation (APF). We found that wing pouch evagination began at pupariation and was clearly evident by 2 h APF (Figure
1A,C). By 6 h APF, wing elongation and hinge region constriction were occurring (Figure
1C), as described previously
. In addition to gross morphological changes such as these, dramatic changes in shape were also observed at the cellular level. Developing wing tissue was stained for DE-cadherin to reveal the apical shape of each cell
. Between 18 and 36 h APF, wing cells adopted a hexagonal shape (a process termed hexagonal packing)
, and vein/intervein shape differences emerged
1B). Each wing cell also formed a wing hair during early pupal stages
, as revealed by F-actin localization (Figure
Between late L3 and 36 h APF, presumptive wing cells also exit the cell cycle. To precisely determine the timing of this event, we used fluorescence-activated cell sorting (FACS) analysis to characterize the temporal changes in DNA content of wing cells (Figure
1C). The wing epithelium remained proliferative through 2 h APF (i.e., high levels of G1-, G2-, and S-phase cells), but by 6 h APF many cells had temporarily arrested cycling with a G2 DNA content. This results in an in vivo semi-synchrony of the final cell cycle. By 18 h APF, most cells in the wing had re-entered the cell cycle and completed their final round of division (i.e., many cells were already in G1, with smaller proportions completing their final S- and G2-phases)
[3, 4]. By 24 h APF, over 95% of the wing was arrested in G1, and only very rare S- or M-phases were observed in the wing blade (see Additional file
1). No further cell cycling was observed in the wing blade between 24 and 36 h APF.
To capture the gene expression dynamics that drive these developmental changes, we performed a microarray analysis. RNA was isolated from late L3 wing discs (the reference sample), and from wings at 2, 6, 18, 24, and 36 h APF. For technical reasons concerning the dissection, samples from the earliest time points (L3, 2 h, and 6 h APF) contained presumptive notum tissue, whereas the later samples (18, 24, and 36 h APF) did not. This array analysis revealed a striking number of transcripts whose level of expression significantly changed during this developmental time course (Figure
1D). For example, comparison between L3 and 36 h APF yielded 6855 genes that were differentially expressed. In total, 8338 genes were differentially expressed at one or more time points.
Gene expression dynamics during differentiation of the wing
With such a large number of differentially expressed genes, we explored computational methods of categorizing and grouping genes to identify patterns within the data. We first used self-organized mapping of genes (ordered by similarity across the five time points) to intuitively display global gene-expression dynamics. This provided a clear view of general changes in gene expression over time, but masked many discrete temporal patterns by grouping them in with the most common trends (see Additional file
We next used hierarchical clustering to organize all 8338 differentially expressed genes. From this analysis, we estimated that the dataset contained between 25 and 30 distinct temporal patterns of gene expression. Based on this information, we used the Genesis program to perform k-means clustering, imposing values of k that ranged from 26–40. We empirically determined that 30 clusters represented the optimal number of temporal expression patterns, without generating too many similar or overlapping clusters. Additional file
3 shows the pattern of expression for each k-means cluster and Supplemental text files 1–30 contain the transcript list and complete data for each cluster. We next hierarchically organized the genes within each k-means cluster, and ordered the clusters to generate a heat map, which represents the global changes in gene expression that occur during metamorphosis of the wing (Figure
Enhancer motifs are enriched in co-regulated clusters
We first examined gene ontology (GO) enrichment within the 30 k-means clusters to see if clusters with similar temporal patterns of expression contained functionally related genes. For this analysis 118 GO terms were used. Sorting clusters based on GO enrichment, however, did not group clusters with similar patterns of expression (see Additional file
4). During differentiation of the wing, therefore, GO terms do not accurately predict the temporal pattern of expression for a given gene.
We next asked whether similarly regulated clusters of genes share known regulatory motifs (i.e., transcription factor binding sites) within upstream enhancer regions. A database of 87 known Drosophila transcription factor-binding motifs (based upon the FLYREG motif set) was used for this analysis
. We also added several motifs to this dataset, including the mammalian E2F motif from Transfac, and custom motifs for Drosophila Myc, Ecdysone Receptor/Ultraspiracle complex (EcR/USP), Trithorax-like (Trl), Medea (Med), and Drosophila E2F. These custom motifs were obtained by performing MEME analysis on previously described target genes
[11–13], and by deriving a consensus from two Drosophila E2F binding sites identified by Yamaguchi and colleagues
. The Med motif was obtained by MEME analysis of confirmed direct Decapentaplegic (Dpp) signaling target genes in the wing (DOK and LAB, manuscript in preparation). Finally, a starvation responsive E-box motif
 was also included. See Additional file
5 for details concerning these custom motifs. Additional file
6 provides the position-specific scoring matrices (PSSM) for all motifs used for this analysis.
For each gene, 1 kb of sequence immediately upstream of the transcriptional start site (TSS) was defined as the DNA-search space. This amount of sequence provided a good compromise between capturing important regulatory elements and minimizing sequence search space, which can decrease enrichment. We empirically tested search spaces of 1, 2, and 3 kb upstream of each TSS, and found the greatest amount of motif enrichment and/or depletion when 1 kb of sequence was analyzed. This finding is consistent with published data concerning Drosophila regulatory elements
[16, 17], but certainly does not identify all binding sites within target-gene enhancers and will omit binding sites outside of the 1-kb window and in introns.
3, the average change in gene expression for each cluster is represented by the heat map at the top, and the most significantly enriched GO term for each cluster is listed (if p < 10-5). Clusters were ordered based on similarities between the temporal patterns of gene expression. We found that a number of motifs known to regulate cell growth, protein synthesis, and the cell cycle (e.g., E2F, Myc, DNA-Replication-related Element Factor (Dref), Mothers against Dpp (Mad), Med, and Brinker binding motifs
[18–24]) were frequently enriched in the upstream sequences of the same clusters. In addition, expression of these clusters generally declined during pupal stages of development and upon cell-cycle exit (Figure
3C). An interesting exception to this rule was cluster 29, which contains many of these same binding sites, but is characterized by increased levels of expression during pupal stages (further discussed below). We also found that motifs associated with tissue differentiation (e.g., Caudal, Crocodile, Engrailed, and Ultrabithorax binding motifs
[25–32]) were significantly depleted in clusters with declining levels of expression (Figure
3C). Instead, these differentiation-associated motifs were enriched in clusters that increased expression during pupal time points. Finally, clusters that exhibited a strong peak of gene expression (log2 > |3|) at 6 h APF (or a dip in the case of cluster 28) were typically associated with Broad binding-site enrichment. This suggests that Broad mediates a transcriptional response to Ecdysone signaling at 6 h APF in the wing (Figure
3C). This likely mediates the G2 arrest that characterizes the 6 h APF wing, as two high-scoring potential Broad binding sites are found upstream of the gene string, whose downregulation mediates the temporary G2 arrest. Our data suggests that in the differentiating wing enhancer motifs are generally enriched in clusters of genes with similar temporal patterns of expression.
Core-promoter types strongly correlate with co-regulated clusters of genes
A previous study used cell sorting coupled with serial analysis of gene expression (SAGE) to reveal a role for Dref in controlling gene-expression changes during differentiation of the Drosophila eye
. Dref is a transcription factor that acts as a component of a transcription-initiation complex containing TBP-Related Factor2 (TRF2)
. The Dref regulatory element (DRE) is found in upstream regulatory regions of many cell-cycle genes (e.g., pcna), but is also one of the most common core-promoter motifs identified in Drosophila. We found DRE enrichment in clusters of genes with declining levels of expression during terminal differentiation, which generally contained genes involved in the cell cycle and ribosomal biogenesis. This finding, combined with recent work demonstrating a role for core-promoter components in regulating gene transcription during development
, prompted us to examine whether core-promoter types are associated with specific temporal patterns of gene expression in the differentiating wing. For each of the 30 k-means clusters, therefore, we asked whether ten core Drosophila promoter motifs (see Materials and Methods) were either enriched or depleted. The 100 bp spanning the TSS (from −60 to +40 bp) was defined as the DNA-search space for this analysis. Certain core-promoter (CP) motifs (i.e., DRE, CP-1, CP-7, and CP-6) were frequently found together in clusters of genes with declining levels of expression during wing differentiation. In contrast, Initiator (INR), Downstream Promoter Element (DPE), and Motif-Ten Element (MTE) motifs were often found together in clusters that increased expression over time. TATA motifs were rarely enriched or depleted in these k-means clusters (Figure
3B). In the developing wing, therefore, core-promoter type strongly correlates with the temporal patterns of gene expression, suggesting a change in core promoter usage upon terminal differentiation.
Clusters enriched for certain motifs contain verified target genes
To ask whether motif enrichment within a cluster indicates that these genes are regulated by a common factor (i.e., the transcription factor that binds the motif), we compared selected clusters to published datasets of verified target genes. For example, cluster 30 is enriched in EcR/USP motifs. To validate that this cluster represents a true EcR/USP regulon, we compared genes within this cluster to previously identified EcR/USP targets
4A). Similarly, clusters 13 and 29 are enriched in the starvation-associated E-box motif, the Myc motif, and the Drosophila E2F motif. Genes from these two clusters were therefore compared to starvation-response genes
, Myc target genes
, and E2F target genes
4B-D). In every case except one, significant enrichment for the independently validated target genes was observed (as measured by hypergeometic probability < 0.05). Observed enrichment was most significant when similar tissues and developmental stages were compared. This likely reflects tissue- and temporal-specific target-gene expression, as expected. Strikingly, genes within clusters 13 and 29 are enriched for both Myc and E2F motifs, yet they are characterized by different core promoter motifs, and are regulated in an opposite fashion (i.e., expression of cluster 13 and 29 decline and rise during wing differentiation, respectively). Nonetheless, the majority of genes in cluster 13 and 29 were independently validated as legitimate Myc and E2F target genes. However, nearly all cluster 13 Myc/E2F targets increase in response to Myc or E2F overexpression (Figure
4C,D, blue in heat maps), whereas nearly all cluster 29 genes decrease in response to Myc or E2F overexpression (Figure
4C,D, yellow in heat maps). Our data suggest that this dramatic switch in gene-expression behavior may depend upon the type of core promoter, although we cannot exclude the contribution of other enhancer motifs to gene expression in these clusters.
Core-promoter type may influence gene-expression behavior
In our analysis, each transcript has a unique identifier, providing the transcription start location (TSS) for each isoform. Thus, all transcripts were treated as separate entities, allowing us to include alternate TSSs in our analysis. This raised the interesting question of whether genes with multiple TSSs that contain different core promoter types, switch core-promoter usage during differentiation.
In the Drosophila genome, 14% of genes are predicted to utilize multiple TSSs. Based on our analysis, a small number, 10.5% of these have identifiable DPE or DRE core promoter motifs. Of this small number of genes, there are 20 genes within the Drosophila genome that have a DRE at one promoter and a DPE (or another promoter motif) at an alternative TSS. In our wing dataset, nine of these genes had at least two temporally regulated transcripts that use fundamentally different core promoters. In all nine cases, the complement of promoter motifs predicted the transcript’s overall expression behavior and cluster assignment (Table
1). Transcripts from the same gene that utilized different core-promoter types did not cluster together, whereas transcripts from different genes with similar promoters frequently clustered together (e.g., clusters 3 and 23; see Table
1). These nine genes represent 22 transcripts. For 20 of the 22 alternate TSS transcripts, the predicted core promoter enrichment for the assigned cluster also correctly predicted the transcript expression pattern. For these genes, core-promoter sequences could exert a greater influence on temporal patterns of expression than more distal regulatory sequences.
Cluster assignments for genes with alternate TSSs
A small number of genes in our data set employ multiple TSSs that are associated with different types of core promoters. Transcripts associated with these genes often exhibited different temporal patterns of expression (i.e., were assigned to different clusters), which was predicted by the type of promoter, in all except for two transcripts (indicated by an asterisk).
Changes in expression of core promoter-associated factors upon terminal differentiation
We next examined the patterns of expression for genes that encode core promoter-binding proteins, general transcription components, and TBP-associated factors (TAFs). Expression data from the wing developmental time course was extracted for these genes, and hierarchical clustering was used to group them according to their temporal pattern of gene expression (Figure
5A, Additional file
7). Significant changes in expression were observed for core promoter element binding proteins (e.g., Dref and Boundary Element-Associated Factor (BEAF)), and TFIID components (e.g., TAFs and testes-specific TAFs (tTAFs)). Expression levels of Dref, several TAFs (2, 5, 7, 8, and 13), and an isoform of BEAF (transcript RA) that binds the DRE
 all declined during terminal differentiation of the wing. In contrast, expression of TRF2 (a Dref partner
), an isoform of BEAF (transcript RB), several TAFs (3, 11, 12, and 10b), and select testes-specific TAFS (males in absentia (mia), cannonball (can), andnohitter (nht)) increased during the developmental time course. These observations support the hypothesis that a change in core promoter usage may be associated with terminal differentiation in the wing.
To verify the expression-array data, western-blot analysis was used to examine protein levels for selected factors in pupal-wing tissue (Figure
5B). By microarray analysis, the tTAF spermatocyte arrest (sa) did not increase expression at the time points shown. This was verified by western blot, and demonstrates selectivity in tTAF regulation within the differentiating wing (i.e., not all tTAFs were co-regulated). Levels of the tTAF Mia increased in the wing at later stages of development (both by RT-PCR and western-blot analyses), whereas total levels of BEAF decreased in the wing over time.
Using inducible RNAi transgenes, we then reduced levels of several TAFs and testes-specific TAFs (tTAFs) and examined the effect on differentiation of the wing. Because many of these genes can affect general transcription, we used the temperature-sensitive Gal80 system, in combination with wing Gal4 drivers, to express RNAi constructs from the mid-larval stage of development (L2) until eclosion (i.e., the adult stage). RNAis to specific TAFs and core promoter-associated genes resulted in disrupted wing phenotypes when expressed via apterous-Gal4. These included bip3 (TAF3) and TFIIAS-2 (Additional file
8). Employing the same experimental protocol for the tTAF can resulted in lethality. When the tTAF nht was inhibited in the posterior wing using engrailed-Gal4 (L2 to eclosion), defects in posterior wing growth and cuticle integrity were observed (Figure
5D). Terminal differentiation in nht RNAi-expressing wings was assessed at the cellular level by staining these wings for F-actin at 34 h APF. Beginning at this stage of development, every cell in the wing blade forms an actin-rich wing hair, which is an early differentiation event. Expression of nht RNAi in the posterior wing led to a cell-autonomous delay in wing-hair formation and terminal differentiation (Figure
5F), revealing an unexpected potential somatic role for nht.
In further support of a function for nht in the wing, we noted that flies hemizygous for the EMS loss of function alleles nhtz5946 or nhtz5947 were sub-viable and 50% of escapers were developmentally delayed by 24 hours during the pupal stage. Nearly all nhtz5946 hemizygous adults exhibited ectopic wing vein phenotypes (Figure
5G) and 10% exhibited multiple wing hair phenotypes (Figure
5G, inset). Fifty percent of nhtz5347/nhtz5946 transheterozygous males exhibited a mild ectopic vein phenotype in the anterior wing and 10% of females raised at 25°C exhibited ectopic vein (Figure
5H-I) and occasional patches of small, thin wing-hairs (Figure
Vein refinement involves the coordinated action of multiple signaling pathways and a developmental delay can affect the vein the refinement and limitation process. Other genotypes exhibiting wing hair developmental delay also exhibit multiple wing hairs or thin wing hair phenotypes
, further suggesting a developmental delay in the wing occurs due to loss of nht. This result, combined with the promoter/motif analysis, suggests that a change in TAFs and core-promoter type preference may play an important role in regulating gene expression during differentiation of the wing. In total, our results indicate that the combinatorial influence of enhancer occupancy and core-promoter usage shapes the spatial and temporal patterns of gene expression that drive epithelial differentiation.
Dramatic gene-expression changes during metamorphosis of the wing
Previous studies with whole animals estimated that 10% of the genome is differentially regulated during metamorphosis
[44–46]. Many tissue-specific changes in gene expression are likely obscured by using whole animals however, as Li and White found that expression levels for 17% of the genes in the genome are significantly affected in at least one of five tissues during the developmental transition into metamorphosis
. Strikingly, our analysis of wing tissue during metamorphosis revealed significant changes in expression for ~50% of the genome (8338 genes). This likely reflects the wide range of developmental and cellular processes that are involved in differentiation of the wing (e.g., cell adhesion, cytoskeletal remodeling, apoptosis, neurogenesis, and cuticle formation). Changes to the cell cycle are particularly dynamic during this developmental time period (late L3 to 36 h APF), as initially proliferative cells temporary arrest in G2, synchronously execute a final cell cycle, and permanently exit the cell cycle in G1 (Figure
1). Our results indicate that such a breadth of developmental changes ultimately depends upon the modulation of several thousand genes. Our data suggests that this dynamic program of gene-expression regulation involves combined inputs from both enhancer sequences and the type of core promoter. Furthermore, this suggests that only a small proportion of genes will truly fit the definition of “constitutive” or “housekeeping” when examining dynamic developmental processes across multiple tissues.
Temporal gene-expression patterns correlate most strongly with core-promoter types
In the differentiating wing, gene expression patterns were not highly correlated with GO term enrichment, i.e., functionally related genes were not co-regulated. Instead, shared upstream-regulatory motifs and the type of core promoter were more accurate predictors of a gene’s temporal pattern of expression. These findings are in contrast to a recent analysis, which suggests that co-regulated genes in Drosophila are more likely to share GO terms than transcription factor binding sites
. These authors based their conclusions on gene-expression data from whole organisms, an approach that may have masked tissue-specific patterns that are necessary to find shared transcription factor binding sites. Instead our analysis is consistent with the findings of FitzGerald et al.,
 which found core promoter type strongly correlated with gene expression in embryos vs. adults. We also found that enrichment for a particular type of core promoter more accurately predicted a gene’s temporal pattern of expression than enrichment for any particular enhancer motif within 1 kb of the TSS. Comparing the same tissue at multiple developmental time points, therefore, has revealed an unappreciated level of transcriptional regulation that depends upon the core-promoter type. As the wing differentiates, we found a shift in expression away from genes with DRE, and certain core promoter types (CP-1, -6, and −7), toward genes that use INR, DPE, and MTE promoter elements. This trend even held true for genes with multiple TSSs containing different core promoters, where in nearly every case we observed a shift in TSS usage during differentiation (Table
1). This shift may reflect changes in TAF and core promoter associated factors during terminal differentiation, and is consistent with studies that have found wing differentiation functions for additional components of the basal transcription machinery (summarized in Additional file
The type of core promoter may dictate target-gene induction or repression in response to Myc and E2F
For general enhancer analysis, promoters such as the heat-shock protein (HSP) 70 TATA promoter are thought to be sufficient. Validation of the reporter-gene expression pattern is typically performed at a limited number of pre-selected time points, and reporters that recapitulate the expected pattern are chosen. Our data suggest however, that such constructs may miss important aspects of temporal regulation that are not captured by an arbitrary choice of core promoter.
Two gene clusters in our dataset (cluster 13 and 29) were enriched for both E2F and Myc binding sites, and regulated known E2F and/or Myc target genes
[15, 40]. However, most genes within cluster 13 were induced by dMyc/dE2F1, whereas most genes within cluster 29 genes were repressed by dMyc/dE2F1 activity (Figure
4). The latter represents a non-canonical response, as these factors primarily activate transcription when overexpressed. We identified a difference in core-promoter sequences between these two clusters, as cluster-13 genes were enriched for DRE, CP-1, -6, and −7 sequences, whereas cluster-29 genes were enriched for DPE, MTE, and INR sequences. This finding, together with our observation that the expression of genes encoding for core promoter-associated proteins were temporally regulated in the differentiating wing, suggests that differences in core-promoter type may underlie the non-canonical behavior of these E2F and Myc target genes. Most studies concerning E2F regulation use either well-characterized endogenous DRE-containing promoters
, cell-cycle gene promoters that contain DREs
, or engineered constructs with an HSP70 TATA-type promoter
. Our data suggest that there may be another level of E2F (and possibly Myc) transcriptional regulation that is missed by using reporters such as these. This idea is supported by findings that an enhancer can be functionally linked to a specific core promoter in Drosophila, and that certain factors can stimulate DPE-dependent transcription and inhibit TATA-dependent transcription
. Future studies that involve reporter constructs with different types of core promoters (like those generated in
) will be necessary to verify and characterize this new level of regulation.
A requirement for a tTAF in wing terminal differentiation
A previous report indicated that expression of nht (a tTAF) is restricted to the male germline
. It was surprising to find therefore, that nht may play a role in the proper timing of terminal differentiation in the pupal wing. However, nht expression had only been examined in embryos (0–24 h after egg deposition (AED)) and adults
. Newer modEncode transcriptome data found that nht transcription is undetectable until pupal stages of development
, whereas the adult wing contains very few viable cells, due to developmentally regulated apoptosis that occurs immediately after eclosion
. We found that nht, as well as several other tTAFs including mia and sa, are expressed in the pupal wing. These results are based on microarray, RT-PCR and western blotting analyses, where possible (Figure
5). Our tTAF expression data do not contradict previous findings, therefore, but instead reveal somatic expression at a different stage of development.
nht mutants are viable but male sterile, consistent with their well-defined role in terminal differentiation of spermatocytes
. Although wing phenotypes associated with nht loss-of-function have not been reported, we used an inducible RNAi construct to inhibit nht function in the pupal wing. This acute knockdown approach revealed a dramatic wing hair differentiation phenotype that may be partially masked in nht mutants by compensatory changes in other TAFs. We also observed wing vein and abnormal wing-hair phenotypes in nht mutants (Figure
5). This genetic data demonstrates an unexpected role for nht in wing development. Our results suggest that non-canonical TAFs may play an important role in terminal differentiation of somatic tissues as well
[62, 63], including TAFs that were previously thought to be germline specific.
The dynamic changes in gene expression during Drosophila wing terminal differentiation encompass approximately 50% of the protein coding genes in the genome, and are directed by combinatorial inputs from both enhancer sequences and the core-promoter type at individual genes. Our results suggest that a change in core promoter preference, likely mediated by a shift in expression of core promoter binding proteins, plays a much more significant role in modulating gene expression during metamorphosis than previously recognized.
For pupal wing microarray hybridizations, 10 pupal wings were dissected from w1118 animals (Bloomington Stock #3605) that had been raised on standard cornmeal/molasses media under uncrowded conditions (less than 50 larvae per vial) at 25°C. Wandering L3 was defined as wandering larvae at the gut half-empty stage (approximately 110 h AED). All pupae were staged from 0 h APF, defined as stiff white pre-pupae (approximately 120 h AED ± 20 min at 25°C). RNA was isolated from dissected tissues using Trizol (Invitrogen, Grand Island, NY), and cleaned using the RNAEasy Kit (Qiagen, Valencia, CA). RNA integrity was confirmed via agarose gel. Using 500–1,000 ng input RNA for each reaction, cDNA synthesis was performed with one subsequent round of T7-dependent linear RNA amplification, using the commercially available Message AmpTM kit from Ambion
. Amplified RNA was quantified via nanodrop, and its integrity confirmed via Bioanalyzer (Agilent, Santa Clara, CA). According to Nimblegen protocols (Madison, WI), 10 mg of RNA was subsequently labeled in a cDNA synthesis reaction and hybridized to Nimblegen 4-plex 60-mer Drosophila expression arrays (
http://www.nimblegen.com). Hybridizations were repeated four times with independently obtained biological replicates to ensure maximal confidence in data reproducibility. NimbleScan software was used for array scanning and quantile normalization
. Gene calls were generated using the Robust Multichip Average (RMA) algorithm
. All arrays in this study were normalized together. MA plots of the array data, showing a linear relationship between intensity and average intensity post-normalization, is provided (see Additional file
10). Statistically significant changes in gene expression (adjusted p < 0.05) were determined using ANOVA. The mean fold change for genes with significant changes was 2.2. Approximately half of the genes exhibited fold changes > 2, and half exhibited less than a 2-fold change. We also compared our results to previously published RT-PCR and microarray data where similar pupal wing stages were possible (Additional file
The data discussed in this publication have been deposited in NCBI’s Gene Expression Omnibus. The GEO Series accession number is GSE36015.
The following lines were used in this study: w1118 (Bloomington Stock, #3605), nhtz5347/nhtz5946 , Df(2 L)A263 (Bloomington Stocks #25159, 6062 and 25160), engrailed-Gal4, UAS-GFP; tub-Gal80TS (engrailed-Gal4, UAS-GFP from
, tub-Gal80TS from
), UAS-nhtRNAi (Vienna Drosophila RNAi Collection, stock #10351). For experiments using the temperature-sensitive (TS) Gal80 system, animals were shifted from 18°C to 28°C at the indicated times.
Gene expression clusters
Genes differentially expressed at one or more time points compared to the wandering L3 reference sample (> 1.5-fold) were grouped into 30 clusters using the k-means clustering algorithm with Euclidean distance in the Genesis program (
http://genome.tugraz.at). Hierarchical clustering was also performed in Genesis to generate heat maps showing transcript changes (log2 ratio of expression compared to L3 reference) for all transcripts with a fold-change ≥ 1.3 (> log2 ± 0.4) at one or more time points. For the gene expression heat map of 30 k-means clusters shown in Figure
3, the average expression for all genes within the indicated cluster at the indicated time point is shown.
Western blotting, RT-PCR, and histology
Western blots using 10 wild-type Drosophila pupal wings at the indicated stages were performed as described
. Blots were incubated with primary antibodies directed against acetylated histone H4 (1:2000; Millipore, Billerica, MA), BEAF (1:200; Developmental Studies Hybridoma Bank (DSHB)), Sa (1:100; kindly provided by Dr. M. Fuller), Mia (1:100; a kind gift from Dr. D. Wassarman). Appropriate horse-radish peroxidase-conjugated secondary antibodies were used at 1:2000 and detected using enhanced chemiluminescence reagents (GE Healthcare, Piscataway, NJ). For RT-PCRs, 10 dissected wings of the indicated stages were homogenized in Trizol for RNA isolation and clean up using the RNAeasy kit (Qiagen) with DNAse treatment to remove genomic DNA. RNA was quantified using a nanodrop, integrity was confirmed using gel electrophoresis, and 500 ng was used for cDNA synthesis via the oligo-dT primed, Superscript III kit (Invitrogen). PCR was carried out using 1/10 of the cDNA reaction using Go Taq (Promega, Madison, WI), with 35 amplification cycles. Primers for mia are available upon request. Mock reactions without RT (−RT in Figure
5) were carried out to ensure PCR products were not from genomic DNA contamination.
Pupal wings of the indicated stages were fixed with 4% paraformaldehyde in phosphate-buffered saline (PBS) and stained using a primary antibody directed against DE-cadherin (DSHB) and an Alexa-fluor 568 secondary antibody as described
. Pupal wings were stained for nuclei using Hoechst 33258 at 0.5 μg/mL in PBS, or F-actin using Rhodamine-labeled phalloidin diluted 1:200 in PBS (Molecular Probes). For F-actin staining, wings were not exposed to detergent. Confocal images of pupal wings were taken using a Zeiss LSM510 confocal microscope at 20× and Zeiss AIM software.
Adult wings were dehydrated in 100% ethanol, placed in methyl salicylate (Sigma, St. Louis, MO) for 10 min, and permanently mounted and flattened in Canada balsam and methyl salicylate (1:1) (Sigma). Adult wings were photographed under brightfield conditions on a Leitz Orthoplan2 at 10-20× magnification, using a Nikon DS-Vi1 color camera and Nikon NIS Elements software.
Twenty wild-type Drosophila wings of the indicated stages were dissected, dissociated in a trypsin/PBS solution and stained for DNA using Hoechst 33342 as described
. Dissociated cells were analyzed using a FACS Vantage Cytometer (BD) with CellQuest software. At least 20,000 cells were measured to generate cell-cycle histograms.
 was used (at p < 0.0001) to identify the locations of relevant DNA-binding motifs (Additional file
6 includes all PSSMs used. Additional file
5 shows custom motifs added to the FLYREG dataset). Core promoter (CP) motif PSSMs have been described
. For all genes, Fimo scans were performed independently at the core promoter (−60 to +40 bp), and across the extended promoter region (−1 kb to −1 bp) of all genes.
To identify motifs that were enriched or depleted in given gene clusters, the sum of all motif occurrences for each cluster was calculated. Permutation tests were then performed to determine the significance of seeing that many motif occurrences at random in a cluster with ‘n’ genes where ‘n’ equals the number of genes for each cluster. For each permutation test, ‘n’ genes were randomly chosen to form a random cluster. The sum of each given motif’s occurrence in a randomly selected set was recorded. One thousand such permutations were performed to calculate the mean and standard deviation of motif occurrences in randomly selected clusters of ‘n’ genes. From the distribution of motif occurrences, a z-score was calculated to express the enrichment or depletion of a given motif in a real cluster of genes compared with randomly-chosen clusters of equal size. A threshold z-score of |3| was chosen as significant for enrichment or depletion of motifs. Permutation tests were performed to calculate z-scores for each cluster and motif combination.
Gene ontology analysis
Gene ontology terms were obtained from FlyBase (
http://www.flybase.org). For each cluster of genes, the enrichment of functional terms associated with that cluster was calculated using the hypergeometric distribution, which calculates the probability of randomly drawing ‘b’ genes with a given functional term from among all ‘N’ genes in the genome, given that there were ‘n’ genes in the cluster (the number of draws from the genome) and ‘B’ total genes with that annotation in the entire genome.
Validation of cluster motif enrichment
Hypergeometric probabilities were calculated using the hypergeometric distribution. In our case: (1) the population size equals the total number of genes on our microarray platform, (2) the number of successes in the population is the total number of genes that are represented in both our dataset and the published dataset under comparison, (3) the sample size is the total number of genes in each cluster (for which there is data in the published dataset under comparison), and (4) the number of successes in the sample equals the total number of genes in that cluster that are validated as significantly changed by the treatment or transcription factor under study in the published dataset under comparison. A p-value less than 0.05 indicated a significant overlap.
Analysis of variance
After puparium formation
DNA replication-related element
Fluorescence activated cell sorting
Find individual motif occurrences
Larval instar stage 3
Multiple em for motif elicitation
Position specific scoring matrix
Serial analysis of gene expression
TATA binding protein
Transcriptional start site
We thank M. Fuller and D. Wassarman for reagents; L. Li for essential advice and analysis of preliminary microarray data; M. Morgan, J. Davison and the FHCRC array facility for help with microarray hybridization, processing, and data analysis. Antibodies were obtained from the Developmental Studies Hybridoma Bank, developed under the auspices of the NICHD and maintained by The University of Iowa, Department of Biology, Iowa City, IA 52242. The indicated Drosophila lines were obtained from the Bloomington Drosophila Stock Center at Indiana University or the Vienna Drosophila RNAi Collection (VDRC). This work was supported by grants from the NIH, GM070887 to BAE and GM086517 to LAB. DDO was supported by NCI 1 U54 CA132381 to the FHCRC.
Division of Basic Sciences, Fred Hutchinson Cancer Research Center
Gladstone Institute for Cardiovascular Disease, J David Gladstone Institutes
German Cancer Research Center (DKFZ)-Center for Molecular Biology Heidelberg (ZMBH) Alliance
Molecular Cellular and Developmental Biology, University of Michigan
Fristrom D, Wilcox M, Fristrom J: The distribution of PS integrins, laminin A and F-actin during key stages in Drosophila wing development.Development 1993,117(2):509–523.PubMed
Fristrom D, Gotwals P, Eaton S, Kornberg TB, Sturtevant M, Bier E, Fristrom JW: Blistered: a gene required for vein/intervein formation in wings of Drosophila.Development 1994,120(9):2661–2671.PubMed
Schubiger M, Palka J: Changing spatial patterns of DNA replication in the developing wing of Drosophila.Dev Biol 1987,123(1):145–153.PubMedView Article
Milan M, Campuzano S, Garcia-Bellido A: Cell cycling and patterned cell proliferation in the Drosophila wing during metamorphosis.Proc Natl Acad Sci U S A 1996,93(21):11687–11692.PubMedView Article
O'Keefe DD, Prober DA, Moyle PS, Rickoll WL, Edgar BA: Egfr/Ras signaling regulates DE-cadherin/Shotgun localization to control vein morphogenesis in the Drosophila wing.Dev Biol 2007,311(1):25–39.PubMedView Article
Classen AK, Anderson KI, Marois E, Eaton S: Hexagonal packing of Drosophila wing epithelial cells by the planar cell polarity pathway.Dev Cell 2005,9(6):805–817.PubMedView Article
Eaton S, Wepf R, Simons K: Roles for Rac1 and Cdc42 in planar polarization and hair outgrowth in the wing of Drosophila.J Cell Biol 1996,135(5):1277–1289.PubMedView Article
Oda H, Uemura T, Harada Y, Iwai Y, Takeichi M: A Drosophila homolog of cadherin associated with armadillo and essential for embryonic cell-cell adhesion.Dev Biol 1994,165(2):716–726.PubMedView Article
Halfon MS, Gallo SM, Bergman CM: REDfly 2.0: an integrated database of cis-regulatory modules and transcription factor binding sites in Drosophila.Nucleic Acids Res 2008,36(Database issue):D594-D598.PubMed
Orian A, van Steensel B, Delrow J, Bussemaker HJ, Li L, Sawado T, Williams E, Loo LW, Cowley SM, Yost C, et al.: Genomic binding by the Drosophila Myc, Max, Mad/Mnt transcription factor network.Genes Dev 2003,17(9):1101–1114.PubMedView Article
Gauhar Z, Sun LV, Hua S, Mason CE, Fuchs F, Li TR, Boutros M, White KP: Genomic mapping of binding regions for the Ecdysone receptor protein complex.Genome Res 2009,19(6):1006–1013.PubMedView Article
van Steensel B, Delrow J, Bussemaker HJ: Genomewide analysis of Drosophila GAGA factor target genes reveals context-dependent DNA binding.Proc Natl Acad Sci U S A 2003,100(5):2580–2585.PubMedView Article
Yamaguchi M, Hayashi Y, Matsukage A: Essential role of E2F recognition sites in regulation of the proliferating cell nuclear antigen gene promoter during Drosophila development.J Biol Chem 1995,270(42):25159–25165.PubMedView Article
Li L, Edgar BA, Grewal SS: Nutritional control of gene expression in Drosophila larvae via TOR. Myc and a novel cis-regulatory element.BMC Cell Biol 2010, 11:7.PubMedView Article
Li XY, Thomas S, Sabo PJ, Eisen MB, Stamatoyannopoulos JA, Biggin MD: The role of chromatin accessibility in directing the widespread, overlapping patterns of Drosophila transcription factor binding.Genome Biol 2011,12(4):R34.PubMedView Article
Thomas S, Li XY, Sabo PJ, Sandstrom R, Thurman RE, Canfield TK, Giste E, Fisher W, Hammonds A, Celniker SE, et al.: Dynamic reprogramming of chromatin accessibility during Drosophila embryo development.Genome Biol 2011,12(5):R43.PubMedView Article
Grewal SS, Li L, Orian A, Eisenman RN, Edgar BA: Myc-dependent regulation of ribosomal RNA synthesis during Drosophila development.Nat Cell Biol 2005,7(3):295–302.PubMedView Article
Hirose F, Yamaguchi M, Kuroda K, Omori A, Hachiya T, Ikeda M, Nishimoto Y, Matsukage A: Isolation and characterization of cDNA for DREF, a promoter-activating factor for Drosophila DNA replication-related genes.J Biol Chem 1996,271(7):3930–3937.PubMedView Article
Hyun J, Jasper H, Bohmann D: DREF is required for efficient growth and cell cycle progression in Drosophila imaginal discs.Mol Cell Biol 2005,25(13):5590–5598.PubMedView Article
Dimova DK, Stevaux O, Frolov MV, Dyson NJ: Cell cycle-dependent and cell cycle-independent control of transcription by the Drosophila E2F/RB pathway.Genes Dev 2003,17(18):2308–2320.PubMedView Article
Martin-Castellanos C, Edgar BA: A characterization of the effects of Dpp signaling on cell growth and proliferation in the Drosophila wing.Development 2002,129(4):1003–1013.PubMed
Rogulja D, Irvine KD: Regulation of cell proliferation by a morphogen gradient.Cell 2005,123(3):449–461.PubMedView Article
Yamashita D, Sano Y, Adachi Y, Okamoto Y, Osada H, Takahashi T, Yamaguchi T, Osumi T, Hirose F: hDREF regulates cell proliferation and expression of ribosomal protein genes.Mol Cell Biol 2007,27(6):2003–2013.PubMedView Article
Pavlopoulos A, Akam M: Hox gene Ultrabithorax regulates distinct sets of target genes at successive stages of Drosophila haltere morphogenesis.Proc Natl Acad Sci U S A 2011,108(7):2855–2860.PubMedView Article
Jo S, Lee H, Kim S, Hwang EM, Park JY, Kang SS, Chung H: Inhibition of PCGF2 enhances granulocytic differentiation of acute promyelocytic leukemia cell line HL-60 via induction of HOXA7.Biochem Biophys Res Commun 2011,416(1–2):86–91.PubMedView Article
Louvi A, Wassef M: Ectopic engrailed 1 expression in the dorsal midline causes cell death, abnormal differentiation of circumventricular organs and errors in axonal pathfinding.Development 2000,127(18):4061–4071.PubMed
Cosgaya JM, Aranda A, Cruces J, Martin-Blanco E: Neuronal differentiation of PC12 cells induced by engrailed homeodomain is DNA-binding specific and independent of MAP kinases.J Cell Sci 1998,111(Pt 16):2377–2384.PubMed
Verzi MP, Shin H, He HH, Sulahian R, Meyer CA, Montgomery RK, Fleet JC, Brown M, Liu XS, Shivdasani RA: Differentiation-specific histone modifications reveal dynamic chromatin interactions and partners for the intestinal transcription factor CDX2.Dev Cell 2010,19(5):713–726.PubMedView Article
Wu LH, Lengyel JA: Role of caudal in hindgut specification and gastrulation suggests homology between Drosophila amnioproctodeal invagination and vertebrate blastopore.Development 1998,125(13):2433–2442.PubMed
Hacker U, Kaufmann E, Hartmann C, Jurgens G, Knochel W, Jackle H: The Drosophila fork head domain protein crocodile is required for the establishment of head structures.EMBO J 1995,14(21):5306–5317.PubMed
Lehmann OJ, Sowden JC, Carlsson P, Jordan T, Bhattacharya SS: Fox's in development and disease.Trends Genet 2003,19(6):339–344.PubMedView Article
Jasper H, Benes V, Atzberger A, Sauer S, Ansorge W, Bohmann D: A genomic switch at the transition from cell proliferation to terminal differentiation in the Drosophila eye.Dev Cell 2002,3(4):511–521.PubMedView Article
Hochheimer A, Zhou S, Zheng S, Holmes MC, Tjian R: TRF2 associates with DREF and directs promoter-selective gene expression in Drosophila.Nature 2002,420(6914):439–445.PubMedView Article
Yamaguchi M, Hayashi Y, Nishimoto Y, Hirose F, Matsukage A: A nucleotide sequence essential for the function of DRE, a common promoter element for Drosophila DNa replication-related genes.J Biol Chem 1995,270(26):15808–15814.PubMedView Article
Ohler U, Liao GC, Niemann H, Rubin GM: Computational analysis of core promoters in the Drosophila genome.Genome Biol 2002,3(12):RESEARCH0087.PubMedView Article
Goodrich JA, Tjian R: Unexpected roles for core promoter recognition factors in cell-type-specific transcription and gene regulation.Nat Rev Genet 2010,11(8):549–558.PubMed
Beckstead RB, Lam G, Thummel CS: The genomic response to 20-hydroxyecdysone at the onset of Drosophila metamorphosis.Genome Biol 2005,6(12):R99.PubMedView Article
Zinke I, Schutz CS, Katzenberger JD, Bauer M, Pankratz MJ: Nutrient control of gene expression in Drosophila: microarray analysis of starvation and sugar-dependent response.EMBO J 2002,21(22):6162–6173.PubMedView Article
Buttitta L, Katzaroff AJ, Edgar B: A robust cell cycle control mechanism limits E2F-induced proliferation of terminally differentiated cells in vivo.J Cell Biol 2010,189(6):981–996.PubMedView Article
Hart CM, Cuvier O, Laemmli UK: Evidence for an antagonistic relationship between the boundary element-associated factor BEAF and the transcription factor DREF.Chromosoma 1999,108(6):375–383.PubMedView Article
Hiller M, Chen X, Pringle MJ, Suchorolski M, Sancak Y, Viswanathan S, Bolival B, Lin TY, Marino S, Fuller MT: Testis-specific TAF homologs collaborate to control a tissue-specific transcription program.Development 2004,131(21):5297–5308.PubMedView Article
Fang X, Adler PN: Regulation of cell shape, wing hair initiation and the actin cytoskeleton by Trc/Fry and Wts/Mats complexes.Dev Biol 2010,341(2):360–374.PubMedView Article
White KP, Rifkin SA, Hurban P, Hogness DS: Microarray analysis of Drosophila development during metamorphosis.Science 1999,286(5447):2179–2184.PubMedView Article
Arbeitman MN, Furlong EE, Imam F, Johnson E, Null BH, Baker BS, Krasnow MA, Scott MP, Davis RW, White KP: Gene expression during the life cycle of Drosophila melanogaster.Science 2002,297(5590):2270–2275.PubMedView Article
Rifkin SA, Kim J, White KP: Evolution of gene expression in the Drosophila melanogaster subgroup.Nat Genet 2003,33(2):138–144.PubMedView Article
Li TR, White KP: Tissue-specific gene expression and ecdysone-regulated genomic networks in Drosophila.Dev Cell 2003,5(1):59–72.PubMedView Article
Marco A, Konikoff C, Karr TL, Kumar S: Relationship between gene co-expression and sharing of transcription factor binding sites in Drosophila melanogaster.Bioinformatics 2009,25(19):2473–2477.PubMedView Article
FitzGerald PC, Sturgill D, Shyakhtenko A, Oliver B, Vinson C: Comparative genomics of Drosophila and human core promoters.Genome Biol 2006,7(7):R53.PubMedView Article
Terriente-Felix A, Lopez-Varea A, de Celis JF: Identification of genes affecting wing patterning through a loss-of-function mutagenesis screen and characterization of med15 function during wing development.Genetics 2010, 185:671–684.PubMedView Article
Aoyagi N, Wassarman DA: Developmental and transcriptional consequences of mutations in Drosophila TAF(II)60.Mol Cell Biol 2001, 21:6808–6819.PubMedView Article
Wassarman DA, Aoyagi N, Pile LA, Schlag EM: TAF250 is required for multiple developmental events in Drosophila.Proc Natl Acad Sci U S A 2000, 97:1154–1159.PubMedView Article
Prince F, Katsuyama T, Oshima Y, Plaza S, Resendez-Perez D, Berry M, Kurata S, Gehring WJ: The YPWM motif links Antennapedia to the basal transcriptional machinery.Development 2008, 135:1669–1679.PubMedView Article
Thacker SA, Bonnette PC, Duronio RJ: The contribution of E2F-regulated transcription to Drosophila PCNA gene function.Curr Biol 2003,13(1):53–58.PubMedView Article
Nicolay BN, Bayarmagnai B, Islam AB, Lopez-Bigas N, Frolov MV: Cooperation between dE2F1 and Yki/Sd defines a distinct transcriptional program necessary to bypass cell cycle exit.Genes Dev 2011,25(4):323–335.PubMedView Article
Lu J, Ruhf ML, Perrimon N, Leder P: A genome-wide RNA interference screen identifies putative chromatin regulators essential for E2F repression.Proc Natl Acad Sci U S A 2007,104(22):9381–9386.PubMedView Article
Butler JE, Kadonaga JT: Enhancer-promoter specificity mediated by DPE or TATA core promoter motifs.Genes Dev 2001,15(19):2515–2519.PubMedView Article
Willy PJ, Kobayashi R, Kadonaga JT: A basal transcription factor that activates or represses transcription.Science 2000,290(5493):982–985.PubMedView Article
Sayal R, Ryu SM, Arnosti DN: Optimization of reporter gene architecture for quantitative measurements of gene expression in the Drosophila embryo.Fly (Austin) 2011,5(1):47–52.View Article
Graveley BR, Brooks AN, Carlson JW, Duff MO, Landolin JM, Yang L, Artieri CG, van Baren MJ, Boley N, Booth BW, et al.: The developmental transcriptome of Drosophila melanogaster.Nature 2011,471(7339):473–479.PubMedView Article
Link N, Chen P, Lu WJ, Pogue K, Chuong A, Mata M, Checketts J, Abrams JM: A collective form of cell death requires homeodomain interacting protein kinase.J Cell Biol 2007,178(4):567–574.PubMedView Article
Deato MD, Tjian R: Switching of the core transcription machinery during myogenesis.Genes Dev 2007,21(17):2137–2149.PubMedView Article
Deato MD, Tjian R: An unexpected role of TAFs and TRFs in skeletal muscle differentiation: switching core promoter complexes.Cold Spring Harb Symp Quant Biol 2008, 73:217–225.PubMedView Article
Bolstad BM, Irizarry RA, Astrand M, Speed TP: A comparison of normalization methods for high density oligonucleotide array data based on variance and bias.Bioinformatics 2003,19(2):185–193.PubMedView Article
Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ, Scherf U, Speed TP: Exploration, normalization, and summaries of high density oligonucleotide array probe level data.Biostatistics 2003,4(2):249–264.PubMedView Article
Ren N, Zhu C, Lee H, Adler PN: Gene expression during Drosophila wing morphogenesis and differentiation.Genetics 2005, 171:625–638.PubMedView Article
Prober DA, Edgar BA: Ras1 promotes cellular growth in the Drosophila wing.Cell 2000,100(4):435–446.PubMedView Article
Buttitta LA, Katzaroff AJ, Perez CL, de la Cruz A, Edgar BA: A double-assurance mechanism controls cell cycle exit upon terminal differentiation in Drosophila.Dev Cell 2007,12(4):631–643.PubMedView Article
de la Cruz AF, Edgar BA: Flow cytometric analysis of Drosophila cells.Methods Mol Biol 2008, 420:373–389.PubMedView Article
Bailey TL, Boden M, Buske FA, Frith M, Grant CE, Clementi L, Ren J, Li WW, Noble WS: MEME SUITE: tools for motif discovery and searching.Nucleic Acids Res 2009,37(Web Server issue):W202-W208.PubMedView Article
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (
http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.