Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

A hyper-dynamic nature of bivalent promoter states underlies coordinated developmental gene expression modules

BMC Genomics volume 15, Article number: 1186 (2014) Cite this article

Abstract

Background

Chromatin remodeling is crucial for proper programing of developmental gene expression. Recent work provides a dynamic view of post-translational histone modifications during differentiation; however there is little insight on the evolution of combinatorial genome-wide patterns of chromatin marks, excluding an essential aspect of developmental gene regulation.

Results

We report here a 15-chromatin state Hidden Markov Model which describes changes in chromatin signatures in relation to transcription profiles during differentiation of human pre-adipocytes into adipocytes. We identify nineteen modules of gene expression reflecting multiple waves of transcriptional up- and down-regulation which characterize adipogenic differentiation. From our model, we developed chromatin state matrices fitting each of these transcription modules to show how the complexity and dynamic nature of chromatin signatures relate to expression patterns. Spatial relationships between chromatin states underlie a high-order chromatin organization in differentiating adipocytes. We show the importance of gene expression level in generating diversity in chromatin signatures, and show that the hyper-dynamic nature of H3K4me2/H3K27me3-marked ‘bivalent’ promoter states underlies many of the gene expression patterns associated with adipogenic differentiation.

Conclusions

Our results reveal the highly dynamic nature of bivalent promoter states within the adipogenic lineage. The data constitute a valuable resource enabling the assessment of possibilities to alter the adipogenic program.

Background

Developmental gene expression entails waves of coordinated transcriptional activation and repression events, which are orchestrated at least in part by the epigenome [1–5], a layer of reversible post-translational modifications on chromatin [6]. These studies invariably indicate that the epigenome is dynamic and that epigenetic modifications are linked to changes in gene expression in a concordant or sometimes seemingly non-concordant manner. However, some features of the epigenome appear to be more static than others: subsets of histone post-translational modifications (hPTMs) form combinatorial associations that are relatively stable and hence can be used to functionally annotate genomic elements [7, 8].

Several stem or progenitor cell-based differentiation models have been used to assess the temporal changes in epigenetic states across the genome [9–11]. A system which has remained under-studied relative to its societal importance is the differentiation of pre-adipocytes into adipocytes. Regulation of the balance between maintenance of a pool of adipocyte progenitors and their differentiation into adipocytes is essential for adipose tissue homeostasis [12, 13]. Adipogenic differentiation is driven by activation of genes encoding transcription factors that synergistically up-regulate target genes involved in adipocyte formation and lipid metabolism [14–16], and by chromatin remodeling notably at regulatory elements essential for transcription factor binding [17]. Genome-wide profiling of hPTMs by chromatin immunoprecipitation and high-throughput sequencing (ChIP-seq) during mouse and human in vitro adipogenic differentiation has led to the identification of thousands of putative adipogenic-specific promoters and enhancers [18]. The current data portray a dynamic view of enrichment in individual hPTMs [17, 18]; however insights on the temporal transitions in combinatorial associations of hPTMs during adipogenic differentiation have been missing, leaving out a key aspect of developmental gene regulation.

A deeper appreciation of the complexity of chromatin signatures can be obtained through Bayesian methods which model combinatorial associations of chromatin marks [19, 20]. Among these, Hidden Markov Modeling (HMM) uses machine learning to discover chromatin states from recurrent combinations of histone modifications, transcription factors and chromatin remodeling factors [19, 21]. From the analysis of a panel of factors in several unrelated cell types [1, 21–25], HMM provides the ability to distinguish functional genomic elements, generating genome-wide profiles of chromatin ‘activity’ [19, 21]. However, the experimental material used in these studies was not chosen to infer a temporal dynamics of chromatin states; thus developmental transitions in chromatin states in a genome-wide context have not been fully explored.

To palliate this gap, a variation of HMM has been proposed as an unsupervised hierarchical model enabling correlations between gene expression patterns and clusters of combinatorial chromatin marks [9]. This model supports the tissue- and cell type-specificity of enhancer activity and associated histone modifications [10, 21, 26], and concurs with recent evidence for distinct mechanisms of gene expression regulation along the genome [27].

Here, we applied ChromHMM, a high-throughput pipeline based on a multivariate HMM [19], using ChIP-seq data for seven chromatin marks [18], in combination with RNA-sequencing (RNA-seq), to discover interrelationships between chromatin states and gene expression patterns during differentiation of human pre-adipocytes. We identify several coordinated gene expression modules and learned a 15-state model which we use to map and quantify temporal transitions in chromatin states across these expression modules. We reveal the importance of gene expression level in generating diversity in chromatin signatures and the hyper-dynamic nature of ‘bivalent’ promoter states during lineage-specific differentiation.

Results

Coordinated gene expression modules characterize adipogenic differentiation

We used deep paired-end RNA-seq to unveil transcriptomic changes at four time points of adipogenic differentiation of human preadipocytes. Numbers of reads per time point ranged from 40 to 102 million, of which numbers of paired alignments were >30 to > 95 million (Additional file 1: Table S1). We analyzed cells two days before adipogenic induction (D-2; undifferentiated proliferating ASCs), immediately prior to adipogenic induction (D0; undifferentiated confluent cells, 48 h after growth factor removal) and on D3 and D9 of differentiation. Principal component analysis shows that the first two components segregate the transcriptome of proliferating and confluent ASCs from that of adipogenic-stimulated cells (Figure 1A), reflecting a major transcriptional switch. This conversion is manifested by the greatest differences both in the numbers of differentially expressed genes (2512 and 2910 up- and down-regulated genes respectively; Additional file 1: Figure S1A) and in the magnitude of differential expression levels (Figure 1B; see Additional file 2: Table S2 for lists of differentially expressed genes). Culturing ASCs to confluency remarkably reduces the overall variability in transcript levels detected in proliferating cells (Figure 1C), consistent with greater individual variations in gene expression patterns in unsynchronized cell populations [28]. This reduced variability is maintained after adipogenic induction (Figure 1C), consistent with the establishment of a coordinated gene expression program. Of note, this does not preclude potentially persistent cell-to-cell variations in transcriptional response to the differentiation stimulus in the populations analyzed [29].

Figure 1
figure 1

Coordinated gene expression modules during adipogenic differentiation. (A) Principal component analysis of differentially expressed genes. (B) Fold change in expression level of each gene plotted against its P-value (expressed in -log10 FPKM). Significance: fold change > 2; α < 0.05. (C) Square coefficient of variance (CV2) of expression as a function of expression level (log FPKM) at each differentiation time point. (D) RNA-seq profiles of inducers and repressors of PPARG expression. Scales (FPKM): KLF15, 0-25; KLF5, 0-10; CEBPA, CEBPB, SREBF1, KLF2: 0-50, GATA2: 0-50, PPARG, 0-250. (E) Hierarchical clustering of genes into 19 expression modules throughout differentiation. Each line represents one transcript. Numbers of genes in each cluster are shown. FPKM scales are constant for each cluster.

Full size image

