Genome-wide analysis of H4K5 acetylation associated with fear memory in mice

Fear memory induces H4K5ac in the hippocampus in a training-dependent manner

To examine the epigenetic and transcriptional profile of genes associated with memory formation in the hippocampus, we trained adult mice on a CFC paradigm (Figure 1A). We chose CFC because it is a robust, long-lasting learning paradigm in which memory for a context can persist for more than one year after a single training session [27, 28]. Mice were exposed to a novel context in which they received a foot-shock, either once (Day 1) or twice on two consecutive days (Day 1 and Day 2), then tested for fear memory 24 hours later (Day 3). After a single foot-shock, the animals expressed a significant freezing response (47.2 ± 16.5%; p < 0.001) compared to control mice (After Shock, FC Day 1; Figure 1B) that was maintained when tested 24 hours later (47.3 ± 10.0%; p < 0.01) (Test Day 2; Figure 1B). However, with a second training session on day 2, the freezing response was increased further by 20% (67.4 ± 14.2%; p <0.001) when tested 24 hours later (Test Day 3; Figure 1B). In control mice, freezing on days 2 and 3 compared to day 1 (Before Shock) was significant (p < 0.01 and p < 0.05, respectively), but was not significant compared to day 1 (After Shock), which is the measure by which we make all comparisons. It is also worth noting that control mice plateau on day 2 while FC mice continue to have higher freezing.

FC has been associated with transcriptional programs that are activated within 1 hour after conditioning, and that persist for up to 6 hours [3, 29]. Subsequent training, however, may increase gene expression, recruit additional genes to reinforce the memory, or prime existing transcriptional programs for rapid induction of genes for synaptic strengthening. Since memory formation has been associated with histone acetylation in the brain, we examined whether memory performance correlates with higher acetylation levels following additional training sessions. We determined the level of H4K5ac, a PTM recently implicated in gene bookmarking, and increased with FC and object recognition memory tasks [4, 13], following one or two days of CFC. Western blots show that H4K5ac was increased approximately 3-fold in the hippocampus 1 hour after one CFC session. With two conditioning sessions, H4K5ac level was increased 4.6-fold over controls following a memory test on day 3 (Figure 1C), suggesting that H4K5ac induction is proportionate to the amount of training. H4K5ac was examined 1 hour after memory test on day 3 because 1) gene expression is activated within 1 hour following fear conditioning and memory retrieval [30–32], 2) memory is consolidated or reconsolidated within 6 hours [3, 29, 32], 3) histone acetylation decreases to baseline levels within 2–4 hours [4, 33], 4) memory for the context is enhanced by an additional training session, and 5) H4K5ac levels are higher at this time point.

Distribution of H4K5ac across the genome and within genes

Previous studies have shown the association of histone acetylation at promoters of a restricted set of canonical genes involved in memory [4, 9, 13], but to date, genome-wide data are limited. Here, we used ChIP-Seq to determine the distribution of H4K5ac across the genome, followed by de novo identification of genes associated with H4K5ac after CFC (after 2 training sessions) in the mouse hippocampus.

Analysis of H4K5ac distribution showed enrichment of reads in the promoter and coding sequence (CDS) of H4K5ac-ChIP samples compared to IgG-IP samples in both FC (Figure 2A and 2B) and controls (Figure 2D and 2E), an increase of 19% (59 million read) and 17.7% (55 million reads), respectively. The targeted enrichment of H4K5ac to gene bodies is consistent with the proposed role of this PTM in transcriptional regulation. Analysis of H4K5ac in genic regions revealed higher acetylation upstream of the transcription start site (TSS), spanning the CDS and extending down to the transcription termination site (TTS) compared to IgG-IP samples (Figure 2C and 2F). Specifically, there was a prominent peak of H4K5ac in the promoter region approximately 800 bp upstream of the TSS, as well as in the CDS 1 kb downstream of the TSS. H4K5ac distribution was similarly enriched in the control group (Figure 2D), suggesting that learning does not change the overall profile of this PTM in the hippocampus. IgG-IP samples showed low coverage in both groups (Figure 2C and 2F) and, thus, are appropriate input controls for H4K5ac-ChIP sequence reads.

