Waves of chromatin modifications in mouse dendritic cells in response to LPS stimulation

Background The importance of transcription factors (TFs) and epigenetic modifications in the control of gene expression is widely accepted. However, causal relationships between changes in TF binding, histone modifications, and gene expression during the response to extracellular stimuli are not well understood. Here, we analyze the ordering of these events on a genome-wide scale in dendritic cells in response to lipopolysaccharide (LPS) stimulation. Results Using a ChIP-seq time series dataset, we find that the LPS-induced accumulation of different histone modifications follows clearly distinct patterns. Increases in H3K4me3 appear to coincide with transcriptional activation. In contrast, H3K9K14ac accumulates early after stimulation, and H3K36me3 at later time points. Integrative analysis with TF binding data reveals potential links between TF activation and dynamics in histone modifications. Especially, LPS-induced increases in H3K9K14ac and H3K4me3 are associated with binding by STAT1/2 and were severely impaired in Stat1−/− cells. Conclusions While the timing of short-term changes of some histone modifications coincides with changes in transcriptional activity, this is not the case for others. In the latter case, dynamics in modifications more likely reflect strict regulation by stimulus-induced TFs and their interactions with chromatin modifiers. Electronic supplementary material The online version of this article (10.1186/s13059-018-1524-z) contains supplementary material, which is available to authorized users.


Background
Epigenetic features, such as histone modifications and DNA methylation, are thought to play a crucial role in controlling the accessibility of DNA to RNA polymerases. Associations have been found between histone modifications and both long-term and short-term cellular processes, including development, heritability of cell type identity, DNA repair, and transcriptional control [1,2]. For cells of the hematopoietic lineage, cell type-defining enhancers are established during differentiation by priming with the H3K4me1 marker [3,4]. After differentiation, signals from the surrounding tissue environment or from pathogens induce changes in histone modifications reflecting the changes in activity of enhancers and promoters, including the de novo establishment of latent enhancers [5][6][7][8][9].
Transcription factors (TFs) are key regulators in the control of epigenetic changes [10,11]. During the long-term process of differentiation, closed chromatin is first bound by pioneer TFs, which results in structural changes that make it accessible to other TFs and RNA polymerase II (Pol2) [6,12]. Similarly, short-term changes in gene expression following stimulation of immune cells are regulated by TFs. This regulation is thought to involve TF binding, induction of changes in histone modifications, and recruitment of Pol2 [13][14][15][16]. However, details of the temporal ordering and causal relationships between these events remain poorly understood [17,18]. Especially, it is unclear whether certain histone modifications are a requirement for, or a result of, TF binding and transcription [19][20][21].
As sentinel cells of the innate immune system, dendritic cells (DCs) are well equipped for detecting the presence of pathogens. Lipopolysaccharide (LPS), a component of the cell wall of Gram-negative bacteria, is recognized by DCs through the membrane-bound Toll-like receptor 4 (TLR4), resulting in the activation of two downstream signaling pathways [22]. One pathway is dependent on the adaptor protein MyD88 and leads to the activation of the TF NF-κB, which induces expression of proinflammatory cytokines. The other pathway involves the receptor protein TRIF, whose activation induces phosphorylation of the TF IRF3 by TBK1 kinase. The activated IRF3 induces expression of type I interferon, which in turn activates the JAK-STAT signaling pathway, by binding to the type I IFN receptor (IFNR) [23].
Here, we present a large-scale study of short-term changes in histone modifications in mouse DCs during the response to LPS. We focused on the timing of increases in histone modifications at promoters and enhancers, relative to the induction of transcription and to TF binding events. We observed that LPS stimulation induced increased levels of H3K9K14ac, H3K27ac, H3K4me3, and H3K36me3 at LPS-induced promoters and enhancers. Surprisingly, we observed clearly distinct patterns: accumulation of H3K9K14ac was early (between 0.5 and 2 h after stimulation), regardless of the timing of transcriptional induction of genes. Accumulation of H3K36me3 was late and spreads from the 3′ end of gene bodies towards the 5′ end, reaching promoters at later time points (between 8 and 24 h). H3K4me3 accumulation was later than that of H3K9K14ac (between 1 and 4 h) and was more correlated with transcriptional induction times. Integrated analysis with genome-wide binding data for 24 TFs revealed possible associations between increases in H3K9K14ac and H3K4me3 and binding by RelA, Irf1, and especially STAT1/2. LPS-induced accumulation of H3K9K14ac and H3K4me3 was severely impaired in Stat1 −/− cells. Together, these results suggest that stimulus-induced dynamics in a subset of histone modifications reflect the timing of activation of stimulus-dependent TFs, while others are more closely associated with transcriptional activity.

Genome-wide measurement of histone modifications at promoter and enhancer regions
To elucidate the temporal ordering of stimulus-induced changes in transcription and chromatin structure, we performed chromatin immunoprecipitation experiments followed by high-throughput sequencing (ChIP-seq) for the following histone modifications in mouse DCs before and after LPS stimulation: H3K4me1, H3K4me3, H3K9K14ac, H3K9me3, H3K27ac, H3K27me3, H3K36me3, and similarly for Pol2 (Additional file 1: Figure S1), for ten time points (0 h, 0.5 h, 1 h, 2 h, 3 h, 4 h, 6 h, 8 h, 16 h, 24 h). We integrated this data with publicly available whole-genome transcription start site (TSS) data (TSS-seq) [24]. All data originated from the same cell type, treated with the same stimulus, and samples taken at the same time points.
Snapshots of the data for a selection of features at four promoters are shown in Additional file 1: Figure S2.
Using this data collection, we defined 24,416 promoters (based on TSS-seq data and Refseq annotations) and 34,079 enhancers (using H3K4me1 high /H3K4me3 low signals) (see the "Methods" section). For this genome-wide set of promoters and enhancers, we estimated the levels of histone modifications, Pol2 binding, and RNA reads over time (see the "Methods" section).