The transition from cell proliferation to confluency (D-2 to D0) is characterized by up-regulation of genes involved in extracellular matrix reorganization, regulation of proliferation, development and signaling, and down-regulation of genes important for cytoskeletal reorganization and signaling functions (P < 10−6 – 10−12; Additional file 2: Table S3), in line with the acquisition of a cell cycle arrest phenotype. Adipogenic stimulation (D0-D3 transition) up-regulates metabolic genes key for adipocyte development including lipid synthesis, metabolism and homeostasis (P < 10−8 – 10−28; Additional file 2: Tables S2 and S3). We importantly confirm expression and up-regulation of positive regulators of PPARγ (peroxisome proliferator-activated receptor γ), a key regulator of adipogenesis [15] (Figure 1D), including Krüppel-like factors (KLF4, KLF5, KLF6, KLF9 and KLF15), STAT proteins (STAT5A, STAT5B), sterol-response element binding protein 1c (SREBP1c; SREBF1) and CCAAT/enhancer-binding proteins (C/EBPs; CEBPA, CEBPB, CEBPD, CEBPE). Conversely, negative regulators of PPARγ (KLF2 and GATA2) are down-regulated (Figure 1D; Additional file 2: Table S2). Our RNA-seq data thus identify massive transcriptional changes leading to the activation of key metabolic genes required for adipocyte formation and function.

To provide a dynamic assessment of transcriptional changes during differentiation, we identified by hierarchical clustering 19 cohorts of genes displaying unique expression profiles throughout differentiation (Figure 1E). We identify cohorts of genes with stable expression levels (clusters 1-4) which regroup 81% of all expressed genes. The remaining genes partition into 15 clusters showing sequential transcriptional induction (that is, up-regulation from zero FPKM) on D0, D3 or D9 (clusters 5-7), up-regulation of already expressed genes (cluster 8), sequential transcriptional down-regulation (to zero FPKM; cluster 9-11), and transient up- or down-regulation (cluster 12-19). Genes sequentially activated or repressed (clusters 5-7, 9-11) or transiently activated or inactivated (clusters 12-19) are involved in signaling and transcription regulation processes important for adipogenic and lipid metabolism functions encoded by cluster 8 (Additional file 1: Figure S1B; Additional file 2: Tables S4 and S5). Our RNA-seq and clustering data reveal therefore the establishment of coherent gene expression modules characterizing the adipogenic differentiation program. Our datasets also constitute a high-depth transcriptome resource mapping the adipogenic process in human primary preadipocytes.

Spatial relationships between chromatin states reveal restricted state transition choices along the genome

We sought to identify a relationship between temporal gene expression changes and enrichment in combinations of chromatin marks during adipogenic differentiation. We used ChromHMM [19] to learn a 15 chromatin state (‘cs’) model from recurrent combinations, in consecutive 200-base pair (bp) bins, of seven chromatin marks (H3K27me3, H3K27ac, H3K4me1, H3K4me2, H3K4me3, H3K36me3 and the CCCTC-binding factor CTCF) profiled by ChIP-seq [18] in preadipocytes on D-2, 0, 3 and 9 of differentiation (Figure 2A). Our model balances information level, interpretability and resolution. We annotated the 15 states into functional genomic elements including active and inactive enhancers, enhancer location (near transcription start sites, within promoters or within gene bodies), active, inactive and ‘bivalent’ promoters and transcribed gene bodies (Figure 2B). ChromHMM also generates transition parameters based on spatial relationships between adjacent genomic segments, representing the sequence of states across the genome (Figure 2C). These parameters show that one state is most commonly followed by the same state or by one other state rather than by many states (Figure 2C). For instance, cs1 (H3K27ac/H3K4me1, annotated as active enhancer sites) and cs2 (H3K27ac/H3K4me1/H3K4me2; active enhancers in promoter regions), or cs3 (H3K4me1, H3K4me2) and cs4 (H3K4me1), frequently follow each other, suggesting embedding of active or inactive enhancers within promoter regions. State 4 can also be followed by ‘blank’ state 7, revealing enhancer sites also within chromatin deserts (Figure 2C).

Figure 2
figure 2

Chromatin states learned from recurrent combinations of seven chromatin marks in proliferating and adipogenic-stimulated preadipocytes. (A) Chromatin modifications and resulting chromatin states (color coded) in a 6 kb region of chromosome 1. (B) ChromHMM emission parameters and functional element annotation from the 15 states. Scale reflects enrichment of a given state in indicated chromatin marks. (C) Transition parameters across the genome. Scale reflects the frequency of occurrence for a give state. (D) Proportions of the genome covered by the 15 chromatin states learned in our model.

Full size image

We find that genome coverage and length of states vary considerably, from 200 bp (a bin size) to tens of kilobases (Additional file 1: Figure S2A). This agrees with the genomic nature of histone modifications defining these states, which mark spatially restricted elements (e.g. promoters or enhancers; cs1) or wider areas such as H3K27me3 domains (cs8) or entire H3K36me3-marked coding regions (cs6) (Figure 2A,D). Our model also enables to estimate that ~16% of the genome contains at least one of the marks examined (Figure 2D). This is in line with recent data from man [21], mouse [9] and Drosophila[22, 30], which based on the marks examined and analysis methods concur in that most of the genome is in a ‘blank’ state. Logically, state 7 is the broadest (Figure 2D; Additional file 1: Figure S2A) and essentially deprived of genes (4.3×10−4 genes per megabase; data not shown). The spatial relationship between chromatin states learned in our HMM demonstrates a limited choice in the the transition from one state to another along the genome, which reflects a higher order chromatin layout. The data also reveal the existence of chromatin state ‘variants’ learned from non-canonical genomic elements, such as canonical enhancer sites (marked by H3K4me1 with or without H3K27ac) within promoter elements or active gene bodies.

Analyzing the load and temporal dynamics of chromatin states

We next established relationships between chromatin state enrichment level or state dynamics in the course of differentiation, and gene expression outputs. We computed, for each expression cluster, the ratios of genes (defined as gene body ± 10 kb) harboring any given state at a given time point. The output was normalized chromatin state heat maps linked to genes at each time point (Figure 3A; see Methods). We also determined time points at which gene ratios significantly differ from the previous or following time point (P < 0.05; t-test with Bonferroni correction; Figure 3A, red dots). Since a difference in gene ratios for any given state between two time points represents a gain or loss of that state for these genes, the maps in effect show the levels of enrichment of genes in a given chromatin state at each time point, and significant changes in chromatin state between time points. Normalization of gene ratios, i.e. chromatin state enrichment levels, within expression clusters allows chromatin state maps to be compared between clusters.

Figure 3
figure 3

Gene expression level correlates with chromatin signature complexity. (A) Heat maps of chromatin state enrichment in gene regions (gene length ± 10 kb) for each expression cluster during adipogenic differentiation. Gene ratios (see text) at each differentiation time point are scaled from low-to-high and reflect chromatin state ‘loads’. Time point transitions at which a difference in gene ratios is significant are marked (dot; P < 0.05; t-test with Bonferroni correction). (B) Matrix of statistically significant differences in global state enrichment between expression clusters; P-values are shown in Additional file 2: Table S6A; Wilcoxon test with Bonferroni correction. (C) Chromatin states and RNA-seq profiles at the inactive EOMES locus and at the expressed adipogenic SREBF1 locus; FPKM ranges are shown. (D) Regression analysis of gene expression level as a function of the number of chromatin states, i.e. complexity of chromatin signatures, in gene regions.

Full size image

Strong gene expression is associated with high chromatin signature complexity