To determine whether the observed profile was specific for H4K5ac, we compared it with H4K12ac, another histone PTM associated with fear memory, from a publicly available dataset [4]. Although H4K5ac and H4K12ac datasets could not be directly compared due to the different CFC training protocols used, the increase of both H4K5ac and H4K12ac immediately following CFC and the higher levels of H4K5ac after two training sessions, suggest that histone acetylation is a consistent marker of memory formation. As with H4K5ac, our analysis of H4K12ac revealed a similar bimodal peak centered at the TSS which was restricted to approximately ± 1 kb relative to the TSS but did not extend into the CDS and TTS as with H4K5ac (Figure 2G). Moreover, H4K12ac had lower enrichment in the promoter than in the CDS, in contrast to H4K5ac, which was largely enriched in the promoter. We were unable to compare H4K12ac controls, as ChIP-Seq controls for sample and experimental conditions for H4K12ac were not available in the public release of this dataset. Together, these data suggest different occupancy and potentially different modes of transcriptional regulation by H4K5ac and H4K12ac following learning [34].

H4K5ac as a marker of actively transcribed genes in the adult hippocampus

We then examined the relationship between H4K5ac and gene transcription using a publicly available whole mouse genome microarray dataset (Agilent) for gene expression immediately after CFC in the mouse hippocampus [4]. We reasoned that because gene expression occurs within 1 hour of both memory consolidation and reconsolidation [3, 29–32, 35], this dataset was appropriate to determine the association between H4K5ac and global gene expression. The 18,023 genes form the expression dataset were ranked by level of expression (from lowest to highest) in FC compared to naïve controls (Figure 3A) and plotted against the average coverage of H4K5ac ± 5 kb relative to the TSS. The level of gene expression was found to correlate to H4K5ac enrichment such that the highest expressed genes had the highest coverage for H4K5ac, while the least expressed genes had the lowest coverage (Figure 3B and 3C; Additional file 1: Figure S3A and 3D). This applied to both groups regardless of training, suggesting that H4K5ac is a general feature of expressed genes. We also confirmed that H4K12ac correlated with the level of gene expression (Figure 3D; Additional file 1: Figure S3B). There was no correlation between gene expression and IgG-IP coverage (Additional file 1: Figure S1A and 1B). These results indicate a clear association between both H4K5ac and H4K12ac and gene expression.

We then identified genes acetylated above average and performed a cross-wise comparison between experimental groups. Based on the average promoter read count of 45 in our dataset, we considered genes with more than 50 reads in the promoter as above average. From a total of 23,235 genes in the dataset, 7,103 genes were identified in the FC group, and 7,708 genes in the control (Figure 3E). Using this criteria, 742 genes (15.1%) were specific for FC, 1,273 genes (21.8%) were specific for control, and 6,029 genes (85% of FC and 78% of control) were common to both groups. We then looked at whether genes with above average H4K5ac after 2 days of CFC were also associated with H4K12ac after one session of CFC. Using an adjusted threshold of 10 reads in promoter due to the lower average coverage, approximately 9 reads in promoter, in the H4K12ac dataset, we identified 4,259 unique genes with above average H4K12ac, of which 2,772 genes (65%) overlapped with genes with above average H4K5ac in FC, and 2,846 genes (67%) with above average H4K5ac in controls (Figure 3E; Additional file 1: Table S1). 2,440 genes overlapped all three groups using this criteria.

The results of these analyses extend our findings that in control conditions most nucleosomes are not only acetylated for H4K5 above the average of all genes, but are also acetylated for H4K12. Interestingly, nearly two-thirds of genes with above average H4K12ac after one session of CFC was found to overlap with above average H4K5ac after 2 days of CFC or context. This suggests that the same set of genes, associated with H4K12ac and induced immediately after CFC, may be upregulated following reinforced training, regardless of the associated histone acetylation used to identify the genes. It also suggests that the same set of genes may be activated after initial learning, during the formation of contextual fear memory, and after memory retrieval, independently of the CFC paradigm.

