Genome Methylation Predicts Age and Longevity of Bats

Bats hold considerable potential for understanding exceptional longevity because some species can live eight times longer than other mammals of similar size [1]. Estimating their age or longevity is difficult because they show few signs of aging. DNA methylation (DNAm) provides a potential solution given its utility for estimating age [2-4] and lifespan [5-7] in humans. Here, we profile DNAm from wing biopsies of nearly 700 individuals representing 26 bat species and demonstrate that DNAm can predict chronological age accurately. Furthermore, the rate DNAm changes at age-informative sites is negatively related to longevity. To identify longevity-informative sites, we compared DNAm rates between three long-lived and two short-lived species. Hypermethylated age and longevity sites are enriched for histone and chromatin features associated with transcriptional regulation and preferentially located in the promoter regions of helix-turn-helix transcription factors (TFs). Predicted TF binding site motifs and enrichment analyses indicate that age-related methylation change is influenced by developmental processes, while longevity-related DNAm change is associated with innate immunity or tumorigenesis genes, suggesting that bat longevity results, in part, from augmented immune response and cancer suppression.

effects. Therefore, to determine if age in bats can be predicted by DNAm we used a custom microarray that assays 37,500 conserved cytosine-phosphate-guanine (i.e. CpG) sites to measure DNAm extracted from wing biopsies of 694 known-aged individuals representing 26 species of bats (Supplementary Table 1). Similar to human epigenetic clocks [4,7], elastic-net regression accurately predicted chronological age from a linear combination of DNAm beta values (henceforth DNAmAge) using 150 CpG sites. Leave-one-out (LOO) cross-validation shows that DNAm predicts age accurately (Fig.1a). Even more accurate age estimation is possible for species (or genera) with sufficient data. For example, the correlation between To assess how well DNAm might predict age in a species not represented by our samples, we conducted a second cross-validation analysis in which data for one species was left out (leave-one-species out; LOSO) and ages were predicted for that species using a clock estimated from the remaining data. This analysis (Fig. 1b) resulted in a correlation between observed and predicted age of 0.8 (MAE=1.3 years). The LOSO analysis also showed that species vary in DNAmAge estimates. For example, Desmodus rotundus (sp. 5, Fig. 1b) exhibit lower values of DNAmAge (suggesting lower aging rates) than Phyllostomus hastatus (sp. 15, Fig. 1b), consistent with the longer lifespan of D. rotundus [1]. Experimental [11,12] and observational data [13] suggest that DNAm rate varies with lifespan in model organisms.
However, interspecific studies have so far used different methods on a few primate, rodent, or canid species [12,13] making it difficult to determine reasons for methylation differences.
Thus, to determine if DNAm rate predicts bat longevity, we incorporated a recent bat phylogeny [14] into a generalized least squares regression (PGLS) to predict the longevity quotient -the ratio of observed to expected maximum lifespan for a mammal of the same body size [1]. We identified a common set of CpG sites for this analysis by first conducting a meta-analysis of all age-DNAm correlations by probe for 19 bat species with 15 or more samples (Methods). We find that mean DNAm rate at 2000 age-associated sites (henceforth, age differentially methylated positions or age DMPs) predicts longevity. Long-lived species have lower rates of change at hypermethylated (Fig. 2a) and hypomethylated sites (Fig. 2b).
Assuming that the rate of change in DNAm reflects epigenetic stability, these results suggest that better epigenetic maintenance is associated with longer maximum lifespan, independent of body size, across bats.
The distribution and function of genomic regions that exhibit age or longevity-related changes in DNAm are not fully understood [7,15]. To identify DMPs associated with longevity, we compared relationships between DNAm and age for three long-lived species and two short-lived species (Fig. 2) from four bat families (Methods). After adjusting for multiple testing (BY FDR = 0.05, [16]), we identified 1491 longevity DMPs, including 694 in which short-lived species gain DNAm faster with age and 797 in which short-lived species lose DNAm faster. In the following, we refer to these sites as longevity DMPs. Both longevity and age DMPs are widely distributed in the genome, but differ in relative abundance across chromosomes (Fig. 3a,b). For example, of the 1228 probes that map to chromosome 1 (syntenic with the human X chromosome) in R. ferrumequinum (the bat genome with the most 5 mapped probes, 34,411, Supplementary Table 3) only 13 are age-associated while 61 are longevity-associated. Not surprisingly, 672 of 838 sites (80.2%) that differ between the sexes in methylation values are on the R. ferrumequinum X. Sex DMPs are independent of age DMPs (6.1% overlap, P = 0.25, Fisher's Exact Test, FET) but not longevity DMPs (5.8% overlap, P = 0.0064, FET). Age DMPs are also not independent of longevity DMPs (22.5% overlap, P < 0.0001, FET, Fig. 3c). Change in DNAm with respect to age correlates with change in DNAm with respect to longevity (r = 0.462, P < 0.0001); when limited to promoter regions, almost all DMPs exhibit hypermethylation and mapped to the same gene in multiple bat genomes (Fig. 3e).
The gene nearest each CpG site often differs between bats and humans (Fig. 4a). In contrast, genomic regions occupied by DMPs are similar among bat species ( Supplementary   Fig. 3). For example, 70% of 3197 probes that map to a promoter region in M. molossus map to a promoter region in R. ferrumequinum (Fig. 4b). Promoter regions in M. molossus (and other bats, Supplementary Fig. 4) are enriched for hypermethylating, but not hypomethylating, age and longevity DMPs (Fig. 4c,d). In bat genomes where CpG islands have been identified, hypermethylating DMPs are much more likely than hypomethylating DMPs to be located in CpG islands for the age (P < 0.0001, FET) and longevity (P < 0.0001, FET) associated sites.
Given that regions near promoters contain more age and longevity DMPs than expected, we evaluated the genes nearest those DMPs for possible functions. In M. molossus, 36% of genes had DMPs associated with age and longevity (Fig. 3d). Other bats show similar patterns ( Supplementary Fig. 4d). Not surprisingly, M. molossus genes with age or longevity DMPs in promoter regions show similar patterns of enrichment among biological processes, 6 i.e. regulation of development, biosynthetic processes, and regulation of transcription are enriched (Fig. 4e). Genes with age DMPs in promoter regions are further enriched for cell differentiation and cellular development. With regard to protein class, gene lists for both age and longevity DMPs are enriched for homeodomain transcription factors containing helixturn-helix motifs (Fig. 4f). These patterns are characteristic of other bat species ( Supplementary Fig. 5), although the gene list composition varies. For example, the same gene was near a hypermethylated site in a promoter region for 92 out of 143 age-associated genes in at least three of the four bat genomes used for identifying longevity DMPs.
Additionally, hypermethylated age genes in bats strongly overlap hypermethylated age genes reported for dogs [18] (e.g. 60 of 316 genes for M. molossus, P = 7.79e-36, FET). In contrast, among hypomethylated age genes, only four of 353 M. molossus genes overlap with dog genes (P = 0.202, FET). Bat age genes include more immunity genes (Fig. 4g) than expected (P = 0.021, FET), but overlap between immunity and longevity genes is far greater (P = 0.001, FET). Bat longevity genes also overlap genes frequently mutated in human tumors Gene promoters poised for transcription can be identified by histone marks and chromatin states [19]. Given that hypermethylation associated with human aging occurs at bivalent chromatin domains [20], we used eFORGEv.2.0 [21] to predict how DMPs likely influence regulatory regions based on epigenetic studies in human or mouse cell lines derived from relevant tissues. We find that sites exhibiting hypermethylation with age and longevity in bat wing tissue are enriched for repressive histone H3 trimethylated at lysine27 7 (H3K27me3) and active H3K4me1 marks in skin, muscle and blood cell lines (Fig. 5a, b).
Hypomethylated age DMPs are enriched in all three tissues for H3K9me3, while hypomethylated longevity DMPs show no enrichment (Fig. 5a, b). Analysis of predicted chromatin states reveals that hypermethylated age DMPs are enriched in all three tissues for repressed polycomb complexes, while hypomethylated age DMPs are enriched for quiescent chromatin states (Fig. 5c). Longevity DMPs, both hypermethylating and hypomethylating, also show enrichment for quiescent states, as well as enrichment for repressive polycomb complexes or enhanced bivalent states in some tissues (Fig. 5d).
Transcription factor (TF) motifs involved in cell cycle regulation and genome stability are enriched among hypermethylating age sites (Fig. 5e) Longevity TF motifs are largely independent of age TF motifs (Fig. 5e), with one exception, c203-transcription factor AP-2 gamma (TFAP2C), which is involved in epidermal cell lineage commitment [27] and regulation of tumor progression [28]. The other longevity TF motifs also have known associations with tumorigenesis. GCM1/3 binds to pleiomorphic adenoma gene-like 1 (Plagl1), which codes for a protein that suppresses cell growth. This gene is often methylated and silenced in cancer cells [29,30]. CNOT3 acts as a tumor suppressor in T-cell acute lymphoblastic leukemia (T-ALL) [31] but can also facilitate development of non-small cell lung cancer [32]. Finally, HIC1, hypermethylated in cancer 1 8 protein, acts as a tumor suppressor and is involved in regulation of p53 DNA damage responses [33,34]. Only a single TF motif, HD/5 in the BARHL2 group [35], was associated with hypomethylated longevity DMPs.
Enrichment analyses [36], using the age and longevity DMPs found in promoter regions to create gene lists for M. molossus, identify several key regulators that are significantly associated with hypermethylated sites, but none with hypomethylated sites ( Fig.   5f and Supplementary Fig. 5c) Orthodenticle homeobox 2 (OTX2) is associated with both age and longevity, whereas other predicted TFs largely differ between age and longevity. The highest-ranked age TF, RE1 silencing transcription factor (REST), is induced during human aging and represses neuronal genes that promote cell death [37]. Note that four of nine transcription regulators predicted to be associated with longevity frequently undergo mutations in human tumors and three are involved in innate immunity (Fig. 5f).
DNAm influences many processes including development [38], gene regulation [39], genomic imprinting [40], X chromosome inactivation [41], transposable element defense [42], and cancer [43]. Over 75% of CpG sites are typically methylated in mammalian cells, but global DNAm declines with age, which can lead to loss of transcriptional control and either cause or contribute to deleterious aging effects [44]. As with other species [18, [45][46][47], agerelated changes in DNAm occur throughout bat genomes. While only 150 CpG sites are sufficient to predict chronological age, these represent only a fraction of sites that correlate with age, because penalized regression excludes highly correlated variables to avoid multicollinearity. Consequently, we carried out a marginal analysis that correlated individual cytosines with age across species to identify age DMPs. At these sites, long-lived species exhibit a lower rate of change in DNAm, while short-lived species exhibit faster increases in 9 DNAm. How those changes contribute to longevity is not entirely clear, but our results suggest several key transcriptional regulators are involved either by binding to DNA at open chromatin or by being enriched among genes with DMPs in promoter regions and modulate the rate at which DNAm changes between short and long-lived species.
Our results are consistent with an epigenetic clock theory of aging that connects beneficial developmental and cell maintenance processes to detrimental processes causing tissue dysfunction [7]. A large body of evidence links age-related hypermethylated sites to genes and genomic regions that influence developmental processes [18,20,48]. The sites that gain DNAm with age also tend to be in CpG islands, consistent with studies in humans [49].
But, in contrast, we find little enrichment for genes associated with hypomethylated sites, and these genes are less likely to be shared across species. We interpret these results to indicate that DNAm loss with age is widespread and not concentrated in particular pathways. These DNAm patterns could differ by tissue, as has been frequently observed [50]. However, bat wing tissue is capable of unusually rapid regeneration [51,52] and consists of multiple tissue types [53], making it particularly useful for measuring age-related changes in DNAm.
DNAm of genes suppressed in stem cells is a hallmark of cancer [48]. Several lines of evidence suggest that bat genes with longevity DMPs are important for cancer suppression and provide enhanced immunity. First, these genes disproportionately include many known to mutate frequently in human cancers or involved in innate immunity. Second, several transcription factors identified by motif analysis act as tumor suppressors, such that if they are silenced by methylation in older individuals, tumor formation should be more likely. Third, among the transcription factors identified from the list of genes with hypermethylated sites in promoter regions, several of them mutate in human cancers. While bats are not immune from cancer [54][55][56][57][58][59], genetic adaptations for tumor suppression have been described for Myotis brandtii [60] and M. myotis [61] to help explain the extreme longevity of those species. Bats also have genetic mechanisms that enable strong antiviral immune responses without inducing damaging inflammatory reactions (e.g. cytokine storm [62]) that may enable them to tolerate high levels of viral exposure [9,10,17,63,64]. The results of this study are consistent with the hypothesis that enhanced epigenetic stability, especially associated with innate immunity and cancer suppression genes, facilitates the exceptional longevity in bats.

Wing tissue samples
The number of tissue samples per species and per sex, as well as the range of ages of individuals of each bat species, are provided in Supplementary Table 1. Wing punches were taken from individually marked animals that were either kept in captivity (15 species) or recaptured as part of long-term field studies (11 species). For 652 samples the individual was marked shortly after birth, so age estimates were exact. For the remainder, age represented a minimum estimate because the individual was not initially banded as a juvenile. We used minimum age estimates when other evidence, such as tooth wear or time since initial capture, indicated that the minimum age estimate was likely to be close to the real age. Below, we provide additional information on when and where samples were taken from either captive or free-ranging animals.
Pallid bats, Antrozous pallidus, were captured between 2005 and 2008 at six sites in central Oregon (44.94°N, 120.38°W) using mist nets over a water source or outside a night roost or with a handnet on an extension pole outside a day-roosting crevice. Each bat was weighed, measured and marked with a numbered band. Adults were distinguished from juveniles by closed epiphyseal gaps [8]. Tissue samples were obtained from wing membranes using 3 mm biopsy punches and stored in 95% ethanol until DNA was extracted using a Qiagen DNeasy Tissue Kit. DNA extracts were stored frozen at -80°C. Live animal procedures conformed to the American Society of Mammalogists guidelines [65]   Samples of velvet free-tailed bats, Molossus molossus, come from a long-term study [71][72][73] in Gamboa, Panama (09°07' N 79°41' W), where the bats roost in crevices in houses.
We captured social groups with mist nets (Ecotone, Gydnia, Poland) at the entrance of roosts during evening emergence and individually marked all bats with a subcutaneous passive integrated transponder (Trovan ID-100, Euro ID, Weilerswist, Germany) at first capture.
Wing tissue samples were taken with a 3 mm biopsy punch and stored in 96% ethanol until Common noctules, Nyctalus noctula, were captured as part of a long-term study [79][80][81]  Wing tissue samples were taken from greater horseshoe bats, Rhinolophus ferrumequinum, by using 3 mm biopsy punches between 2016 and 2018 from wild female bats as part of a long-term study at a maternity colony in Gloucestershire, UK [86][87][88]. Bats were captured at the roost with hand nets, and all individuals were weighed, ringed with stainless-steel rings and morphometric data such as forearm length recorded. All bats studied were first marked as infants, so we could be certain of their age. Bats were aged between 1-21 years, with the 40 individuals selected in a fairly even manner across this age span. All  [86][87][88]. Bats were mist-netted in the vicinity of their roosts. Wing tissue was sampled with a 4 mm biopsy punch, individuals were marked with colored plastic bands, sexed, measured and age class determined (juvenile: 0-4 months, subadult: 5-10 months, or adult>10 months) [89]. Age was determined exactly for individuals that were banded as pups and recaptured as adults. Ethanol (80%) was used to preserve tissue samples, and a salt-chloroform procedure or Qiagen BioSprint 96 DNA Blood Kit was used for DNA isolation [89,90].  [91,92]. Bats were captured with mist nets when entering or leaving their day roosts, individually banded with two coloured plastic bands on their forearms, and a wing tissue biopsy sample (4mm) preserved in 80% ethanol was taken.
Age was determined exactly for individuals that were banded as pups and recaptured as adults. DNA was extracted with a salt-chloroform procedure or with the Qiagen BioSprint 96 DNA Blood Kit [91,93].
Mexican free-tailed bats, Tadarida brasiliensis, are housed at Bat World Sanctuary, a licensed non-profit bat rehabilitation facility and accredited sanctuary in Weatherford, Texas.
Most individuals sampled were rescued as pups, although some were rescued as adults, making their exact age unknown. Individuals are group housed in large indoor enclosures.
Wing membrane biopsies (4 mm) were collected in August 2019 and stored in Zymo DNA Shield until DNA was extracted using a Zymo miniprep plus kit.
After extraction DNA concentration was estimated with a QuBit and samples were concentrated, if necessary, to reach a minimum of 10 ng/µl in 20 µl. To estimate rates of methylation we limit analyses to the 23 species for which we had more than 10 samples from known-aged individuals. Maximum lifespan for each species was obtained from [1] or from captivity records (Supplementary Table 1).

DNA methylation profiling
DNA methylation was quantified using a custom Infinium methylation array, "HorvathMammalMethylChip40" that was designed using an alignment of 62 mammalian

Probe mapping and annotation
We used sequences and annotations for ten bat genomes (Supplementary Table 3), which include six recently published reference assemblies [19], to locate each 50 bp probe on the array. The alignment was done using the QUASR package [95] with the assumption for bisulfite conversion treatment of the genomic DNA. For each species' genome sequence, QUASR creates an in-silico-bisulfite-treated version of the genome. The set of nucleotide sequences of the designed probes, which includes degenerate base positions due to the bisulfite conversion, was expanded into a larger set of nucleotide sequences representing 20 every possible combination of degenerate bases. We then ran QUASR (a wrapper for Bowtie2) with parameters -k 2 --strata --best -v 3 and bisulfite = "undir" to align the enlarged set of probe sequences to each prepared genome. From these files, we collected only alignments where the entire length of the probe perfectly matched to the genome sequence (i.e. the CIGAR string 50M and flag XM=0).
Following the alignment, the CpGs were annotated based on the distance to the closest transcriptional start site using the Chipseeker package [96]. A gff file with these was created using these positions, sorted by scaffold and position, and compared to the location of each probe in BAM format. We report probes whose variants only mapped to one unique locus in a

Creation of epigenetic clocks using penalized regression
We developed epigenetic clocks for bat wing tissue by regressing chronological age on all CpGs that map to at least one of the ten bat genomes. To improve linear fit we transformed age to sqrt(age+1). Penalized regression models were created in the R package "glmnet" [97].
We investigated models produced by "elastic net" regression (alpha=0.5). The optimal penalty parameters in all cases were determined automatically by using a 10-fold internal cross-validation (cv.glmnet) on the training set. By definition, the alpha value for the elastic net regression was set to 0.5 (midpoint between Ridge and Lasso-type regression) and was not optimized for model performance. We performed two cross-validation schemes for arriving at unbiased estimates of the accuracy of the different DNAme based age estimators. One type consisted of leaving out a single sample (LOO) from the regression, predicting an age for that sample by regressing an elastic net on the methylation profiles of all other samples, and iterating over all samples. We conducted LOO analyses using all samples from all species, using all samples from each species, and using all samples from several species in the same genus (see Supplement for results of the latter two analyses). The second type consisted of leaving out a single species (LOSO) from the regression, thereby predicting the age of each sample using the data for all other species.

Identification of differentially methylated positions for age and longevity
To identify differentially methylated positions (DMPs) associated with age, we computed the Pearson correlation coefficient between methylation level (β) and transformed chronological age (sqrt(age +1)) for each of the 37,500 sites for the 19 species with 15 or more samples Table 1). The significance of each site across species was then evaluated using Stouffer's unweighted z-test [98]. CpG sites were ranked by significance and the top 2000 sites were selected for subsequent analyses and are referred to as age DMPs. For probes with contrasting patterns in different species, direction was assigned based on the most frequent direction across species to ensure mean methylation rates are comprised of the same set of sites in each species. Because we used all sites on the array, some sites do not map to a unique position in all available bat genomes. The supplement indicates how many sites map to each of the four species with a genome that were used for identifying longevity DMPs.

(Supplementary
To identify DMPs associated with longevity we compared methylation rates between three long-lived species (R. ferrumequinum, D. rotundus, and M. myotis) and two short-lived species (M. molossus and L. yerbabuenae). We chose these five species because they represent three independent lineages of increased longevity [1] and because high-quality genome assemblies are available for four of them [17]. We used a linear mixed-effects model (nlme) to model methylation level (β) in response to transformed chronological age (sqrt(age +1)), longevity category, and their interaction, with species included as a random effect. We defined probes as longevity-associated if the p-value of the interaction term was less than 0.05 after Benjamini-Yekutieli (BY) false discovery rate correction [16]. In this analysis, a positive interaction means a steeper positive slope for the short-lived species relative to the long-lived species. If the main effect of age is positive (hypermethylation) and the interaction is positive, then short-lived species are gaining methylation faster. If the main effect is negative and the interaction is negative, then short-lived species are losing methylation faster.

Phylogenetic analysis of bat longevity
Using phylogenetic generalized least squares regression (PGLS) we tested the effect of mean rate of methylation change on longevity using both the longevity quotient (LQ) and maximum longevity (log-transformed). LQ is the ratio of the observed species maximum lifespan to the maximum lifespan predicted for a nonflying placental mammal of the same body mass [1].
We present results for LQ in the text and for a model containing both log(maximum longevity) in Supplementary Table 2. For each species with at least 10 known-age samples, we calculated the mean rates of hypermethylation and hypomethylation using the top 2000 age-associated DMPs as described above. Hypermethylation and hypomethylation rates were tested separately. Log-transformed body mass [from 1] was included as a covariate in each log(maximum longevity) model. Phylogenetic relationships among bats are based on a recent maximum likelihood tree [14]. Models were fit via maximum likelihood using the "gls" function of the "nlme" package [99] and assume a Brownian model of trait evolution.

Probe and gene enrichment analyses
To determine how changes in methylation influence age and longevity, we conducted enrichment analyses on the CpG probes and on the genes nearest to them. We used eFORGE 2.0 [23, 98] to test for enrichment among age or longevity DMPs that either increase or decrease in methylation in comparison with five histone marks and 15 chromatin states mapped in cell lines by the Epigenomics Roadmap Consortium (http://www.ncbi.nlm.nih.gov/epigenomics). Bat wing tissue is unusual in that it contains epithelial skin, muscle, blood and elastin [53]. Consequently, we limited enrichment analyses to data from cell lines derived either from skin, blood or muscle. We also restricted the 24 analysis to probes mapped in a bat genome at least 1 kb apart. We used Demodus rotundus to provide a background probe set but obtained very similar results by using other bat genomes, e.g. Eptesicus fuscus or Pteropus vampyrus, available in eFORGE as backgrounds for the 37K array. We present enrichment values for each DMP set as the -log10 p binomial value and consider those outside the 95th percentile of the binomial distribution after correction for multiple testing [16] as significant.
We identified putative transcription factors that could utilize open chromatin and bind to the DNA by testing for enrichment in each DMP set for predicted binding sites among the probes on the 37K array. Binding sites were included if their FIMO (Find Individual Motif Occurrence) [100] p-value was less then 10e-5. We then used a hypergeometric test (phyper) to evaluate overlap between probe sets and transcription factor motifs obtained from four transcription factor databases: TRANSFAC [101], UniPROBE [102], HT-Selex [103], and JASPAR [104]. Redundant transcription factor motifs were then consolidated into clusters [cf.
105] to identify distinct transcription factors. Function was inferred using information derived primarily from studies in mouse and humans [106].
We used several approaches for determining the type and function of genes associated with age and longevity DMPs. First, we identified the gene (using human orthologs) with the nearest transcription start site to every mapped probe for each of the four species, R.  Table 3). Given that the probes were designed to align to regions conserved across all mammals, we suspect some of the differences among species in gene associations reflect variation in genome assembly or annotation. In addition, an important caveat to keep in mind is that the CpGs on the array do not randomly sample the genome. Thus, even when we use mapped probes or the genes near them as background for enrichment tests, there is potential for bias given that the probes are in highly conserved regions. We assumed genes were associated with hypermethylated DMPs if they had more hypermethylated than hypomethylated sites nearest their transcription start site. We present results in the text for DMP-gene associations for M. molossus because it was the only short-lived species with a genome, but we summarize the DMP-gene associations for the other three species in Supplementary figures. Because we anticipated the mechanisms responsible for causing increases in methylation over time likely differs from those causing decreases, we conducted separate enrichment tests for genes with hypermethylated and hypomethylated sites associated with age and longevity using Panther v.14 [107,108] in relation to biological process, molecular function, and protein class. We carried out enrichment tests using genes with DMPs in promoter regions because promoter regions showed enrichment for hypermethylated sites.
We also used the significant age and longevity gene promoter lists to predict possible transcription factor regulators using BART, Binding Analysis for Regulation of Transcription [36], which correlates the cis-regulatory profile derived from a gene set to the genomic binding profiles of 918 transcription regulators using over 7000 human ChIP-seq datasets. We report the Irwin-Hall P value, which indicates significance of the rank integrated over three test statistics [36].
In addition, we carried out three additional analyses to assess gene function. The first utilized a list of 394 genes associated with changes in methylation over the lifespan of dogs