We first examined the relationship between gene expression patterns and global chromatin state enrichment level (or ‘load’). From the expression cluster-based chromatin state maps, we computed a statistics matrix of significant differences in global chromatin state enrichment between expression clusters (Figure 3B). To this end, we averaged all gene ratios within each cluster for all states (excluding state 7) and compared ratios between clusters. We find that strongly expressed genes (cluster 3, 4, 8) display greater chromatin state enrichment than weakly expressed genes (clusters 1, 2) (Figures 3A,B; red box 1; P < 0.01, Wilcoxon test; see Additional file 2: Table S6A for P-values). This is exemplified for individual genes: EOMES, essentially not expressed in ASCs or adipocytes (cluster 1) is marked by bank state 7 and the repressive states cs8 and cs9, whereas the pro-adipogenic factor SREBF1, expressed at all stages (cluster 8), is heavily enriched in states (Figure 3C; see also Additional file 1: Figure S2B). Clusters of strongly expressed genes (clusters 3, 4, 8) also display greater state enrichment than all other clusters (clusters 9-19; Figure 3A,B), which can be related to the lower expression level in the latter (see Figure 1E). Linear regression analysis shows that gene expression level correlates with global enrichment of chromatin in hPTMs (Pearson correlation 0.52; Figure 3D). This would be expected from the nature of the hPTMs examined, which for the most part characterize active chromatin domains. Lastly, genes that are stably induced or repressed (clusters 5-7 and 9-11) or transiently up- or down-regulated (clusters 12-19) show similar chromatin state enrichment level regardless of the time point at which gene expression changes are detected (Figure 3B, blue boxes 1 and 2). Thus, gene up- or down-regulation, and timing thereof, do not affect the global chromatin state load in the gene regions implicated.

Our results indicate that based on our HMM, chromatin signature complexity positively correlates with gene expression level, more so than with a switch in expression pattern such as activation or repression. This implies that timing of gene activation or repression is not determined by global epigenetic load; chromatin state enrichment is rather related to gene expression level.

Temporal dynamics of chromatin states

The temporal chromatin state maps we generated enable a dynamic analysis of chromatin state enrichment in the course of differentiation. We first find that the total number of significant state alterations between differentiation time points (Figure 3A, dots) is greater after adipogenic induction (D0-D3 transition) than during progression towards an adipogenic phenotype (D3-D9 transition; Figure 4A; P ≤ 0.001, t-test with Bonferroni correction). This suggests that once adipogenic commitment is initiated, patterns of chromatin modifications tend to stabilize.

Figure 4
figure 4

Temporal and dynamic changes in chromatin states. (A) Numbers of significant chromatin state changes (all clusters confounded) at the D0/D3 and D3/D9 transitions (t-test with Bonferroni correction). (B) Dynamics of chromatin states in each expression cluster. Graphs represent the sum of changes in P-values for each cluster. (C) Matrix of significant differences in chromatin state dynamics between clusters: clusters in columns (top) show greater state dynamics than clusters in rows (left). P-values for (B) and (C) are shown in Additional file 2: Table S6B (one-sided Wilcoxon tests). (D) Browser view of chromatin state changes and RNA-seq profiles at the PPARG locus during differentiation; the red boxed area is enlarged (zoom, right).

Full size image

Next, we compared global chromatin state dynamics between expression clusters by calculating the P-value of differences in state enrichment levels between time points in each cluster, averaging these P-values within clusters, and comparing cluster-averaged P-values. We first find that strongly expressed genes, in addition to being heavily epigenetically marked, display highly dynamic chromatin; in other words, states are locally dynamically changing between differentiation time points (Figure 3A, dots; Figure 4B, clusters 2-4). Second, chromatin states at genes that are expressed and further up-regulated during differentiation (cluster 8) are by far the most dynamic (Figure 4B,C; P = 2×10−3 to 4×10−8, one-sided Wilcoxon test; Additional file 2: Table S6B). For instance at the PPARG locus two dominant states, blank cs7 and weak enhancer state cs4 (H3K4me1) are substituted by active enhancer states (cs1, 2, 3, 12) and logically a state annotating transcribed gene bodies (cs6; H3K36me3; Figure 4D). Third, short temporal gene activation events, detected at one time point only, are associated with the lowest state dynamics (Figure 4B, clusters 12, 13). It is notably lower than the chromatin dynamics of genes showing longer lasting changes or several up- and down-regulation events changes (Figure 4B,C; clusters 5, 6, 14, 15 and 19). These observations suggest that chromatin is less prone to remodeling when genes are only transiently induced than when a more durable change occurs; this is again consistent with transcription being associated with dynamic chromatin states.

Bivalent promoter states are the most dynamic during differentiation, while weak enhancers retain their chromatin signature

Analysis of significant temporal changes in individual chromatin states during differentiation reveals the dynamic nature of specific promoter states while in contrast, enhancer states are more stable. We find that, together with cs6 (H3K36me3), cs5 (H3K4me1/H3K36me3) is the least dynamic: it shows the lowest frequency of change throughout clusters (see Figure 3A) and changes of least significance between time points (Figure 5A). State 5 is annotated as ‘enhancer in active gene bodies’ (see Figure 2A), and reflects the variable localization of enhancers relative to the promoters they regulate [3]. Stability of cs5 during adipogenesis is exemplified by the steadily expressed DNTTIP2 locus (Figure 5B) and is evidenced for several gene expression modules (Figure 5C). Additional enhancer states marked by H3K4me1, such as cs1, cs2 and cs4 (a ‘weak’ enhancer state) [21], or cs12 and cs13, are also among the temporally least dynamic (Figure 5A,B). State 1, an active enhancer state marked by H3K4me1 and H3K27ac, is more dynamic than weak enhancer states (Figure 5A; P = 0.049, two-sample Wilcoxon test): this may be explained by the connection between H3K27ac and enhancer activity, i.e. expression of the gene(s) it regulates [31]. We infer from these data that enhancers retain epigenetic identity through H3K4me1 marking in the adipogenic lineage.

Figure 5
figure 5

Bivalent promoter states are the most dynamic during differentiation. (A) Significance level of chromatin state change in all clusters during differentiation. (B) Browser views of chromatin state stability (cs5, 12, 13; top panel, stars) and dynamics (cs9, bottom panel, stars) throughout differentiation. Note the acquisition of cs9 (H3K4me2/H3K27me3) on the SORT1 promoter on D0. RNA-seq profiles are also shown (scale: 0-200 FKPM). Window sizes shown: top, 8 kb; bottom, 6 kb. (C) Chromatin state enrichment profiles during adipogenic differentiation, illustrating the low and high dynamics of cs5 (left) and cs9 (right), respectively, for genes in the various clusters shown.

Full size image

In contrast to enhancer states, the ‘bivalent’ promoter state cs9 (H3K4me2/H3K27me3) is the most dynamic during differentiation (Figure 5A; P = 0.02 to 2.6x10−6). This is exemplified by the SORT1 promoter (Figure 5B), controlling developmental expression of sortilin-1, involved in lipoprotein metabolism [32]. The dynamics of cs9 is consistent with enrichment of this apparently bivalent state on developmentally-regulated promoters in embryonic stem cells [33, 34] and in ASCs (data not shown), and its resolution into a H3K4 or K27 methylated state coincident with promoter activation or repression, respectively.