H4K5ac is associated with both promoter and coding regions

Nucleosome occupancy studies have shown that acetylated and methylated histones are enriched in the promoter of highly expressed genes, but subsequently removed or replaced in the CDS [36–38]. To investigate the positional effect of nucleosomes with H4K5ac on transcription, we clustered genes based on their acetylation profile ± 2 kb relative to the TSS. Five H4K5ac clusters were identified in FC: one in the CDS (cluster 1), one with relatively no enrichment (cluster 2), and three in the promoter (clusters 3, 4, and 5) (Figure 4A). Genes with H4K5ac that feature in either the promoter or the CDS (clusters 1 and 3–5) constituted a larger proportion of highly expressed genes, while genes with relatively no enrichment (cluster 2) accounted for the largest proportion of genes with low expression (Figure 4B). Genes clustered for H4K5ac in controls had profiles and cluster contributions relative to expression comparable to FC (Additional file 1: Figure S2A and 2B). For H4K12ac-clustered genes, we obtained two in the promoter (clusters 2 and 5) and two in the CDS (clusters 1 and 3), which contributed to a greater proportion of highly expressed genes compared to the non-enriched cluster (cluster 4) (Figure 4C and 4D). In contrast, IgG-IP-clustered genes, which were not enriched for H4K5ac, had equal distribution in low, moderate, and highly expressed genes, regardless of training or the histone mark (Figure 4E and 4F; Additional file 1: Figure S2C and 2D). Promoter, CDS, and 3’-UTR-associated genes correlated with H4K5ac and H4K12ac, with and without CFC, but did not correlate with IgG-IP clusters (Additional file 1: S3A and 3E).

These findings suggest that H4K5ac in the promoter and/or CDS may be a feature of highly expressed genes. To validate this observation, we examined the profile of H4K5ac in Sfi1 and Phactr3, two representative genes differentially acetylated for H4K5ac in CFC and involved in cell division in mitotic cells and in memory processes [39, 40], respectively (Additional file 1: Table S2). In Sfi1, Phactr3, and Phactr3 splice variants, H4K5ac was targeted specifically to the CDS (green) (Figure 4G and 4H). For Sfi1, H4K5ac was also highly enriched in the adjacent CDS of Pisd-ps1/3 (blue; Figure 4G), and downstream of the TTS in an intergenic region preceding the CDS of Eif4enif1 (pink; Figure 4G). In contrast, the CDS of Eif4enif1 and Drg1 showed dramatically lower H4K5ac. The overlap of H4K5ac in the CDS of Sfi1 and Pisd-ps1/3 translated to similar expression levels for Sfi1 (15.19; shown as log2 expression and hereafter) and Pisd-ps1/3 (14.72) but not for Eif4enif1 (11.48) or Drg1 (12.44), which had lower enrichment for H4K5ac. For Phactr3, H4K5ac coverage was lower in intergenic and CDS of neighboring genes Zfp931, Sycp2, and Ppp1r3d (pink; Figure 4H). The effect of H4K5ac on gene expression was also clearly evident for Phactr3 (15.07) and neighboring genes, Zfp931 (11.42), Sycp2 (3.97), and Ppp1r3d (11.51), which show lower expression levels. This provides further evidence that the level of H4K5ac enrichment in the CDS is directly proportional to the level of gene transcription.

TF binding sites proximal to the TSS increase the statistical probability of H4K5ac-nucleosome occupancy in the promoter

