Body weight and high‐fat diet are associated with epigenetic aging in female members of the BXD murine family

Abstract DNA methylation (DNAm) is shaped by genetic and environmental factors and modulated by aging. Here, we examine interrelations between epigenetic aging, body weight (BW), and life span in 12 isogenic strains from the BXD family of mice that exhibit over twofold variation in longevity. Genome‐wide DNAm was assayed in 70 liver specimens from predominantly female cases, 6–25 months old, that were maintained on normal chow or high‐fat diet (HFD). We defined subsets of CpG regions associated with age, BW at young adulthood, and strain‐by‐diet‐dependent life span. These age‐associated differentially methylated CpG regions (age‐DMRs) featured distinct genomic characteristics, with DNAm gains over time occurring in sites such as promoters and exons that have high CpG density and low average methylation. CpG regions associated with BW were enriched in introns, tended to have lower methylation in mice with higher BW, and were inversely correlated with gene expression (i.e., higher mRNA levels in mice with higher BW). CpG regions associated with life span were linked to genes involved in life span modulation, including the telomerase reverse transcriptase gene, Tert, which had both lower methylation and higher expression in long‐lived strains. An epigenetic clock defined from age‐DMRs revealed accelerated aging in mice belonging to strains with shorter life spans. Both higher BW and the HFD were associated with accelerated epigenetic aging. Our results highlight the age‐accelerating effect of heavier BW. Furthermore, we demonstrate that the measure of epigenetic aging derived from age‐DMRs can predict genotype and diet‐induced differences in life span among female BXD members.


| INTRODUC TI ON
The "epigenetic clock" based on DNA methylation (DNAm) has emerged as a widely used biomarker of aging that surpasses telomere length assays in its accuracy and utility (Breitling et al., 2016;Marioni et al., 2018). Often referred to as DNAm age (DNAmAge), the CpGbased estimator of biological age comes in a few different versions for both humans and mice (Hannum et al., 2013;Horvath, 2013;Levine et al., 2018;Petkovich et al., 2017;Stubbs et al., 2017;Thompson et al., 2018;Wang et al., 2017). All these clocks share a common feature-they rely on the methylation status of preselected subsets of CpGs that are each assigned weights and are used collectively to estimate age. A critical question has been: are these DNAmAge clocks detecting changes that are purely a function of time, and therefore, correlates of chronological age? Or are they providing a measure of the intrinsic pace of biological aging that can be related to health, fitness, and life expectancy? Evidence from human epidemiological studies indicates that some versions of the clock perform well at predicting life expectancy, and a younger DNAmAge relative to chronological age implies decelerated biological aging, and is associated with a lower risk of disease and increased longevity (Chen et al., 2016;Levine et al., 2018;Lu et al., 2019;Marioni, Shah, McRae, Chen, et al., 2015;Marioni, Shah, McRae, Ritchie, et al., 2015).
However, to date, most DNAmAge estimations in mice have been modeled on the canonical C57BL/6J (B6) background strain (Cole et al., 2017;Petkovich et al., 2017;Stubbs et al., 2017;Thompson et al., 2018;Wang et al., 2017). Similar to interindividual variation in humans, aging trajectories vary considerably among different mouse strains, and common DNA variants contribute to the pace of normal aging and longevity (Yuan et al., 2009). Whether the differential rates of epigenetic aging can discern normative life span differences between mouse strains remains an open question. Furthermore, if body weight reduction due to CR or dwarfing mutations can slow down the clock, then another question is: can subtle differences in normal body weight also have an impact on DNAmAge?
The BXD recombinant inbred strains are the most deeply phenotyped mouse genetic reference panel, and have a long history in aging and longevity research (De Haan & Van Zant, 1999;Gelman, Watson, Bronson, & Yunis, 1988;Haan, Gelman, Watson, Yunis, & Van Zant, 1998). Notably, the different BXD sibling strains exhibit wide variation in life expectancy. Some strains have a mean life expectancy of less than 15 months (e.g., BXD13, BXD5), while others typically live well over two or even 3 years (e.g., BXD19, BXD65) (Gelman et al., 1988;Haan et al., 1998;Lang et al., 2010;Roy et al., 2019). Longevity data in the BXDs have been collected since the 1980s. Life span data continue to be collected from an enlarged family that now consists of 150 sets of isogenic siblings (Ashbrook et al., 2019;Gelman et al., 1988;Roy et al., 2019). The BXDs were derived by crossing two parental strains, B6 and DBA/2J (D2), and then inbreeding the progeny (Ashbrook et al., 2019). Genomes of the BXDs therefore represent random recombinations of the B6 and D2 genomes, and each strain is a unique mosaic of homozygous B6 or D2 genotypes. The D2 strain is considered to have a more accelerated aging profile, and has consistently shorter life span than B6 (Yuan, Peters, & Paigen, 2011;Yuan et al., 2009). Other age-associated parameters include rapid thymic involution (Hsu, Li, Zhang, & Mountz, 2005), quicker replicative senescence of hematopoietic stem cells (De Haan & Van Zant, 1999), and increased tail tendon breakage in D2 compared to B6 (Sloane et al., 2011). Due to random assortment of gene variants, the progeny BXDs have a greater range of variation in life span and aging traits (De Haan & Van Zant, 1999;Roy et al., 2019), and provide a unique resource with which to dissect the interrelations between epigenetic aging and longevity.
Here, we leveraged the extensive longevity data generated for the BXDs  to evaluate associations between body weight, epigenetic aging, and life span. We used affinity-capture enrichment with the methyl-CpG-binding domain protein (MBD), followed by deep sequencing (MBD-seq) to profile the aging liver methylome of 12 BXD strains (Aberg, Chan, Xie, Shabalin, & van den Oord, 2018). To examine the impact of a common metabolic stressor, we also quantified the methylome on a subset of BXD cases maintained on a high-fat diet (HFD), a diet that decreased life span by as much as 12% . Our goal was to chart differentially methylated regions related to aging, and to strain differences in body weight and life span, and to examine correlation with gene expression. We computed a simple DNAm clock using the age-dependent CpG regions, and this revealed strain differences in rates of epigenetic aging. Overall, our results show that both higher body weight and HFD can have an accelerating effect on the epigenetic age, and this can then be related to life span.