The temporal relationship between cs9 and transcript levels during differentiation is complex (Figure 5C) and reflects both a regulatory and priming role of this state on developmental gene expression. For instance, increase in cs9 from D0 to D3 on genes activated on D9, when cs9 load declines (cluster 7; Figures 2A and 5A), suggests that bivalent promoters are marked early during differentiation, prior to gene activation. In contrast, cs9 may mark expressed genes for repression (Figure 5C, cluster 11) or, consistent with a recent report [35], it may target genes after transcriptional inactivation is initiated (Figure 5C, cluster 10). cs9 load can also parallel transcript levels during differentiation (e.g. cluster 15; Figure 5C). Our data suggest that temporal dynamics of cs9 is linked to the H3K27me3 mark, as cs11, 12 and 13, promoter states deprived of H3K27me3, are less dynamic than cs9 (Figure 5A; P = 2.6×10−5 to 3.9×10−6). The temporal dynamics of the bivalent promoter state cs9 identified here infers its sensitivity to fluctuating changes in response to environmental stimuli.

Discussion

We report a 15-chromatin state HMM which describes temporal changes in chromatin signatures in relation to gene expression patterns throughout adipogenic differentiation of human primary ASCs. Identification of coherent gene expression modules and of chromatin states fitting these modules demonstrates the complexity and dynamic nature of hPTM combinations during adipogenesis. These expression modules reveal distinct patterns including sustained expression, up- and down-regulation transitions, single-pulse patterns and oscillatory patterns. Interestingly, single-pulse expression patterns prevail in responses to stimuli [36], suggesting an adaptive response of ASCs to environmental changes such as cell cycle exit and adipogenic induction. These pulses initiate a downstream cascade of expression changes with sequential offsets which may involve feed-forward loops [37].

Active chromatin regions display tremendous diversity in epigenetic signatures (Figure 6A), suggestive of plasticity and amenability to remodeling. This complexity would be anticipated from the modifications examined, which mostly mark active chromatin. Our findings nonetheless also show the remarkable dynamics of these states at a given locus in the cell populations analyzed. Clearly, the diversity of chromatin state patterns may reflect complexity in a population of cells inasmuch as complexity within each cell in the population. Support to this interpretation is lent by a recent single-cell RNA-seq study of cell differentiation [29]. Single-cell quantitative epigenetic analyses would help clarify this issue but remain currently inaccessible at the level of resolution required. Gene expression may regardless play an important role in generating this chromatin diversity, perhaps by providing a chromatin template accessible for histone modifiers. The functional interplay between chromatin modifiers, hPTMs and transcriptional outputs has been difficult to underpin because this relationship is development- and locus-dependent [3, 10, 38, 39]). Chromatin engineering strategies may prove useful in dissecting the effect of histone modifications on gene expression [40].

Figure 6
figure 6

A model of temporal patterns of chromatin state changes during differentiation. (A) Increased gene expression level is associated with enhanced enrichment in chromatin states (color bars). Promoter regions (P) are particularly prone to modifications. G, gene body; E, enhancer. (B) Patterns of chromatin state dynamics during adipogenic differentiation. A blank state is commonly temporally followed by loading of chromatin states (color bars). States can then be replaced by another (substitution), become a site of emergence of one or more other states (insertion), or may shift in position (shift).

Full size image

Gene activation and inactivation is also associated with a remodeling of hPTMs, though not always concomitant to expression changes. Both the extent of chromatin modifications and chromatin state dynamics are lower in gene regions that are transiently up- or down-regulated than in highly active domains. This suggests that temporal gene expression is under control of factors other than hPTMs in a pre-disposed chromatin environment [16]. Supporting this idea, ‘stand-by’ occupancy of regulatory regions by pioneer transcription factors precedes activation of developmentally-controlled genes, conferring transcriptional competence [41–43]. Cooperative binding of transcription factors also plays a key role in the induction of adipogenic genes [15, 44]. Further, the correlation between transcription factor binding and target gene activation timing [45], together with differences in transcription factor affinity [46], could underlie the gene activation offsets within the adipogenic program.

Our temporal HMM reveals several patterns of chromatin state changes during adipogenic differentiation (Figure 6B). (i) The simplest is the formation of one or more state from a blank state, which is often detected at loci that become activated. (ii) Conversely, upon or following gene inactivation, one or more state can be lost, giving rise to a ‘blank’ state. (iii) One state can also be substituted by another: this occurs in the form of replacement of one state by another or by insertion of a state into a larger pre-existing domain. This is illustrated by the introduction of enhancer sites within gene bodies concomitant with gene up-regulation, such as our observations on PPARG. (iv) States may also spatially shift over a limited area, leading to the perception of local replacement; this may represent a regulated mechanism or a stochastic event of functional importance for proper gene regulation [47]. Restricted positional shifts of chromatin states are consistent with nucleosome sliding and histone replacement (along with their PTMs) events that occur on promoters and bodies of active genes [48–51].

The finding that the bivalent promoter state cs9 is the most dynamic during adipogenesis is surprising because promoter states tend to be conserved between cell types [21]. cs9 dynamics seems to be linked to the Polycomb-associated H3K27me3 mark, as promoter states with no H3K27 methylation are significantly less prone to change. The dynamic nature of Polycomb repressor complex 2 and of ensuing H3K27 trimethylation has been documented [52, 53]. Notably, an important caveat in the identification of a ‘bivalent’ state is that it remains unclear whether co-methylation on K4 and K27 occurs at the same locus or whether this state reflects heterogeneity in the cell populations. HMM data are compatible with both interpretations because input data into the model result from epigenomic information from not only all cells in a given population but also from all cell populations (i.e. here, all time points) examined [19]. Gain of cs9 on developmentally-regulated genes before transcriptional activation, or conversely loss of cs9 (i.e. loss of K4 or K27 methylation) reinforces the importance of Polycomb-mediated marking for proper adipogenic differentiation [54].

Inasmuch as promoter bivalency may precede developmental gene activation, de novo enhancer establishment by H3K4me1 marking emerges as a predictor of lineage-specific enhancer usage and downstream expression of the gene(s) it regulates [3, 10]. Thousands of H3K4me1-marked enhancers found in undifferentiated ASCs are retained through differentiation, in line with the adipogenic commitment of ASCs [55, 56]. In contrast, active enhancers (H3K4me1/H3K27ac; cs1 and cs2) are more dynamic, in keeping with the cell type specificity of H3K27ac [3, 10, 23]. A functional advantage in maintaining enhancer identity within a lineage may be by favoring the formation of ‘hubs’ of signal integration and relay to efficiently remodel chromatin and activate lineage-specific gene networks [17, 57, 58].

A remaining question is the origin of the multiple histone modifications detected on promoters and enhancers, including some known to commonly mark either element. This may reflect the embedding of enhancers in distinct domains such as promoters or transcribed exons, illustrating the variable localization of enhancer elements relative to promoters and genes they modulate [3]. Recent evidence indicates that enhancers interact with other elements, particularly promoters, through 3-dimensional chromatin looping [59–62]. Future studies will be important to determine the extent to which three-dimensional conformation of the genome impacts chromatin states, their developmental transitions and gene expression.

Conclusions

Genome-scale modeling of chromatin into 15 states in the course of adipogenic differentiation demonstrates a hyper-dynamic nature of bivalent promoter states, which underlies distinct and coherent gene expression cohorts. Our results take our understanding of transitions in chromatin organization through the genome and in a developmental context to a new level, and constitute a resource enabling the assessment of possibilities to manipulate the adipogenic differentiation program.

Methods

Cells and adipogenic induction