We next examined whether high levels of gene expression associated with H4K5ac is linked to permissible TF binding. We scanned the promoter region 2 kb upstream of the TSS for conserved TFBS, and computed the percentage of expressed genes with H4K5ac at that position (Figure 5A). For expressed genes, the percentage of acetylated genes was significantly lower across all positions with a consensus TFBS compared to positions without a known TFBS. Unexpressed genes accounted for approximately 20% of genes with H4K5ac. Our assumption is that having a TFBS at a specific position, on average, increases the probability that TF binding occurs at that position relative to a random sequence position in the presence of H4K5ac. To refine our search and identify regions in the promoter where TF binding may affect H4K5ac occupancy, we profiled the coverage of H4K5ac on all genes, on genes with a TFBS at 500 bp, 800 bp or 1100 bp upstream of the TSS, and on genes with no TFBS 100 bp upstream of the TSS (Figure 5B). Using the average coverage of H4K5ac of all genes as baseline, we observed that the presence of a TFBS at position −500 bp or −800 bp, and −1100 bp resulted in modest a reduction in H4K5ac relative to baseline coverage at that position. However, genes with no TFBS upstream of 100 bp resulted in significantly higher H4K5ac in both the promoter and CDS, approximately ±1 kb relative to the TSS.

Based on the increase of H4K5ac coverage in the absence of a TFBS upstream of 100 bp, we focused our analysis in this region, proximal to the TSS. We compared the contribution of acetylated gene clusters (Figure 4A, 4C, 4E; Additional file 1: Figure S2A and 2C) in the presence or absence of a TFBS relative to 150 bp of the TSS: either no TFBS present in the promoter (within 2 kb of the TSS) or no TFBS, one TFBS, or multiple TFBS 150 bp upstream of the TSS (Figure 5C-E; Additional file 1: Figure S4A and 4B). Gene clusters with relatively no enrichment for H4K5ac or H4K12ac constituted a larger proportion of genes regardless of whether a TFBS was present or not (cluster 2, yellow; Figure 5C; cluster 4, blue; Figure 5D). However, in the presence of at least one TFBS within 150 bp of the TSS (no TFBS > 150), the contribution of cluster 4 for H4K5ac in FC (blue, nearest the TSS in the promoter; Figure 4A and 5C), cluster 3 for H4K5ac in control (green, nearest the TSS in the promoter; Additional file 1: Figure S2A and 4A), and cluster 1 for H4K12ac after CFC (red, nearest the TSS in the CDS; Figure 4C and 5D) increased from approximately 10% to 20%, compared to the same clusters when no TFBS was present. To a lesser extent, cluster contribution was also increased in the presence of one TFBS 150 bp upstream of the TSS, but was diminished in the presence of multiple TFBS. These observations provide novel insight into H4K5ac-mediated regulation of gene transcription and support the notion that TF binding and acetylation are mutually exclusive in the promoter [41]. However, H4K5ac is increased when TF binding occurs proximal to the TSS.

The observed increase in acetylation and transcription at proximal TFBS may be attributed to the recruitment of transcriptional machinery including TFs and RNA polymerase II, which is also known to occupy positions near the TSS in actively transcribed genes [42]. Additionally, recent ENCODE studies have shown that a set of TFs is strongly associated to positions proximal to the TSS and that transcriptional initiation is determined by stereotyped TF binding in this region, approximately 100 to 200 bp upstream of the TSS [43, 44]. Acetylated nucleosomes further away in the promoter, greater than 1 kb from the TSS, may either be more strongly bound and less easily displaced by TF binding, or they may be regulatory regions which do not depend on the presence or acetylation of nucleosomes [45]. As expected, IgG-IP control clusters were uniformly proportioned in the presence or absence of a TFBS (Figure 5E; Additional file 1: Figure S4B). Together, these data suggest that since H4K5ac is associated with increased gene expression, enrichment of H4K5ac proximal to the TSS may be a reliable marker of actively transcribed genes.

Genes differentially acetylated for H4K5 are associated with fear memory in the hippocampus

