Time makes histone H3 modifications drift in mouse liver

Global patterns of histone mark profiles identified age and diet groups

Livers were obtained from C57BL/6 female mice of 3, 6 and 12 months of age fed with standard diet (SD) and 12 months of age fed with 30% caloric restriction (CR). We immunoprecipitated and obtained histone profiles for H3K4me3, H3K27me3, H3K27ac and H3K4me1 for a total amount of 108 sequences. In particular, we produced 30 profiles for the group SD 3m (7 H3K4me3, 7 H3K27me3, 8 H3K27ac, 8 H3K4me1), 20 for SD 6m (3 H3K4me3, 6 H3K27me3, 6 H3K27ac, 5 H3K4me1), 29 for SD 12m (8 H3K4me3, 8 H3K27me3, 7 H3K27ac, 6 H3K4me1) and 29 for CR 12m (5 H3K4me3, 6 H3K27me3, 9 H3K27ac, 9 H3K4me1) as summarized in the Supplementary Table 1.

With the initial purpose to globally characterize similarities and differences between all the 108 histone modification profiles, considering age and diet as combined factors, a dimensionality reduction technique (Uniform Manifold Approximation and Projection, UMAP) was used to analyze the normalized, genome-wide signal of each respective histone mark in consecutive, non-overlapping bins of 10 kb. The UMAP resulting from that matrix is shown in Figure 1A and 1B. Each dot represents one of the 108 ChIP-seq samples and its relative position in the map considers the signal of all the 10 kb bins of the corresponding animal and histone mark. At a first comparison, including all profiles (Figure 1A), samples clustered largely by histone mark. The heterogeneity among individuals of each group, estimated by the sum of standard deviations in both dimensions, is greater in H3K27me3 (1.65) and H3K27ac (1.37) compared to the other two histone marks, with the H3K4me3 (0.79) profiles being the most homogeneous group (Figure 1A). The comparison between profiles obtained by a single histone mark clearly distinguished younger (blue dots) from older groups (yellow and red dots) in the cases of the H3K4me3 and H3K27me3, but not for H3K4me1 or H3K27ac samples (Figure 1B). Notably, the 6 months old groups (green dots) tend to be located in between the young and old samples, supporting the notion that a temporal trend may exist. From this global analysis, CR profiles (yellow dots) are not distinguishable from the aged-matched SD counterparts (red dots) except for the case of H3K27ac which shows samples separated by diet.

Figure 1. Similarity between all genome-wide histone modifications profiles represented through dimensional reduction using the UMAP algorithm (A), samples are colored by histone modification and aggregate primarily by histone mark. (B), amplification of previous UMAP with samples split by histone mark and colored by the experimental group they belong to (different age and diet).

Then, we calculated Spearman’s correlation coefficients between each pair of samples, based on the same genome-wide signal in consecutive, non-overlapping bins of 10 kb, and produced the correlation heatmap between each pair-wise comparison of samples (Figure 2A and Supplementary Figure 1).

Figure 2. Genome-wide correlation patterns between histone modifications by group. (A) Heatmap showing the Spearman’s correlation coefficient for each combination of experimental groups after merging the replicates. (B) Box plots of Spearman’s correlation coefficients of each replicate between different groups of histone modifications. Two-sided Wilcoxon test p-value are shown.

Consistent with the proximity seen in first component of the UMAP, the correlation analyses demonstrate that H3K4me1 and H3K27ac groups of profiles were highly similar and are hierarchically clustered into a separate group. Remarkably, unsupervised hierarchical clustering separated the young and old samples (3 months old mice vs. older mice) of both H3K4me3 and H3K27me3 samples. The numbers above the horizontal lines indicate p-values from a two-sided Wilcoxon test between the respective groups. Correlation between H3K4me3 (Figure 2B, left panel) and H3K27me3 (Figure 2B, second panel from left) with H3K4me1/H3K27ac profiles decreases with age, unaffected by diet. Correlation between H3K4me3 and H3K27me3 (Figure 2B, third panel from left) increases with age, with a smaller effect in the CR 12m group, while correlation between enhancer-related histone modifications (Figure 2B, right panel) remains high in all groups and does not present an age- or diet-related trend (p-value, two-sided Wilcoxon test, are shown above each comparison).

Overall, these findings show that changes occur with time in marking histone H3, in particular with K4me3 and K27me3, in the liver of mice. On a genome-wide scale, promoter-related histone mark profiles and enhancer-related histone mark profiles tend to become more similar with age. Moreover, the promoter-activating and -repressive histone mark profiles become more dissimilar with age, and this phenomenon is partially prevented by caloric restriction (Figure 2B first two panels on the left).

H3K4me3 and H3K27me3 signal at TSS changes with time