Human adipose tissue stromal cells ASCs [63] were cultured under proliferative conditions in DMEM/F12 (Life Technologies) containing 10% FBS, 20 ng/ml basic fibroblast growth factor and 10 ng/ml epidermal growth factor (Sigma-Aldrich). Cells at passage 6-7 were used for differentiation in two independent experiments. Two days prior to adipogenic induction, cells were harvested and either used for analysis (D-2 time point), or reseeded to confluency in DMEM/F12/10% FBS without growth factors. Adipogenesis was induced on D0 by adding 0.5 μM 1-methyl-3 isobutylxanthine (Dumex Alpharma), 1 μM dexamethasone (Dumex Alpharma), 10 μg/ml insulin (Sigma-Aldrich) and 200 μM indomethacin (Dumex Alpharma) for up to 9 days [63].

RNA isolation and RNA-sequencing

Total RNA was isolated on D-2, D0, D3 and D9 of adipogenic stimulation using the Ambion TRIzol® Reagent RNA extraction kit (Life Technologies). Libraries were prepared according to the Illumina protocol and sequenced to generate 100 base pair paired-end reads on an Illumina HiSeq2500. RNA-seq reads were processed using the Tuxedo pipeline [64]. TopHat [65] was used to align reads with no mismatch against human genome UCSC hg19 with default settings, applying the bowtie2 [66] preset ‘-very sensitive’. Cufflinks and cuffdiff were run using default settings and bias correction. Results were analyzed and visualized in R through cummeRbund. Genes with a log (fold change) > 2 and α < 0.05 were considered as differentially expressed. Coefficient of variation (CV2) was plotted against log (FPKM) (fragments per kilobase of exon per million fragments mapped) to identify differences in distribution of expression patterns. Principal component analysis was done using the first two principal components. Additional scripting was done in Perl or R.

Clustering genes based on expression patterns

Genes used for clustering were obtained from the Tuxedo protocol with FPKM > 0 at at least one time point examined. Genes with similar expression patterns were clustered through hierarchical clustering using hclust in R. Gene Ontology enrichment analysis was done in Gorilla using default parameters [67].

Attributing chromatin states to genes and gene expression clusters

ChIP-seq datasets of hPTMs and CTCF were obtained from a previous study [18] using a similar cell source, differentiation protocol and time line. ChIP-seq reads were re-mapped using bowtie and ‘-best’ default settings. hPTM and CTCF enrichment data were converted into chromatin states in consecutive 200-base pair bins using ChromHMM [19]. Options were selected to learn a 15-state model using the Baum-Welch training algorithm. States were linked to genes by computing the presence of a given state in the gene body plus a 10 kb extension upstream and downstream (‘gene ± 10 kb’) to take into account regulatory regions. This 10 kb value was set to 50% of the upper quantile of the distance of a changing chromatin state to the nearest gene/upstream or downstream) during differentiation. We found the mean distance to be 20 kb, calculated from the following data: 22.31 kb at the D-2/D0 transition, 22.52 kb at the D0/D3 transition, and 22.81 kb at the D3/D9 transition (data not shown). Note that increasing the extension window size from 10 kb to 22 kb did not significantly alter the number of chromatin states attributed to genes (data not shown). Chromatin state counts were normalized between gene expression clusters by dividing the number of genes (±10 kb) containing a given chromatin state at each differentiation time point by the total number of genes in each entire expression cluster. This generated ‘gene ratios’, which were used in statistical analyses. Normalization enabled a comparison of chromatin state levels within and between gene expression clusters.

Statistical analyses

Statistical analyses were done in R. Regression and correlation analysis of chromatin state enrichment as a function of gene expression level was done using the Pearson method from numbers of chromatin states in gene regions (gene ± 10 kb; see above) and log (FPKM + 1) values. P-values for differences in enrichment of a chromatin state between differentiation time points were calculated using parametric one sample t-tests and Bonferroni correction to remove false positives, comparing each gene ratio difference to the overall difference in ratios in the cluster. Differences between clusters for chromatin state enrichment were identified using a non-parametric two sample Wilcoxon test and Bonferroni correction to remove false positives, as samples showed a skewed distribution; to this end, we used gene ratios at each time point for one cluster and compared them to all other clusters in a pair-wise manner. Differences in overall significance of chromatin state changes between clusters were tested by a non-parametric one-sided Wilcoxon test since distribution of the P-values was skewed. Graphs identifying the overall differences in the significance of P-values (chromatin state changes) were generated by computing the -log10 value of the sum of the P-values.

Data viewing

Browser views of gene tracks, ChIP-seq data and chromatin states are shown using Integrated Genomics Viewer (IGV; broadinstitute.org/igv) [68]. Unless otherwise stated genes considered in the analyses are from the Illumina iGenomes gene annotation with UCSC data source for hg19 (support.illumina.com/sequencing/sequencing_software/igenome.ilmn).

Data access

Our RNA-seq data are available from NCBI GEO accession number GSE60237. Published ChIP-seq data [18] were from NCBI GEO accession number GSE20752.

Abbreviations

ASC:

Adipose tissue stromal cell

ChIP:

Chromatin immunoprecipitation

ChIP-seq:

ChIP and high-throughput sequencing

cs:

chromatin state

RNA-seq:

RNA-sequencing

TSS:

Transcription start site.