The high percentage of genes with above average H4K5ac in both FC and controls suggest that this modification is important and that it is subject to tight regulation in the context of transcription-dependent memory formation. Using a criteria-based approach, we found that ~15% of genes were uniquely acetylated for H4K5 with CFC (Figure 3E), however, this did not account for differentially acetylated genes. We also found that H4K5ac correlates to global gene expression levels. Thus, to identify specific genes induced by learning and increased H4K5ac levels in the hippocampus, we used a top-down approach – rather than identifying specific genes activated by learning through differential gene expression, we identified highly expressed genes through differential acetylation of H4K5 in FC compared to controls. We used a peak-calling algorithm to scan the genome at 300 bp intervals for differentially acetylated regions between FC and controls. Using model-based analysis of ChIP-Seq (MACS), we obtained consensus coverage of H4K5ac-enriched regions across the mouse genome [46]. Out of 20,238 peaks identified for H4K5ac in FC by MACS, 708 peaks were found −4000 to −2000 bp relative to the TSS, 3,370 peaks were found in the promoter (−2000 to 0 bp), and 1,340 peaks were found in the CDS (0 to +2000 bp). Of these, we identified 241 regions significantly acetylated for H4K5 in FC, 115 of which were associated with gene bodies representing 114 unique genes, and 126 within intergenic regions (Additional file 1: Table S2).

To validate the results obtained with MACS, we repeated the analysis with three other published algorithms for ChIP-Seq analysis, including SICER, EpiChip, and Genomatix NGS analyzer (Additional file 1: Figure S5A – 5D) [47–49]. We performed a cross-wise comparison of genes identified with the algorithms to genes identified using pre-defined criteria, including genes with more than 50 reads in the promoter (Raw H4K5), previously defined as above average, or genes with more than 50 reads in the promoter with CFC but 40 reads or less in controls (Diff H4K5), analogous to algorithm-based differential acetylation (Figure 3E; Additional file 1: Figure S5E; see Methods). Of all genes identified by MACS, approximately 70% overlapped with SICER, the other most widely used algorithm for differential peak finding. Thus, we considered the genes identified by MACS as a reliable and representative gene set to evaluate further.

Genes differentially acetylated for H4K5 in FC are associated with memory processes

In contrast, the 47 genes differentially acetylated for H4K5 in controls were classified into brain processes such as negative regulation of axogenesis, of neurogenesis, and of cell development, but also contributed to normal brain development and neuronal differentiation (Table 1; Additional file 1: Table S3 and 5). Pathway analysis for genes identified in controls showed enrichment for normal neuronal processes such as axon guidance, but also for genes associated with long-term depression, a form of synaptic plasticity typically associated with synaptic weakening (Table 2). The repressive functional categories and pathways enriched in controls suggest that training counteracts these pathways for memory formation. Alternatively, pathways upregulated in controls may be those that are needed to maintain homeostatic processes and basal neuronal functions in the absence of learning.

Animals and contextual fear conditioning

Experiments were conducted using adult C57Bl6/J males (4–7 months old). Mice were housed under standard conditions with a 12 hour reversed light–dark cycle and access to food and water ad libitum. All animals were maintained in accordance with the Federation of Swiss Cantonal Veterinary Office and European Community Council Directive (86/609/EEC) guidelines.

Mice were habituated to the testing room and handled for three days prior to training and testing. They were then trained in a contextual fear-conditioning paradigm using a TSE Fear Conditioning System. The training consisted of a 3 min. exposure to the conditioning context followed by a brief electric shock (0.7 mA for 1s), then left for an additional 3 min. in the conditioning context. Mice that were not re-conditioned were euthanized 1 hour after the initial fear-conditioning session. Mice that were to be further fear-conditioned were trained on the second day and the memory test performed 24 hours later on the third day. Single trial CFC is known to produce a robust, long-lasting memory, however subsequent training has been shown to strengthen the memory and prevent random association of shock with re-exposure [6–10]. Furthermore, as re-exposure to the context on day 3 increased freezing, euthanasia was performed within one hour of the memory test on day 3, but before the 6-hour reconsolidation window and before extinction could take place [3, 29]. The control group was handled and trained in the same manner but without a foot shock.