Histone modifications H3K4me3 and H3K27me3 are enriched at the transcription start site (TSS), mostly in a mutually exclusive manner (except for bivalent promoters).

To further describe the genomic regions targeted by the time-dependent H3K4me3 and H3K27me3 changes, we focused the analysis on the TSS, comparing the shape of the mean signal of these two histone marks in all 24,418 annotated TSS regions in groups of mice of different age. This analysis revealed consistently stronger (higher amplitude) and wider peaks of H3K4me3 and H3K27me3 in the TSS regions of young mice compared to the older ones (Figure 3A).

Figure 3. Analysis of promoter-related histone marks in transcription start sites. (A) H3K4me3 and H3K27me3 profiles around all transcription start sites. (B) Density plot of the H3K4me3 and H3K27me3 signal in all TSS, SD 12m over SD 3m. (C) Distribution of H3K4me3-to-H3K27me3 signal ratio across all TSS.

Comparing the distribution of H3K4me3 signal between mice aged 3 and 12 months, a slightly decreased signal was observed across most TSS (Figure 3B, left panel). In contrast, the H3K27me3 signal distribution (Figure 3B, right panel) differed more among the two age groups with a subset of TSS having less H3K27me3 in old compared to young mice. We also investigated bivalency of H3K4me3 and H3K27me3 across all annotated TSS regions in young and older samples. As shown in Figure 3C, the ratio of H3K4me3-to-H3K27me3 signal revealed an almost unimodal distribution in 3 months old animals which, with age, assumes a more bimodal distribution in 12 months old mice and an intermediate shape in the 6 months old mice. This transition indicates that the proportion of putative bivalent promoters, having a flexible state, based on the presence of both H3K4me3 and H3K27me3 signal, decreases with time.

Both age and diet affect chromatin state transitions

To determine the functional elements which are most affected by the observed time-dependent changes in histone modifications, we performed an integrated chromatin state analysis using all histone mark profiles, taking advantage of the multivariate model approach offered by the ChromHMM tool [27]. In this analysis, chromatin states are defined as recurring combinations of the four studied histone marks in the collected ChIP-seq profiles. Each region of the genome in each of the experimental groups (in consecutive bins of 200 bp) is then assigned one of the chromatin states. As expected, given the histone modifications we investigated, multiple TSS- (state 1 and 2) and enhancer-related (states 3–5) states were readily identified (Figure 4A, left panel). These states are characterized by different combinations of H3K4me3, H3K27ac, H3K4me1, and the absence of H3K27me3. Instead, chromatin state 6 is represented only by the repressive histone modification H3K27me3. State 8 contained none of the histone marks analyzed in this study, and, peculiarly, state 7 presents a new signature of weak H3K4me3 and minimal H3K27me3 signal in the absence of the other two histone marks. Annotation of the identified chromatin states was done based on the combinations of histone marks using related ENCODE studies as a reference [27]. With the exception of chromatin state 7, all other states had been previously detected and annotated in the ENCODE studies.

Figure 4. Chromatin state definition and transitions between states with age. (A) Chromatin state definition, labels and enrichment of selected groups of functional elements in young samples. Number in the heatmap on the left are probabilities of observing a given histone modification in the respective chromatin state. Number in the heatmap on the right are enrichment values obtained by ChromHMM for each combination of chromatin state and feature. (B) alluvial plot showing the transition between chromatin states in the SD 3m (center), SD 12m (left), and CR 12m (right) groups. Bar height represents genome coverage in percent (all states together are 100%).

We validated the learned model by checking which known genomic regions and functional elements are enriched in each chromatin state (Figure 4A, right panel and Supplementary Figure 2). For example, chromatin state 1 presents the classic signature of active TSS and is coherently enriched for TSS regions (RefSeq) as well as CPG islands. Chromatin state 2, previously identified in the ENCODE project as “flanking TSS regions” is enriched in TSS regions +/− 2 kb, transcription end sites (TES), and RefSeq gene bodies. As expected, the known histone modification signature of enhancers (state 3) is enriched in enhancers of the adult liver [28]. Consistently, enhancers identified as active only in the embryonic liver [28] show a stronger enrichment for weak and poised enhancer signatures (states 4 and 5). Apart from embryonic enhancers, also simple repeats show an enrichment in state 5. Chromatin state 6, which is characterized by a distinct H3K27me3 signal, is coherently enriched for Polycomb-repressed chromatin and is the most enriched state on chromosome X. As TSS regions are not enriched in chromatin state 7, even though it contains faint signals of both H3K4me3 and H3K27me3, it is likely not the signature of bivalent promoters. Instead, it is enriched in long terminal repeat retrotransposons (LTR), long interspersed nuclear elements (LINE), and constitutively lamina-associated domains (cLAD) [29]. These functional regions also show enrichment in chromatin state 8 which is characterized by the absence of all four histone marks and covers most of the genome (roughly 50–70%). This finding is consistent with the previously mentioned model from the ENCODE project where the empty chromatin state represents the largest portion of the genome. Notably, the genome column in Figure 4A (right heatmap) represents the genome coverage of each state. As expected, TSS-associated chromatin states 1 and 2 are found only in a very small fraction of the genome. Instead, about 20% of the genome is covered by enhancer-related chromatin states 3, 4 and 5. While state 6 accounts for roughly 10% of the genome, the large majority (more than 60%) is in state 8.

