BACKGROUND
transcriptional control of gene expression via localisation and binding of transcription factors at the vicinity of gene loci lies at the heart of metazoan development.
with the advances in sequencing protocols, chip sequencing is rapidly becoming the preferred tool to find genome-wide binding patterns  of a transcription factor in a given cell type.
as chip sequencing protocols mature, peaks of multiple transcription factors have been characterized in single cell types to study combinatorial control  <cit> .
one of the striking observations of these multi-tf studies was detection of so-called ‘transcriptional hotspots’.
in  <dig>  moorman et al.
<cit>  generated genome-wide binding profiles of seven transcription factors in d. melanogaster and identified a subset of peaks bound by all seven tfs .
of these  <dig> hotspots, when tested using transgenic assays, 94% acted as enhancers strongly activating the neighbouring developmentally important genes  <cit> .
the presence of hotspot regions was further confirmed in other species such as worm  <cit>  and humans  <cit> .
hotspots in c. elegans and d. melanogaster were enriched for sequence binding motifs of many tfs including signature motifs such as gaga and zelda  <cit> .
however, analysis of human encode data did not support these observations.
in humans, hotspots were specifically deprived of cis-regulatory motifs and no signature motif similar to the ‘gaga’ motif was found  <cit> .
at the other end of the spectrum, most tf binding events were ‘singletons’  accounting for nearly one third of the binding events in a cell type.
these regions when studied in d. melanogaster transgenic assays did not drive patterned reporter gene expression leading to the conclusion that they do not act as strong developmental enhancers  <cit> .
this leads to a series of questions about properties of both hotspots and singletons.
for example, are they truly distinct genomic regions, do they have a characteristic sequence or chromatin signature?
in order to answer the above questions, we collected genome-wide binding patterns of multiple transcription-related factors in ten murine cell types.
for each cell type, the peaks were classified into three groups: singleton genomic regions occupied by only one tf , hotspot genomic regions occupied by most tfs under study  and combinatorial genomic regions occupied by a combination of tfs .
the genomic regions largely clustered according to the group suggesting distinct genomic regions marked for occupancy independent of cell type.
the singletons and combinatorials were bound neighbouring similar genes in all cell types while hotspot peaks occurred near a distinct set of genes in each cell type and showed functional enrichment for cell type specific genes.
though hotspots were enriched for many tf sequence motifs, no signature motif such as gaga motif was found in murine cell types.
we identified h3k4me <dig> chromatin modification distinguishing hotspots from singletons where hotspots showed a bimodal h3k4me <dig> peak whereas singletons were mono-modal.
finally, we collected the differentially regulated genes after perturbations of multiple transcription-related factors in es cells, to show that hotspots were preferentially bound by activators and not repressors.
as genes differentially expressed after setdb <dig> knockout but not dnmt knockout were enriched for es hotspot genes we suggest that setdb <dig> might be involved in silencing hotspots.
RESULTS
combinatorial binding events overlap across multiple cell types
a chip sequencing experiment typically identifies thousands to hundreds of thousands of genome-wide binding sites of a tf in a cell type.
in order to investigate if distinct characteristics of peaks grouped based on the number of binding factors at a binding location, we collected genome-wide binding events of transcription-related factors  across  <dig> different normal and cancer murine cell types.
in each cell type, chip sequencing data sets  for  <dig> to  <dig> tfs  were collected.
peaks were called in each sample  with the number of peaks ranging from  <dig>  to  <dig>  in  <dig> samples .
the number of tfs at binding sites followed an exponential curve such that about half of the regions were bound by only one tf in any given cell type and less than  <dig> % occupied by all tfs studied.
we divided all peaks in a cell type into three groups: singletons or low-occupancy, peaks bound by only one transcription-related factor , combinatorials or mid-occupancy, peaks bound by a combination of transcription-related factors  and hotspots or high-occupancy, peaks bound by more than five tfs studied in a given cell type .
on average about 50% of all the peaks in a cell type were singletons and combinatorials while only  <dig> -2% was classified as hotspots .
we assigned peaks in these three groups to genomic locations: promoter, 5′ utr, 3′ utr, exon, intron or inter-genic using homer  <cit> .
the location of peaks followed a uni-modal distribution with a peak at the tss .
singleton peaks were specifically under-represented in promoter regions  and 5′ utr regions  .figure  <dig>
combinatorial binding events overlap across multiple cell types.
a. schematic diagram of a genomic region with genome tracks of six chip sequencing samples in a cell type marking singletons , combinatorials  and hotspots .
b. a table listing all  <dig> cell types used along with the number of transcription factors c. bar graph of fraction of peaks in hotspots across  <dig> cell types.
d. heatmap of pearson’s correlation coefficients of peaks in  <dig> sets  with a hierarchical clustering tree showing that combinatorials clustered according to the groups.
e. box plots for peak heights for the three groups  for t cells and macrophages.
in order to investigate the overlap of the three groups across cell types, we constructed a hierarchical tree of  <dig> peak sets .
to construct hierarchical tree, we first integrated of all  <dig> peak sets across  <dig> cell types resulted in  <dig>  unique peaks.
we then constructed a binary matrix of  <dig>  rows and  <dig> columns recording the presence or absence of peaks in each group for each cell type and built hierarchical tree using pearson’s correlation coefficient as a distance measure.
many samples clustered according to the group across multiple cell types .
the peaks in three groups in a given cell type are mutually exclusive by definition but the fact that combinatorials across many cell types cluster according to the groups is a notable observation.
this suggests that each genomic location is marked for occupancy i.e.
mid-occupied loci in one cell type were more likely to be bound by multiple factors rather than one factor in another cell type.
hotspots of related cell types cluster into two tight clusters.
the two embryonic stem  cell samples were the most similar to each other as expected.
the four related blood cell types – hpc, erythroid, mk progenitors and b cells formed second separate dense cluster of hotspot peaks.
in d. melanogaster, peaks occupied at high levels by transcription factors in vivo drive patterned gene expression, whereas those occupied only at lower levels mostly do not  <cit> .
we calculated average peak height for the three groups where peak height was defined as the total number of tags in a 400 bp window around the peak summit .
singletons consistently showed significantly lower peak height compared to combinatorials and hotspots .
for example, the average peak height of fli <dig> bound singleton regions in t cells is significantly lower than combinatorials  as were the average peak heights of cebpβ bound regions in macrophages where singletons were significantly lower  than combinatorials .
finally, we collected chip sequencing data for oct <dig> in es cells generated by three independent labs  <cit>  and found about  <dig> oct <dig> peaks overlapping in all three experiments.
the  <dig> peaks were enriched for combinatorials and hotspots but were depleted for singleton regions .
combinatorials and hotspot binding events were found consistently across different experiments.
taken together, singletons were depleted in promoter and 5’ utr regions, had lower peak height and were not consistently found across multiple experiments.
hotspots preferentially bind in cell type-specific gene neighbourhoods
the hierarchical clustering of  <dig> peak sets clustered them largely according to occupancy.
to investigate if the group-wise clustering of peaks is also reflected at a gene level, we mapped all peaks to the nearest gene resulting into a total of  <dig>  genes.
we created a binary gene matrix of  <dig>  rows and  <dig> columns and built a hierarchical tree using pearson’s correlation coefficient as a distance measure.
in contrast to the peak overlap, the gene overlap clustered all singleton and combinatorial peaks together forming one cluster while hotspot genes showed a lower overlap across cell types .
by considering peaks only in promoters and performing a hierarchical clustering, the three groups clustered largely together and combinatorials formed one tight cluster .
thus, combinatorials were located in the vicinity of similar genes across cell types but hotspots were present near a distinct set of genes in each sample.
the functional enrichment of hotspot genes highlighted enrichment for cell type specific functions , such as b cell hotspots which were over-represented for the b cell receptor signalling pathway  and b cell activation , while stem cell hotspots were enriched for stem cell differentiation  and stem cell development .
though binding events of all peaks were enriched for cell type specific functions, hotspot genes had a much higher cell type specific functional enrichment than all peaks in a cell type.figure  <dig>
hotspots preferentially bind in cell type-specific gene neighbourhoods.
a. a heatmap of pearson’s correlation coefficients of genes in  <dig> sets  with a hierarchical clustering tree showing that singletons and combinatorials clustered together in a tight cluster while hotspot gene sets were very distinct from each other.
b. fraction of es super enhancers defined by whyte et al.
<cit>  in three groups showing hotspots significantly overlapped with super enhancers.
c. fraction of high oct <dig> bound regions in oct4+/− compared to oct4+/+  <cit>  es cells in three groups showing hotspots enriched for high oct <dig> bound regions in oct4+/− and also pluripotent genes d. functional enrichments of hotspot genes in  <dig> cell types along with p values showing that hotspot genes were enriched for cell type specific genes.
recently, three groups generated a list of cell type specific enhancers in mouse es cells using diverse experiments.
whyte et al.
<cit>  showed that master transcription factors such as oct, sox, and nanog form ‘super-enhancers’ at key cell identity genes that span large domains and drive cell-type-specific gene expression program.
we calculated the overlap of es cell peaks in all three groups with super-enhancers and noted that about 30% of super-enhancers overlapped with hotspots in es cells .
karwacki-neisius et al.
<cit>  showed that oct4+/−  es cells have increased oct <dig> binding at pluripotency enhancers.
only the enhancer regions strongly bound by oct <dig> in oct4+/− compared to oct4+/+ had a preferential bias towards hotspot regions .
as oct <dig> binding in oct4+/− was shown to be specifically enriched near pluripotent genes, this confirmed that hotspots occupied cell type specific gene loci.
stadler et al.
<cit>  generated a methylome map in murine es cells and noted that low methylated regions mark multi-tf bound enhancer regions in the vicinity of cell type specific genes.
the overlap of high- and low-methylated regions with three groups in es cells confirmed that hotspots were specifically depleted for high-methylated and enriched for low-methylated regions in es cells .
taken together, by comparing the three groups against three complementary experimental datasets we concluded that hotspots lie preferentially in the vicinity of cell type specific genes.
hotspots are enriched for sequence motifs of multiple transcription factors and are preferentially bound by pioneering factors
in order to search for properties distinguishing the three groups, we calculated the average mammalian conservation score for the genome sequences of peaks in all three groups.
there was no statistically significant difference in the fraction of conserved peaks between the three groups .
in c. elegans, hotspot-associated genes were more likely to be essential and were enriched for a variety of functions such as growth, reproduction, and larval and embryonic development  <cit> .
to check if this observation holds true in the mammalian system, we downloaded essential genes in mouse  <cit>  and investigated whether any of the three groups were preferentially bound near essential genes.
though there was a relatively smaller fraction of essential genes in singletons compared to hotspots, the difference was not statistically significant .
thus, genomic conservation and gene essentiality did not distinguish the three groups.
in d. melanogaster, hotspots were enriched for sequence motifs of many but not all tfs under study and those without sequence motif enrichment were shown to bind via protein-protein interaction  <cit> .
they were enriched for sequence motifs of two global activators zelda and gaga and this sequence signature alone is sufficient to distinguish the hotspot regions  <cit> .
they were also enriched for beaf- <dig> and trl/gaf insulator motifs  <cit> .
in c. elegans, cis-regulatory sequences of tfs are not always present in hotspots, the binding was associated with open chromatin  <cit> .
in line with this, hotspots were enriched for ‘gaga’ sequence motifs  <cit> .
on the other hand, the depletion of cis-regulatory motifs was observed in human hotspots where the sequence motif for tfs in ‘cold’ regions were more than two- fold more enriched than in hotspot regions and no ‘gaga’ or similar motif was found  <cit> , but they mentioned the possible presence of other sequence-specific motifs.
indeed foley et al.
found that on average, any given hot promoter contained about four different types of motifs  <cit> .
in summary there is no evidence so far regarding a specific cis-regulatory signature in hotspot peaks in mammals, nor is there a consensus whether hotspots are enriched or depleted for tf sequence motifs.
to check if any sequence motif signature discriminates the three groups, we calculated the number of occurrences of mammalian cis-regulatory sequence motifs and compared them between the three groups.
in most cell types, the motifs associated with the transcription factors known to be key players for that cell type were significantly enriched in hotspot regions compared to singletons.
for example, gata  and runx  motifs were enriched in b cell hotspots, e-box  and gfi <dig>  motifs were enriched in erythroid hotspots while meis  and scl-gata  motifs were enriched in hpc hotspots .
in contrast, there was a significant depletion of some motifs in hotspots compared to singletons.
in particular, the ctcf motif was less abundant in hotspots across multiple cell types .
in summary, hotspots were enriched for many known transcription factor motifs  including the one not characterized by chip sequencing such as hox motif in hpcs .
as gc content can bias the motif enrichment, we confirmed that the gc content across three groups was not significantly different across ten cell types .
moreover, we validated that the motif enrichment was consistent only when the promoter peaks were considered in three groups .figure  <dig>
hotspots are enriched for sequence motifs of multiple transcription factors.
a. bar graphs of fraction of peaks with a given sequence motif in the three groups  with the name of the sequence motif along with the sequence in iupaac format for three representative motifs for each of b cells, erythroid and hpcs .
b. a table of cell type, transcription-related factors analysed by chip sequencing, transcription-related factors binding to more than 90% of hotspot peaks and statistically significant known sequence motifs found in hotspots for b cells, erythroid and hpcs .
siersbaek et al.
identified ~ <dig>  tf hotspots during adipogenesis and demonstrated that cebpβ marks a large number of these hotspots and is required for chromatin remodeling during differentiation  <cit> .
this hypothesizes that pioneering tfs in a cell type mark hotspots.
we identified tfs preferentially binding in hotspots  and found pu.
<dig> marking b cell hotspots and cebps marking macrophage hotspots .
med <dig> was bound at more than 99% of es hotspots which proposes hotspots to be under control of mediator complex  <cit> .
thus hotspots regions were indeed bound by pioneering factors across multiple cell types.
on the other hand, only around 1% of peaks of suz <dig> overlapped with hotspots .
as suz <dig> predominantly acts as a transcriptional repressor this suggests that hotspots might act as transcriptional hubs mainly for activation.
altogether, murine hotspots did not contain a signature sequence motif but were enriched for motifs of many transcription factors including pioneering factors which mainly act as transcription activators.
hotspots lie near highly expressed genes and are enriched for enhancers in nucleosome flanked regions
in order to investigated whether each group shows a distinct chromatin signature, we collected data for three activating chromatin modifications in enhancers and promoters , ccctc-binding factor , rna polymerase ii  binding as well as rna sequencing data in murine erythroleukemia  cells from the encode resource.
hotspot peaks showed significantly higher binding in all six datasets compared to the other two groups.
in d. melanogaster, hotspot genes were expressed at higher levels during early embryogenesis but the difference is reduced later in development  <cit> .
in humans, using encode data, it was shown that genes associated with ‘hot’ regions were expressed at higher levels  <cit> .
in line with this, h3k4me <dig>  polii binding as well as gene expression measured in rpkm values using rna sequencing data  in mel cells were significantly higher in hotspots .figure  <dig>
hotspots lie near highly expressed genes and are enriched for enhancers in nucleosome flanked regions.
a. boxplots for the peak height for h3k4me <dig>  rna polymerase ii using chip sequencing as well as average rpkm value using rna sequencing data in mel cells for the three groups  showing higher peak height for all chromatin modifications in hotspots.
b. bar graph of fraction of vista enhancers  <cit>  and boxplots for the peak height for h3k4me <dig>  h3k27ac using chip sequencing in the three groups in mel cells concluding hotspots were enriched for enhancers.
c. the average density plot of mel h3k4me <dig> in the  <dig> kb window for each of the three groups showing hotspots had a bi-modal h3k4me <dig> peak while singletons have a mono-modal peak.
pennachio et al.
<cit>  predicted about  <dig>  high-confidence murine tissue-specific enhancers for  <dig> tissue types by integrating tissue-specific expression data, conservation information and cis-regulatory motifs and further experimentally validated them using transgenic assays in mouse .
mel hotspots were enriched in vista enhancers compared to the singletons and combinatorials .
hotspots in other cell types were also enriched for vista enhancers .
mel hotspots had significantly higher h3k4me <dig> as well as h3k27ac .
histone h3k4me <dig> marks for enhancer regions whereas histone h3k27ac distinguishes active enhancers from inactive/poised enhancer elements containing h3k4me <dig> alone  <cit> .
hotspots therefore mark active enhancer regions in a given cell type.
hoffman et al.
<cit>  showed high nucleosome positioning at the center of bimodal h3k4me <dig> peaks compared to mono-modal in mouse islets and liver.
the average density plots of h3k4me <dig> showed a very strong bimodal distribution at hotspots while singletons showed a mono-modal distribution .
this suggests that multiple transcription-related factors can indeed displace nucleosomes for strong enhancer activity at hotspot regions.
hotspots were therefore enriched for nucleosome flanked multiple tf-bound loci.
as expected, genes associated with flanked tf-bound loci were more abundantly expressed than those associated with nucleosomal loci, consistent with flanked sites being active enhancer elements  <cit> .
in order to validate these observations in another cell type, we collected genome-wide chromatin modification data sets in es cells.
similar to mel cells  <cit> , all activating chromatin modifications as well as expression values in rpkm using rna sequencing data showed significant enrichment in hotspots .
the repressive chromatin modifications such as h3k27me <dig> were depleted in hotspots .
taken together, we integrated genome wide chromatin binding patterns with tf chip sequencing data to show that hotspot regions were enriched for enhancer regions which are located near highly expressed cell type specific genes.
furthermore, the h3k4me <dig> profile showed a strong bimodal peak at hotspots suggesting nucleosome displacement for transcription activation.
es hotspots are preferentially occupied by activators and hotspots regions and might be under setdb <dig> control for silencing
in order to establish the functional relevance of the binding in three groups, we collected a range of differentially expressed gene sets after tf perturbations in es cells.
firstly we collected differentially expressed genes after oct <dig> knockdown from two independent experiments  <cit> .
es hotspot genes were specifically enriched for down-regulated genes  after oct <dig> knockdown but not up-regulated genes .
this suggests that oct <dig> binds to hotspot regions to act as a transcriptional activator.figure  <dig>
es hotspots are preferentially occupied by activators.
a. overlap of genes in three groups in es ng dataset with up- and down-regulated genes upon oct <dig> knockdown in es cells showing hotspots were enriched for only down-regulated genes.
b. fraction of up- and down-regulated genes upon dnmt and set1db knock out in three groups in es ng dataset showing hotspot enrichment only in set1db ko down-regulated genes.
nishiyama et al.
<cit>  identified potential functional targets by identifying differentially expressed genes after overexpressing  <dig> tfs in es cells.
es hotspot genes were over-represented for differentially expressed regulators such as etv <dig> .
however, es hotspot genes were under-represented  in differentially expressed genes after overexpression of repressors such as atf <dig> , gadd45a , mybl <dig> , rhox <dig>  and tcf <dig> .
this corroborates the previous observation that hotspot regions were predominantly occupied by activators of transcription .
not only transcription factors but also other chromatin modifiers such as baf <dig>  <cit> , a subunit of the baf chromatin remodelling complex, as well as hdac <dig>  <cit> , a histone deacetylase involved in the removal of acetyl groups from core histones, were preferentially bound at es hotspots .
enrichment in hotspot loci suggests hdac <dig> might act as a transcriptional activator.
though hdacs are generally considered repressors of transcription, it has been noted recently that hdac <dig> can indeed act as an activator  <cit> .
as hotspots mark highly active enhancers driving cell type specific gene expression, these regions need to be silenced in other cell types.
karimi et al.
<cit>  showed dna methylation and setdb <dig> methylation work on different sets of genes, where only 7% of the genes up-regulated in a setdb <dig> knockout also were up-regulated in dnmt knock out in es cells.
no significant difference was observed between singletons, combinatorial or hotspot genes for differentially expressed genes after dnmt knock out .
challen et al.
<cit>  found ~ <dig>  loci under dnmt3a control in hpcs.
there was no bias in dnmt <dig> occupied loci in hpc, neither towards cpg islands in the three groups.
however, genes differentially expressed after setdb <dig> knock out were enriched for hotspots  suggesting that hotspots might be regulated by setdb <dig>
altogether, hotspots might be used by activators to control expression of neighbouring genes as they were enriched for differentially expressed genes after perturbation of activators in es cells.
furthermore, we suggested that es hotspots might be under setdb <dig> control for silencing.
CONCLUSIONS
the combinatorial action of multiple tfs controls regulation of gene expression where tfs activate or repress gene expression via binding to cis-regulatory sequence elements in gene loci.
identification of functional target genes regulated by specific transcription factors is critical for understanding of the molecular mechanisms behind transcriptional control.
multi-tf chip sequencing studies have led to a remarkable observation that there were regions in the genome bound by almost all of the factors, the so called ‘transcriptional hotspots’  <cit> .
these regions acted as active enhancers of developmentally essential genes  <cit> .
on the other hand, low occupancy regions did not drive reporter gene expression  <cit>  questioning the functional significance of these peaks.
in order to study whether there are indeed classes of peaks with distinct characteristics based on their occupancy, we collected  <dig> tf binding chip sequencing datasets in  <dig> murine cell types.
all peaks in a cell type were classified into three groups based on the number of factors bound: singletons , combinatorials  and hotspots .
singletons were depleted in promoters and 3’ utr, had a low peak height and were less conserved across multiple biological replicates putting in doubt their functional relevance.
on the other hand, conserved cis-regulatory sequences were thought be more functionally relevant but the three groups could not be distinguished with respect to sequence conservation.
indeed, it has been shown in yeast that the low affinity protein-dna interactions were conserved as well as being functionally relevant  <cit> .
the hierarchical clustering of  <dig> peak sets  clustered combinatorials according to occupancy.
putative target gene sets of singletons and combinatorial peaks clustered into one tight cluster whereas each hotspot gene set consisted of unique set of genes with cell type specific functional gene enrichment.
the presence of hotspot regions near cell type specific genes has previously been observed in many cell types across species .
we have shown previously that most tf peaks cluster according to the cell type  <cit> .
gaga and zelda signature sequence motifs were associated with hotspots in c. elegans and d. melanogaster  <cit> .
we did not find a signature sequence motif for any of the groups but hotspots were enriched for sequence motifs of many key transcription factors as well as for pioneering factors in a cell type.
hotspots were specifically depleted for ctcf sequence motif across all cell types.
though depleted in ctcf sequence motifs compared to singletons, the fraction of ctcf peaks and their mean peak height at hotspots is higher .
we further intersected three groups with diverse chromatin marks and concluded that hotspots were enriched for enhancers present in the vicinity of highly expressed genes.
moreover, we identified a chromatin modification that can distinguish hotspots from singletons.
bimodal h3k3me <dig> peaks were specifically associated with hotspots, and mono-modal peaks with singletons.
this suggests that three events, binding by multiple factors, creation of a nucleosome-depleted region at the binding site and establishment of bivalent h3k4me <dig> marks go hand in hand.
we so far confirmed that hotspots had strong occupancy  and were enriched for cell type specific enhancers but their functional relevance is still a puzzle.
there are various theories such as it has been suggested that they might serve as sinks for tfs  <cit> .
in humans, as hotspots were depleted for tf motifs, hotspots were thought to be highly open chromatin marks where tfs bind non-specifically.
contrary to this, we find hotspots enriched for tf motifs especially for pioneering cell type specific factors such as pu.
<dig> in b cells.
pu.
<dig> binding has been shown to initiate nucleosome modeling followed by h3k <dig> mono-methylation  <cit> .
our hypothesis therefore is that pioneering factors establish stable open chromatin regions at hotspots by recruiting specific chromatin marks which in turn facilitates binding of additional tfs at their motifs or non-specifically.
finally, as hotspots are cell type-specific active enhancer regions, they need to be silenced in other cell types.
hotspot genes in es cells significantly overlapped with genes perturbed after setdb <dig> knockout and not dnmt knockout suggesting that setdb <dig> might play a role in es hotspot silencing.

 transcriptional control of gene expression via localisation and binding of transcription factors at the vicinity of gene loci lies at the heart of metazoan development. with the advances in sequencing protocols, chip sequencing is rapidly becoming the preferred tool to find genome-wide binding patterns of a transcription factor in a given cell type. as chip sequencing protocols mature, peaks of multiple transcription factors have been characterized in single cell types to study combinatorial control. one of the striking observations of these multi-tf studies was detection of so-called ‘transcriptional hotspots’. in moorman et al. generated genome-wide binding profiles of seven transcription factors in d. melanogaster and identified a subset of peaks bound by all seven tfs. of these hotspots, when tested using transgenic assays, 94% acted as enhancers strongly activating the neighbouring developmentally important genes. the presence of hotspot regions was further confirmed in other species such as worm and humans. hotspots in c. elegans and d. melanogaster were enriched for sequence binding motifs of many tfs including signature motifs such as gaga and zelda. however, analysis of human encode data did not support these observations. in humans, hotspots were specifically deprived of cis-regulatory motifs and no signature motif similar to the ‘gaga’ motif was found. at the other end of the spectrum, most tf binding events were ‘singletons’ accounting for nearly one third of the binding events in a cell type. these regions when studied in d. melanogaster transgenic assays did not drive patterned reporter gene expression leading to the conclusion that they do not act as strong developmental enhancers. this leads to a series of questions about properties of both hotspots and singletons. for example, are they truly distinct genomic regions, do they have a characteristic sequence or chromatin signature? in order to answer the above questions, we collected genome-wide binding patterns of multiple transcription-related factors in ten murine cell types. for each cell type, the peaks were classified into three groups: singleton genomic regions occupied by only one tf, hotspot genomic regions occupied by most tfs under study and combinatorial genomic regions occupied by a combination of tfs. the genomic regions largely clustered according to the group suggesting distinct genomic regions marked for occupancy independent of cell type. the singletons and combinatorials were bound neighbouring similar genes in all cell types while hotspot peaks occurred near a distinct set of genes in each cell type and showed functional enrichment for cell type specific genes. though hotspots were enriched for many tf sequence motifs, no signature motif such as gaga motif was found in murine cell types. we identified h3k4me chromatin modification distinguishing hotspots from singletons where hotspots showed a bimodal h3k4me peak whereas singletons were mono-modal. finally, we collected the differentially regulated genes after perturbations of multiple transcription-related factors in es cells, to show that hotspots were preferentially bound by activators and not repressors. as genes differentially expressed after setdb knockout but not dnmt knockout were enriched for es hotspot genes we suggest that setdb might be involved in silencing hotspots.  combinatorial binding events overlap across multiple cell types a chip sequencing experiment typically identifies thousands to hundreds of thousands of genome-wide binding sites of a tf in a cell type. in order to investigate if distinct characteristics of peaks grouped based on the number of binding factors at a binding location, we collected genome-wide binding events of transcription-related factors across different normal and cancer murine cell types. in each cell type, chip sequencing data sets for to tfs were collected. peaks were called in each sample with the number of peaks ranging from to in samples. the number of tfs at binding sites followed an exponential curve such that about half of the regions were bound by only one tf in any given cell type and less than % occupied by all tfs studied. we divided all peaks in a cell type into three groups: singletons or low-occupancy, peaks bound by only one transcription-related factor, combinatorials or mid-occupancy, peaks bound by a combination of transcription-related factors and hotspots or high-occupancy, peaks bound by more than five tfs studied in a given cell type. on average about 50% of all the peaks in a cell type were singletons and combinatorials while only -2% was classified as hotspots. we assigned peaks in these three groups to genomic locations: promoter, 5′ utr, 3′ utr, exon, intron or inter-genic using homer. the location of peaks followed a uni-modal distribution with a peak at the tss. singleton peaks were specifically under-represented in promoter regions and 5′ utr regions.figure combinatorial binding events overlap across multiple cell types. a. schematic diagram of a genomic region with genome tracks of six chip sequencing samples in a cell type marking singletons, combinatorials and hotspots. b. a table listing all cell types used along with the number of transcription factors c. bar graph of fraction of peaks in hotspots across cell types. d. heatmap of pearson’s correlation coefficients of peaks in sets with a hierarchical clustering tree showing that combinatorials clustered according to the groups. e. box plots for peak heights for the three groups for t cells and macrophages. in order to investigate the overlap of the three groups across cell types, we constructed a hierarchical tree of peak sets. to construct hierarchical tree, we first integrated of all peak sets across cell types resulted in unique peaks. we then constructed a binary matrix of rows and columns recording the presence or absence of peaks in each group for each cell type and built hierarchical tree using pearson’s correlation coefficient as a distance measure. many samples clustered according to the group across multiple cell types. the peaks in three groups in a given cell type are mutually exclusive by definition but the fact that combinatorials across many cell types cluster according to the groups is a notable observation. this suggests that each genomic location is marked for occupancy i.e. mid-occupied loci in one cell type were more likely to be bound by multiple factors rather than one factor in another cell type. hotspots of related cell types cluster into two tight clusters. the two embryonic stem cell samples were the most similar to each other as expected. the four related blood cell types – hpc, erythroid, mk progenitors and b cells formed second separate dense cluster of hotspot peaks. in d. melanogaster, peaks occupied at high levels by transcription factors in vivo drive patterned gene expression, whereas those occupied only at lower levels mostly do not. we calculated average peak height for the three groups where peak height was defined as the total number of tags in a 400 bp window around the peak summit. singletons consistently showed significantly lower peak height compared to combinatorials and hotspots. for example, the average peak height of fli bound singleton regions in t cells is significantly lower than combinatorials as were the average peak heights of cebpβ bound regions in macrophages where singletons were significantly lower than combinatorials. finally, we collected chip sequencing data for oct in es cells generated by three independent labs and found about oct peaks overlapping in all three experiments. the peaks were enriched for combinatorials and hotspots but were depleted for singleton regions. combinatorials and hotspot binding events were found consistently across different experiments. taken together, singletons were depleted in promoter and 5’ utr regions, had lower peak height and were not consistently found across multiple experiments. hotspots preferentially bind in cell type-specific gene neighbourhoods the hierarchical clustering of peak sets clustered them largely according to occupancy. to investigate if the group-wise clustering of peaks is also reflected at a gene level, we mapped all peaks to the nearest gene resulting into a total of genes. we created a binary gene matrix of rows and columns and built a hierarchical tree using pearson’s correlation coefficient as a distance measure. in contrast to the peak overlap, the gene overlap clustered all singleton and combinatorial peaks together forming one cluster while hotspot genes showed a lower overlap across cell types. by considering peaks only in promoters and performing a hierarchical clustering, the three groups clustered largely together and combinatorials formed one tight cluster. thus, combinatorials were located in the vicinity of similar genes across cell types but hotspots were present near a distinct set of genes in each sample. the functional enrichment of hotspot genes highlighted enrichment for cell type specific functions, such as b cell hotspots which were over-represented for the b cell receptor signalling pathway and b cell activation, while stem cell hotspots were enriched for stem cell differentiation and stem cell development. though binding events of all peaks were enriched for cell type specific functions, hotspot genes had a much higher cell type specific functional enrichment than all peaks in a cell type.figure hotspots preferentially bind in cell type-specific gene neighbourhoods. a. a heatmap of pearson’s correlation coefficients of genes in sets with a hierarchical clustering tree showing that singletons and combinatorials clustered together in a tight cluster while hotspot gene sets were very distinct from each other. b. fraction of es super enhancers defined by whyte et al. in three groups showing hotspots significantly overlapped with super enhancers. c. fraction of high oct bound regions in oct4+/− compared to oct4+/+ es cells in three groups showing hotspots enriched for high oct bound regions in oct4+/− and also pluripotent genes d. functional enrichments of hotspot genes in cell types along with p values showing that hotspot genes were enriched for cell type specific genes. recently, three groups generated a list of cell type specific enhancers in mouse es cells using diverse experiments. whyte et al. showed that master transcription factors such as oct, sox, and nanog form ‘super-enhancers’ at key cell identity genes that span large domains and drive cell-type-specific gene expression program. we calculated the overlap of es cell peaks in all three groups with super-enhancers and noted that about 30% of super-enhancers overlapped with hotspots in es cells. karwacki-neisius et al. showed that oct4+/− es cells have increased oct binding at pluripotency enhancers. only the enhancer regions strongly bound by oct in oct4+/− compared to oct4+/+ had a preferential bias towards hotspot regions. as oct binding in oct4+/− was shown to be specifically enriched near pluripotent genes, this confirmed that hotspots occupied cell type specific gene loci. stadler et al. generated a methylome map in murine es cells and noted that low methylated regions mark multi-tf bound enhancer regions in the vicinity of cell type specific genes. the overlap of high- and low-methylated regions with three groups in es cells confirmed that hotspots were specifically depleted for high-methylated and enriched for low-methylated regions in es cells. taken together, by comparing the three groups against three complementary experimental datasets we concluded that hotspots lie preferentially in the vicinity of cell type specific genes. hotspots are enriched for sequence motifs of multiple transcription factors and are preferentially bound by pioneering factors in order to search for properties distinguishing the three groups, we calculated the average mammalian conservation score for the genome sequences of peaks in all three groups. there was no statistically significant difference in the fraction of conserved peaks between the three groups. in c. elegans, hotspot-associated genes were more likely to be essential and were enriched for a variety of functions such as growth, reproduction, and larval and embryonic development. to check if this observation holds true in the mammalian system, we downloaded essential genes in mouse and investigated whether any of the three groups were preferentially bound near essential genes. though there was a relatively smaller fraction of essential genes in singletons compared to hotspots, the difference was not statistically significant. thus, genomic conservation and gene essentiality did not distinguish the three groups. in d. melanogaster, hotspots were enriched for sequence motifs of many but not all tfs under study and those without sequence motif enrichment were shown to bind via protein-protein interaction. they were enriched for sequence motifs of two global activators zelda and gaga and this sequence signature alone is sufficient to distinguish the hotspot regions. they were also enriched for beaf- and trl/gaf insulator motifs. in c. elegans, cis-regulatory sequences of tfs are not always present in hotspots, the binding was associated with open chromatin. in line with this, hotspots were enriched for ‘gaga’ sequence motifs. on the other hand, the depletion of cis-regulatory motifs was observed in human hotspots where the sequence motif for tfs in ‘cold’ regions were more than two- fold more enriched than in hotspot regions and no ‘gaga’ or similar motif was found, but they mentioned the possible presence of other sequence-specific motifs. indeed foley et al. found that on average, any given hot promoter contained about four different types of motifs. in summary there is no evidence so far regarding a specific cis-regulatory signature in hotspot peaks in mammals, nor is there a consensus whether hotspots are enriched or depleted for tf sequence motifs. to check if any sequence motif signature discriminates the three groups, we calculated the number of occurrences of mammalian cis-regulatory sequence motifs and compared them between the three groups. in most cell types, the motifs associated with the transcription factors known to be key players for that cell type were significantly enriched in hotspot regions compared to singletons. for example, gata and runx motifs were enriched in b cell hotspots, e-box and gfi motifs were enriched in erythroid hotspots while meis and scl-gata motifs were enriched in hpc hotspots. in contrast, there was a significant depletion of some motifs in hotspots compared to singletons. in particular, the ctcf motif was less abundant in hotspots across multiple cell types. in summary, hotspots were enriched for many known transcription factor motifs including the one not characterized by chip sequencing such as hox motif in hpcs. as gc content can bias the motif enrichment, we confirmed that the gc content across three groups was not significantly different across ten cell types. moreover, we validated that the motif enrichment was consistent only when the promoter peaks were considered in three groups.figure hotspots are enriched for sequence motifs of multiple transcription factors. a. bar graphs of fraction of peaks with a given sequence motif in the three groups with the name of the sequence motif along with the sequence in iupaac format for three representative motifs for each of b cells, erythroid and hpcs. b. a table of cell type, transcription-related factors analysed by chip sequencing, transcription-related factors binding to more than 90% of hotspot peaks and statistically significant known sequence motifs found in hotspots for b cells, erythroid and hpcs. siersbaek et al. identified ~ tf hotspots during adipogenesis and demonstrated that cebpβ marks a large number of these hotspots and is required for chromatin remodeling during differentiation. this hypothesizes that pioneering tfs in a cell type mark hotspots. we identified tfs preferentially binding in hotspots and found pu. marking b cell hotspots and cebps marking macrophage hotspots. med was bound at more than 99% of es hotspots which proposes hotspots to be under control of mediator complex. thus hotspots regions were indeed bound by pioneering factors across multiple cell types. on the other hand, only around 1% of peaks of suz overlapped with hotspots. as suz predominantly acts as a transcriptional repressor this suggests that hotspots might act as transcriptional hubs mainly for activation. altogether, murine hotspots did not contain a signature sequence motif but were enriched for motifs of many transcription factors including pioneering factors which mainly act as transcription activators. hotspots lie near highly expressed genes and are enriched for enhancers in nucleosome flanked regions in order to investigated whether each group shows a distinct chromatin signature, we collected data for three activating chromatin modifications in enhancers and promoters, ccctc-binding factor, rna polymerase ii binding as well as rna sequencing data in murine erythroleukemia cells from the encode resource. hotspot peaks showed significantly higher binding in all six datasets compared to the other two groups. in d. melanogaster, hotspot genes were expressed at higher levels during early embryogenesis but the difference is reduced later in development. in humans, using encode data, it was shown that genes associated with ‘hot’ regions were expressed at higher levels. in line with this, h3k4me polii binding as well as gene expression measured in rpkm values using rna sequencing data in mel cells were significantly higher in hotspots.figure hotspots lie near highly expressed genes and are enriched for enhancers in nucleosome flanked regions. a. boxplots for the peak height for h3k4me rna polymerase ii using chip sequencing as well as average rpkm value using rna sequencing data in mel cells for the three groups showing higher peak height for all chromatin modifications in hotspots. b. bar graph of fraction of vista enhancers and boxplots for the peak height for h3k4me h3k27ac using chip sequencing in the three groups in mel cells concluding hotspots were enriched for enhancers. c. the average density plot of mel h3k4me in the kb window for each of the three groups showing hotspots had a bi-modal h3k4me peak while singletons have a mono-modal peak. pennachio et al. predicted about high-confidence murine tissue-specific enhancers for tissue types by integrating tissue-specific expression data, conservation information and cis-regulatory motifs and further experimentally validated them using transgenic assays in mouse. mel hotspots were enriched in vista enhancers compared to the singletons and combinatorials. hotspots in other cell types were also enriched for vista enhancers. mel hotspots had significantly higher h3k4me as well as h3k27ac. histone h3k4me marks for enhancer regions whereas histone h3k27ac distinguishes active enhancers from inactive/poised enhancer elements containing h3k4me alone. hotspots therefore mark active enhancer regions in a given cell type. hoffman et al. showed high nucleosome positioning at the center of bimodal h3k4me peaks compared to mono-modal in mouse islets and liver. the average density plots of h3k4me showed a very strong bimodal distribution at hotspots while singletons showed a mono-modal distribution. this suggests that multiple transcription-related factors can indeed displace nucleosomes for strong enhancer activity at hotspot regions. hotspots were therefore enriched for nucleosome flanked multiple tf-bound loci. as expected, genes associated with flanked tf-bound loci were more abundantly expressed than those associated with nucleosomal loci, consistent with flanked sites being active enhancer elements. in order to validate these observations in another cell type, we collected genome-wide chromatin modification data sets in es cells. similar to mel cells, all activating chromatin modifications as well as expression values in rpkm using rna sequencing data showed significant enrichment in hotspots. the repressive chromatin modifications such as h3k27me were depleted in hotspots. taken together, we integrated genome wide chromatin binding patterns with tf chip sequencing data to show that hotspot regions were enriched for enhancer regions which are located near highly expressed cell type specific genes. furthermore, the h3k4me profile showed a strong bimodal peak at hotspots suggesting nucleosome displacement for transcription activation. es hotspots are preferentially occupied by activators and hotspots regions and might be under setdb control for silencing in order to establish the functional relevance of the binding in three groups, we collected a range of differentially expressed gene sets after tf perturbations in es cells. firstly we collected differentially expressed genes after oct knockdown from two independent experiments. es hotspot genes were specifically enriched for down-regulated genes after oct knockdown but not up-regulated genes. this suggests that oct binds to hotspot regions to act as a transcriptional activator.figure es hotspots are preferentially occupied by activators. a. overlap of genes in three groups in es ng dataset with up- and down-regulated genes upon oct knockdown in es cells showing hotspots were enriched for only down-regulated genes. b. fraction of up- and down-regulated genes upon dnmt and set1db knock out in three groups in es ng dataset showing hotspot enrichment only in set1db ko down-regulated genes. nishiyama et al. identified potential functional targets by identifying differentially expressed genes after overexpressing tfs in es cells. es hotspot genes were over-represented for differentially expressed regulators such as etv. however, es hotspot genes were under-represented in differentially expressed genes after overexpression of repressors such as atf, gadd45a, mybl, rhox and tcf. this corroborates the previous observation that hotspot regions were predominantly occupied by activators of transcription. not only transcription factors but also other chromatin modifiers such as baf, a subunit of the baf chromatin remodelling complex, as well as hdac, a histone deacetylase involved in the removal of acetyl groups from core histones, were preferentially bound at es hotspots. enrichment in hotspot loci suggests hdac might act as a transcriptional activator. though hdacs are generally considered repressors of transcription, it has been noted recently that hdac can indeed act as an activator. as hotspots mark highly active enhancers driving cell type specific gene expression, these regions need to be silenced in other cell types. karimi et al. showed dna methylation and setdb methylation work on different sets of genes, where only 7% of the genes up-regulated in a setdb knockout also were up-regulated in dnmt knock out in es cells. no significant difference was observed between singletons, combinatorial or hotspot genes for differentially expressed genes after dnmt knock out. challen et al. found ~ loci under dnmt3a control in hpcs. there was no bias in dnmt occupied loci in hpc, neither towards cpg islands in the three groups. however, genes differentially expressed after setdb knock out were enriched for hotspots suggesting that hotspots might be regulated by setdb altogether, hotspots might be used by activators to control expression of neighbouring genes as they were enriched for differentially expressed genes after perturbation of activators in es cells. furthermore, we suggested that es hotspots might be under setdb control for silencing.  the combinatorial action of multiple tfs controls regulation of gene expression where tfs activate or repress gene expression via binding to cis-regulatory sequence elements in gene loci. identification of functional target genes regulated by specific transcription factors is critical for understanding of the molecular mechanisms behind transcriptional control. multi-tf chip sequencing studies have led to a remarkable observation that there were regions in the genome bound by almost all of the factors, the so called ‘transcriptional hotspots’. these regions acted as active enhancers of developmentally essential genes. on the other hand, low occupancy regions did not drive reporter gene expression questioning the functional significance of these peaks. in order to study whether there are indeed classes of peaks with distinct characteristics based on their occupancy, we collected tf binding chip sequencing datasets in murine cell types. all peaks in a cell type were classified into three groups based on the number of factors bound: singletons, combinatorials and hotspots. singletons were depleted in promoters and 3’ utr, had a low peak height and were less conserved across multiple biological replicates putting in doubt their functional relevance. on the other hand, conserved cis-regulatory sequences were thought be more functionally relevant but the three groups could not be distinguished with respect to sequence conservation. indeed, it has been shown in yeast that the low affinity protein-dna interactions were conserved as well as being functionally relevant. the hierarchical clustering of peak sets clustered combinatorials according to occupancy. putative target gene sets of singletons and combinatorial peaks clustered into one tight cluster whereas each hotspot gene set consisted of unique set of genes with cell type specific functional gene enrichment. the presence of hotspot regions near cell type specific genes has previously been observed in many cell types across species. we have shown previously that most tf peaks cluster according to the cell type. gaga and zelda signature sequence motifs were associated with hotspots in c. elegans and d. melanogaster. we did not find a signature sequence motif for any of the groups but hotspots were enriched for sequence motifs of many key transcription factors as well as for pioneering factors in a cell type. hotspots were specifically depleted for ctcf sequence motif across all cell types. though depleted in ctcf sequence motifs compared to singletons, the fraction of ctcf peaks and their mean peak height at hotspots is higher. we further intersected three groups with diverse chromatin marks and concluded that hotspots were enriched for enhancers present in the vicinity of highly expressed genes. moreover, we identified a chromatin modification that can distinguish hotspots from singletons. bimodal h3k3me peaks were specifically associated with hotspots, and mono-modal peaks with singletons. this suggests that three events, binding by multiple factors, creation of a nucleosome-depleted region at the binding site and establishment of bivalent h3k4me marks go hand in hand. we so far confirmed that hotspots had strong occupancy and were enriched for cell type specific enhancers but their functional relevance is still a puzzle. there are various theories such as it has been suggested that they might serve as sinks for tfs. in humans, as hotspots were depleted for tf motifs, hotspots were thought to be highly open chromatin marks where tfs bind non-specifically. contrary to this, we find hotspots enriched for tf motifs especially for pioneering cell type specific factors such as pu. in b cells. pu. binding has been shown to initiate nucleosome modeling followed by h3k mono-methylation. our hypothesis therefore is that pioneering factors establish stable open chromatin regions at hotspots by recruiting specific chromatin marks which in turn facilitates binding of additional tfs at their motifs or non-specifically. finally, as hotspots are cell type-specific active enhancer regions, they need to be silenced in other cell types. hotspot genes in es cells significantly overlapped with genes perturbed after setdb knockout and not dnmt knockout suggesting that setdb might play a role in es hotspot silencing. 