Comparisons between groups were analyzed by paired student’s t-test or one-way ANOVA with Tukey post hoc analysis, where appropriate. GraphPad Prism was used for statistical analysis and significance was set at *p ≤ 0.05, **p ≤ 0.01, and ***p ≤ 0.001. All data are shown as mean ± SEM.

Nuclear extraction and Western blots

Nuclear protein extraction was performed as previously described with the following modifications [13]. Hippocampi were dissected and homogenized in 100 μl nuclear inhibition buffer (NIB) at pH 7.4 containing 3.75 mM Tris–HCl, 15mM KCl, 3.75 mM NaCl, 250 μM EDTA, 50 μM EGTA, 30% (v/v) glycerol, and 15 mM β-mercaptoethanol, with the addition of 1:200 proteinase inhibitor cocktail (Sigma-Aldrich), 1:500 PMSF (Sigma-Aldrich) and 1:100 phosphatase inhibitor cocktail (Sigma-Aldrich). The structures were then uniformly homogenized with a 22G syringe and centrifuged at 14,000 rpm for 30 min. The supernatant and pellet, containing cytoplasmic and nuclear material, respectively, was separated and resuspended in another 100 μl NIB with appropriate inhibitors. The pellet was re-homogenized with a 26G syringe and centrifuged at 14,000 rpm for 30 min.

15 μg of proteins from nuclear extracts was mixed with 4× LDS sample buffer (NuPage; Invitrogen) and 10% β-Mercaptoethanol to a final volume of 20 μl and loaded on a Novex 4-12% Bis-Tris Gel (NuPage; Invitrogen). Proteins were then transferred onto a nitrocellulose membrane (Bio-Rad), blocked (Rockland IR blocking buffer; Rockland Immunochemicals), and incubated with primary and secondary antibodies. Whole purified histones were run in parallel to confirm histone subunits (ImmunoVision) and Precision Plus protein dual color standards were used to determine molecular weights (Bio-Rad Laboratories). Bands were identified and quantified using an Odyssey IR scanner (LI-COR Biosciences) and the H4K5ac (11 kDa) signal was normalized to β-actin (42 kDa). Primary antibodies used were anti-acetyl H4K5 [72] (1:1000; Millipore) and monoclonal β-actin (1:1000; Sigma-Aldrich); secondary antibodies used were goat-anti-rabbit (IRDye 680 nm; 1:10,000) and goat-anti-mouse (IRDye 800 nm; 1:10,000; LI-COR Biosciences).

Quantitative real-time PCR

Total RNA was extracted from hippocampus using TRIzol reagent and 1 μg of RNA was reverse-transcribed using the SuperScript First Strand Synthesis II system (Invitrogen). Equal amounts of cDNA from each sample were run in duplicate along with an endogenous control, Gapdh, on a Light Cycler 480 (Roche AG). Crossing point (Cp) values, which are more reliable and reproducible than Ct values, were obtained using the second derivative maximum method (Roche AG). Comparative analysis on Cp values was performed and expressed as fold change over the average of controls [73]. Mean and SEM values were obtained for each and analyzed using two-tailed paired t-tests to determine statistical significance (p <0.05). Oligonucleotides used for quantitative real-time PCR are listed in (Additional file 1: Table S6).

Chromatin immunoprecipitation