Then, we studied the dynamics by which chromatin states transformed into each other at different age or diverse diet regimen. An alluvial plot helped us to follow these transitions between chromatin states in the SD 3m, SD 12m, and CR 12m groups (Figure 4B). The bar height represents genome coverage in percent (all states together are 100%). Many, but not all, chromatin state transitions observed with age, especially the loss of quiescent chromatin (state 8), are not present or present at a smaller degree in the comparison with the CR 12m group. State 7 shows the biggest fold change between young (SD 3m) and old (SD 12m) mice kept in standard diet. Therefore, we labeled it “age-related state”. Mice kept in caloric restriction show only a slightly elevated genome coverage of state 7 compared to young mice. A similar trend can be observed for state 6, however with smaller changes between young and old SD mice and instead a bigger change between young and old CR samples.

In Supplementary Table 2 we report the percentage of the genome of older groups (6 months and 12 months in both SD and CR) which was found to change chromatin state compared to the youngest group (SD 3m). In all groups, most of the genome (around 68–79%) does not undergo state changes. However, in standard diet, we observed a progressive increase of the portion of the genome which has changed chromatin state, from about 20% in the first three months to more than 30% at 12 months. Remarkably, mice kept in caloric restriction do not show the same extent of state transition (24.2%). Furthermore, we found that genomic regions which were annotated as quiescent chromatin (state 8) in 3 months old mice are assigned either state 4 (active enhancer 2), 6 (repressed Polycomb) or 7 (Repeats/cLAD) in 12 months aged mice depending on the diet (Figure 4B). This means that going from the absence of all four histone marks, the old samples show slight presence of either H3K27me3, H3K4me3 or H3K27ac in the respective genomic regions.

To validate the potential role in aging of the genes whose histones H3 appeared targeted by the observed time- or diet-dependent changes, we used a peak calling tool and differential binding analysis. Ingenuity Pathways Analysis of the differentially H3K4me4-marked genes revealed several pathways involving these gene sets (Supplementary Figure 3). However, no clear gene signatures typical of aging pathways emerged from these gene sets. On the opposite, for the H3K27me3, much less loci showed ages or diet effects (Supplementary Figure 4). We have also investigated presumed active (positive for H3K27ac) or primed (barely positive for H3K27ac and strong for H3K4me1) enhancer regions and enhancer regions identified as previously characterized in the mouse [28] without revealing any different temporal trend or CR effect.

Mice

All aspects of animal care and experimentation were performed in accordance with the Guide for the Care and Use of Laboratory Animals published by the US National Institutes of Health (NIH Publication No. 85-23, revised 1996) and the Italian Laws (D.L.vo 116/92and following additions), which enforces EU86/609 Directive. Experiments were approved by the local Ethical Committee, Organismo Preposto al Benessere degli Animali (OPBA), and notified to the Ministry of Health.

Chromatin was extracted from livers collected from inbred C57BL/6 female mice of 3, 6 and 12 months of age fed standard diet (SD) and 12 months of age fed 30% caloric restriction (CR) (Supplementary Table 1). At 3 months of age, the mice were randomly divided in the CR and SD groups and started the CR treatment. To calculate the amount of food in order to provide 30% less calories to the mice with respect to the ad libitum SD, food intake was measured daily in a subset of mice from 3 to 6 months of age [37].

All 3 months old females were fertile and plateaued at above 20 grams of body weight. Survival of these females was consistently 100% in the first year. The rare occurrence of injuries, deviant behavior and any external or internal, as determined upon necroscopy, signs of abnormalities including minimal alopecia, dermatitis or histiocytosis, determined the exclusion of the individual from the study. Histological examination of the livers used to extract chromatin did not reveal structural differences in infiltrating cells or residue of hematopoiesis. Altogether, 40 mice were bred and 31 livers collected from 3 months (n = 8), 6 months (n = 6), 12 months old mice fed SD (n = 8) or CR (n = 9), and fixed and embedded in paraffin blocks.

Chromatin immunoprecipitation and sequencing