References

  1. Roy S, Ernst J, Kharchenko PV, Kheradpour P, Negre N, Eaton ML, Landolin JM, Bristow CA, Ma L, Lin MF, Washietl S, Arshinoff BI, Ay F, Meyer PE, Robine N, Washington NL, Di SL, Berezikov E, Brown CD, Candeias R, Carlson JW, Carr A, Jungreis I, Marbach D, Sealfon R, Tolstorukov MY, Will S, Alekseyenko AA, Artieri C, Booth BW, et al: Identification of functional elements and regulatory circuits by Drosophila modENCODE. Science. 2010, 330: 1787-1797.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  2. Lindeman LC, Andersen IS, Reiner AH, Li N, Aanes H, Østrup O, Winata CL, Mathavan S, Müller F, Aleström P, Collas P: Pre-patterning of developmental gene expression by modified histones before zygotic genome activation. Dev Cell. 2011, 21: 993-1004. 10.1016/j.devcel.2011.10.008.

    Article  CAS  PubMed  Google Scholar 

  3. Bonn S, Zinzen RP, Girardot C, Gustafson EH, Perez-Gonzalez A, Delhomme N, Ghavi-Helm Y, Wilczynski B, Riddell A, Furlong EE: Tissue-specific analysis of chromatin state identifies temporal signatures of enhancer activity during embryonic development. Nat Genet. 2012, 44: 148-156. 10.1038/ng.1064.

    Article  CAS  PubMed  Google Scholar 

  4. Ghavi-Helm Y, Furlong EE: Analyzing transcription factor occupancy during embryo development using ChIP-seq. Methods Mol Biol. 2012, 786: 229-245. 10.1007/978-1-61779-292-2_14. 229-245

    Article  CAS  PubMed  Google Scholar 

  5. Wilczynski B, Liu YH, Yeo ZX, Furlong EE: Predicting spatial and temporal gene expression using an integrative model of transcription factor occupancy and chromatin state. PLoS Comput Biol. 2012, 8: e1002798-10.1371/journal.pcbi.1002798.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  6. Bernstein BE, Meissner A, Lander ES: The mammalian epigenome. Cell. 2007, 128: 669-681. 10.1016/j.cell.2007.01.033.

    Article  CAS  PubMed  Google Scholar 

  7. Bernstein BE, Stamatoyannopoulos JA, Costello JF, Ren B, Milosavljevic A, Meissner A, Kellis M, Marra MA, Beaudet AL, Ecker JR, Farnham PJ, Hirst M, Lander ES, Mikkelsen TS, Thomson JA: The NIH roadmap epigenomics mapping consortium. Nat Biotechnol. 2010, 28: 1045-1048. 10.1038/nbt1010-1045.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  8. Zhou VW, Goren A, Bernstein BE: Charting histone modifications and the functional organization of mammalian genomes. Nat Rev Genet. 2011, 12: 7-18.

    Article  PubMed  Google Scholar 

  9. Yu P, Xiao S, Xin X, Song CX, Huang W, McDee D, Tanaka T, Wang T, He C, Zhong S: Spatiotemporal clustering of the epigenome reveals rules of dynamic gene regulation. Genome Res. 2013, 23: 352-364. 10.1101/gr.144949.112.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  10. Lara-Astiaso D, Weiner A, Lorenzo-Vivas E, Zaretsky I, Jaitin DA, David E, Keren-Shaul H, Mildner A, Winter D, Jung S, Friedman N, Amit I: Chromatin state dynamics during blood formation. Science. 2014, 345: 943-949. 10.1126/science.1256271.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  11. Saeed S, Quintin J, Kerstens HH, Rao NA, Aghajanirefah A, Matarese F, Cheng SC, Ratter J, Berentsen K, van der Ent MA, Sharifi N, Janssen-Megens EM, Ter HM, Mandoli A, Van ST, Ng A, Burden F, Downes K, Frontini M, Kumar V, Giamarellos-Bourboulis EJ, Ouwehand WH, van der Meer JW, Joosten LA, Wijmenga C, Martens JH, Xavier RJ, Logie C, Netea MG, Stunnenberg HG: Epigenetic programming of monocyte-to-macrophage differentiation and trained innate immunity. Science. 2014, 345: 1251086-10.1126/science.1251086.

    Article  PubMed Central  PubMed  Google Scholar 

  12. Rosen ED, Spiegelman BM: Adipocytes as regulators of energy balance and glucose homeostasis. Nature. 2006, 444: 847-853. 10.1038/nature05483.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  13. Cristancho AG, Lazar MA: Forming functional fat: a growing understanding of adipocyte differentiation. Nat Rev Mol Cell Biol. 2011, 12: 722-734. 10.1038/nrm3198.

    Article  CAS  PubMed  Google Scholar 

  14. Lefterova MI, Haakonsson AK, Lazar MA, Mandrup S: PPARgamma and the global map of adipogenesis and beyond. Trends Endocrinol Metab. 2014, 25: 293-302. 10.1016/j.tem.2014.04.001.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  15. Siersbaek R, Nielsen R, Mandrup S: Transcriptional networks and chromatin remodeling controlling adipogenesis. Trends Endocrinol Metab. 2012, 23: 56-64. 10.1016/j.tem.2011.10.001.

    Article  CAS  PubMed  Google Scholar 

  16. Siersbaek R, Baek S, Rabiee A, Nielsen R, Traynor S, Clark N, Sandelin A, Jensen ON, Sung MH, Hager GL, Mandrup S: Molecular architecture of transcription factor hotspots in early adipogenesis. Cell Rep. 2014, 7: 1434-1442. 10.1016/j.celrep.2014.04.043.

    Article  CAS  PubMed  Google Scholar 

  17. Siersbaek R, Rabiee A, Nielsen R, Sidoli S, Traynor S, Loft A, Poulsen LL, Rogowska-Wrzesinska A, Jensen ON, Mandrup S: Transcription factor cooperativity in early adipogenic hotspots and super-enhancers. Cell Rep. 2014, 7: 1443-1455. 10.1016/j.celrep.2014.04.042.

    Article  CAS  PubMed  Google Scholar 

  18. Mikkelsen TS, Xu Z, Zhang X, Wang L, Gimble JM, Lander ES, Rosen ED: Comparative epigenomic analysis of murine and human adipogenesis. Cell. 2010, 143: 156-169. 10.1016/j.cell.2010.09.006.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  19. Ernst J, Kellis M: ChromHMM: automating chromatin-state discovery and characterization. Nat Methods. 2012, 9: 215-216. 10.1038/nmeth.1906.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  20. Hoffman MM, Buske OJ, Wang J, Weng Z, Bilmes JA, Noble WS: Unsupervised pattern discovery in human chromatin structure through genomic segmentation. Nat Methods. 2012, 9: 473-476. 10.1038/nmeth.1937.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  21. Ernst J, Kheradpour P, Mikkelsen TS, Shoresh N, Ward LD, Epstein CB, Zhang X, Wang L, Issner R, Coyne M, Ku M, Durham T, Kellis M, Bernstein BE: Mapping and analysis of chromatin state dynamics in nine human cell types. Nature. 2011, 473: 43-49. 10.1038/nature09906.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  22. Kharchenko PV, Alekseyenko AA, Schwartz YB, Minoda A, Riddle NC, Ernst J, Sabo PJ, Larschan E, Gorchakov AA, Gu T, Linder-Basso D, Plachetka A, Shanower G, Tolstorukov MY, Luquette LJ, Xi R, Jung YL, Park RW, Bishop EP, Canfield TK, Sandstrom R, Thurman RE, MacAlpine DM, Stamatoyannopoulos JA, Kellis M, Elgin SC, Kuroda MI, Pirrotta V, Karpen GH, Park PJ: Comprehensive analysis of the chromatin landscape in Drosophila melanogaster. Nature. 2011, 471: 480-485. 10.1038/nature09725.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  23. Negre N, Brown CD, Ma L, Bristow CA, Miller SW, Wagner U, Kheradpour P, Eaton ML, Loriaux P, Sealfon R, Li Z, Ishii H, Spokony RF, Chen J, Hwang L, Cheng C, Auburn RP, Davis MB, Domanus M, Shah PK, Morrison CA, Zieba J, Suchy S, Senderowicz L, Victorsen A, Bild NA, Grundstad AJ, Hanley D, MacAlpine DM, Mannervik M, et al: A cis-regulatory map of the Drosophila genome. Nature. 2011, 471: 527-531. 10.1038/nature09990.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  24. Ram O, Goren A, Amit I, Shoresh N, Yosef N, Ernst J, Kellis M, Gymrek M, Issner R, Coyne M, Durham T, Zhang X, Donaghey J, Epstein CB, Regev A, Bernstein BE: Combinatorial patterning of chromatin regulators uncovered by genome-wide location analysis in human cells. Cell. 2011, 147: 1628-1639. 10.1016/j.cell.2011.09.057.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  25. Zhu J, Adli M, Zou JY, Verstappen G, Coyne M, Zhang X, Durham T, Miri M, Deshpande V, De Jager PL, Bennett DA, Houmard JA, Muoio DM, Onder TT, Camahort R, Cowan CA, Meissner A, Epstein CB, Shoresh N, Bernstein BE: Genome-wide chromatin state transitions associated with developmental and environmental cues. Cell. 2013, 152: 642-654. 10.1016/j.cell.2012.12.033.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  26. Heintzman ND, Hon GC, Hawkins RD, Kheradpour P, Stark A, Harp LF, Ye Z, Lee LK, Stuart RK, Ching CW, Ching KA, Ntosiewicz-Bourget JE, Liu H, Zhang X, Green RD, Lobanenkov VV, Stewart R, Thomson JA, Crawford GE, Kellis M, Ren B: Histone modifications at human enhancers reflect global cell-type-specific gene expression. Nature. 2009, 459: 108-112. 10.1038/nature07829.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  27. Lee Y, Ghosh D, Zhang Y: Regression hidden Markov modeling reveals heterogeneous gene expression regulation: a case study in mouse embryonic stem cells. BMC Genomics. 2014, 15: 360-10.1186/1471-2164-15-360.

    Article  PubMed Central  PubMed  Google Scholar 

  28. Marinov GK, Williams BA, McCue K, Schroth GP, Gertz J, Myers RM, Wold BJ: From single-cell to cell-pool transcriptomes: stochasticity in gene expression and RNA splicing. Genome Res. 2014, 24: 496-510. 10.1101/gr.161034.113.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  29. Trapnell C, Cacchiarelli D, Grimsby J, Pokharel P, Li S, Morse M, Lennon NJ, Livak KJ, Mikkelsen TS, Rinn JL: The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol. 2014, 32: 381-386. 10.1038/nbt.2859.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  30. Filion GJ, van Bemmel JG, Braunschweig U, Talhout W, Kind J, Ward LD, Brugman W, de Castro IJ, Kerkhoven RM, Bussemaker HJ, van Steensel B: Systematic protein location mapping reveals five principal chromatin types in Drosophila cells. Cell. 2010, 143: 212-224. 10.1016/j.cell.2010.09.009.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  31. Rada-Iglesias A, Bajpai R, Swigut T, Brugmann SA, Flynn RA, Wysocka J: A unique chromatin signature uncovers early developmental enhancers in humans. Nature. 2011, 470: 279-283. 10.1038/nature09692.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  32. Strong A, Patel K, Rader DJ: Sortilin and lipoprotein metabolism: making sense out of complexity. Curr Opin Lipidol. 2014, Aug 12. [Epub ahead of print]

    Google Scholar 

  33. Azuara V, Perry P, Sauer S, Spivakov M, Jorgensen HF, John RM, Gouti M, Casanova M, Warnes G, Merkenschlager M, Fisher AG: Chromatin signatures of pluripotent cell lines. Nat Cell Biol. 2006, 8: 532-538. 10.1038/ncb1403.

    Article  CAS  PubMed  Google Scholar 

  34. Bernstein BE, Mikkelsen TS, Xie X, Kamal M, Huebert DJ, Cuff J, Fry B, Meissner A, Wernig M, Plath K, Jaenisch R, Wagschal A, Feil R, Schreiber SL, Lander ES: A bivalent chromatin structure marks key developmental genes in embryonic stem cells. Cell. 2006, 125: 315-326. 10.1016/j.cell.2006.02.041.

    Article  CAS  PubMed  Google Scholar 

  35. Riising EM, Comet I, Leblanc B, Wu X, Johansen JV, Helin K: Gene silencing triggers polycomb repressive complex 2 recruitment to CpG islands genome wide. Mol Cell. 2014, 55: 347-360. 10.1016/j.molcel.2014.06.005.

    Article  CAS  PubMed  Google Scholar 

  36. Chechik G, Oh E, Rando O, Weissman J, Regev A, Koller D: Activity motifs reveal principles of timing in transcriptional control of the yeast metabolic network. Nat Biotechnol. 2008, 26: 1251-1259. 10.1038/nbt.1499.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  37. Yosef N, Regev A: Impulse control: temporal dynamics in gene transcription. Cell. 2011, 144: 886-896. 10.1016/j.cell.2011.02.015.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  38. Batta K, Zhang Z, Yen K, Goffman DB, Pugh BF: Genome-wide function of H2B ubiquitylation in promoter and genic regions. Genes Dev. 2011, 25: 2254-2265. 10.1101/gad.177238.111.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  39. Wu L, Li L, Zhou B, Qin Z, Dou Y: H2B ubiquitylation promotes RNA Pol II processivity via PAF1 and pTEFb. Mol Cell. 2014, 54: 920-931. 10.1016/j.molcel.2014.04.013.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  40. Pick H, Kilic S, Fierz B: Engineering chromatin states: Chemical and synthetic biology approaches to investigate histone modification function. Biochim Biophys Acta. 1839, 2014: 644-656.

    Google Scholar 

  41. Blythe SA, Cha SW, Tadjuidje E, Heasman J, Klein PS: Beta-Catenin primes organizer gene expression by recruiting a histone H3 arginine 8 methyltransferase, Prmt2. Dev Cell. 2010, 19: 220-231. 10.1016/j.devcel.2010.07.007.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  42. Liber D, Domaschenz R, Holmqvist PH, Mazzarella L, Georgiou A, Leleu M, Fisher AG, Labosky PA, Dillon N: Epigenetic priming of a pre-B cell-specific enhancer through binding of Sox2 and Foxd3 at the ESC stage. Cell Stem Cell. 2010, 7: 114-126. 10.1016/j.stem.2010.05.020.

    Article  CAS  PubMed  Google Scholar 

  43. Zaret KS, Carroll JS: Pioneer transcription factors: establishing competence for gene expression. Genes Dev. 2011, 25: 2227-2241. 10.1101/gad.176826.111.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  44. Lefterova MI, Lazar MA: New developments in adipogenesis. Trends Endocrinol Metab. 2009, 20: 107-114. 10.1016/j.tem.2008.11.005.

    Article  CAS  PubMed  Google Scholar 

  45. Chechik G, Koller D: Timing of gene expression responses to environmental changes. J Comput Biol. 2009, 16: 279-290. 10.1089/cmb.2008.13TT.

    Article  CAS  PubMed  Google Scholar 

  46. Lam FH, Steger DJ, O'Shea EK: Chromatin decouples promoter threshold from dynamic range. Nature. 2008, 453: 246-250. 10.1038/nature06867.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  47. Luger K, Dechassa ML, Tremethick DJ: New insights into nucleosome and chromatin structure: an ordered state or a disordered affair?. Nat Rev Mol Cell Biol. 2012, 13: 436-447. 10.1038/nrm3382.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  48. Henikoff S: Nucleosome destabilization in the epigenetic regulation of gene expression. Nat Rev Genet. 2008, 9: 15-26. 10.1038/nrg2206.

    Article  CAS  PubMed  Google Scholar 

  49. Conerly ML, Teves SS, Diolaiti D, Ulrich M, Eisenman RN, Henikoff S: Changes in H2A.Z occupancy and DNA methylation during B-cell lymphomagenesis. Genome Res. 2010, 20: 1383-1390. 10.1101/gr.106542.110.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  50. Soboleva TA, Nekrasov M, Ryan DP, Tremethick DJ: Histone variants at the transcription start-site. Trends Genet. 2014, 30: 199-209. 10.1016/j.tig.2014.03.002.

    Article  CAS  PubMed  Google Scholar 

  51. Weber CM, Ramachandran S, Henikoff S: Nucleosomes Are Context-Specific, H2A.Z-Modulated Barriers to RNA Polymerase. Mol Cell. 2014, 53: 819-830. 10.1016/j.molcel.2014.02.014.

    Article  CAS  PubMed  Google Scholar 

  52. Bracken AP, Dietrich N, Pasini D, Hansen KH, Helin K: Genome-wide mapping of Polycomb target genes unravels their roles in cell fate transitions. Genes Dev. 2006, 20: 1123-1136. 10.1101/gad.381706.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  53. Lee TI, Jenner RG, Boyer LA, Guenther MG, Levine SS, Kumar RM, Chevalier B, Johnstone SE, Cole MF, Isono K, Koseki H, Fuchikami T, Abe K, Murray HL, Zucker JP, Yuan B, Bell GW, Herbolsheimer E, Hannett NM, Sun K, Odom DT, Otte AP, Volkert TL, Bartel DP, Melton DA, Gifford DK, Jaenisch R, Young RA: Control of developmental regulators by Polycomb in human embryonic stem cells. Cell. 2006, 125: 301-313. 10.1016/j.cell.2006.02.043.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  54. Wang L, Jin Q, Lee JE, Su IH, Ge K: Histone H3K27 methyltransferase Ezh2 represses Wnt genes to facilitate adipogenesis. Proc Natl Acad Sci U S A. 2010, 107: 7317-7322. 10.1073/pnas.1000031107.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  55. Boquest AC, Noer A, Collas P: Epigenetic programming of mesenchymal stem cells from human adipose tissue. Stem Cell Rev. 2006, 2: 319-329. 10.1007/BF02698059.

    Article  CAS  PubMed  Google Scholar 

  56. Sørensen AL, Jacobsen BM, Reiner AH, Andersen IS, Collas P: Promoter DNA methylation patterns of differentiated cells are largely programmed at the progenitor stage. Mol Biol Cell. 2010, 21: 2066-2077. 10.1091/mbc.E10-01-0018.

    Article  PubMed Central  PubMed  Google Scholar 

  57. Steger DJ, Grant GR, Schupp M, Tomaru T, Lefterova MI, Schug J, Manduchi E, Stoeckert CJ, Lazar MA: Propagation of adipogenic signals through an epigenomic transition state. Genes Dev. 2010, 24: 1035-1044. 10.1101/gad.1907110.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  58. Buecker C, Wysocka J: Enhancers as information integration hubs in development: lessons from genomics. Trends Genet. 2012, 28: 276-284. 10.1016/j.tig.2012.02.008.

    Article  CAS  PubMed  Google Scholar 

  59. Chepelev I, Wei G, Wangsa D, Tang Q, Zhao K: Characterization of genome-wide enhancer-promoter interactions reveals co-expression of interacting genes and modes of higher order chromatin organization. Cell Res. 2012, 22: 490-503. 10.1038/cr.2012.15.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  60. Li G, Ruan X, Auerbach RK, Sandhu KS, Zheng M, Wang P, Poh HM, Goh Y, Lim J, Zhang J, Sim HS, Peh SQ, Mulawadi FH, Ong CT, Orlov YL, Hong S, Zhang Z, Landt S, Raha D, Euskirchen G, Wei CL, Ge W, Wang H, Davis C, Fisher-Aylor KI, Mortazavi A, Gerstein M, Gingeras T, Wold B, Sun Y, et al: Extensive promoter-centered chromatin interactions provide a topological basis for transcription regulation. Cell. 2012, 148: 84-98. 10.1016/j.cell.2011.12.014.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  61. Sanyal A, Lajoie BR, Jain G, Dekker J: The long-range interaction landscape of gene promoters. Nature. 2012, 489: 109-113. 10.1038/nature11279.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  62. Ghavi-Helm Y, Klein FA, Pakozdi T, Ciglar L, Noordermeer D, Huber W, Furlong EE: Enhancer loops appear stable during development and are associated with paused polymerase. Nature. 2014, 512: 96-100.

    CAS  PubMed  Google Scholar 

  63. Boquest AC, Shahdadfar A, Fronsdal K, Sigurjonsson O, Tunheim SH, Collas P, Brinchmann JE: Isolation and transcription profiling of purified uncultured human stromal stem cells: alteration of gene expression after in vitro cell culture. Mol Biol Cell. 2005, 16: 1131-1141. 10.1091/mbc.E04-10-0949.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  64. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L: Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010, 28: 511-515. 10.1038/nbt.1621.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  65. Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L: Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 2012, 7: 562-578. 10.1038/nprot.2012.016.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  66. Langmead B, Salzberg SL: Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012, 9: 357-359. 10.1038/nmeth.1923.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

  67. Eden E, Navon R, Steinfeld I, Lipson D, Yakhini Z: GOrilla: a tool for discovery and visualization of enriched GO terms in ranked gene lists. BMC Bioinformatics. 2009, 10: 48-10-

    Article  Google Scholar 

  68. Robinson JT, Thorvaldsdottir H, Winckler W, Guttman M, Lander ES, Getz G, Mesirov JP: Integrative genomics viewer. Nat Biotechnol. 2011, 29: 24-26. 10.1038/nbt.1754.

    Article  CAS  PubMed Central  PubMed  Google Scholar 

