Dynamic regulation of epigenomic landscapes during hematopoiesis

Background Human blood develops from self-renewing hematopoietic stem cells to terminal lineages and necessitates regulator and effector gene expression changes; each cell type specifically expresses a subset of genes to carry out a specific function. Gene expression changes coincide with histone modification, histone variant deposition, and recruitment of transcription-related enzymes to specific genetic loci. Transcriptional regulation has been mostly studied using in vitro systems while epigenetic changes occurring during in vivo development remain poorly understood. Results By integrating previously published and novel global expression profiles from human CD34+/CD133+ hematopoietic stem and progenitor cells (HSPCs), in vivo differentiated human CD4+ T-cells and CD19+ B-cells, and in vitro differentiated CD36+ erythrocyte precursors, we identified hundreds of transcripts specifically expressed in each cell type. To relate concurrent epigenomic changes to expression, we examined genome-wide distributions of H3K4me1, H3K4me3, H3K27me1, H3K27me3, histone variant H2A.Z, ATP-dependent chromatin remodeler BRG1, and RNA Polymerase II in these cell types, as well as embryonic stem cells. These datasets revealed that numerous differentiation genes are primed for subsequent downstream expression by BRG1 and PolII binding in HSPCs, as well as the bivalent H3K4me3 and H3K27me3 modifications in the HSPCs prior to their expression in downstream, differentiated cell types; much HSPC bivalency is retained from embryonic stem cells. After differentiation, bivalency resolves to active chromatin configuration in the specific lineage, while it remains in parallel differentiated lineages. PolII and BRG1 are lost in closer lineages; bivalency resolves to silent monovalency in more distant lineages. Correlation of expression with epigenomic changes predicts tens of thousands of potential common and tissue-specific enhancers, which may contribute to expression patterns and differentiation pathways. Conclusions Several crucial lineage factors are bivalently prepared for their eventual expression or repression. Bivalency is not only resolved during differentiation but is also established in a step-wise manner in differentiated cell types. We note a progressive, specific silencing of alternate lineage genes in certain cell types coinciding with H3K27me3 enrichment, though expression silencing is maintained in its absence. Globally, the expression of type-specific genes across many cell types correlates strongly with their epigenetic profiles. These epigenomic data appear useful for further understanding mechanisms of differentiation and function of human blood lineages.


Background
As every cell type within one organism presumably contains the same genes, differential usage of the genome likely modulates functional differences between cell types. Human blood cells derive from self-renewing hematopoietic stem and progenitor cells (HSPCs) [1]. Differentiating HSPCs undergo progressive loss of differentiation capability in response to environmental conditions, resulting in terminally differentiated blood cells ( Figure 1A). Red blood cells derive from HSPCs via a common myeloid progenitor (CMP) and have a nucleated erythrocyte precursor (pRBC) [2][3][4][5][6]. T and B lymphocytes are related immune cell types, both arising through a common lymphoid progenitor (CLP) [7].
As many previous results were obtained using cell lines cultured in vitro, they may not reflect dynamic changes of histone modifications occurring during differentiation in vivo. Particularly, although it has been suggested that differentiation genes are bivalently modified in ES cells, it is unclear how bivalency resolves during differentiation in vivo. We therefore took advantage of the well-characterized human hematopoietic system and tested the relationship of tissue-specific gene expression to epigenomic changes, and identified common and tissue-specific regulatory elements that may confer tissue-specific transcriptional regulation. Our data indicate that key differentiation genes are primed in progenitor cells by bivalent modification and RNA PolII binding, some even in embryonic stem cells, prior to their expression in the downstream differentiated lineages. Bivalent modification is retained in closely related parallel cell types but is resolved to silent chromatin structure in more remote lineages. Using these epigenomic data, we have identified thousands of common and tissue-specific putative enhancers that might have critical roles in controlling the cell fate decisions during hematopoietic differentiation.