Chromatin immunoprecipitation (ChIP) was performed as previously described [13], with the following modifications. Briefly, three hippocampal samples for each group were individually cross-linked with 1% formaldehyde, quenched with 0.125 M glycine, and spun down at 1500 rpm for 5 min at 4°C. To isolate chromatin, samples were washed and homogenized in 2 ml cell lysis buffer containing proteinase and phosphatase inhibitors with a Dounce homogenizer. Samples were centrifuged at 4000 rpm for 5 min. and homogenized again in 1 ml nuclear lysis buffer with inhibitors. DNA was sheared using a Baendelin Sono Plus to a fragment length of 600–800 bp. Total genomic DNA (input) was quantified and 80 μg of chromatin from each sample was immunoprecipitated overnight at 4°C with either 5 μl of anti-acetyl-H4K5 (07–327; Millipore) or 5 μl of IgG (17–685; Millipore) as a negative control. After incubation, nucleosome complexes were isolated with 60 μl of protein A agarose/salmon sperm DNA slurry (Millipore) for 1 h at 4°C. The complexes were washed and dissociated from the beads by incubation in 1% SDS in TE and nuclear lysis buffer at 65°C for 10 min. Histones were then digested with proteinase K for 1 h at 45°C and the DNA was finally extracted with phenol/chloroform/isoamyl alcohol and ethanol precipitation. DNA concentrations were measured on a Nanodrop (Thermo Fisher Scientific) and further verified on a Qubit fluorometer (Invitrogen). Uniformity of fragment size and quality control was validated on a 2100 BioAnalyzer (Agilent Technologies).

ChIP-Seq library preparation

Library preparation was according to recommended guidelines (Life Technologies). From both ChIP and input control samples, 200 ng of DNA was further sonicated at 4°C to a mean fragment size of between 100 to 150 bp using the Covaris S2 sonicator. The DNA was then end-repaired using end-polishing enzymes such that damaged DNA with protruding 5’ or 3’ ends were blunt-ended and phosphorylated. Following repair, the samples were purified using a column purification kit and the blunt ends were ligated with 1 μl of multiplex adaptors. The ligated samples were then nick translated and amplified according to the SOLiD Fragment Library Barcoding protocol and column purified separately. The libraries were then quantitated using a Qubit fluorometer. 20 μl of each library was size-selected for ligation products of 170–230 bp using 2% E-gels and pooled following gel purification. Finally, equimolar amounts of each barcoded library were mixed together before ePCR followed by sequencing.

SOLiD sequencing and mapping statistics

Sequencing was performed on an Applied Biosystems SOLiD 3 platform. Image acquisition and base calling was automated on the SOLiD Instrument Control Software system. The color space reads were mapped and aligned to the current assembly of the mouse genome (mm9, NCBI Build 37 dating July 2007; UCSC genome browser) using the mapping tool of the Bioscope v1.2.1 software suite (Applied Biosystems). Only reads with a maximum of 4 failed color calls and quality values larger than 8 were considered for contiguous mapping. The reads were mapped allowing a maximum of 6 color mismatches and reads with up to 10 mappings on the genome were reported in a SAM file. This file was used for subsequent identification of enriched regions. Sequence data from this study has been submitted to NCBI Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo/) and assigned the identifier (accession no. GSE30325).

From a total of 309 million (309,346,614) 50-bp ChIP-seq reads, 230 million (229,838,436; 74.3%) were uniquely mapped to the current mouse reference genome with a mismatch allowance of ≤ 6 per 50 consecutive bases (Additional file 1: Table S1). The total number of sequenced reads was equivalent to ~6.2 complete mouse genomes (15.5 Gb), while the mappable reads were equivalent to ~4.6 genomes (11.5 Gb). We obtained an average of ~45 reads per promoter region, 783 and 894 reads per CDS for FC and control, respectively, with lower read counts for mock IgG-immunoprecipitated (IgG-IP) control samples (12 and 10 reads per promoter, 287 and 245 reads per CDS for FC and control, respectively) (Additional file 1: Table S1). An equivalent H4K12ac ChIP-seq dataset from Peleg et al. was obtained from Galaxy-Central (sm1186088) at <main.g2.bx.psu.edu/u/fischerlab/h/sm1186088> and re-analyzed using our workflow. With the H4K12ac dataset, we obtained 5.53 million total reads, of which 4.04 million were unique reads (73.1%) with an average coverage of 8.7 reads per promoter and 123 reads per CDS (Additional file 1: Table S1). The higher sequence coverage of H4K5ac in control, ~13.3% more mapped reads (4.9 million reads) compared to FC, may account for the larger number of genes identified in control with our exclusion criteria (>50 reads). The lower coverage (5.53 million reads in total) in H4K12ac may also explain the smaller percentage of genes found to overlap with H4K5ac.

