Co-analysis of methylation platforms for signatures of biological aging in the domestic dog reveals previously unexplored confounding factors

A framework to detect biological age in the dog

We define biological age to be the quotient between the age of a sample and its inferred lifespan based on breed averages (Figure 1). To account for datasets containing mixed-breed samples, where no average breed lifespan could be inferred, we additionally defined biological age to be the product between the sample age and its individual weight, as lifespan and weight are generally anticorrelated in the dog [15] (Figure 2). These definitions can be statistically formalized as an interaction term between chronological age and lifespan or weight (hereby also referred to as moderators or moderator variables), and have some key differences with the model containing only the moderator variables as a main effect [8, 28] (Figure 1). The interaction term implies a changing rate of methylation as it correlates with age and longevity or weight (Figure 1). Conversely, the presence of significant moderators likely indicates that the methylation values at a site are associated with lifespan or weight, but independent of age. Furthermore, the known non-linear relationships between methylation and age may also mask interaction effects and produce inaccurate summary statistics (Figure 1). At a single site level, nonlinearities are overcome by using power transformations [29–31] or nonparametric approaches such as BayesAge [32] (Methods), while at a whole methylome level the epigenetic pacemaker theory developed by Snir et al. [33] ensures that a non-linearity trend is estimated globally.

Dog breed standard lifespans and weights display high phylogenetic signals

Our analysis makes use of three previously published dog methylation datasets, each with differing breed compositions and sample size, featuring three distinct methylation platforms and two tissue types. The three studies have a well-balanced composition of males and females and a similar distribution of ages and weights (inter quantile ranges 2.49 ± 0.50 – 8.41 ± 1.26 years (yrs) and 10.27 ± 1.14 - 29.47 ± 0.10 kg, respectively, and Figure 2). Importantly, the capture sequencing dataset was composed of both purebred and mixed-breed dogs, and therefore lifespan values could not be assigned to all dogs. However, empirical weight measurements were available for all mixed-breed dogs in the dataset. Only 11 mixed-breed dogs of the185 had not reached maturity at the time indicated by the study metadata. The weights of all purebred dogs, including those that had not reached adulthood, were replaced with their respective breed standard adult weights [3]. We replicate the previously observed anticorrelation [15] between breed weights and lifespans for the breeds present in all datasets (ρ=-0.59 and F-test p-value=4.29e-77), and note that the relationship holds for each dataset individually (Figure 2). The use of standard breed weights and lifespans can obscure environmental effects such as overfeeding or under-exercising, which affect the weight and longevity of individual dogs. Both lifespan and weight displayed highly significant phylogenetic signals compared to the canonical dog breed phylogeny (Pagel’s [34] λ = 0.65 and 1 respectively). This highlights the difficulty encountered in decoupling not only lifespan from weight but also the effects of population stratification from both traits.

We also sought to determine how many sites were shared across the different datasets. However, the sparsity of array sites and capture experiments in addition to the stochastic dropout of lower depth sequencing could lead to an underestimation of shared sites. To make the comparison between arrays and sequencing data more less strict, we extended the genomic range of every site in every dataset 50 bps upstream and downstream and merged any overlapping loci. We postulate that neighboring methylation sites generally belong to larger scale elements, such as CpG islands or SINEs, and therefore have dependent methylation values. Although even a single bp overlap between datasets was considered a match, we observed very little overlap between platforms, suggesting that fully concordant cross-platform results are not necessarily expected (Table 1 and Figure 3).

All three studies were able to independently derive chronological age clocks and determine sex from methylation data. However, consistent with the absence of shared markers among datasets, there were notable differences in major axes of variation. The mammalian methylation array was still able to predict sex even after removing any sites located in sex chromosomes but performed poorly at breed determination based on a methylation distance matrix (Mantel test p-value 0.155) (Figure 3). Conversely, capture sequencing offered a much more accurate recapitulation of phylogeny (Mantel test p-value <1e-5), which could be validated using SNPs extracted from the same dataset. Average RRBS methylation distances did not recapitulate phylogeny (Mantel test p-value 0.195), but genetic variation could still be extracted from the dataset using non-methylated bases (Methods).

Sparse, biological age epigenetic clocks are confounded by phylogeny

We next created biological age clocks moderated by lifespan and weight and corrected for phylogenetic relationships between samples. In order to enhance comparability across datasets we limited the number of markers that make up any clock to ≤ 25, as this value lies in the interval between the number of markers corresponding to the minimum cross validating penalty value and the penalty value one standard error above it (Figure 4 and Supplementary Figure 1, and Methods). We posit that marker sparsity in biological age clocks can prevent overfitting and promote the inclusion of sites with large biological age effects (Supplementary Figure 1). Additionally, the ability to create sparse, cross-validating chronological age clocks in all datasets (Figure 4) motivated the assessment of such clocks for biological age. Of note, all biological and chronological age epigenetic clocks passed ten-fold cross-validation both under ordinary and generalized least squares error models except for RRBS moderated by weight (Figure 4C, 4D and Supplementary Figure 1). Finally, we applied the epigenetic pacemaker and BayesAge clocks [32] to the methylation data and regressed the output methylation state curves against biological age.