| BXD strain selection and life span and body weight characteristics
The present work is based on data collected from two separate cohorts of BXD mice. Life span data were collected from a group of female BXDs, referred to as the "longevity cohort," that were kept either on ad libitum standard chow (control diet or CD) or on HFD, and allowed to age until mortality. The BXD strains exhibited a wide range in natural life span, and HFD reduced overall longevity. Details on this longevity cohort are reported in Roy et al. (2019). Based on the life span data, 12 members of the BXD panel (including F1 hybrids) were selected for DNAm assays as they were representative of the wide variation in life span. This included related sub-strains with differences in life span (BXD65/BXD65b, and BXD73/BXD73b; Table 1 and Figure 1a). Figure 1a plots the ages at natural death for mice in the longevity cohort, and sample sizes for life span determination for these selected strains range from 6 to 22. For five of these, we included mice maintained on HFD as these strains showed significant reduction in life span (specifically, BXD65 on HFD; Figure 1a) or were related sub-strains with variable response to HFD (e.g., BXD48 had slight life span reduction, but BXD48a appeared unaffected by HFD; Table 1). Based on life span averages, each strain-by-diet group was classified as short-lived (mean life span <600 days), mediumlived (600-750 days), and long-lived (>800 days; no strain on HFD were in this group). Note that a strain classified as long-lived on CD may be classified as short-lived on HFD (Table 1).
A parallel, identically treated and strain-matched group, the "biospecimen cohort," was maintained for tissue collection at different ages, and liver samples were obtained from this. We selected 70 liver specimens for the corresponding 17 strain-by-diet groups.
These were chosen so that distribution of age at time of tissue collection was closely matched across the life span groups (Figure 1b; individual-level sample information in Table S1). We note that aside from three male cases for BXD102, B6D2F1, and D2B6F1 (Table S1), all liver specimens were from females. Mice were initially weighed at young adulthood (mean age of 134 ± 81 days), and at this stage, age was a significant correlate of weight ( Figure 1c). We refer to this as baseline body weight or BW0. The mice assigned to HFD were introduced to the diet at the time of initial weighing, with the exception of four that were introduced 90 days later (see Figure 1d). Mice spent 73-614 days on HFD (Table S1) before final weighing on the day of sample collection, by which time, the group on HFD had become significantly heavier to the matched strains on CD (HFD = 41 ± 12 g vs. CD = 26 ± 6 g, p < 0.0001; n = 34). For the HFD group, the final weight was uncorrelated with the number of days on HFD. The weight of the liver on HFD was slightly heavier but the effect was not statistically significant, likely due to the modest sample size for HFD (HFD = 1.29 ± 0.23 g vs. CD = 1.22 ± 0.23 g; p = 0.37, n = 34).
There was significant strain variation in BW0 even after adjusting for age, and the long-lived F1s in particular, had higher body weights, and this hybrid vigor was apparent with or without the male cases ( Figure 1e). The final weight of mice also showed significant strain variation ( Figure 1e). Liver weight appeared fairly consistent across strains. By final weighing, age was no longer a significant correlate of body weight (r = 0.01) or liver weight (r = 0.12). Instead, BW0 remained a significant predictor of final weight in both the CD and HFD mice (Figure 1f). We performed a multivariable regression to evaluate the degree to which the final body weight was predicted by BW0, diet, strain, and chronological age. This showed that the strongest predictors of final weight were diet (F (1,55) = 61, p = <0.0001), followed by baseline weight (F (1,55) = 7, p = 0.009), and strain (F (11,55) = 3, p = 0.0095), but not the final age of mice.    (Table S3).