Results
Characterization of cell type-specific expression in HSPCs, erythrocyte precursors, B, and T-cells  Figure 1 Many transcripts show cell type-specific expression in hematopoietic subsets. (A) All blood cells derive from hematopoietic stem and progenitor cells (HSPCs) through common progenitors, including the common lymphoid and common myeloid progenitors (CLP, CMP). T and B-cells arise from CLPs, whereas red blood cell precursors (pRBCs) differentiate from CMPs. (B) Pairwise comparison of expression profiles from four cell types results in many differentially expressed transcripts but relatively few transcripts with cell type-specific expression. Total: the number of differentially expressed genes between two cell types. The numbers of cell type-specific genes are indicated below the panel.
(TS) gene expression therein. We calculated gene expression profiles using RNA-Seq [32,38] to identify pairwise differential expression ( Figure 1B). E.g., 794 genes were more highly expressed in HSPCs than in B-cells and 1,622 were more highly expressed in B-cells than in HSPCs. Since B-cells derive from HSPCs, this indicates that 794 and 1,622 genes are repressed and activated, respectively, during differentiation. To identify TS genes, i.e. those expressed in only one cell type, we intersected those highly expressed in each of the pair-wise comparisons ( Figure 1B). Indeed, their overall expression was significantly higher in their specific cell type (Additional file 1: Figure S1A-S1D). These gene lists included the markers used for isolating or defining the populations. E.g., CD133 and CD34 were HSPC-specific; CD4 was T-cell-specific; CD19 was B-cell -specific; and CD36 was pRBC-specific.
In addition to cell surface markers, key transcription regulators were also TS (Additional file 1: Table S1). Among these genes were PAX5, the master regulator of B-cell development; GATA3, an important regulator of T-cell differentiation; FOXP3, a regulator of regulatory T-cells; and GATA1, a master regulator of RBC development. That these genes of known type-specific importance were called type-specific shows that our tests extract high-confidence type-specific transcripts.
KEGG pathway analyses of TS genes showed associations with expected functions. The B-cell receptor signaling and IgA production pathways were enriched in B-cell-specific genes (Additional file 1: Figure S2C); T-cell-specific genes were similarly overrepresented in the T-cell receptor pathway, as well as in immunodeficiency and cell adhesion molecules (Additional file 1: Figure S2B). Although erythrocyte transcripts, including many versions of hemoglobin, were specifically expressed in pRBCs, these genes were enriched in KEGG pathways of more general biochemical function (Additional file 1: Figure S2D).

