Hotspots of transcription factor colocalization in the genome of Drosophila melanogaster

We began our study in Drosophila by mapping the in vivo genomic-binding sites of seven transcription factors from different classes and with diverse physiological functions: Bicoid (Bcd), GAGA factor (Gaf), Jun-related antigen (Jra), Max, odd paired (opa), Ecdysone receptor (EcR) isoform B1 (EcRB1), and its heterodimerization partner Ultraspiracle (USP) (Table 1, which is published as supporting information on the PNAS web site). Except for heterodimerization of EcRB1 and USP, no physical or functional interactions have been reported between any of these proteins (www.flybase.org).

The localization patterns of these seven factors were determined in the embryonic Kc cell line, which provides a homogeneous cell population. We used the DamID technology, which has been used to map genomic binding of a variety of proteins (14, 16–20). DamID involves the in vivo expression of a chimeric protein consisting of a transcription factor of interest fused to Escherichia coli DNA adenine methyltransferase (Dam) (15). The expression levels of the Dam-fusion proteins are kept very low to avoid mistargeting due to overexpression (14, 15). DNA in the close vicinity of the natural binding sites of the transcription factor is methylated preferentially by the tethered Dam. Methylated DNA fragments are isolated by using a PCR-based method, and microarrays then are used to detect the pattern of targeted adenine methylation, from which the localization pattern of the transcription factor can be deduced (14, 16, 18, 20). We used genomic tiling arrays containing 3,648 genomic fragments of 430–920 bp, together covering 2.9 Mbp of the Adh-cactus region of chromosome 2 of Drosophila (16).

The results for all seven transcription factors are displayed graphically in a detailed map of the Adh-cactus region (Figs. 1A and B and 2A; see also Fig. 7, which is published as supporting information on the PNAS web site). The fraction of probed regions targeted by individual proteins ranged between 3–8%. Surprisingly, we found that there is a high degree of overlap between the localization profiles of these seven factors. Each factor shares at least 67% of its binding sites with one or more other factors, and many sites are targeted by all seven factors (see below). Correlation analysis showed that the similarity of the profiles is in all cases statistically highly significant (Fig. 8, which is published as supporting information on the PNAS web site). The high degree of overlap is not restricted to the 2.9 Mbp Adh-cactus region; we also observed it for three of the transcription factors (Bcd, Gaf, and Jra) on arrays containing ≈12,000 cDNA fragments, corresponding to ≈60% of the Drosophila genes (Fig. 9, which is published as supporting information on the PNAS web site). MatrixREDUCE (21, 22) analysis of the DamID profiles of Bcd, EcRB1, Gaf, Max, and Jra shows that they correlate strongly with the in vitro sequence specificity predicted by using their respective known binding motifs (Table 2, which is published as supporting information on the PNAS web site; also see Materials and Methods). USP correlates weakly with both its own binding motif and the binding motif of its heterodimerization partner, EcR (P = 0.057 and 0.035, respectively.). Using a reporter gene containing an ecdysone response element, we could also show that Dam does not interfere with the function of EcRB1 (L.V.S. and K.P.W., unpublished results). These results suggest that the Dam-fusion proteins are functional. Taken together, our results reveal a high degree of colocalization for a wide range of transcription factors that (with the exception of EcRB1 and USP) were thought to be functionally unrelated.

To further verify our DamID results, we performed ChIP, a method for detecting protein–DNA interactions that is fundamentally different from DamID (4, 5, 23). The ChIP profile of endogenous Gaf (i.e., in untransfected cells) is in good agreement with the Gaf DamID profile (Fig. 3). This agreement also was reported recently in an independent comparison (24). We also confirmed the DamID microarray data for a subset of transcription factors and genomic locations by using quantitative PCR (Fig. 10, which is published as supporting information on the PNAS web site), ruling out that our results are caused by a microarray artifact.

To investigate whether the high degree of colocalization is caused by nonspecific protein–DNA interactions, we determined the DamID profile of the DNA-binding domain (DBD) of the sequence-specific yeast transcription factor Gal4p (Fig. 1C). We found that Gal4p DBD binds significantly to 16 fragments on the array (corresponding to 11 loci), but there is no overlap with the localization profile of the other proteins (Figs. 1 and 8). The observed DamID profile of Gal4 DBD correlates strongly with the Gal4p consensus motif (Table 2), which confirms that the Gal4 DBD-Dam fusion protein is functional. We conclude that the colocalization of the seven Drosophila transcription factors is not simply due to nonspecific DNA binding.