| Strain-dependent patterns in global features of the methylome
We started with a principal component analysis (PCA). A plot of the top principal components, PC1 and PC2 (captured 19% and 13% F I G U R E 1 Age distribution, and life span and body weight characteristics. (a) Each point depicts a longevity cohort mouse that was allowed to age till mortality either on control diet (CD; black circles) or high-fat diet (HFD; red crosses). There was a total of 225 mice, and mean sample size was 13 (ranging from 6 to 22) per strain-by-diet. (b) Each point depicts a biospecimen cohort mouse used for tissue harvest and methylome assay. Age distribution (y-axis) is uniform across the three life span groups. (c) At baseline, before introduction to HFD, age was significantly correlated with body weight. (d) Individual trajectory of change in body weight from baseline to final age is plotted with each mouse colored by strain. Mice in the HFD group were introduced to the diet at the time of initial weight measurement, with the exception of 4 mice belonging to BXD48a, BXD65, and BXD73b (black arrows) that were placed on HFD at age ~145 (entry to HFD marked by black rectangles for the 4 mice). (e) The bar plots show average body weight at young adulthood (original baseline weight, and residual values after adjustment for age), and final body weight. Error bars are standard error. (f) The baseline body weight (x-axis) was a significant predictor of individual body weight at older age (y-axis); r = 0.53 (p < 0.0001) for all mice; r = 0.77 (p < 0.0001) CD group; r = 0.53 (p = 0.04) for 15 HFD group. (g) Baseline weight was negatively correlated with mean life span (y-axis) for the strain-by-diet groups. However, when separated by diet group, this inverse correlation was significant only for the CD group; r = -0.27 (p = 0.04) for all mice; r = -0.31 (p = 0.04) for the CD group; ns for HFD groups of the variance, respectively), showed clustering of samples by strain identity, irrespective of diet ( Figure 2a). The one exception was a BXD65 on HFD that plotted away from the BXD65 cluster, and due to its questionable identity, it was excluded from downstream analyses, and remaining statistical tests were done in 69 samples.
Sub-strains (e.g., BXD73/BXD73b; BXD65/BXD65b) also clustered in close proximity. Unsupervised hierarchical clustering confirmed the clustering of samples by strain identity rather than age or diet groups (QC plots in Figure S1).
The top 5 PCs collectively explained 58% of the variance (Table   S1), and we examined whether these were associated with age, diet, body and liver weights, and strain life span. Age was not a significant correlate of any of the 5 PCs. For strains with matched CD and HFD samples, the PCs did not differentiate between the diets. For the weight measurements, PC1 was a significant negative correlate of BW0 (Figure 2b), and final body and liver weights (correlations with and without F1s in Table S4). When partitioned by diet, the negative correlations remained significant only for the CD group. None of the other PCs were associated with the weight measurements. For the strain longevity data, PC4 (8% of variance) was the strongest correlate of mean, median, minimum, and maximum life span ( Figure 2c; Table S4). PC2 and PC3 (11% of variance) were also significantly correlated with maximum life span (Table S4). These correlations between PCs and life span remained significant when restricted to the CD group.
We computed the overall mean methylation for genic regions (i.e., CpG regions that overlap annotated gene features), and inter-  (Table S4).
Taken together, the methylome-wide analysis showed that the PCs capture strain-dependent differences in overall mean methylation. These PCs were also associated with body weight, and strainlevel life span, but not age and diet.