Epigenetic changes at inducible promoters and their enhancers
Recent studies using the same cell type and stimulus showed that most changes in gene expression patterns were controlled at the transcriptional level, without widespread changes in RNA degradation rates [25,26]. We therefore defined 1413 LPS-induced promoters based on increases in TSS-seq reads after LPS stimulation. Similarly, for both promoters and enhancers, we defined significant increases in histone modifications and Pol2 binding by comparison to pre-stimulation levels. Our analysis suggested that changes were in general rare; only 0.7 to 5.3% of all promoters (Fig. 1a) and 0.2 to 11.0% of all enhancers (Fig. 1b) experienced significant increases in histone modifications and Pol2 binding. However, changes were frequent at LPS-induced promoters, especially for markers of activity such as Pol2 binding, H3K4me3, H3K27ac, and H3K9K14ac, as well as for H3K36me3 (Fig. 1a). For example, while only 957 promoters (out of a total of 24,416 promoters; 3.9%) experienced significant increases in H3K9K14ac, this included 27.6% of the LPS-induced promoters (390 out of 1413 promoters). To a lesser extent, we observed the same tendency at associated enhancers (Fig. 1b). The smaller differences at enhancers are likely to be caused by imperfect assignments of enhancers to LPS-induced promoters (i.e., we naively assigned enhancers to their most proximal promoter). Analysis of an independent ChIP-seq dataset originating from LPS-treated macrophages [6] revealed a high consistency between DCs and macrophages in LPS-induced increases in Pol2 binding, H3K27ac, and H3K4me3 at promoters and enhancers (see Additional file 1: section "Analysis of histone modification changes in LPS-treated macrophages" and Figure S3). The overlap in increases in H3K4me1 at enhancers was lower, though still statistically significant (p < 1e−4, based on 10,000 randomizations), possibly reflecting differences between DC-and macrophage-specific enhancers and the molecular processes that define these cell types.
LPS-induced promoters were less frequently associated with CpG islands (57%) than stably expressed promoters (87%, Additional file 1: Figure S4A) [27]. Non-CpG promoters more frequently had lower basal levels (i.e., levels at 0 h, before stimulation) of activation-associated histone modifications, such as H3K27ac, H3K9K14ac, and H3K4me3, and similarly lower levels of Pol2 binding and pre-stimulation gene expression (Additional file 1: Figure  S4B). This partly explains the higher frequency of significant increases in histone modifications at LPS-induced promoters (Fig. 1a) and the higher fold-induction of genes associated with non-CpG promoters (Additional file 1: Figure S4C).
Previous studies have reported only limited combinatorial complexity between histone modifications, i.e., subsets of modifications are highly correlated in their occurrence [28,29]. In our data too, basal levels of activation markers at promoters and, to a lesser degree, at enhancers were highly correlated (Additional file 1: Figure S5). Stimulus-induced accumulations of histone modifications and Pol2 binding at promoters and enhancers further support this view: increases in H3K9K14ac, H3K4me3, H3K36me3, H3K27ac, Pol2 binding, and transcription often occurred at the same promoters (Fig. 1c). Similarly, increases in H3K9K14ac, H3K27ac, Pol2 binding, and transcription often coincided at enhancer regions (Fig. 1d). In general, activated regions experienced increases in several activation markers.