Download references

Acknowledgements

Andrew Reiner is thanked for advice on bioinformatics. Kristin Vekterud and Anita Sørensen are thanks for technical assistance. This work was supported by the Research Council of Norway and the Norwegian Stem Cell Center, South East Health Norway (ARO) and the Molecular Life Science program of the University of Oslo (AS).

Author information

Authors and Affiliations

  1. Department of Molecular Medicine, Institute of Basic Medical Sciences, Faculty of Medicine, University of Oslo, 0317, Oslo, Norway

    Akshay Shah, Anja Oldenburg & Philippe Collas

  2. Norwegian Center for Stem Cell Research, Oslo University Hospital, Oslo, Norway

    Akshay Shah, Anja Oldenburg & Philippe Collas

Authors
  1. Akshay Shah
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Anja Oldenburg
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Philippe Collas
    View author publications

    You can also search for this author in PubMed Google Scholar

Corresponding author

Correspondence to Philippe Collas.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

AS designed the study, analyzed data, made figures and wrote the manuscript. ARO contributed to study design, carried out differentiation experiments, reviewed the manuscript. PC designed the study, made figures, wrote the manuscript and supervised the work. All authors read and approved the final manuscript.

Electronic supplementary material

12864_2014_6925_MOESM1_ESM.pdf