Differential peak calling and data mining analysis

Peak finding was performed using a Model-based Analysis of ChIP-Seq (MACS, version 1.3.7.1 - Oktoberfest) algorithm [46]. To determine genes differentially enriched for H4K5ac in the respective groups, we ran MACS on fear-conditioned against non-fear-conditioned control and vice versa. H4K5ac peaks were identified in MACS with the following parameters: effective genome size = 1.87e+09, tag size = 50, bandwidth = 300, m-fold = 4, and P-value cutoff = 1.00e-5. We also used the Statistical model for the Identification of chip-Enriched Regions (SICER, version 1.1) to call differentially acetylated peaks between groups [47]. We used the following parameters for SICER: redundancy threshold = 1, window size = 200, fragment size = 150, effective genome fraction = 0.7, gap size = 400, FDR = 1.00e-3, and filtered post-analysis for genes with P-value = 1.00e-5. We further compared results to the Genomatix NGS analyzer with Auto-Claverie algorithm with the following parameters: window size = 100 and P-value = 0.05, filtered post-analysis for genes with P-value = 1.00e-5 [49]. EpiChip analysis was performed according to standard protocols (epichip.sourceforge.net), except gene scoring was performed ± 5000 from the 5’ start position [48]. H4K12ac ChIP-Seq data, by CFC in young mice, was obtained from the public repository at Galaxy-Central (sm1186088 at <main.g2.bx.psu.edu/u/fischerlab/h/sm1186088>). Control ChIP-Seq data for H4K12ac, for sample or experimental condition, was not available [4].

Gene ontology and pathway analysis

To determine functional gene enrichment and interaction networks of genes differentially acetylated in fear-conditioned compared to non-fear-conditioned controls, we used the genes identified in MACS for functional annotation. From the 241 differentially acetylated regions identified in fear-conditioned over control, 115 unique peaks were associated in the promoter or coding region of genes. From the 77 differentially acetylated regions identified in control over fear-conditioned, 42 unique peaks were associated with gene bodies. We used The Database for Annotation, Visualization and Integrated Discovery (DAVID, v6.7; http://david.abcc.ncifcrf.gov) for the analysis of functionally-enriched genes in our respective gene lists [74]. Settings were set at a count threshold of 2 and EASE score of 0.1, a more conservative test than Fisher’s Exact test. We also used Web-based Gene Set Analysis Toolkit V2 (WebGestalt; http://bioinfo.vanderbilt.edu/webgestalt) for the analysis of functionally-enriched genes in our respective gene lists [75]. Genes were analyzed using a hypergeometric test with multiple adjustment using the method of Benjamini & Hochberg and categorized into their respective classes or pathway associations based on the Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.genome.jp/kegg/kegg2.html).

Gene expression analysis

Gene expression data was obtained from the Gene Expression Omnibus repository at NCBI (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE20270), and processed and analyzed with R/Bioconductor. Spot intensities were normalized using quantile normalization. For the comparison with the acetylation levels, the genes were subdivided in ten equally sized bins according to the average expression.

Transcription factor binding site analysis

Known transcription factor binding sites (TFBS) were downloaded from the CisRED database [76]. A total of 223,000 binding sites was used to analyze whether the presence of a known TFBS at a given position in the promoter determines the acetylation level at that position. Genes were divided into expressed and unexpressed genes, and expressed genes were further subdivided into two groups depending on whether a TFBS was annotated at that position. For each group we computed the percentage of genes acetylated at position × in step widths of 10, from 0 to 2000 bp upstream of the TSS.

Acetylation profile clustering

We computed acetylation profiles in the ± 2 kb region around the transcription start site (TSS) and used k-means clustering to subdivide the profiles into 5 clusters. We subsequently built cross-tabulation tables to check whether cluster membership correlates with the expression level and/or with the presence of a known TFBS in certain regions. Clusters were generated in an unsupervised fashion and correlation between acetylation scores and gene expression was computed using Spearman's rank correlation.