Several histone modifications are induced at a specific time after stimulation
Previous studies have reported considerable dynamics in histone modifications in response to environmental stimuli (see the "Background" section) based on the analysis of small numbers of time points. Our dataset, however, allows the analysis of the order and timing of changes over an extended time period after stimulation. To this end, we analyzed the induction times of transcription activity, Pol2 binding, and histone modifications.
First, we inferred the transcriptional induction time of the 1413 LPS-induced genes (see the "Methods" section and Fig. 2a). In addition, we defined a set of 772 promoters with highly stable activity over the entire time course.
As a proof of concept, using an independent time series of RNA-seq samples, we confirmed that significant increases in RNAs are seen at LPS-induced promoters in a consistent temporal order (Fig. 2b). For example, at promoters with early induction of transcription initiation (TSS-seq), there was an early induction of mapped RNA-seq reads, while those with later induction have later induction of mapped reads. Stably expressed genes lack induction of mapped RNA reads at their promoter. Significant increases in Pol2 binding were less frequent but followed a similar pattern (Fig. 2c).
However, the accumulation of histone modifications showed more varying patterns ( Fig. 2d-g). Increases in H3K9K14ac were in general early, between 0.5 and 2 h after stimulation (Fig. 2e), although promoters with early induction of transcription (0.5 h, 1 h, 2 h) tended to have early increases in H3K9K14ac (at 0.5 h). Even genes with transcriptional induction between 3 and 6 h had increases in H3K9K14ac between 0.5 and 2 h after stimulation. Therefore, the increases in acetylation for these promoters preceded induction of transcription. shown for the genome-wide set of promoters (green) and for the LPS-induced promoters (orange). Increases in H3K4me3, H3K9K14ac, H3K27ac, H3K36me3, RNA, and Pol2 binding are observed frequently at LPS-induced promoters. Significance of differences was estimated using Fisher's exact test; *p < 1e−4; **p < 1e−6; ***p < 1e−10. b Same as a, for enhancers. c, d Heatmaps indicating the overlap in induction of pairs of features. Colors represent p values (− log 10 ) of Fisher's exact test. White, low overlap; red, high overlap. Plots are shown for promoters (c) and enhancers (d) Significant increases later than 3 h after stimulation were rare. In addition, increases were rare at promoters with late induction (16 h, 24 h) or at stably active promoters.
Increases in H3K4me3 were concentrated between 1 and 4 h after stimulation (Fig. 2f ). In contrast with H3K9K14ac, increases in H3K4me3 were rare at time point 0.5 h. Accumulation of H3K4me3 was frequent at promoters with transcriptional induction between 1 and 4 h, but-in contrast with H3K9K14ac-it was rare at immediate-early promoters.
Finally, H3K36me3 was only induced at later time points (between 8 and 24 h), regardless of transcriptional induction times of promoters (Fig. 2g). In contrast with H3K9K14ac and H3K4me3, H3K36me3 is located within gene bodies and peaks towards their 3′ end (Additional file 1: Figure S6) [30]. Upon stimulation, H3K36me3 gradually accumulated within the gene bodies of LPS-induced genes, spreading towards the 5′ end, and reached the promoter region at the later time points in our time series (Additional file 1: Figure   S6A). Stably expressed genes had on average high basal levels of H3K36me3, with only limited changes over time. However, interestingly, even for stably expressed genes, an accumulation of H3K36me3 was observed towards their 5′ end at time points 16-24 h (Additional file 1: Figure S6B).
Remarkably, the induction times of H3K9K14ac, H3K4me3, and H3K36me3 at promoters did not change depending on their basal levels (Additional file 1: Figure  S7); regardless of their pre-stimulus levels, increases in H3K9K14ac were early, followed by H3K4me3, and H3K36me3 accumulation was late. This might indicate that a common mechanism is regulating these accumulations, regardless of basal levels. No differences in the accumulation times were observed between non-CpG promoters and CpG island-associated promoters (Additional file 1: Figure S8).
Compared to H3K9K14ac, H3K4me3, and H3K36me3, significant increases in H3K27ac appeared to be less frequent at promoters (Fig. 1a), and their timing coincided with the induction of transcription (Fig. 2d). Increases in H3K9me3, H3K27me3, and H3K4me1 were rare at promoters (Additional file 1: Figure S9). The early accumulation of H3K9K14ac, followed by H3K4me3, was confirmed using an independent replicate TSS-seq data and ChIP-seq time series dataset with lower resolution (time points 0 h, 1 h, 2 h, and 4 h; Additional file 1: Figure S10). Although accumulation of both modifications was earlier in the duplicate data than in the original time series, their relative ordering was preserved: Increases in H3K9K14ac at 1 h were more frequent than H3K4me3 in the original data (58% vs 30% of loci), as well as in the replicate data (89% vs 71% of loci). Additional replication was performed using RT-qPCR and ChIP-qPCR measuring RNA, H3K9K14ac, H3K4me3 (see wild-type (WT) data in Additional file 1: Figure S11), and H3K36me3 (Additional file 1: Figure  S12) at the promoters of nine LPS-induced genes. Here too, accumulation of H3K9K14ac and H3K4me3 occurred early, with H3K9K14ac preceding H3K4me3, while accumulation of H3K36me3 occurred later. The early (between 0 and 4 h) timing of LPS-induced increases in H3K27ac and H3K4me3 were further supported by the analysis of an independent ChIP-seq time series dataset (0, 4, and 24 h) originating from LPS-treated macrophages [6] (see Additional file 1: section "Analysis of histone modification changes in LPS-treated macrophages" and Figure S13).