Using phylogenetically corrected error models while enforcing sparsity offered key insights into the interpretation of biological age across the three datasets (Figure 4). A qualitative comparison of ordinary least squares sparse models to their phylogenetically corrected analogs suggests a stark suppression regarding significance of the interaction terms in both the mammalian methylation array and capture sequencing datasets (Figure 4B, 4D, 4F, second column). In the case of RRBS, phylogenetic errors hampered cross validation (Supplementary Figure 1). This demonstrates that we cannot construct reliable, sparse biological age clocks that are uncounfounded by phylogeny using weight or lifespan as a moderator. Non-sparse biological age clocks generated using minimum cross-validating penalty values (Methods and Supplementary Figures 1, 2) included a variable number of markers and therefore were much less comparable. Notably, the mammalian methylation array required up to 264 markers to achieve its optimal biological age prediction power, and visual examination of those markers revealed poor correlation between individual sites and biological age (Supplementary Figure 3). In contrast, the capture sequencing epigenetic clock required 45 markers and retained a significant yet attenuated biological age signal even under the phylogenetically corrected error model.

Biological age clocks constructed using breed standard weight outperformed those constructed using lifespan across all datasets, even under sparse conditions and after correcting for phylogeny (Figure 4). Similarly, ordinary least squares clocks qualitatively outperformed those constructed using phylogenetic least squares, suggesting that the relationship between breed standard weights and lifespans and aging rates is often concordant with phylogenetic tree relationships. The epigenetic pacemaker and BayesAge models revealed small but significant biological age effects in almost all datasets. However, these effects could not be decoupled from phylogeny and a decorrelating transformation of the methylation matrices akin to that used in our generalized, penalized regression model yielded no valid results. As such, the epigenetic pacemaker and BayesAge predictions were deemed more comparable to penalized, ordinary least squares regression models than those that were phylogenetically corrected.

Biological age methylation markers are scarce and display low effect sizes

A site-by-site examination of biological age signals revealed few biological age candidate sites both before and after correcting for phylogeny (Figure 5). Because sites more tightly associated with chronological age produce smaller interaction effects, the p-values for age and the interaction term between age and the moderator variable are not independent. Nevertheless, we report a general depletion of sites that are simultaneously significant for age and the interaction term between age and the moderator, indicating that very few sites display true potential biological age effects. Even then, the sites that make up the biological age epigenetic clocks tend to be the best candidates for biological age effects (Figure 5), highlighting that the phylogenetic and ordinary penalized regression models are selecting for the correct features of their component sites. Only the RRBS biological age epigenetic clock displayed an excess of sites independently associated with weight (Figure 5B, second panel). Indeed, RRBS presents the largest inflation of weight-associated p-values of all datasets, especially given its small sample size (Figure 5), highlighting that it is potentially composed of more phylogenetically concordant sites than the two other datasets.

Gene ontology enrichment in biological age loci

We performed gene ontology (GO) enrichment analysis on the set of genes 5kb upstream or downstream of the 50 top biological age candidate loci. As an insufficient number of significant sites were found to be associated with biological age, the top candidate loci were arbitrarily picked based on a linear combination of p- values of age and the interaction between age and the moderator, relaxing significance to include 50 sites (Figure 5). As expected, due to the small number of shared sites between datasets, the top ten loci with the smallest biological age p-values in each dataset did not colocalize (Supplementary Figure 4). No GO terms were found to be significant in the capture sequencing dataset. Both the RRBS and the mammalian methylation array datasets were significantly enriched for GO terms related to skeletal and nervous system morphogenesis such as “GO:0007389 pattern specification process”, “GO:0048598 embryonic organ morphogenesis” and “GO:0021675 nerve development”. Enrichment for the same or similar terms was reported in the respective studies for chronological age. Of note, the main loci driving these terms in both the RRBS and mammalian methylation array datasets are the HOX gene clusters (in particular HOXA1-10), which have previously been identified as major markers of biological age in mice [36]. Albeit consisting of different positions, the HOXA1-10 region is also represented in the capture sequencing experiment, but it was not found to be significant for biological age effects.

Data reprocessing and metadata regularization