they have been identified recently as cell type specific enhancers regulating developmentally essential genes in many species such as worm, fly and humans. transcriptional hotspots are enriched for tissue specific enhancers near cell type specific highly expressed genes. the genes neighbouring hotspots had a small overlap with hotspot genes in other cell types and were highly enriched for cell type specific function. in es cells, differentially expressed genes after perturbation of activators were enriched for hotspot genes suggesting hotspots primarily act as transcriptional activator hubs. in es cells, they are predicted to act as transcriptional activator hubs and might be under setdb control for silencing. the in-depth analysis of hotspots across multiple cell types in same species still remains to be explored and can bring new biological insights. we therefore collected transcription-related factor chip sequencing data sets in ten murine cell types and classified the peaks in each cell type in three groups according to binding occupancy as singletons, combinatorials and hotspots. transcriptional hotspots are defined as genomic regions bound by multiple factors. hotspots were enriched for sequence motifs of key tfs in that cell type and more than 90% of hotspots were occupied by pioneering factors. finally, we proposed that es hotspots might be under control of setdb and not dnmt for silencing. the peaks in the three groups clustered largely according to the occupancy, suggesting priming of genomic loci for mid occupancy irrespective of cell type. though we did not find any sequence signature in the three groups, the h3k4me binding profile had bimodal peaks at hotspots, distinguishing hotspots from mono-modal h3k4me singletons.
 transcriptional hotspots are defined as genomic regions bound by multiple factors. they have been identified recently as cell type specific enhancers regulating developmentally essential genes in many species such as worm, fly and humans. the in-depth analysis of hotspots across multiple cell types in same species still remains to be explored and can bring new biological insights.  we therefore collected transcription-related factor chip sequencing data sets in ten murine cell types and classified the peaks in each cell type in three groups according to binding occupancy as singletons, combinatorials and hotspots. the peaks in the three groups clustered largely according to the occupancy, suggesting priming of genomic loci for mid occupancy irrespective of cell type. we then characterized hotspots for diverse structural functional properties. the genes neighbouring hotspots had a small overlap with hotspot genes in other cell types and were highly enriched for cell type specific function. hotspots were enriched for sequence motifs of key tfs in that cell type and more than 90% of hotspots were occupied by pioneering factors. though we did not find any sequence signature in the three groups, the h3k4me binding profile had bimodal peaks at hotspots, distinguishing hotspots from mono-modal h3k4me singletons. in es cells, differentially expressed genes after perturbation of activators were enriched for hotspot genes suggesting hotspots primarily act as transcriptional activator hubs. finally, we proposed that es hotspots might be under control of setdb and not dnmt for silencing.  transcriptional hotspots are enriched for tissue specific enhancers near cell type specific highly expressed genes. in es cells, they are predicted to act as transcriptional activator hubs and might be under setdb control for silencing. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords chip sequencingtranscription hotspotstranscription regulationdata analysisissue-copyright-statement© the author 2014 