| Characterizing differentially methylated CpG regions
We next applied an epigenome-wide association study (EWAS) approach to identify CpG regions that were associated with age, BW0, and strain median life span using a multiple linear mixed model. Age was a significant predictor of site-specific DNAm (Figure 3a), and F I G U R E 2 Global features of the methylome. (a) Scatter plot between the top 2 principal components-PC1 (19% of variance) and PC2 (13% of variance)shows a strong population structure with mice clustering by strain identity (color coded). For strains with cases on both standard chow (CD; solid circles) and high-fat diet (HFD; squares), the HFD and CD samples co-cluster. (b) Body weight at young adulthood has a significant negative correlation with PC1 (r = -0.3, p = 0.02, n = 69). When partitioned by diet, the negative correlation is significant only for the CD (r = -0.30, p = 0.05). For HFD group, the negative correlation does not reach statistical significance (r = -0.35, p = 0.22). (c) The methylome in turn may be predictive of life span, and PC4 is strongly correlated with life span (r = 0.85, p = <0.0001, n = 69). When partitioned by diet group, the correlation is significant for the CD (r = 0.90; p < 0.0001), but not for the HFD group (r = 0.35; 26 regions, covering 319 CpGs, were above the 10% Bonferroni threshold (unadjusted p ≤ 2.7 × 10 −7 ). These strong age-dependent differentially methylation regions (age-DMRs) were mostly located within genes, and 25 of the 26 bins were associated with increased methylation with age (age hypermethylation; Table S5). Although

BW0 and the life span traits had strong associations with global
DNAm patterns, after partly accounting for the strain-dependent effects with the mixed model, only seven CpG regions were significantly associated with BW0 at the 10% Bonferroni threshold ( Figure 3b). All these BW0 associated differentially methylation regions (BW0-DMRs) had lower methylation among mice with higher BW0 (negative regression estimates; Table S6). At the same Bonferroni threshold, only three CpG regions were associated with strain-level median life span and are referred to as LS-DMRs ( Figure 3c; Table S7).
Given the non-independence of adjacent CpG regions, we then applied a lenient threshold of uncorrected p ≤ 1.0 × 10 −4 to define the general characteristics of the sites associated with age, BW0, and life span. In total, 306 CpG regions were age-DMRs at this suggestive threshold (Figure 3a; Table S5). Of these, 57% were age-hypermethylated (Table 2). Compared to the background set of 368,300 CpG bins, the age-DMRs were highly enriched in genic regions such as promoters and exons, and CpG island, and depleted in intergenic regions (enrichment and depletion p-values in Table S3). For each age-DMR, we computed the average methylation and CpG density, and compared these to the age regression coefficients, which convey the change as a function of age. The regression estimates appeared to be highly dependent on local genomic characteristics, and the most pronounced changes involved age hypermethylation (positive regression coefficients) in bins with high CpG density ( Figure 3d) and F I G U R E 3 Features of differentially methylated CpGs regions (DMRs). Each point in the Manhattan plot represents the location of a CpG region (x-axis: autosomal chromosomes 1-19, and chromosome X as 20), and the association -log 10 p (y-axis) for (a) age effect, (b) body weight at young adulthood (BW0), and (c) median life span (LS). The genome-wide significant threshold is set at -log 10 (2.7e−7) (red line; 10% Bonferroni threshold for 368,300 tests) and the suggestive threshold at -log 10 (1.0e−4) (blue line). For the age-DMRs, the regression coefficients for age (i.e., change in DNA methylation per unit change in age in days, log 10 scale), and whether a site gained (positive coefficients, burgundy) or lost methylation (negative coefficients; sandy brown) was highly dependent on (d) the CpG density, and (e) mean methylation. For BW0-DMRs, the body weight associated coefficients also showed correlations with the CpG counts within the bin (f), and sites that were negatively associated with BW0 (blue) had lower mean methylation than the sites that were positively associated with BW0 (red) (g). For the LS-DMRs, whether a site was positively (red) or negatively (blue) associated with strain median life span was only modestly dependent on (h) the CpG counts, and not dependent on (i) the mean methylation levels lower average methylation (Figure 3e). In contrast, age-DMRs that lost methylation with age (age-hypomethylated; negative regression coefficients) featured lower CpG density, and higher average methylation (Figure 3d,e). Functional annotations of these regions identified significant enrichment in gene sets involved in establishment of cellular polarity (e.g., Wnt5a, Lrrd1, Ptk7) (Table S8). We consulted the human GWAS catalog to identify genes represented by the age-DMRs that have been significantly associated with human aging and longevity (Table S5), and this identified only one gene, Cux2 (Buniello et al., 2019).
At the suggestive threshold, 689 CpG regions were classified as BW0-DMRs, and the majority of these (517 or 75%) were negatively correlated with BW0 (Table 2). Introns were the most enriched gene feature (Tables S3 and S6). Compared to the age-DMRs, the BW0-DMRs were in regions with relatively lower CpG density and had only 3 regions within CpG islands (Figure 3f). The regions that were negatively correlated with BW0 in particular featured lower CpG densities and also lower mean methylation levels ( Figure 3g).
For life span, there were 124 LS-DMRs (Table S7), and a slight majority of these (59%, mostly intergenic sites) were positively correlated with higher median life span ( Table 2)

| Relating differentially methylated regions to gene expression
Fifty-two of the samples with liver MBD-seq also had matched liver RNA-seq data, and we used this group to examine whether the DMRs were associated with gene expression. We linked the age-, BW0-, and LS-DMRs to the corresponding transcript. Of the 306 age-DMRs, 279 were paired with mRNA. For these pairs, we tested how many of the DMRs were cis-correlated with gene expression, and how many of the transcripts were also correlated with age of mice (Table S9). At an uncorrected p ≤ 0.05 (|r| ≥ 0.27), 79 age-DMRs were correlated with the expression of cognate genes (Figure 4a).
The age-hypermethylated DMRs mostly showed positive correlations with gene expression such that the transcripts also showed increased expression with age (e.g., Jak3, Amn, Tradd). The age-hypomethylated DMRs were more likely to be negatively correlated with gene expression, and this set of transcripts also showed increased expression with age (Nfkbia, Slit3).
For the BW0-DMRs, 614 of the CpG regions were paired to corresponding transcripts. At p ≤ 0.05 (|r| ≥ 0.27), 121 BW0-DMRs were correlated with gene expression (Table S10). The overall pattern indicated that the DMRs that had lower methylation in mice with higher BW0 (i.e., negative regression estimates) tended to be negatively correlated with gene expression, and there was higher expression of these genes in mice with higher BW0 (Figure 4b). This included Mc5r, Nfkb1, and Tcf4. The few DMRs that had positive associations with BW0 were more likely to be positively correlated with gene expression (e.g., Aldh18a1, Madd, Rap2a).  lived group had significantly lower methylation levels compared to both the medium-and short-lived groups (F 2,66 = 7.67, p = 0.001), and the long-lived group had significantly higher gene expression compared to the short-lived group (F 2,49 = 3,31, p = 0.04; Figure 4d).

| Measure of age acceleration from ageassociated CpG regions
We next evaluated whether we can derive a measure of differential rates of aging from the age-DMRs that could be predictive of life span differences. We first summarized the age-dependent changes by computing the weighted averages for the set of 306 age-DMRs.
The weighted averages, as expected, had a strong correlation with the chronological age of mice, and for a more direct comparison, the values were scaled to the age range for the 69 samples (Table S1).
We refer to this DMR-based age estimate as DMRmAge, and this showed a nearly linear correlation with chronological age at r = 0.88 (p < 0.0001, n = 69; Figure 5a). Age acceleration was derived as the residuals from the regression of DMRmAge on chronological age (Horvath, 2013;Thompson et al., 2018), with positive values indicating an older or accelerated biological age, and negative values indicating a decelerated biological age, and we refer to this as DMRmAge-acc (Table S1).
We then treated the DNAmAge-acc as the outcome variable, and used multivariable regression to evaluate the association with body weight (BW0 or final weight), diet, and the strain maximum life span.
In the model with BW0 as an explanatory variable, the DNAmAge-acc To verify that these associations are robust, we repeated the EWAS using a randomly selected subset of 55 female methylomes (these are identified in Table S1). With the reduced sample size, 237 CpG regions were age-DMRs at p ≤ 1.0 × 10 −4 . Using these age-DMRs and respective regression coefficients from the subsample analysis, we re-calculate the DMRmAge and DMRmAge-acc (Table   S1), and evaluated whether the age acceleration will show consistent associations with life span, body weight, and diet in the 55 subsamples, and also in the 14 excluded samples (the 3 male samples were assigned to this test set). The DMRmAge was significantly correlated with chronological age in both sets with r = 0.83 in the 55 subsamples, and r = 0.92 in the test set of 14 samples ( Figure S3a). The DMRmAge-acc defined from this continued to show inverse correlation with life span, and this was highly significant in the test set ( Figure S3b). Notably, the two BXD102 (long-lived strain) samples in (a) (d) the test set, including the male BXD102, had the most decelerated clocks among the 14 (Table S1). Both BW0 and the final body weight were positively correlated with the DMRmAge-acc, although this did not reach statistical significance in either sets. Four HFD mice had been randomly assigned to the test set, and all four had positive DMRmAge-acc. In the subsampled set, the HFD group had significantly higher age acceleration compared to the strain-matched CD mice (60.49 ± 33.05 for HFD, 3.87 ± 55.40 for CD, p = 0.008, n = 25; Figure S3c). Similarly, HFD was associated with accelerated epigenetic aging in the strain-matched samples in the test set (43.57 ± 19.83 for HFD, -12.02 ± 19.98 for CD, p = 0.008, n = 8; Figure S2c).
Taken together, the analysis demonstrates three points. First, that the age-DMR-based estimates of age acceleration are predictive of strain-dependent differences in life span among female BXDs; second, that higher body weight, even before introduction to HFD, is associated with more accelerated aging; and third, HFD, which results in strain-dependent increase in body weight, significantly accelerates aging, as measured by DMRmAge-acc.

| Genomic features of differentially methylated regions
In the present study, we parsed the variance in the methylome and defined site-specific methylation differences that may be attributed to a strain-level phenotype, median life span, and two individuallevel variables-age and BW0. For the age-DMRs, the time-dependent patterns were consistent with previous reports (Ciccarone, Tagliatesta, Caiafa, & Zampieri, 2018;Sziraki et al., 2018). As in Sziraki et al. (2018), methylation loss over time occurred mostly in regions with higher average DNAm, and methylation gains occurred mostly in regions with lower average DNAm. Since DNAm was quantified over 150 bp non-overlapping bins, we were also able to relate the differential methylation patterns to the local CpG density. In particular, for the age-hypermethylated DMRs, the increase in DNAm with aging had a strong positive correlation with CpG density. This too is consistent with reports that CpG dense regions-a feature of CpG islands, which typically remain unmethylated-are the sites that tend to gain methylation with age (Rakyan et al., 2010). The genes represented by the age-DMRs included a few notable members such as Cyp46a1 and Abca7, which are involved in cholesterol metabolism and implicated in Alzheimer's disease (Carter, 2007), and few members of the WNT signaling and mesenchyme developmental pathways (e.g., Fzd1, Fzd8, Wnt5a, Jak3, Ptk7, Nrp2). The current data replicated the CpG islands in C1ql3 and Ptk7, which we previously reported as age-hypermethylated sites in the BXD parental strains, B6 and D2 (Mozhui & Pandey, 2017). A region-based functional annotation analysis revealed a significant enrichment in genes involved in cell polarity, an aspect of cells that is established during development through interaction with the WNT polarity signaling pathway, and which becomes dysregulated during aging (Berger, Wodarz, & Borchers, 2017;Budovsky, Fraifeld, & Aronov, 2011).
Body weight, even at young adulthood and before introduction to HFD, was significantly associated with the global methylome patterns and was also predictive of strain longevity. We therefore included BW0 as a predictor variable in the multiple regression model.
The BW0-DMRs included a few genes that have been previously associated with body weight in human GWAS at p < 5 × 10 −8 (Buniello et al., 2019). This included an intron DMR in the fat mass and obesity associated Fto gene, which plays a key role in energy homeostasis, and consistently shown to influence body weight in humans (Zhou, Simmons, Lai, Hambly, & McLachlan, 2017). The BW0-DMRs were mostly located within genes, particularly introns, and a striking aspect of the BW0-DMRs was that 75% were negatively associated with BW0. For the BW0-DMRs that were linked to the corresponding transcript, the overall pattern indicated a negative correlation between DNAm and gene expression. This means that while DNAm at these sites were, in general, inversely correlated with body weight, the transcript levels of the corresponding genes generally had higher expression in mice that were heavier. been linked to aging and longevity by human GWAS Teumer et al., 2016).