We leveraged three previously published methylation studies in dogs that used different experimental designs: RRBS in dog blood (SRA: SRP065666), an application of the mammalian methylation array in dog blood (GEO: GSE223748 and https://mydata.clockfoundation.org/app/mammalian-consortium-data-browser), and capture sequencing in dog saliva (https://www.tandfonline.com/action/downloadSupplement?doi=10.1080%2F15592294.2022.2069385&file=kepi_a_2069385_sm0108.zip). SNP data for the capture sequencing experiment were requested from the authors. The three studies offered pre-processed and filtered beta values readily available for analysis. We chose, however, to re-process the RRBS dataset to recover more methylation sites and improve comparability. After adapter trimming [42], the raw reads were aligned against the CanFam3.1 Tasha reference [43] and beta values were calculated using the Bismark [44] suite, excluding the first and last three nucleotides, as they are part of the restriction site and reflect biased nucleotide compositions.

bismark -N 1 –output_dir ${output_dir} ${input_fastq} && bismark_methylation_extractor –ignore 3 –ignore_3prime 3 –bedGraph –output ${output} ${input_bam}.

Methylation sites missing in at most four samples were tabulated and resulted in a total of 244,334 sites. We found over 98% Pearson correlation between our methylation calls and the original study. SNPs were called from aligned reads using hard heterozygous allele balance (AB) thresholds of 0.1 0.1 to a total of 158,989. All C →T and A →G substitutions were discarded.

In addition to the metadata provided by each study, we complemented the information tables with breed-based expected lifespan using, if possible, the information published in Horvath et al. [3] and complementing the missing breeds with AKC [45] records.

Construction of phylogenetic variance-covariance and distance matrices

In the analyses where SNPs could be extracted from the corresponding platform (RRBS and capture sequencing), we used Gower distances between samples to create global phylogenies which were then converted into variance-covariance matrices using the vcv function in the R ape package. In the case of the mammalian methylation array, where no SNP panel was available, we inferred population stratification based on breed. First, the average distance between all samples belonging to all possible breed pairs within a reference breed topology was calculated. Then, the resulting tree was subsampled to contain the breeds represented in the mammalian methylation array and the distance value of each breed broadcasted to all samples belonging to that breed (Supplementary Figure 6). Three out of ninety-four breeds were not represented in the reference topology, so a proxy breed from the same putative clade was assigned in their stead (English Setter → Irish Setter, Jack Russel terrier → Glen of Imaal terrier, Weimaraner → Wirehaired Pointing Griffon). Pagel’s lambda was calculated on this canonical phylogeny using the phytools [46] R package. Mantel tests between methylation and genetic variance-covariance matrices were calculated with the mantel.test function provided in the R ape package using 100,000 permutations. To correct for degeneracy in the resulting matrix, we raised each symmetrical entry of the matrix to a random power very close to one of the form F~N(1, 0.01) so that $${D^{\prime }}_{ij}={D^{\prime }}_{ji}=D_{ij}^F$$ where N is a normal distribution draw, D is the original matrix and D′ is the resulting jittered matrix. Distance matrices were converted to variance-covariance matrices using the vcv function in the R ape package. Note that this approach allows for an arbitrary shrinkage of the determinant of the jittered matrix, which will ultimately influence the magnitude of the regression likelihood. Finally, to ensure positive semi-definiteness, we applied the nearPD R base [47] function to the resulting matrix (Supplementary Figure 6).

Phylogenetic penalized regression and epigenome-wide association study

Before applying penalized regression, we performed two corrections to the methylation matrix based on the two previous methods sections: first, we linearized the methylation matrix, site by site, using the Box-Cox transformation implemented in the R MASS [48] package. Second, we performed a decorrelation transformation by multiplying the linearized methylation matrix and response variable with the inverse Cholesky-factored variance-covariance matrix to account for phylogenetic dependence between samples. To validate this procedure, we proved that it is possible to use the glmnet proximal gradient descent framework with decorrelated, generalized errors (Supplementary Text 1 and Supplementary Figures 9, 10). We limited the number of markers that make up any clock to ≤ 25, as this value lies in the interval between the number of markers corresponding to the minimum cross validating penalty value and the penalty value one standard error above it (Supplementary Figure 1). All biological and chronological age epigenetic clocks, except for RRBS moderated by weight, passed ten-fold cross-validation both under ordinary and generalized least squares error models (Figure 4C, 4D and Supplementary Figure 1). Generalized-penalized regression models were then deployed using the R glmnet [49] package with fifteen-fold cross-validation, using age as the dependent variable. Given the diverse breed composition of the three datasets, ten-fold cross-validation separated all the components of more than one breed into a single cross-validation fold, ensuring non-redundancy between training and test dataset partitions. Our results were robust to permutation of the folds.

We applied independent, site-wise generalized least squares (GLS) regression models using methylation as the response variable and sex as a covariate. We implemented an exact solution to GLS regression to bypass the calculation of the inverse variance-covariance matrix for every site, which is standard in R packages such as gls. The interaction term was explicitly modeled to be a quotient in the case of lifespan moderation and a product for weight. Specifically, we tested whether the addition of an interaction term such as lifespan or weight, improves the relationship between methylation and chronological age. For any site i let M1, M2 and M3 be three separate linear models containing the following terms:

$$M_1:=BC(methylatio{n}_i)=Age+Sex+{\epsilon }_i$$

$$\begin{array}{l}{M}_2:=BC(methylatio{n}_i)=Age+Moderator\\ +Moderator:Age+Sex+{\xi }_i\end{array}$$

$$M_3:=BC(methylatio{n}_i)=Moderator+Sex+{\delta }_i$$

$$\forall i\in \{1,2,\dots ,N\},N:numberofCpGsites$$

$$\begin{array}{l}BC(X)~Box-Coxtransformationofvariable\\ X(methylation)\end{array}$$

$$\begin{array}{l}Residuals{\epsilon }_i,{\xi }_i,{\delta }_i~N(0,\Sigma ),\\ 0:vectorof0s;\Sigma :phylogenyvariance\\ -covariancematrixa\end{array}$$

The increase in likelihood between the model M2 and M1 as well as the p-value of the interaction term in M2 allowed us to query the linear regression model for biological age signals. Additionally, the interaction term in M2 was compared against the isolated lifespan term in model M3 to rule out an independent effect of lifespan in the regression.

Epigenetic pacemaker and bayesAge

We re-implemented the epigenetic pacemaker model described by Snir et al. [33] in R including cross-validation. Additionally, we implemented a version of the BayesAge [50] inspired age determination algorithm described by Mboning et al. [32]. In brief, for every relevant site, a non-parametric curve relating methylation and age is built in a training cross-validation fold using locally estimated scatterplot smoothing (LOESS) [51]. These curves are queried for the probability of a test sample being any age given its methylation status. This algorithm can be interpreted as a realization of a leave-one-out cross-validation of the methylation-age LOESS regression fit and can be readily extended to n-fold cross-validation by doing bulk prediction in test fold n instead of a single-sample. The final age likelihood of a sample is calculated by multiplying the probabilities for any given age over all sites.

As we do not possess a breakdown of methylated and unmethylated reads, we modified the algorithm to handle M/beta values. Instead of using a (negative) binomial link function between age and methylation, we simply use the LOESS prediction and standard deviation values for a given methylation value to create a pseudo-probability distribution of age (Supplementary Figure 5). Particularly, we model P(Age = x)~N(m, s) where m is the predicted age given a methylation value and s is the standard deviation of the LOESS curve at point m. We note that, contrary to the binomial distribution, this approach accounts for the heteroskedasticity of the LOESS curve but does not account for the within-sample coverage variability as this information is unobtainable from M/beta values.

Cell type deconvolution and gene ontology analysis

We first attempted to use the human blood methylation atlas to perform cell fraction deconvolution on the RRBS and mammalian methylation array datasets. We observed that standard deconvolution algorithms such as DeconRNASeq in R returned similar outputs if a proper methylation matrix or a matrix of random numbers were used as an input (Supplementary Figure 7). To assess the possibility to perform deconvolution using the shared sites between the mammalian methylation array and the human methylation atlas, we reframed the deconvolution as a constrained regression problem. Then, we derived multiple regression p-values which assess whether the complete regression model implies an improvement from an only-intercept model. We used the R quadprog library to estimate the regression coefficients enforcing the restrictions that all coefficients should be positive and sum to one. We observe that deconvolution regression models for blood are not significant in dogs and barely reach significance in humans when applied to the subset of CpGs present in the mammalian methylation array (Supplementary Figure 8). Sheep methylation data [4] and data from a blood capture sequencing methylation assay in Labrador retrievers [12] were added for confirmation.

We performed gene ontology (GO) enrichment analysis using the web server webGestalt [52] on the set of significant genes derived from each dataset. Each ontology search was performed against a background of all genes represented in the corresponding dataset. We used human gene ontology networks in place of the dog, as they are more densely populated with curated terms.

Data and materials availability

All data needed to evaluate the conclusions in this paper are present in the paper, code and/or the Supplementary Materials. The raw data for the three datasets can also be obtained from: SRA: SRP065666, GEO: GSE223748, https://mydata.clockfoundation.org/app/mammalian-consortium-data-browser, https://www.tandfonline.com/action/downloadSupplement?doi=10.1080%2F15592294.2022.2069385&file=kepi_a_2069385_sm0108.zip. The code used to analyze the processed data is available at https://github.com/aserresarmero/bioage_methylation.