H3K4me3, H3K27me3, H3K4me1 and H3K27ac profiles were obtained from the aforementioned samples using a modified version of the protocol for extraction of chromatin from formalin-fixed paraffin-embedded (FFPE) tissue samples [38]. Quality control of all ChIPs were performed before sequencing by real-time QPCR, measuring the enrichment observed of bound with respect to the input DNA, on liver transcriptionally active and inactive promoter regions or enhancer regions. Sequencing was performed on an Illumina HiSeq 2000, 50 bp read length, single-stranded.

Processing of FASTQ files and data analysis

Reads were aligned to the mm10 reference genome using “bwa aln” (v0.6.2-r126). Unmapped reads (reads with a MAPQ smaller than 1), reads that mapped outside of chr1-19 and chrX, as well as duplicate reads were removed using samtools (v0.1.18). Data, if not specified otherwise, was processed inside a Singularity (v2.6.0) container using R (v3.5.1). For some analysis, the BAM files of the replicates of each group were merged using bamtools (v2.5.1), and subsequently indexed using samtools (v1.7). BigWig files for each sample and group were generated from the respective BAM files using deepTools bamCoverage (v3.1.3) with a bin size of 10 bp, BPM-normalization (chrX was ignored for normalization), and reads extended to 200 bp.

Genome-wide signal for UMAP and correlation: for each sample and group, the genome-wide signal was retrieved for consecutive, non-overlapping bins of 10 kb using deepTools multiBigwigSummary (v3.1.3). Uniform Manifold Approximation and Projection (UMAP) was calculated using the genome-wide signal of each sample in bins of 10 kb using the ”uwot” R package. Spearman’s correlation coefficients between samples and groups were calculated using deepTools plotCorrelation (v3.1.3) based on the genome-wide signal of 10 kb bins.

The heatmaps were generated using the “ComplexHeatmap” (v1.20.0) and “circlize” (v0.4.5) R packages. Statistical significance shown in the box plots were generated using the “ggsignif” (v0.4.0) R package using a Wilcoxon test.

TSS profiles, signal and ratio distribution and enhancer signal

Annotation of transcription start sites (TSS) was taken from the curated RefSeq set, downloaded from UCSC Genome Browser on 01.03.2018. The list of TSS includes the TSS of the longest transcript for each gene. The signal around TSS was calculated as mean signal in bins of 50 bp, with a range of +/− 5 kb around the TSS, using deepTools “computeMatrix” (v3.1.3). Missing data was treated as zero. The output was then plotted (for profile) or summed up by TSS to generate the overall signal. For the ratio, summed signal values for H3K4me3 signal was divided by the H3K27me3 signal for each individual TSS. Relative levels of H3K4me3 signal on selected TSS regions obtained upon ChIP-seq from comparing 3 months with respect to 12 months old samples were validated by investigating gene expression fold by reverse transcriptase-QPCR on same samples (Supplementary Figure 5).

Chromatin state analysis

Chromatin state analysis was performed using ChromHMM (v1.17). Annotation of TSS represents the same list of TSS as above with the range adjusted to +/− 2 kb. Similarly, annotation of transcription end sites (TES) was taken from the curated RefSeq set, downloaded from UCSC Genome Browser on 01.03.2018. The list of TES includes the TES of the longest transcript for each gene. Annotation of chrX was taken from the mm10.txt file provided with ChromHMM. Long terminal repeats and long interspersed nuclear elements annotation was downloaded from the UCSC Genome Browser for mm9 and then converted to mm10 using the UCSC Genome Browser LiftOver feature. Annotation of constitutively lamina-associated domains (cLAD) was taken from the report of Peric-Hupkes et al. [39] and converted from mm9 to mm10 using the UCSC Genome Browser LiftOver feature (mm9 files can be found in the GSE17051 data set in the GEO database). Merged BAM files of each group were binarized using the “BinarizeBam” command. Then, the chromatin state model was trained using the “LearnModel” command for 8 outcome states. States were then reordered using the “Reorder” command, then segmentation and enrichment with sets of genomic locations were generated using the “MakeSegmentation” and “OverlapEnrichment” commands. For the transition analysis, BED files of all groups containing the chromatin state for every location in the genome were cut into bins of 200 bp using awk (v4.0.2) and bedtools (v2.26.0). Alluvial plot was generated with the “ggalluvial” R package (v0.9.1).

More information about chromatin states can be found at http://compbio.mit.edu/ChromHMM/ and https://www.encodeproject.org/publications/aa02932b-39c5-43d8-9ee1-2059fb604a81/.

Data availability

ChIP-seq data are accessible on GEO repository (https://www.ncbi.nlm.nih.gov/geo/) with accession number GSE179090. All the data that support the figures and the other findings are available from the authors upon request.