| Building clocks from age-dependent CpG regions
Several different versions of the DNAmAge estimator are now available for both mice and humans (Hannum et al., 2013;Horvath, 2013;Petkovich et al., 2017;Stubbs et al., 2017;Thompson et al., 2018;Wang et al., 2017). The standard protocol for developing DNAmAge clocks starts by applying a regression algorithm in a training dataset, followed by age estimation in validation cohorts to gauge the accuracy of the clock. Our goal in this study was not to develop another DNAmAge clock. Given the sample size of the present study, that would not have been a feasible pursuit. Instead, our goal was to test whether the age-dependent CpG regions could provide an estimate of epigenetic aging that could discern life span differences between the mouse strains, and that could be related to body weight and diet. For this, we simply summarized the age-DMRs by computing the weighted averages for each sample. The weighted averages, unsurprisingly, correlated strongly with the chronological age of mice. More importantly, the age acceleration derived from the age-DMR-based clock was (a) inversely correlated with the strain life span phenotype, (b) was significantly more accelerated in the HFD group, and (c) was positively correlated with body weight. This conveys that the age-DMRs can estimate genetically modulated differences in rates of biological aging and life span, and is a modifiable outcome. Furthermore, the results highlight the interdependence between body weight, diet, and health and aging, and our observations agree with the well-known influence of body mass on longevity, and the more favorable health profile associated with lower body mass within a species (Bartke, 2012).
The strain with the most decelerated clock, and presumably