Correlation between LPS-induced TF binding and increases in epigenetic features
To reveal potential regulatory mechanisms underlying the epigenetic changes induced by LPS, we performed an integrative analysis of our histone modification data with TF binding data. For this, we used a publicly available ChIP-seq dataset for 24 TFs with high expression in mouse DCs [31], before and after treatment with LPS (typical time points include 0 h, 0.5 h, 1 h, and 2 h; see the "Methods" section).
Initial analysis confirmed the known widespread binding of promoters by PU.1 and C/EBPβ, and to a lesser degree by IRF4, JUNB, and ATF3 [31] (Additional file 1: Figure  S14A), and the known association between H3K4me1 and binding by PU.1 and C/EBPβ (Additional file 1: Figure  S15A,B) [12,15]. LPS-induced promoters were frequently bound by TFs controlling the response to LPS, such as NF-κB (subunits NFKB1, REL, and RELA) and STAT family members (Additional file 1: Figure S14B).
Focusing on the overlap between LPS-induced TF binding at promoters and enhancers, and accumulation of epigenetic features, we found that binding of promoters by RelA, IRF1, STAT1, and STAT2 was especially associated with increases in H3K9K14ac, H3K4me3, H3K36me3, transcription, and to a lesser degree Pol2 binding and H3K27ac (Fig. 3, left; Fisher's exact test). For example, of the 418 promoter regions that become newly bound by STAT1 after stimulation, 223 (53.3%) experience increases in H3K9K14ac (vs 3.0% of promoters not bound by STAT1; p: 8.3E-205). LPS-induced binding by the same four TFs was also strongly associated with increases in H3K9K14ac and H3K27ac at enhancers (Fig. 3, right). Combinations of these four TFs often bind to the same promoters and enhancers (Additional file 1: Figure S14C,D), and STAT1 functions both as a homodimer or as a heterodimer with STAT2 [32]. LPS-induced TFs, including NF-κB and STAT family members, have been shown to bind preferentially at loci that are pre-bound by PU.1, C/EBPβ, IRF4, JUNB, and ATF3 [31]. Accordingly, histone modifications were also more frequently observed at regions that were pre-bound by these five TFs (Additional file 1: Figure S16).
Together, these results suggest a strong correlation between increases in activation-associated histone modifications and LPS-induced binding by RelA, IRF1, STAT1, and STAT2.

STAT1 and STAT2 binding coincides with accumulation of H3K9K14ac and precedes accumulation of H3K4me3
The relative timing of LPS-induced TF binding events and increases in histone modifications can reflect potential causal relationships. Particularly, many LPS-induced promoters show increases in H3K9K14ac between 0.5 and 2 h after LPS stimulation (Fig. 2e), and we found a strong overlap between increases in H3K9K14ac and binding by STAT1 (Fig. 3). STAT1 is not active before stimulation, and its activity is only induced about 2 h after LPS stimulation [34], resulting in a strong increase in STAT1-bound loci (from 56 STAT1-bound loci at 0 h to 1740 loci at 2 h; Additional file 1: Figure S14B).
We observed a particularly strong coincidence in timing between STAT1 binding and increases in H3K9K14ac (Fig. 4a): genomic regions that become bound by STAT1 at 2 h show a coinciding sharp increase in H3K9K14ac around the STAT1 binding sites. At promoters and enhancers that became bound by STAT1 at 2 h, the induction of H3K9K14ac was particularly frequent (Fig. 4b, c). Of the 407 promoters and 378 enhancers that become bound by STAT1 at 2 h after stimulation, 222 (54%) and 214 (57%) have an increase in H3K9K14ac (vs only 3.0% of promoters and 3.3% of enhancers lacking STAT1 binding). These increases were especially frequent at the 2-h time point (Fig. 4b, c).
Similar to H3K9K14ac, we observed a general increase in H3K4me3 around STAT1 binding sites (Fig. 4d), between 2 to 4 h after stimulation. Accordingly, only 21 STAT1-bound promoters (out of 409; 5.1%) had significant increases between 0.5 to 1 h, but an additional 140 promoters (34%) experienced increases at the following time points (2-4 h; Fig. 4e). As noted above, H3K4me3 was in general absent at enhancers. Similar patterns were observed for enhancers and promoters bound by STAT2 2 h after stimulation (Additional file 1: Figure S17). In contrast, regions bound by RelA (Additional file 1: Figure S18) and IRF1 (Additional file 1: Figure S19) showed increased levels of H3K27ac and to a lesser degree H3K9K14ac at earlier time points. Associations with H3K9K14ac induction after 2 h were weak compared to STAT1/2. Average increases in H3K4me3 at RelA-and IRF1-bound regions were only modest (Additional file 1: Figure S18G-I and S19G-I), suggesting that the association between RelA-and IRF1-binding and H3K4me3 as seen in Fig. 3 is mostly through co-binding at STAT1/2-bound regions. Associations between histone modifications and binding by other TFs were in general weak (not shown; see Fig. 3). No changes were observed in H3K4me1 at STAT1/2-bound regions (Additional file 1: Figure S20A). Although there was a tendency for STAT1/2-bound loci to have increases in H3K27ac, binding seemed to slightly lag behind H3K27ac induction (Additional file 1: Figure S20B). Finally, although STAT1/2-bound regions tended to experience increases in H3K36me3, there was a large time lag between binding and increases in this modification (Additional file 1: Figure S20C). This is also true for other TFs, such as RelA and IRF1, and even PU.1 and C/EBPβ, regardless of the timing of TF binding (Additional file 1: Figure S15C-F).
These results suggest possible causal relationships between STAT1/2 binding and the accumulation of H3K9K14ac and H3K4me3. The specific timing of increases in these modifications might reflect the timing of activation of these TFs, resulting in the recruitment of acetyl transferases and methyl transferases to specific promoter and enhancer regions.  Figure S21A) and confirmed that they were frequently bound by STAT1/2 (Additional file 1: Figure S21B). We observed that promoters of TRIF-dependent and STAT1/2-bound genes frequently had LPS-induced increases in H3K9K14ac and H3K4me3 (Additional file 1: Figure S21C TRIF-dependent and STAT1/2-bound genes (in particular Ifit1 and Rsad2) showed increases in H3K9K14ac and H3K4me3 in WT but not in KO cells (Additional file 1: Figure S11 and section "A Subset of STAT1/2 Target Genes lack Induction of H3K9K14ac and H3K4me3 in Trif −/− , Irf3 −/− , and Ifnar1 −/− cells").
Furthermore, stimulation of WT cells using IFN-β induced expression of Ifit1 and Rsad2, and accumulation of H3K9K14ac and H3K4me3 at their promoters (Additional file 1: Figure S22). In this system, the activation of the IFNR signaling pathway and of STAT1/2, is independent of TRIF. Accordingly, this IFN-β-induced accumulation of H3K9K14ac and H3K4me3 was not affected in Trif −/− cells, further supporting a role for STAT1/ 2 in the control of these modifications at these genes.
Finally, we performed new ChIP-seq analysis of H3K9K14ac and H3K4me3 in WT and Stat1 −/− DCs (Fig. 5). Genomic regions that are bound by STAT1 showed a sharp increase in H3K9K14ac (Fig. 5a) and H3K4me3 (Fig. 5d) in WT cells, reproducing our observations from our first time series data (Fig. 4a, d). However, this increase was completely abrogated in Stat1 −/− cells (Fig. 5a, d). Focusing on promoter sequences, we noted 321 promoters that had increases in H3K9K14ac in WT but not in KO (Fig. 5b). These promoters were frequently bound by STAT1/2, IRF1, and NF-κB in WT cells (Fig. 5c). On the other hand, 184 promoters had increases in H3K9K14ac in the Stat1 −/− cells but not in WT (Fig. 5b). These promoters lack binding by STAT1/2 in WT (Fig. 5c), and the KO-specific increase in H3K9K14ac might be the result of a different set of TFs recruiting histone modifiers to these promoters, in the absence of functional STAT1. One such TF might be HIF1A, which binds a subset of these promoters but not promoters with H3K9K14ac increases in WT (Fig. 5c) and has been reported to be repressed by STAT1 [35]. Similar observations were made for H3K4me3 induction in WT and Stat1 −/− cells (Fig. 5e, f).

Discussion
The concept of active genes being in an open chromatin conformation was introduced several decades ago [36], but the contribution of histone modifications to the control of gene activity remains controversial [17]. In contrast, the contribution of TFs to regulating gene expression is widely recognized [37], and several studies have identified important crosstalk between TFs and histone modifiers in the regulation of the response to immune stimuli [6,7,[38][39][40][41][42]. Nevertheless, our understanding about causal relationships between TF binding, changes in histone modifications, and changes in transcriptional activity of genes in response to stimuli is still lacking.
Analysis of the ordering of events over time can reveal insights into possible causal relationships or independence between them. Here, we have presented an integrative study of the timing and ordering of changes in histone modifications, in function of transcriptional induction in response to an immune stimulus. Our results suggest that, while the dynamics of some histone modifications are closely associated with transcriptional activity, other modifications appear to be induced at specific time frames after stimulation. For a subset of modifications (e.g., H3K9K14ac and H3K36me3), these time frames appear to be non-coincident with the timing of induction of transcription.
In our dataset, we roughly observed three patterns of modifications. The first was early induction of H3K9K14ac, which occurs mainly in the first 2 h after stimulation. A second pattern consisted of increases in H3K4me3 and H3K27ac, roughly coinciding with induction of transcription. Finally, a third pattern consisted of changes in H3K36me3, occurring only around 8-24 h after stimulation. Although H3K4me3 is widely used as a marker for active genes, the functional role of this modification is still unclear. For example, the deletion of Set1, the only H3K4 methyltransferase in yeast, resulted in slower growth than in wild type but otherwise appears to have only limited effects on transcription [19]. Other studies too have reported a lack of a direct effect of H3K4me3 on transcription [20,21]. Several experiments by Cano-Rodriguez and colleagues illustrate that transcription can be transiently induced in the absence of H3K4me3 and that loci-specific induction of H3K4me3 had no or limited effect on transcription [43]. Another study showed that H3K4 methyltransferase Wbp7/MLL4 controls expression of only a small fraction of genes directly [44]. In contrast, fluorescence microscopy experiments have shown that H3K27ac levels can alter Pol2 kinetics by up to 50% [21].
Since the induction of remodeling appears to occur specifically at LPS-induced genes, it is likely that histone modifiers are recruited by one or more LPS-activated TFs to specific target regions in the genome defined by the binding specificity of the TFs. In this model, primary response regulators could control immediate stimulus-induced changes in transcription and histone modifications, while the later "waves" could depend to different degrees on (1) the process of transcription itself, (2) subsequent activation of secondary regulators, and (3) the presence of other histone modifications (Fig. 6). This fits well with our observations for STAT1/2 and the induction of H3K9K14ac and H3K4me3, within specific time frames and mostly restricted to LPS-induced promoters, and the later establishment of H3K36me3. Other studies have reported associations between STAT1 binding and changes in epigenetic markers following environmental stimulation, including the activation of latent enhancers [6] and histone acetylation [40,45]. Moreover, epigenetic priming by histone acetylation through STAT1 binding to promoters and enhancers of Tnf, Il6, and Il12b has been reported, resulting in enhanced TF and Pol2 recruitment after subsequent TLR4 activation [46]. These primed regions were reported to have sustained binding by STAT1 and IRF1 and prolonged associations with CBP/p300 and constitute a stable, stimulus-induced chromatin state. The step-by-step establishment of histone modifications could reflect one way of regulating this process, with combinations of regulators deciding whether a locus will reach a stably active/poised state or whether it will return to the basal inactive state (Fig. 6). Since TFs such as STAT1 are also known to induce gene expression, one might expect the timing of increases in histone modifications to co-occur with induction of expression. However, as we described here, and as supported by the above studies, this is not necessarily the case. Gene expression is known to be regulated by combinations of TFs, and in this study too, we noticed that LPS-activated TFs such as NF-κB, IRF1, and STATs often bound to the same loci (Additional file 1: Figure S14), which were moreover often pre-bound by several other TFs, including PU.1 and C/EBPβ. Discrepancies between timing of expression induction and accumulation of histone modifications could be caused by different requirements for combinatorial binding. This could also explain widely reported "non-functional" TF binding, where TF binding does not seem to affect the activity of nearby genes [47]. Such "non-functional" TF binding might instead trigger changes in histone modifications that remain unnoticed and affect gene activity in more subtle ways.
Although many studies have compared histone modifications before and after stimulation, most lack sufficient time points and resolution to allow analysis of temporal ordering of changes. One recent study in yeast reported results that are partly similar to ours [48]: specific modifications (especially, but not only, acetylation) occur at earlier time frames during the response of yeast to diamide stress and others at later time points. Another study in yeast showed that H3K9ac deposition appeared before the passing of the replication fork during DNA replication, while tri-methylations took more time to be established [49]. Interestingly, in these studies, typical time frames for changes in histone modifications (including H3K36me3) are less than 1 h after stimulation or replication. In contrast, changes in H3K36me3 in our data appeared 8-24 h after stimulation. Thus, time scales of stimulus-induced epigenetic changes in multicellular, higher mammalian systems might be considerably longer. Interestingly, increases in H3K36me3 around 16-24 h often coincide with a decrease in histone acetylation towards pre-stimulation levels at LPS-induced promoters. A study in yeast suggested that H3K36me3 plays a role in the activation of a histone deacetylase [50] and might therefore play a role in the return to a basal state of histone modifications and terminating the response to stimulus.

Conclusions
Our time series ChIP-seq data and analysis present a first genome-wide view of the timing and order of accumulation of histone modifications during a stress response in mammalian immune cells. The stimulus-induced accumulation of H3K9K14ac, H3K4me3, H3K27ac, and H3K36me3 followed distinct patterns over time. Integrative analysis suggests a role for STAT1/2 in triggering increases in H3K9K14ac and H3K4me3 at stimulus-dependent promoters and enhancers. Differences in interactions between histone modifiers, TFs, and the transcriptional machinery are likely causes for the different patterns of dynamics in histone modifications.

Reagents, cells, and mice
Bone marrow cells were prepared from C57BL/6 female mice and were cultured in RPMI 1640 supplemented with 10% of fetal bovine serum under the presence of murine granulocyte/monocyte colony-stimulating factor (GM-CSF, purchased from Peprotech) at the concentration of 10 ng/ mL. Floating cells were harvested as bone marrow-derived dendritic cells (BM-DCs) after 6 days of culture with changing medium every 2 days. The cells were stimulated with LPS (Salmonella minnesota Re595, purchased from Sigma) at the concentration of 100 ng/mL for 0, 0.5, 1, 2, 3, 4, 6, 8, 16, and 24 h and were subjected to RNA extraction or fixation. Murine IFN-β was purchased from Pestka Biomedical Laboratories and was used to stimulate the cells at the concentration of 1 × 10^2 unit/mL. TRIF-, IRF3-, or IFNR-deficient mice have been described previously [51][52][53]. Stat1-deficient mouse was described previously [54].

RNA extraction and RT-qPCR
One million BM-DCs were stimulated with LPS for indicated times and subjected to RNA extraction by using TRIzol (Invitrogen) according to the manufacturer's instruction. RNAs were reverse transcribed by using Reva-Tra Ace (Toyobo). The resulting cDNAs were used for qPCR by using Thunderbird SYBR master mix (Toyobo) and custom primer sets (Additional file 1: Table S1). qPCR was performed by using LightCycler Nano (Roche).

ChIP-qPCR
ChIP was done as above, except 4 × 10^6 cells were used. The resulting ChIP-DNAs were subjected to qPCR similar to the RT-qPCR procedure, using custom primer sets (Additional file 1: Table S2).

Peak calling and processing of ChIP-seq data
For each histone modification and for Pol2 binding data, we aligned reads to the genome and conducted peak calling and further processing as follows.
We mapped sequenced reads of ChIP-seq IP and control (WCE) samples using Bowtie2 (version 2.0.2), using the parameter "very-sensitive," against the mm10 version of the mouse genome [55]. Processing of alignment results, including filtering out low MAPQ alignments (MAPQ score < 30), was performed using samtools [56].
We predicted peaks for each time point using MACS (version 1.4.2) [57], using each IP sample as input and its corresponding WCE sample as control. To improve the detection of both narrow and broad peaks, peak calling was performed using default settings and also using the "nomodel" parameter with "shiftsize" set to 73. Negative control peaks were also predicted in the control sample using the IP sample as reference. Using the predicted peaks and negative control peaks, we set a threshold score corresponding to a false discovery rate (FDR) of 0.01 (number of negative control peaks vs true peaks), for each time point separately. All genomic regions with predicted peaks were collected over all 10 time points, and overlapping peak regions between time points were merged together. Moreover, we merged together peak regions separated by less than 500 bps. This gave us a collection of all genomic regions associated with a peak region in at least one sample of the time series.
In a next step, we counted the number of reads mapped to each region at each time point for both the IP samples and WCE control samples. Using these counts, we performed a read count correction, as described by Lee et al. [58]. Briefly, this method subtracts from the number of IP sample reads aligned to each peak region the expected number of non-specific reads given the number of reads aligned to the region in the corresponding WCE sample. The resulting corrected read count is an estimate of the number of IP reads in a region that would remain if no WCE reads are present [58]. This correction is necessary for the quantitative comparison of ChIP signals over time in the downstream analysis.
Finally, the corrected read counts were converted to reads per kilobase per million reads (RPKM) values (using read counts and the lengths of each region) and normalized using quantile normalization, under the assumption that their genome-wide distribution does not change substantially during each time series. The normalized RPKM values were converted to reads per million reads (ppm) values.

TSS-seq data processing and promoter definition
TSS-seq data for BM-DCs before and after stimulation with LPS was obtained from the study by Liang et al. [24] (DDBJ accession number DRA001234). TSS-seq data reflects transcriptional activity but also allows for the detection of TSSs on a genome-wide scale at a 1 base resolution [59]. Mapping of TSS-seq samples was done using Bowtie2, as for ChIP-seq data. The location (5′ base) of the alignment of TSS-seq reads to the genome indicates the nucleotide at which transcription was started. In many promoters, transcription is initiated preferably at one or a few bases. Because of this particular distribution of TSS-seq reads mapped to the genome, default peak calling approaches cannot be applied. Instead, we used the following scanning window approach for defining regions with significantly high number of aligned TSS-seq reads.
The number of TSS-seq reads mapped to the genome in windows of size 1, 10, 50, 100, 500, and 1000 bases were counted in a strand-specific way, in steps of 1, 1, 5, 10, 50, and 100 bases. As a control, a large number of sequences was randomly selected from the mouse genome and mapped using the same strategy, until an identical number of alignments as in the true data was obtained. For these random regions too, the number of reads was counted using the same scanning window approach. The distribution of actual read counts and control read counts was used to define a FDR-based threshold (FDR: 0.001) for each window size. For overlapping regions with significantly high read counts, the region with the lowest associated FDR was retained.
In order to remove potentially noisy TSSs, we removed TSSs that were located within 3′ UTRs and TSSs located > 50 kb upstream of any known gene. For the remaining TSSs, we used a simple model (see Additional file 1: section Supplementary Methods) (1) to decide the representative TSS location in case a promoter region contained several candidate main TSSs, (2) to remove TSS-seq hits lacking typical features of promoters (e.g., presence of only TSS-seq reads in absence of histone modifications and Pol2 binding), and (3) to decide the main promoter of a gene in case there were multiple candidates. Finally, we obtained 9964 remaining high-confidence TSSs, each assigned to 1 single Refseq gene.
These TSS-seq-based TSSs were supplemented with 14,453 non-overlapping Refseq-based TSSs for all Refseq genes which did not have an assigned high-confidence TSS-seq-based TSS. Most of the genes associated with these TSSs had lower expression in our RNA-seq data (mostly RPKM is 0 or < 1; not shown). Together, TSS-seq-based TSSs and Refseq-based TSSs resulted in a total of 24,416 promoter regions.
CpG-associated promoters were defined as those having a predicted CpG island (from the UCSC Genome Browser Database) in the region − 1 to + 1 kb surrounding the TSS [60]. Other promoters were considered to be non-CpG promoters.

Definition of enhancers
Enhancers were defined based on the signals of H3K4me1 and H3K4me3. First, we collected all genomic regions with significantly high levels of H3K4me1 (see the "Peak calling and processing of ChIP-seq data" section) in at least one of the ten time points. Regions located proximally (< 2 kb distance) to promoter regions and exons were removed, because they are likely to be weak H3K4me1 peaks observed around promoters, as were H3K4me1-positive regions of excessively large size (> 10 kb). Finally, we removed regions with H3K4me1 < H3K4me3 * 5, resulting in 34,072 remaining enhancers.
Enhancers were naively assigned to the nearest promoter (TSS-seq based or Refseq-based) that was < 150 kb separated from it (center-to-center). For 30,448 enhancers (89%), a promoter could be assigned.

Public ChIP-seq data for TFs
Genome-wide binding data (ChIP-seq) is available for mouse DCs before and after stimulation with LPS, for a set of 24 TFs with a known role of importance and/or high expression in DCs [31] (GEO accession number GSE36104). TFs (or TF subunits) included in this dataset are Ahr, ATF3, C/EBPβ, CTCF, E2F1, E2F4, EGR1, EGR2, ETS2, HIF1a, IRF1, IRF2, IRF4, JUNB, MafF, NFKB1, PU.1, Rel, RelA, RelB, RUNX1, STAT1, STAT2, and STAT3. Typically, time points in this data are 0 h, 0.5 h, 1 h, and 2 h following LPS stimulation (some TFs lack one or more time points). We used the ChIP-seq-based peak scores and score threshold as provided by the original study as an indicator of significant TF binding.
Promoters (region − 1 to + 1 kb around TSS) and enhancers (entire enhancer region or region − 1 to + 1 kb around the enhancer center for enhancers < 2 kb in size) were considered to be bound by a TF if they overlapped a ChIP-seq peak with a significantly high peak score. New binding events by a TF at a region were defined as time points with a significantly high score where all previous time points lacked significant binding.

Definition of induction of histone modifications and Pol2 binding
In order to analyze induction times of increases in histone modifications and Pol2 binding, we defined the induction time of a feature as the first time point at which a significant increase was observed compared to its original basal levels (at 0 h). Significant increases were defined using an approach similar to methods such as by DESeq and voom [61,62], which evaluate changes between samples taking into account the expected variance or dispersion in read counts in function of mean read counts. This approach is necessary because regions with low read counts typically experience high fold-changes because of statistical noise in the data. Here, we modified this approach to be applicable to our data (10 time points without replicates; ppm values per promoter/enhancer region).
The values of all histone modifications, Pol2, RNA-seq, and TSS-seq reads (ppms, for each time point) were collected for all promoters (region − 1 to + 1 kb) and enhancers (entire enhancer region or region − 1 to + 1 kb around the enhancer center for enhancers < 2 kb in size). For each feature (all histone modifications and Pol2 binding), we calculated the median and standard deviation in ppm values for each region, over the 10 time points. Dispersion was defined as follow: where d x, f , s x, f , and m x, f represent the dispersion, standard deviation, and median of feature f in region x over the 10 time points of the time series, respectively. Fitting a second-order polynomial function on the log(d x, f ) as a function of log(m x, f ) for all promoter and enhancer regions, we obtained expected dispersion values in function of median ppm value (see for example Additional file 1: Figure S23 for H3K9K14ac). From fitted dispersion values, fitted standard deviation values s x, f, fitted were calculated (see Eq. 1), and 0-h-based Z-scores were calculated as follows: where Z x, f, t is the Z-score of feature f in region x at time point t, and ppm x, f, t is the ppm value of feature f in region x at time point t. A region x was defined to have a significant induction of feature f if there was at least 1 time point t where Z x, f, t ≥ 4. To further exclude low-signal regions, we added this additional threshold: the region should have a ppm value ≥ the 25th percentile of non-0 values in at least 1 time point. For the regions with a significant induction, the induction time was defined as the first time point t where Z x, f, t ≥ 2. We used a similar approach to define LPS-induced promoters using TSS-seq data (see below). For the analysis of induction times of H3K9K14ac, H3K4me3, and H3K36me3 at enhancers in function of their pre-stimulation basal levels (Additional file 1: Figure  S7), we divided promoters into three classes according to their basal levels of each modifications as follows: Promoters lacking a modifications altogether (0 tag reads after correction described above) were considered as one class ("absent"). The remaining promoters were sorted according to their basal level of the modification and were divided into two classes ("low basal level" and "high basal level") containing the same number of promoters.

Definition of LPS-induced promoters and unchanged promoters
LPS-induced promoters were defined using TSS-seq ppm values. LPS-induced promoters should have Z x, TSS − seq, t ≥ 4 for at least 1 time point and have TSS-seq ppm ≥ 1 at at least 1 time point. Only TSS-seq reads aligned in the sense orientation were considered for this (e.g., they should fit the orientation of the associated gene). For each of the thus obtained 1413 LPS-induced promoters, the transcription induction time was defined as the first time point for which Z x, TSS − seq, t ≥ 2 was observed. Unchanged promoters were defined as those promoters having absolute values of Z x, TSS − seq, t < 1 for all time points, leading to 772 promoters.
RNA-seq data processing for wild-type, Trif −/− , and Myd88 −/− cells RNA-seq data for mouse BM-DCs treated with LPS were obtained from the study by Patil et al. [63] (DDBJ accession number DRA001131). This data includes time series data for WT, as well as Trif −/− mice and Myd88 −/− mice.
Mapping of RNA-seq data was performed using TopHat (version 2.0.6) and Bowtie2 (version 2.0.2) [55,64]. Mapped reads were converted to RPKM values [65] using gene annotation data provided by TopHat. RNA-seq data obtained from the Myd88 −/− and Trif −/− mice was processed in the same way. RPKM values were subjected to quantile normalization over all 10 time points.
For genes corresponding to the LPS-induced promoters, the maximum fold-induction was calculated in the WT RNA-seq data. The same was done in the Trif −/− RNA-seq data and in the Myd88 −/− RNA-seq data. TRIF-dependent genes were defined as genes for which the fold-induction was more than five times lower in the Trif −/− data than in WT, leading to 141 TRIF-dependent genes (Additional file 1: Figure S21A). Similarly, 66 MyD88-dependent genes (not shown) were defined as having more than five times lower induction in the Myd88 −/− than in WT.
Duplicate ChIP-seq data and Stat1 −/− data We generated an independent duplicate time series for dendritic cells (DCs) treated with LPS (0, 1, 2, and 4 h), including TSS-seq, and ChIP-seq (H3K9K14ac and H3K4me3) data as described above. Data was processed in the same way as the original time series dataset (see the "Methods" section), and the induction times of H3K9K14ac and H3K4me3 at LPS-induced promoters were estimated. To facilitate the comparison between the duplicate data and the original (longer) time series, we also re-analyzed the original data using only time points 0, 1, 2, and 4 h (the same time points as the duplicate data). Stat1 KO DCs were treated with LPS for 0 or 4 h along with wild-type DCs, and ChIP-seq (H3K9K14ac and H3K4me3) data were obtained as described above.

Fisher's exact test
We used Fisher's exact test to evaluate the significance of differences between induced and non-induced promoters and enhancers (Fig. 1a, b), the significance of associations between changes of pairs of features (Fig. 1c, d), and the association between TF binding and increases in histone modifications, Pol2 binding, and transcription ( Fig. 3 and Additional file 1: Figure S16).