Time makes histone H 3 modifications drift

To disclose the epigenetic drift of time passing, we determined the genome-wide distributions of monoand tri-methylated lysine 4 and acetylated and tri-methylated lysine 27 of histone H3 in the livers of healthy 3, 6 and 12 months old C57BL/6 mice. The comparison of different age profiles of histone H3 marks revealed global redistribution of histone H3 modifications with time, in particular in intergenic regions and near transcription start sites, as well as altered correlation between the profiles of different histone modifications. Moreover, feeding mice with caloric restriction diet, a treatment known to retard aging, preserved younger state of histone H3 in these genomic regions.

signal of each respective histone mark in consecutive, non-overlapping bins of 10 kb. The UMAP resulting from that matrix is shown in Fig. 1A and B. Each dot represents one of the 108 ChIP-seq samples and its relative position in the map takes into account the signal of all the 10 kb bins of the corresponding animal and histone mark. At a first comparison, including all profiles (Fig. 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 (Fig. 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 H3K27m3, but not for H3K4me1 or H3K27ac samples (Fig. 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. Similar results were obtained when using other bin sizes (2 kb, 50 kb), as well as principal component analysis dimensional reduction algorithm (results not shown). a b 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, samples are split by histone mark and colored by the experimental group they belong to (different age and diet).
Then, we calculated Spearman 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 ( Fig. 2a and Fig. S1). 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, unbiased hierarchical clustering separated the young and old samples (  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. 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 27,517 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 (Fig. 3a).
Comparing the distribution of H3K4me3 signal between mice aged 3 and 12 months, a slightly decreased signal was observed across most TSS (Fig. 3b, left panel). In contrast, the H3K27me3 signal distribution (Fig. 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 Fig. 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.
a b c 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.

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 25 . 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 (Fig. 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 25 . With the exception of chromatin state 7, all other states had been previously detected and annotated in the ENCODE studies.
We validated the learned model by checking which known genomic regions and functional elements are enriched in each chromatin state (Fig. 4a, right panel). 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 (not shown). 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 (not shown). As expected, the known histone modification signature of enhancers (state 3) is enriched in enhancers of the adult liver 26 . Consistently, enhancers identified as active only in the embryonic liver 26 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 (not shown). 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) 27 . 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 Fig. 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 (Fig. 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 "agerelated 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 table. S2 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 (23.6%). 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). 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.

Discussion
Changes in covalent modifications of the histone tails as well as in histone variant incorporation and in DNA methylation have been extensively studied during the development and aging of different organisms 28 . Based on these studies, chromatin remodeling, which enables the epigenetic reprogramming of tissue functions, has been suggested to occur throughout lifespan 22 . However, purely time-dependent changes (i.e. in the absence of developmental processes or aging-associated compliances) of histone modifications have not yet been systematically assessed in mammals.
Here, we have characterized the dynamic patterns of four key histone H3 modifications in the liver of healthy C57BL/6 mice within a year of age to determine whether, and if so, to what extent, time impacts chromatin state. In the C57BL/6 strain of laboratory mice, the time period following maturity, starting at 2-3 months of age, is characterized by approximately 6-8 months (depending on the individual) of healthy conditions showing no major signs of frailty. Different approaches were followed to describe the changes in histone modifications within each experimental group and between the groups to identify age-and CR-specific effects. To compare genome-wide histone modifications profiles, we employed dimensional reduction and correlation methods revealing global differences between mice aged 3 and 12 months, most strikingly for the tri-methylation of H3K4 and H3K27. Moreover, correlation between transcription-related (H3K4me3, H3K27me3) and enhancer-related (H3K27ac, H3K4me1) histone modification profiles significantly decreases while correlation between the two transcription-related histone modifications increases with age. The comparisons we implemented of the different histone H3 modification profiles obtained at 3, 6 and 12 months of age indicate that few months of life were sufficient to identify global changes of K4 and K27 histone H3 trimethylation. These observations highlight the potential to use genome-wide histone modification profiles to estimate the age of a mouse, similar to the DNA methylation clock proposed by Steve Horvath in 2013 29 . However, these changes were not affected by CR, which instead was found to prevent the loss of H3K4me3 around transcription start sites occurring with increasing age. Consistent with what has been reported about the effect of CR on proteins 30 and histone acetylation 31,32 in mammals, our analysis revealed that CR particularly affects the genomewide acetylation patterns of H3K27.
The average signal at the TSS of both H3K4 and H3K27 trimethylation showed an almost linear temporal decay and the width of the H3K4me3 peak at TSS collapsed from 3 to 6 months of age. Moreover, at the TSS, the H3K27me3 distribution changed from unimodal at 3 months to a bimodal pattern at 12 months, suggesting an ongoing process of activating/deactivating genes beyond the developmental period. The mechanisms underlying the vast changes of H3 modifications we observed are not clear. The width of H3K4me3 peaks around the TSS region has been linked to transcriptional consistency 33 , while changes in H3K4 and H3K27 tri-methylation bivalency may affect the stability of the promoter landscape. Regardless of their potential as marker of epigenetic aging, some of the changes observed over time do not occur in mice kept in CR. The annotation of genomic regions to different chromatin states and investigation of transition patterns between chromatin states among the different age and diet groups revealed that a consistent fraction (20-30%) of the genome changed chromatin state over time. At domain levels, no major differences characterized the groups except for a decreased fraction of poised enhancers. On the contrary, states assigned to quiescent chromatin underwent a massive transition in aged mice which was not observed in mice following the CR diet. The unmarked regions in 3 months old mice were actually lost in older mice kept in SD diet, revealing, together with increased poised enhancer domain, a trend of increased repression and H3K4 tri-methylation on repetitive regions.
In conclusion, while our data do not establish that the observed changes in H3 modification are causally involved in aging, they shed light on the age-related epigenetic changes in mammals. Furthermore, they unveil chromatin remodeling in regions of the genome, most of which remain to be functionally understood.

Methods
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) (Tab. S1). 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 34 . 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 genomewide signal of each sample in bins of 10 kb using the "uwot" R package. Spearman 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.
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. 35 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).