Transcription factors rely on interactions with other proteins, such as corepressors, coactivators, and chromatin proteins, to regulate gene expression. We analyzed whether these proteins also associate with transcription factor-binding sites. Indeed, the corepressors Groucho (Gro), Rpd3, and Sin3 and the coactivator Brahma (Brm) colocalize strongly with the transcription factors and with each other (Figs. 2, 7, and 8). The heterochromatin proteins HP1a and Su(var)3-9 do not associate with the same loci as the transcription factors (Figs. 2, 7, and 8), but they do target transposable elements and a few other loci, consistent with previously reported data (16, 18). Heterochromatin protein 1c (HP1c) a euchromatic homologue of HP1a of unknown function (18, 25), shows a similar profile as the cofactors (Figs. 2, 7, and 8), suggesting that HP1c may act as a transcriptional regulator.

Visual inspection of the DamID patterns indicated that the Adh-cactus region contains many sites where all transcription factors and cofactors associate together. To identify such sites in a systematic and unbiased way, we applied machine learning methods to classify protein composition along the Adh-cactus region of Drosophila into a number of distinct “chromatin types.” We first used a self-organizing map (SOM) algorithm (26) to generate a 2D representation of the DamID data for all 14 proteins, in which different areas correspond to genomic regions with different protein compositions (Fig. 4A). Next, we used adaptive K means clustering (29) to assign each probed genomic fragment to one of eight distinct chromatin types (Fig. 4B). This number of chromatin types was determined in an unsupervised manner (29). Although this classification is likely to represent an oversimplification of chromatin diversity, it provides a useful framework for further analysis. Chromatin types that can be distinguished are characterized by, for instance, binding by HP1a and Su(var)3-9 only (“heterochromatin”; fragments of type 5) and no binding of any protein (fragments of type 2). Fragments of type 1 show high binding by all proteins except HP1a and Su(var)3-9 and include 166 fragments congregating in 61 loci. We will refer to loci of type 1 as “colocalization hotspots” (Figs. 2B and 7). We also visualized the DamID data of GAL4-DBD and the ChIP data of Gaf in the SOM plane (Fig. 4C). This visualization shows clearly that GAL4-DBD as detected by DamID is completely absent from hotspots, whereas binding of Gaf as detected by ChIP is strongly enriched in hotspots, in agreement with the Gaf DamID results.

We computationally screened a collection of transcription factors of known sequence specificity (30, 31) for preferential binding to colocalization hotspots (see Materials and Methods). Fig. 5 shows that the hotspot regions on average have higher predicted affinity for three of the proteins that we mapped in this study, Gaf, Jra, and Max. In addition, we found motifs of several other proteins to be enriched in the colocalization hotspots (Table 3, which is published as supporting information on the PNAS web site), predicting that these proteins also bind to hotspots. Surprisingly, the motifs for Bcd, EcR, and USP are not enriched in hotspots (Fig. 5). This result suggests that these latter factors are not recruited to hotspots via protein–DNA interactions mediated by their own DNA-binding domain but rather via protein–protein interactions with one or more other DNA-bound proteins.

We investigated the mechanism of regulator targeting to hotspots more directly for the transcription factor Bcd. For this purpose, we generated DamID profiles of Bcd mutants that either lacked a functional DBD (Bcd^K50A; refs. 32 and 33) or consisted of the DBD alone (Fig. 7). Strikingly, Bcd^K50A localizes to the hotspots, whereas the DBD alone is not detected in most hotspots (Fig. 4D). The Bcd DBD binds significantly to 69 fragments on the array, only ≈25% of which also are bound by WT Bcd. Its DamID signal correlates significantly better with the predicted affinity for Bcd than the DamID profile of WT Bcd; by contrast and as expected, the Bcd^K50A DamID signal entirely lacks correlation with the Bcd consensus (Table 2). Taken together, these results suggest that protein–protein interactions help target Bcd to the hotspots, possibly in competition with the direct recruitment of Bcd to its DNA-binding sites in the genome.

We noticed that colocalization hotspots are often located in predicted genes: 41 of the 61 hotspot regions overlap with 37 genes (Fig. 2C). To test whether hotspots occur preferentially in actively transcribed regions, we first mapped the binding of the 18-kDa subunit of RNA polymerase II (Pol II) in the Adh-cactus region. The DamID profile of Pol II shows strong overlap with the binding profiles of the transcription factors (Figs. 1, 7, and 8), and hotspots were significantly enriched in Pol II (Fig. 4E). Essentially the same result was obtained when we analyzed a previously reported Pol II-binding profile in Kc cells that was obtained by ChIP (27), again indicating that the overlap is not caused by a bias in the DamID method (Figs 4E and 8). Furthermore, comparison of our data to chromosomewide expression profiling in Kc cells (27) also shows that mRNA levels correlate with the binding of the transcription factors (Fig. 8) and are significantly elevated in colocalization hotspots (Fig. 4E).