Chromatin environment at transcription start sites
To investigate relationships between gene expression and chromatin environment in progenitor and differentiated cells, we profiled genome-wide distributions of several histone modifications (H3K4me1, H3K4me3, H3K27me1, and H3K27me3), histone variant H2A.Z, chromatin remodeler BRG1, and RNA PolII using ChIP-Seq (see methods). We then plotted tag densities surrounding transcription start sites (TSSs) of HSPC-specific genes (Figure 2A -2E). As expected for expressed genes, PolII was highly enriched in these promoters in HSPCs ( Figure 2B) while no appreciable binding of PolII to these promoters occurred in pRBCs ( Figure 2C), B-cells ( Figure 2D) and T-cells ( Figure 2E). Similarly, BRG1 was enriched in these promoters only in HSPCs. H3K4me3 and H2A.Z were enriched in promoters of HSPC-specific genes not only in HSPCs but also in T-and B-cells, but and less so in pRBCs. While nominal promoter H3K27me3 signals were detected in HSPCs, they were highly enriched in the other three cell types, confirming that gaining H3K27me3 at these HSPC-specific genes during differentiation is associated with silencing (compare Figure 2B to 2C, D and E). To investigate priming of HSPCspecific genes, we investigated their histone modification status in embryonic stem cells (ESCs) (Figure 2A) using published datasets [39]. In ESCs, these promoters were enriched in H3K4me3, indicating preparation for their expression in HSPCs. However, the promoters were also marked with H3K27me3, which was lost at the HSPC stage, indicating further that their expression is associated with loss of H3K27me3.
To expand these observations to further differentiated cells, we profiled chromatin proteins at B-cell-specific genes ( Figure 2F-J). We note that, although these genes were not yet expressed in HSPCs, they had been associated with high levels of H3K4me3, H2A.Z and BRG1 ( Figure 2G), and with high levels of H3K4me3 in ESCs ( Figure 2F), suggesting that priming of these genes occurred before their expression in the differentiated cells. Binding of PolII was also detected on 17% of these genes in HSPCs, suggesting some were poised for expression. T-cells and B-cells descend from similar progenitors. B-cell-specific genes retained substantial signals of H3K4me3, H2A.Z and BRG1 in T-cells ( Figure 2J), suggesting that these genes may have not been fully silenced and may retain expression potential.
The above analyses of cell type-specific genes indicate that, while H3K4me3 signals are generally correlated with gene expression though are sometimes present even in non-expressing cells, changes in H3K27me3 are more closely related to inactive expression states in different cell types. To extend this observation, we sorted all genes into 200 groups according to their RPKM in HSPCs and averaged ChIP-Seq read densities of these two modifications across a 6kb region surrounding TSSs ( Figure 3A and 3D). Indeed, average profiles of H3K4me3 at genes highly expressed in HPSCs were consistent across all cell types, which may be related to the H3K4me3 signals at housekeeping genes, and levels of H3K4me3 positively associated with transcription ( Figure 3A). DNA sequence analysis indicated that promoters associated with constitutive H3K4me3 were enriched in CpG islands ( Figure 3B). Of 17,261 promoters consistently enriched in H3K4me3 in the four cell types, 93% also contained a UCSC-defined CpG island [40]; merely 37% of promoters with no or inconsistent H3K4me3 enrichment contained a CpG island. The correlation between H3K4me3 and CpG islands was observed previously [31], however the additional correlation with H3K4me3 consistency was not. We also noted that all cell types showed enrichment of H3K4me3 upstream and downstream of a depletion at the TSS, and that H3K4me3 extended further downstream of the TSS than upstream. This pattern was similar to PolII at the same TSS groups, sorted and calculated in the same manner ( Figure 3C). Despite fewer PolII-bound than H3K4me3-associated promoters, the majority of TSS groups showed some presence of PolII positively correlating with expression. H3K27me3 signals were mainly enriched in gene groups with low expression and occupied much of the interrogated region, contrasting with H3K4me3 and PolII ( Figure 3D).
The widespread signals of H3K27me3 suggest that it occupies broad domains, consistent with previous reports [32,41]. To examine the size of H3K27me3 domains during differentiation, we compared the sizes of islands across the four cell types. The distribution of H3K27me3 island sizes    peaked at around 1.5kb in ESCs, 4.2kb in HSPCs, and increased substantially in pRBC, T and B-cells ( Figure 3E), demonstrating that differentiation results in larger regions of H3K27me3 enrichment. We also found that H3K27me3 marks larger genomic regions in downstream cell types ( Figure 3F). To further investigate changes and growth of H3K27me3modified regions during differentiation, we performed a clustering analysis of regions enriched in H3K27me3 in these cell types ( Figure 3G). We evenly split each united region into ≤ 2kbp sections, creating 222,914 regions enriched in H3K27me3 in at least one cell type. We found that most of the H3K27me3-enriched regions changed during differentiation, underscoring the variability of H3K27me3 targeting. In fact, most regions enriched in H3K27me3 in downstream cell types were not enriched in HSPCs, and there was no cluster enriched in H3K27me3 in all five cell types. There were many regions that specifically gained H3K27me3 in either the T, B, or pRBC stages, but relatively few that became enriched in all three. After assigning regions to their nearest gene within 20kb, we noted that regions enriched in H3K27me3 in a single cell type coincided with genes with non-hematopoietic function (DATA NOT SHOWN). Genes nearest to H3K27me3-enriched regions in pRBCs alone were enriched in cytokine-cytokine receptor interaction genes (DATA NOT SHOWN), indicating repression of immunological pathways in this myeloid lineage. Several type-specific genes coincided with regions enriched in H3K27me3 in one cell type, as highlighted in Figure 3G, indicating silencing of signature gene expression. Curiously, 80% of genes associated with the cluster of regions enriched solely in HSPCs (upper leftmost cluster) remained silent in downstream cell types, even in the absence of H3K27me3. This indicates that, while H3K27me3 enrichment is indicative of repression, its absence is not necessarily indicative of activation as subsequent factors may maintain silencing and/or activation factors may not be present.
The dynamic regulation of H3K27me3 during differentiation is exemplified at the HOXA and HOXB developmental gene loci, which show initial enrichment followed by substantial loss then reconstitution at select genes, corresponding with their expression [31,32] (Additional file 1: Figure S3). This supports widespread H3K27me3 signals together with loss of active histone modifications stabilizing a silent chromatin conformation after differentiation.