| Technical considerations and limitations
Before concluding, we should address a few caveats. To our knowledge, almost all existing DNAm clocks have relied on bisulfite-based assays, which have the advantage of providing single CpG resolution (Hannum et al., 2013;Horvath, 2013;Levine et al., 2018;Petkovich et al., 2017;Stubbs et al., 2017;Thompson et al., 2018;Wang et al., 2017). For the mouse model, the reduced representation data (RRBS) generated by Petkovich et al. (2017) from ~141 B6 mice have been used in different studies to define tens of thousands of age-dependent CpG sites (over 43,000 CpGs in Lowe et al. (2018), and over 146,000 in Sziraki et al. (2018)). Compared to these, the present work identified only 306 age-DMRs at the suggestive threshold that covered 2691 CpG sites. This is still a relatively modest number of age-dependent CpGs, and this is likely due to the small sample size and the genotype heterogeneity of the present cohort.
Another contributing factor may be the methylome assay we used, as MBD-seq provides lower quantitative sensitivity (Gujar, Liang, Wong, & Mozhui, 2018). Nonetheless, MBD-sequencing still delivers highly sensitive and replicable quantification of genome-wide methylation (Aberg et al., 2018). Using alternate assay methods also demonstrates that the DNAm age estimators are robust to the techniques used to quantify DNAm.
We note that the present study is specific to females. We included a few male cases in the methylome assay, and the age-DMRbased clock closely estimated the chronological ages of the three males as well. Longevity is reported to be highly strain-dependent in the BXD panel, with significant correlation in median life span between male and female BXDs, but also with notable sex differences in the underlying genetics (Lang et al., 2010). The decelerated DMRmAge-acc in the one male BXD102 also suggests strain similarity in epigenetic aging rate between a male and the females of a longlived strain. However, with such a limited sampling of male cases, no conclusion can be drawn, and their inclusion may have added more heterogeneity than information. However, excluding the male cases, as was done for the subset analysis in the 55 female samples, did not alter the main findings of the study.