Recent studies in Drosophila cells have shown that the histone variant H3.3 is specifically deposited in transcribed genes and their flanking regions (28). Analysis of these data revealed that colocalization hotspots are strongly enriched in H3.3 (Fig. 4E). We found that 95% of the hotspot probes have H3.3/H3 ratios that are above the mean of the entire Adh-cactus region (P < 2e−16; Student t test), indicating that colocalization hotspots virtually are restricted to chromatin regions containing H3.3. Taken together, our analyses indicate that hotspots occur preferentially in regions of active transcription.

We wondered whether colocalization hotspots, which we identified in Kc cells, also may be linked to sites of transcription during fly development. We therefore analyzed the developmental expression profiles of genes overlapping with hotspots by using published genomewide expression profiles from six different developmental stages (34). Strikingly, this analysis (Fig. 6) showed that during embryogenesis, hotspot-containing genes are expressed at significantly higher levels than hotspot-free genes. This difference disappears later in development. A modest increase of hotspot-associated gene expression also is observed in adult females, which may be due to the fact that female ovaries contain large amounts of oocytes with maternally loaded mRNA. These results indicate that colocalization hotspots occur preferentially in regions with increased transcription activity during early development.

We report here that the target loci of several unrelated transcription factors and cofactors strongly overlap with each other in the genome of Drosophila. Colocalization hotspots, which are targeted by all seven transcription factors and six of eight other proteins tested, represent extreme cases of this phenomenon. Because Drosophila contains ≈700 transcription factors (35), we predict that many more factors will show strong colocalization and that hotspots may recruit up to a few hundred different proteins. This observation contrasts with results in yeast where overlap between transcription factors was relatively rare and hotspots were not apparent (7, 8).

The mechanism by which such a large number of proteins are recruited to specific genomic loci remains to be elucidated. Our observation that Bcd does not require its DBD for hotspot association, together with the lack of enrichment of the consensus motifs for USP, Jra, and Bcd in hotspots, indicates that at least some transcription factors do not directly bind the DNA in hotspots, but instead are recruited through protein–protein interactions. We note that the DamID maps represent the average protein occupancy over time and over many cells, and, therefore, it is quite possible that the composition of hotspot-associated complexes is variable and dynamic (36).

The function of hotspots is unclear, but at least three models can be envisioned that are not mutually exclusive. First, hotspots may represent general “sinks” or “buffers” that sequester many regulator molecules. Thus, hotspots may regulate the effective concentration of free transcription factor molecules in the nucleus. This model is in agreement with our observation that wild-type Bcd fails to bind to several loci containing the Bcd consensus motif, whereas the Bcd DBD (that is not tethered to hotspots) does bind to these loci. Second, hotspots may be directly involved in the regulation of nearby genes, similar to enhancers. Unlike classical enhancers, however, hotspots may associate with hundreds of proteins, many of which may have only a minor contribution to the overall effect on gene expression. Third, hotspots could mediate physical interactions between distant loci in the genome, as has been reported for several transcription factor-rich loci (37, 38). Regardless of what could be the regulatory function of hotspots, our results have important implications for our understanding of the mechanisms that determine the target specificity of transcription factors and for the interpretation of genomewide protein occupancy data in higher eukaryotes.

We thank Susan Parkhurst and Amir Orian (Fred Hutchinson Cancer Research Center, Seattle, WA) and David Wassarman (University of Wisconsin, Madison, WI) for plasmids; Nicolas Nègre and Giacomo Cavalli (Institute of Human Genetics, Centre National de la Recherche Scientifique, Montpellier, France) for ChIP protocol and Gaf antibody; Marcel van Batenburg and Barrett Foat for technical advice; and Melissa Davis for preparing microarrays. Some materials were received through the Drosophila Genomics Resource Center. This work was supported by a Human Frontier Science Program Grant (to B.v.S., K.P.W., and H.J.B.); National Institutes of Health Grants HG003008, CA12152, and R24 GM074105 (to H.J.B.) and R01 HG03231 (to K.P.W); the W.M. Keck Foundation (K.P.W.), the Beckman Foundation (K.P.W.), and a European Young Investigators Award (to B.v.S.).