Bivalent marking of promoters in HSPC and resolutions upon differentiation
The above indicates that genes specifically expressed in downstream cell types are associated with active chromatin marks, e.g. H3K4me3, in upstream cell types, although they are not expressed and may be previously marked with silencing H3K27me3. The coexistence of H3K4me3 and H3K27me3, termed bivalent modification, was discovered in ES and T-cells and proposed as a preparation for genes to be expressed in response to environmental cues [18,[31][32][33][34]. We sought to understand differentiation-coupled bivalency resolution by constructing a heatmap of promoter bivalency status ( Figure 4A). Most of the 5,345 promoters showing bivalency in any of our five cell types were bivalent in ESCs and most of these lost bivalency in downstream cell types. Stem cells had more bivalent genes than the more committed cell types, but several genes developed bivalency in downstream cell types.
To find examples of genes showing bivalency and resolution, we investigated the bivalent promoters in our cell types (Additional file 1: Figure S4). A minor fraction of these was bound by PolII. The Venn diagram shows different partitioning of bivalent promoters in various cell types (Additional file 1: Figure S4). Several important genes with TS function were bivalently marked in progenitor or parallel non-expressing cells. The zinc finger transcription factor GATA3 is essential for development of T-cells [42]. Our data indicate that, although the Gata3 promoter was associated with H3K27me3 in HSPCs and ESCs, consistent with its silent state, it also was marked by H3K4me3, suggesting that the gene is primed for expression in progenitor cells ( Figure 4C, Additional file 1: Figure S5). When expressed in T-cells, H3K27me3 disappeared from the GATA3 gene, accompanied by appearance of PolII. We note that both the H3K4me3 and H3K27me3 signals remained in B-cells ( Figure   Pol II  Similarly, the gene for paired box transcription factor PAX5, the master regulator of B-cells, was associated with H3K4me3 and PolII, but not H3K27me3 in B-cells ( Figure 4D, Additional file 1: Figure S5). Our data indicate that the PAX5 promoter is bivalent and bound by PolII in HSPCs, and bivalent in ESCs, suggesting that the gene has been prepared for expression ( Figure 4D in grey). Interestingly, PAX5 was still bivalent in T-cells, although the PolII signal disappeared (Figure 4D in green). In contrast, PAX5 contained only H3K27me3 in pRBCs, indicating that the gene is fully silenced therein ( Figure 4D in red).

RNA-Seq
In addition to the regulators discussed above, several other factors were bivalent in different combinations of the cell types (Additional file 1: Figure S4). CD8A, a defining marker of CD8+ T-cells, was bivalent in all except the CD4+ T-cells, where it was monovalently active, indicating a potential for plasticity among T-cell subsets. Another factor, TCF7, which is essential for T-cell development as well as regulation of the self-renewal/differentiation switch [43], was bivalent in HSPCs and ESCs, as was cardiac development regulator GATA5. In addition to cell typerelated genes, several chromatin genes were bivalent. Notably, all genes highlighted in Additional file 1: Figure S4, except the histone genes and HDAC8, were also bivalent in ESCs.
Globally, we found 2,811 promoters bivalent in HSPCs and calculated their bivalency status in each differentiated cell type ( Figure 4B). Curiously, 20-40% of these promoters remained bivalent in differentiated cells, some of which may be important for further differentiation as previously seen for T helper cells [44]. Generally, more HSPC bivalent promoters lost H3K27me3 than lost H3K4me3, except in pRBCs. We also investigated promoters that became bivalent in downstream cell types and analyzed their status in HSPCs. Strikingly few promoters had had neither mark in HSPCs. We noted that pRBCs, which differ from the other cell types as they were differentiated ex vivo from HSPCs [32], showed slight differences in ratios of bivalent HSPC promoter resolutions-many more lost H3K4me3 than remained bivalent, and fewer newly bivalent promoters had had H3K27me3 enrichment in HSPCs. This may reflect the degree of commitment of the cell type, as B and T-cells differentiate further. Although many bivalent promoters lost H3K27me3 enrichment in downstream cell types, most did not gain appreciable expression as a result (DATA NOT SHOWN). This further indicates a role for other factors maintaining silencing in the absence of H3K27me3. Many more genes were bivalently marked in ESCs than HSPCs. Nearly 40% of promoters bivalent in HSPCs were also bivalent in ESCs. Surprisingly few promoters bivalent in HSPCs were marked exclusively by H3K27me3 in ESCs, further indicating that genes are primed for expression notably early in development.
We also examined bivalent priming of TS genes ( Figure 4E). These results indicate that many genes expressed in downstream cell types are primed for expression in progenitor cells. Several TS genes with transcription factor activity were bivalent in HSPCs (Additional file 1: Figure S6), underscoring the importance of bivalency in controlling expression of lineage regulators. However, priming of these genes with PolII was not extensive in HSPCs. Several of these genes had multiple isoforms, which showed differences in bivalency at their discrete promoters, while all isoforms were typespecific by their expression. This suggests a role for chromatin in mediating isoform-specific expression.

Prediction and analysis of distal regulatory elements from chromatin profiles
Distal regulatory elements play critical roles in lineage fate decisions and may contain cell type-specific chromatin modification patterns [29,45,46]. Chromatin modifications and binding of chromatin-modifying enzymes indicate the presence of potential regulatory elements [47,48]. Enrichment of H3K27ac and binding of p300 are detected at active enhancers [45,46,49]. In order to completely survey both active and poised enhancers, we used enrichment of H2A.Z and/or H3K4me1 to predict genomic regions outside of a promoter (see methods) as potential enhancers, since previous studies reported that enhancers are associated with H3K4me1 [50] and the histone variant H2A.Z [27,51]. Among enriched sites in each cell type ( Figure 5A second column), we identified 11,283 regions enriched in these marks across all four cell types (core potential enhancers or CPEs). We also found regions enriched in these marks in only one cell type (cell-specific potential enhancers or CSPEs) ( Figure 5A third column). Among the regions enriched in the marks were several from the hemoglobin beta chain gene locus, which have been reported as functional enhancers (Additional file 1: Table 2). Other known enhancers predicted by our method included the first intron of the RUNX1 gene, which contains a functional enhancer, a PAX5 enhancer, and the PLAT enhancer. These results indicate that our method is successful in finding regions with enhancer activity.
To examine chromatin modifications at these predicted enhancers, we plotted ChIP-Seq reads at the 11,283 CPEs ( Figure 5B-E). CPEs were consistently enriched in our defining marks (H2A.Z and H3K4me1) and active marks (H3K4me3, PolII, BRG1), but were depleted of H3K27me3. In contrast to CPEs, HSPC CSPEs exhibited enhancer-like chromatin only in HSPCs ( Figure 6A). In other cell types, HSPC CSPEs lost active H2A.Z and H3K4me1 marks, and gained H3K27me3. Curiously, H3K27me1, whose function is poorly understood, was present at HSPC CSPEs in all cell types ( Figure 6B, C, D).
Similarly, B-cell CSPEs showed elevated active marks in B-cells ( Figure 6G), and lost active marks but gained the repressive H3K27me3 mark in other cell types ( Figure 6E, F, H). We note that B-cell CSPEs were already associated with elevated levels of BRG1 binding in HSPCs ( Figure 6E), suggesting that chromatin remodeling by BRG1 may be required for subsequent establishment of B-cell-specific enhancers, consistent with previous observations in erythrocyte differentiation [8]. The presence of poorly understood H3K27me1 is puzzling, as previous analyses have shown that H3K27me1 in gene bodies positively correlates with expression [16], and that it may show some enrichment in enhancers [52]. That these TS elements contain repressive marks in the other cell types is logical, since, combined with the lack of active marks, this could result in silencing of target genes.
To investigate how predicted enhancers correlate with the expression of nearby genes, we compared the counts of CSPEs associated with TS genes and versus CSPEs associated with all genes with RPKM greater than 0 ( Figure 7A). We found that type-specific expression tended to associate with more CSPEs, indicating that these CSPEs tend to be expression activators.
To associate histone marks at CSPEs with expression, we compared the distributions of pRBC expression of genes associated with pRBC-specific CSPEs enriched in different marks ( Figure 7B). Genes associated with a CSPE were more highly expressed than all genes (p < 2.2 × 10 -16 ), further indicating that these elements positively correlate with expression. Genes associated with a PolII-bound CSPE tended to show higher expression (vs. all genes p < 2.2 × 10 -16 , vs. CSPE-containing genes p < 2.2 × 10 -16 ) ( Figure 7B). PolII at these potential enhancers may result from direct binding of PolII to enhancers or indirectly from enhancer-promoter interactions. Several histone modifications at CSPEs were tested for their impact on expression of their nearest genes (for full discussion, see Additional file 1), showing strong evidence that chromatin environments of these elements correlates strongly with expression.
Here we compared expression and epigenomic differences among human hematopoietic progenitors and  downstream lineages. Our results indicate that genes activated in downstream lineages are marked by bivalent histone modifications and some are also bound by RNA PolII in progenitor cells. Bivalent modifications resolve in specific lineages but may remain in parallel differentiated lineages, depending on the relative distance of the cell types. We also predict enhancer-like elements that may contribute to the observed tissuespecific expression. Co-existence of active and inactive chromatin modifications at promoters has been detected in various cell types [18,[31][32][33]54]. Genes with such bivalency preferentially show tissue-specific expression in downstream cell types [55]. These promoters generally resolve to monovalency during differentiation [31]; most retain H3K4me3 in their specific cell type and retain H3K27me3 in alternate cell fates [31,55]. Although fewer bivalent genes have been observed in differentiated cell types than in stem cells, it is unclear how bivalency resolves during differentiation of progenitor cells to their direct downstream lineages in vivo. Our analysis of hematopoietic stem and progenitor cells (HSPCs) and in vitro differentiated erythrocyte precursors (pRBCs), in vivo differentiated CD4+ T-cells, and CD19+ B-cells in the blood compartment indicated that many genes critical for HSPC differentiation are primed by bivalent modification and may also be bound by PolII in HSPCs prior to differentiation. E.g., the T-cell regulator GATA3 and the B-cell master regulator PAX5 were both bivalent in HSPCs; they became uniquely expressed and resolved to monovalently active in T and B-cells, respectively. PolII bound PAX5 in HSPCs. Bivalency is not limited to promoters; we noted several enhancer-like regions showing the presence of both active and repressive marks in specific cell types (DATA NOT SHOWN). As suggested by previous publications, bivalency primes these critical genes for expression during differentiation [18,31]. Obviously, bivalent genes are heterogeneous in both histone modification enrichments and expression behavior. They can be separated into different groups based on the ratio between H3K4me3 and H3K27me3 enrichment, and the enormous variation in this ratio translates into important differences in gene expression behavior [18,56,57]. Curiously, GATA3 remained bivalent in B-cells and PAX5 remained bivalent in T-cells, indicating that these closely related lymphocyte types might possess some potential for plasticity between them. This is consistent with a previous report that bivalent modification of key regulators in T helper cells are linked to their plasticity [58]. Although PAX5 is bivalent in both the progenitor HSPCs and T-cells, the chromatin remodeling factor BRG1 and PolII are associated with the promoter in HSPCs but not in T-cells, suggesting that changes in chromatin and expression potential have occurred after differentiation to T-cells. In contrast, both the GATA3 and PAX5 genes lost H3K4me3 in pRBCs, silencing them in this distantly related cell type, indicating little likelihood of direct transdifferentiation between these lymphoid and myeloid cell types. Indeed, conversion of Pax5-deleted B-cells into T-cells has been performed, albeit through de-differentiation into progenitor cell types [59]. Factors implicated in chromatin structure and modification, C/EBPα and C/EBPβ, can transdifferentiate B-cells into macrophage-like cells [60,61]. GATA3 itself was shown to mediate histone modifications at great distances from its binding sites in mouse [44]. Together, these results indicate substantial opportunities for epigeneticmediated transdifferentiation of blood lineages, which is furthered by our analysis of incomplete silencing of TS factors in parallel cell types Establishment of novel bivalency in partially differentiated progenitors is an underaddressed trait of development. Several genes that were not bivalent in HSPCs became bivalent downstream. Strikingly few genes had neither mark in HSPCs and subsequently become bivalent. Bivalency seems to thus be both established and resolved in a stepwise manner, as relatively few genes lose both active and repressive marks between differentiation stages ( Figure 4B). For example, while H3K27me3  at the HOX loci changes drastically from ESC to HSPC and beyond (Additional file 1: Figure S3), most of the genes in these clusters retained H3K4me3 in the promoter regions (DATA NOT SHOWN). This may be explained as H3K27me3 targeting genes for selective repression initially, while H3K4me3 changes are secondary or occur on final commitment. In fact, we see that many more changes in H3K4me3 occur between HSPCs and the strongly committed pRBCs, while most other notable changes between cell types are H3K27me3-related. The stability of promoter H3K4me3 may be related to the presence of CpG islands within the promoter. This establishment of bivalency in non stem-like populations may require further attention to fully investigate its importance. We acknowledge that there are other silencing mechanisms such as H3K9me3 or DNA methylation in the cells, which may or may not overlap with the H3K27me3 pathway. More analyses of these different chromatin modifications are required for a more complete understanding of gene priming, repression and activation during differentiation of these cells. ChIP-Seq is performed on a population of cells, and homogeneity of this population is crucial to analyses that describe subsets [16]. Although the cells used in this study were defined and purity-assessed by the cell surface markers, bivalent modification detection may have resulted from cellular heterogeneity [62]. We have previously shown [18] via sequential ChIP that both H3K4me3 and H3K27me3 can exist at the same promoter, but whether or not these marks exist on the same histone, or are functional in a single-nucleosome context is still unknown.

Conclusion
In this study we identified tissue-specific genes in the human blood compartment and compared gene expression status with chromatin modification patterns. Based on the gene expression and chromatin modification patterns, we have predicted tens of thousands of potential core regulatory elements shared by all cell types and potential tissuespecific regulatory elements and show that a combination of both H3K4me1 and H2A.Z is a better predictor of enhancer activity than either alone. We have shown that bivalent chromatin modification in the well-characterized human hematopoietic system is not only resolved but also established during differentiation in a generally stepwise manner. While much of H3K4me3 enrichment is stable across cell types, H3K27me3 varies more frequently and grows to cover more of the genome during differentiation. These epigenomic data and identification of potential regulatory elements will be useful for further understanding the mechanisms governing decision making and differentiation of hematopoietic cell types.

Data sources
Data for CD34+/CD133+ hematopoietic stem and progenitor cells, CD36+ red blood cell precursors, and CD4+ T cells were downloaded from the NCBI Short Read Archive. A full summary of data sources is available as Additional file 1: Table 3. Counts of uniquely mapped reads in library, non-redundant reads, and reads falling in statistically enriched islands are Additional file 1: Table 4. New data can be downloaded via NCBI GEO using accession GSE39229.
Naïve B cells were purified from human blood using human Naïve B cell isolation kit II kits (Miltenyi, #130-091-150). The cells were digested with MNase to generate mainly mononucleosomes with minor fraction of dinucleosomes for histone modification mapping. For mapping enzyme target sites, the cells were crosslinked with formaldehyde treatment and chromatin fragmented to 200 to 500 bp by sonication. Chromatin from 5 × 10 6 cells was used for each ChIP experiment.

ChIP-Seq and analysis
ChIP-Seq was performed as previously described [16]. Illumina reads were mapped to the hg18 genome using Bowtie [63], allowing only one position per read (−m 1), and filtered to allow only one read per position. Y chromosome reads were disregarded.
We used MACS version 1.4.0 RC 2 [64] with input control libraries from corresponding cell types, a p-value threshold of 1 × 10 -8 , and -g hs to detect highconfidence Pol II-binding sites. Genes were considered Pol II-bound if they contained a peak between 5 kbp upstream of the TSS to 3 kbp downstream of the TES.
Regions enriched in Brg1 or modified or variant histones were detected using SICER version 1.03 [65]. The percent of hg18 uniquely mappable by 25 bp reads (68%) was retrieved from [66]. We used a window size of 200 for Brg1, H2A.Z, H3K4me1, and H3K4me3. Window sizes for H3K27me1 and H3K27me3 were predicted using an unpublished version of SICER [IN PREPARATION vs In Preparation]. We used gap sizes of 0 windows for H3K27me1 and H3K27me3, 1 window for H3K4me1 and H3K4me3, 2 windows for H2A.Z, and 3 windows for Brg1. Fragment size was 150, and the FDR cutoff for statistical enrichment was 1 × 10 -5 . Due to the potentially large size of windows, we nibbled the edges of outermost windows of each island to maximize the difference in read density between the remaining portion and the removed portion of the outermost windows. This feature will be implemented in the next version of SICER [IN PREPARATION vs In Preparation]. UCSC browser tracks were created using 200 bp window sizes, a +/− 75 bp shift by strand, and were normalized to the 10 7 reads in the library.
Reads used for alignments with respect to TSSs and enhancers were required to come from statistically enriched islands. For enhancers, reads were shifted +/− 75 stranddependent bps, sorted by their starts into 50 equally sized bins, counted, and normalized by bin size, TSS number, and sequencing read count. For TSS profiles, reads were aligned relative to RefSeq TSSs separated by strand, counted in 10 bp bins, inverted for negative strand genes, summed, smoothed over 4 windows on either side, and normalized by sequencing read count. Both sets of profiles were plotted using GNUPLOT [67]. We used MeV [68,69] to display heatmaps in Figure 3 A, C, and D, without correction of zeroes. Figure 3B area-proportional Venn diagrams were created using the VennDiagram package in R [70]. CpG islands were downloaded from the UCSC Genome Browser [71] Figure 3F shows a heatmap of H3K27me3 reads. Islands of H3K27me3 from all five cell types were united, and fragmented equally into ≤ 2 kbp fragments. Read counts of H3K27me3 in these fragments were normalized by sequencing library size, clustered using k means (k = 20) clustering, and displayed in a heatmap using R, sorted by cluster sum of H3K27me3 reads [70].
Promoters were considered to be bivalent if they had a region statistically enriched in both H3K4me3 and H3K27me3 within 500 bp of the TSS. Enhancers were considered bivalent if they had any H3K27me3 enrichment in their locus. Type-specific transcription factors bivalently prepared in HSPCs were predicted using GeneCards version 3 [72] to extract Gene Ontology terms [73], and required to have the term "transcription factor." We predicted enhancers by taking the union of H2A.Zand H3K4me1-enriched islands using SICER. Enhancers, and H3K27me3 regions in Figure 3F were associated with transcripts if they were located between 20 kbp upstream of the TSS and 0.5 kbp upstream of the TSS, or between TES and 20 kbp downstream of the TES. P-values in Figure 7A were calculated using a two-tailed Kolmogorof-Smirnoff test in R version 2.6.2 [70]. The highest 5% of RPKM values were removed from distribution analysis in Figure 7B.

RNA-Seq and analysis
RNA-Seq was performed as previously described [8]. Illumina reads were aligned to the hg18 genome using TopHat [74] with standard parameters. Mapped reads were then converted to BED format using SAMtools [75]. UCSC browser tracks [71] were created from BED reads using no shift and a window size of 20bp. We calculated RPKM values for all RefSeq transcripts not on the Y chromosome using a previously described method [8]. Y chromsome reads were disregarded as some subjects were female. Pairwise differential expression was calculated using EdgeR [76] with log(fold-change) ≥ 5 and FDR ≤ 1 × 10 -5 thresholds. We found type-specific genes by taking the intersection of the three lists of significantly more highly expressed genes from the three pairwise comparisons per cell type. KEGG pathway analysis was performed using DAVID [77] KEGG pathway enrichment with standard settings.

Additional file
Additional file 1. Methods. Figure S1A-S1D: Expression of typespecific genes in four cell typesrelated to Figure 1. Figure S2A-S2D. Type-specific genes are enriched in functional pathwaysrelated to Figure 1. Figure S3A-S3B. Enrichment of H3K27me3 at HOXA and HOXB locirelated to Figure 3. Figure S4. Bivalency of promoters across multiple cell types with examplesrelated to Figure 3. Figure S5A-B. Density of chromatin proteins at Gata3 and Pax5 genesrelated to Figure 4. Figure S6. Resolution and preparation of transcription factor genes bivalent in HSPCrelated to Figure 3. Table S1. Excel document of lists of type-specific genes. Table S2. Known enhancers and their enrichments in H3K4me1 and/or H2A.Z in four cell types. Table S3. Files used in analysis and their sources. Table S4. Mapped sequencing reads, unique reads, and unique reads in enriched islands. Discussion of Figure 7B.