| Statistical analyses
Descriptions of animal protocols, sample processing for MBD-seq, sequence alignment, initial bioinformatics, data filtering, and quality checks are provided in Appendix S1 and Figure S1. Additional description of statistics and transcriptome analysis are also described in Appendix S1.
For the EWAS to detect DMRs, we applied the following model using the lme4 R package (v4_1.1-21) (Bates, Mächler, Bolker, & Walker, 2015): lmer(logRPKM ~age + BW0 + median-LifeSpan + (1|StrainDiet)). This was first applied to the full set of 69 samples. Manhattan plots were generated using the qqman R package (Turner, 2018). For the DMRs, we evaluated relative enrichment in genomic features (i.e., introns, exons, and CpG islands) compared to the background set using the hypergeometric test in R (R codes in Table S3). Following that, we carried out functional annotation and enrichment analysis against the whole background genome using the GREAT application (version 4.0.4) (McLean et al., 2010). Enrichment was calculated using both a binomial test and a hypergeometric test, and we report only categories that were significant by both methods (FDR < 0.05). To search for human GWAS hits, we referred to the GWAS catalog (Buniello et al., 2019), and searched for the terms "body weight", "aging", and "longevity". For variants with reported association p ≤ 5 × 10 −8 , the mapped human gene symbols were then matched to the corresponding mouse genes associated with the DMR.  Thompson et al. 2018(Thompson et al., 2018, the "age acceleration" was computed as the residuals after fitting the predicted age to chronological ages: residuals(lm(DMRmAge ~Age)).

| Estimating epigenetic age and age acceleration
The DMR EWAS was repeated using the same mixed model (lmer(logRPKM ~age + BW0 + medianLifeSpan + (1|StrainDiet))) in a randomly selected subset of 55 female samples. Random subsampling was done using the R function, sample() with n = 55 and replace = FALSE, after excluding the 3 males. From this, age-DMRs (age-associated at p ≤ 1.0 × 10 −4 , 237 CpG regions) and the respective regression coefficients were then used to re-compute the DMRmAge using the same weighted averaging method described above. We then used the 14 samples excluded from the EWAS as a test set to examine the association of DMRmAge-acc with life span, body weight, and diet.

| CON CLUS ION
Our results demonstrate that the epigenetic clock defined from age-DMRs is sensitive to subtle differences in natural life span among female mice that arise from common genetic variants, and is modifiable by environmental interventions, such as diet. The intercorrelations between epigenetic aging, body weight, and longevity also provide evidence that the methylome could provide a mechanistic link between the well-known effect of body mass on aging and life span.

ACK N OWLED G EM ENTS
The present work relied on resources from the BXD Aging Colony, and we thank the entire team, particularly Dr. Lu

CO N FLI C T O F I NTE R E S T
No competing interests.

AUTH O R CO NTR I B UTI O N S
KM supervised the study, contributed to conception of study design and analysis, and wrote the manuscript. RWW is the principal investigator of the BXD Aging Colony and contributed to the conception of the aging study, and provided access to the biorepository resource. JVSS, AHBH, and EGW contributed to the laboratory work, and JVSS performed the primary bioinformatics and initial data processing. DA and SR contributed to analysis and data acquisition. All authors contributed to and approved the final version of the manuscript.

E TH I C A L A PPROVA L
All animal procedures were in accordance to protocol approved by the Institutional Animal Care and Use Committee (IACUC) at the University of Tennessee Health Sciencne Center.

DATA AVA I L A B I L I T Y S TAT E M E N T
The normalized MBD-seq data for the 368,300 CpG bins (including chromosomal coordinates, region annotations, CpG density, and variant counts) are available from the NCBI NIH Gene Expression Omnibus. The full raw fastq files are available from the NCBI NIH Sequence Repository Archive. GEO accession ID is GSE137277, and the linked SRA accession ID is SRP221380.