Additional file 1: A PDF file containing two figures including. Figure S1. Transcriptional changes elicited by adipogenic differentiation of human ASCs segregate undifferentiated cells from adipogenic stimulated cells; Figure S2. Genome coverage of chromatin states identified and relation to gene expression. Table S1. Numbers of RNA-seq reads mapped. (PDF 361 KB)

12864_2014_6925_MOESM2_ESM.zip

Additional file 2: A zip file containing five Excel tables. Table S2. Genes up- or down-regulated at each stage of differentiation; Table S3. GO terms enriched among genes up- or down-regulated at each stage of differentiation; Table S4. GO terms enriched among genes included in the 19 expression clusters; Table S5: Genes included in the 19 expression clusters; Table S6: Comparisons of enrichment in chromatin states and of chromatin state dynamics between expression clusters. Tables show P-values of (A) comparisons of enrichment in chromatin states (gene ratios; Wilcoxon test with Bonferroni correction and (B) comparisons of chromatin state dynamics (gene ratio differences; one-sided Wilcoxon tests). (ZIP 1 MB)

Authors’ original submitted files for images

Rights and permissions

Open Access  This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.

The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

To view a copy of this licence, visit https://creativecommons.org/licenses/by/4.0/.

The Creative Commons Public Domain Dedication waiver (https://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Shah, A., Oldenburg, A. & Collas, P. A hyper-dynamic nature of bivalent promoter states underlies coordinated developmental gene expression modules. BMC Genomics 15, 1186 (2014). https://doi.org/10.1186/1471-2164-15-1186

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/1471-2164-15-1186

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us