On the comparative biology of mammalian telomeres: Telomere length co‐evolves with body mass, lifespan and cancer risk

Telomeres, the short repetitive DNA sequences that cap the ends of linear chromosomes, shorten during cell division and are implicated in senescence in most species. Telomerase can rebuild telomeres but is repressed in many mammals that exhibit replicative senescence, presumably as a tumour suppression mechanism. It is therefore important to understand the co‐evolution of telomere biology and life‐history traits that has shaped the diversity of senescence patterns across species. Gomes et al. previously produced a large data set on telomere length (TL), telomerase activity, body mass and lifespan among 57 mammal species. We re‐analysed their data using the same phylogenetic multiple regressions and with several additional analyses to test the robustness of the findings. We found substantial inconsistencies in our results compared to Gomes et al.'s. Consistent with Gomes et al. we found an inverse association between TL and lifespan. Contrary to the analyses in Gomes et al., we found a generally robust inverse association between TL and mass, and only weak nonrobust evidence for an association between telomerase activity and mass. These results suggest that shorter TL may have been selected for in larger and longer lived species, probably as a mechanism to suppress cancer. We support this hypothesis by showing that longer telomeres predict higher cancer risk across 22 species. Furthermore, we find that domesticated species have longer telomeres. Our results call into question past interpretations of the co‐evolution of telomere biology and life‐history traits and stress the need for careful attention to model construction.

Cross-species analyses have furthered understanding of evolutionary trade-offs shaping telomere length (TL) and telomerase activity (TA). Haussmann et al., (2003) showed that telomere loss rate, but not TL, was inversely correlated with maximum lifespan in birds, which has since been corroborated across tetrapod species (Dantzer & Fletcher, 2015;Pepke & Eisenberg, 2020;Sudyka et al., 2016;Tricola et al., 2018;Whittemore et al., 2019). Seluanov et al., (2007) found that TA was negatively correlated with body mass, but not lifespan in rodents. Since, all else being equal, increased TA is expected to lead to slower telomere erosion and thus longer lifespan (Bodnar et al., 1998), Seluanov et al.'s results apparently contrast with the results on telomere shortening rates reported above. Seluanov et al. also found no significant associations of TL with body mass or lifespan. In an analysis of 57 mammalian species, Gomes et al., (2011) reported a negative relationship between TA and mass, but not lifespan, across mammals. Furthermore, Gomes et al. found that lifespan was inversely correlated with TL while accounting for body mass, but there was no independent association between TL and mass. Here we build on these past findings, first trying to replicate the analyses of Gomes et al., then proceeding to further examine how telomere biology may influence cross-species variation in cancer risk and finally how telomere biology may have been shaped by selection related to domestication.
Telomere biology is widely assumed to influence cancer risk and thus for different species to evolve different telomere dynamics in response to differing fitness costs of cancer. Cancer occurs in almost all vertebrates, but among wild animals, cancer incidence rates have been little studied, but seem to vary considerably (Albuquerque et al., 2018;Boddy et al., 2020;Effron et al., 1977;Pesavento et al., 2018). For instance, cancer is rarely reported in long-lived whales, bats or naked mole-rats, while rats and mice are known to be much more prone to cancer (Albuquerque et al., 2018). Different species have evolved different anticancer strategies and mechanisms; for example, long-lived species that reproduce late in life are expected to be under stronger selection pressure to evolve efficient cancer defences compared to short-lived species (Seluanov et al., 2008Weinstein & Ciszek, 2002). Despite larger and longer lived species having more cells that exist for a longer duration, they do not show higher cancer rates than smaller, short-lived species (Peto et al., 1975). This phenomenon, known as "Peto's Paradox," is generally assumed to be explained by species with larger body sizes and longer lives having increased fitness costs from cancer development.
Thus, longer telomeres and higher TA are predicted to increase the risk of cancer (e.gArtandi & DePinho, 2010;Aviv et al., 2017).
However, there are reasons to question the associations between TL and cancer risk (reviewed by Eisenberg & Kuzawa, 2018). Shorter telomeres can cause chromosomal instability that promotes cancer.
Longer telomeres improve immune function, which is important in combatting many infectious diseases which cause cancers and in combatting the development of malignancies. To better characterize the association between telomere biology and cancer, using crossspecies data sets we examine whether longer telomeres and TA predict increased cancer incidence, as is often assumed.
In both mice (Mus musculus and Peromyscus leucopus) and pearl millet (Pennisetum glaucum), domesticated strains have been observed to have longer telomeres than their wild counterparts (Bickle et al., 1998;Hemann & Greider, 2000;Kotrschal et al., 2007;Manning et al., 2002;Seluanov et al., 2008;Sridevi et al., 2002;Weinstein & Ciszek, 2002). Multiple explanations for this putative association between domestication and longer TL have been suggested (Eisenberg, 2011;Manning et al., 2002;Weinstein & Ciszek, 2002). The data set generated by Gomes et al., (2011) included nine domesticated species out of 57 species, which may affect the potential associations between TL, body mass and lifespan (Pepke & Eisenberg, 2020). To better characterize how general the association between domestication and TL is, we examine this association across the cross-species data set.
Our analyses primarily build on a cross-species data set generated by Gomes et al., (2011 Gomes et al., (2011) measured mean TL in kilobases in adult individuals using the telomere restriction fragment (TRF) method. TA was measured in cultured fibroblasts relative to an adenocarcinoma cell line (%H1299). Following Gomes et al., (2011) we used the best species-level mammal supertree from Bininda-Emonds et al., (2007), which is available within the r package "PhyloOrchard" (O'Meara et al., 2013). The phylogeny was pruned to match the data set using the r package "ape" (Paradis & Schliep, 2018). We used the phylogenetic generalized least squares method (Grafen, 1989) implemented as the "pgls" function in the package "caper" (Orme et al., 2018) to perform the phylogenetic regressions. The phylogenetic signal, Pagel's λ, was optimized numerically within the default bounds (0-1) using maximum likelihood. We used estimates of adult body mass (kg) and maximum lifespan (years) reported in Gomes et al., (2011), which were derived from the AnAge: The Animal Ageing and Longevity Database (Tacutu et al., 2018). Because some of these were different from current estimates in AnAge, we compiled new values for a separate re-analysis (Table S2). Following Gomes et al., (2011), TL, mass and lifespan were all log 10 -transformed. We performed two multivariate phylogenetic regressions (outlined in figure   2 in Gomes et al., 2011) of body mass (response variable) predicted by TA and TL, and of maximum lifespan (response variable) predicted by TA, TL and body mass. The phylogeny and TL evolution were visualized using the "phytools" package (Revell, 2012) in which ancestral states were estimated following Felsenstein (1985) but using maximum likelihood with the function "fastAnc."

| Associations between neoplasia incidence and telomere biology
We compiled neoplasia incidence rates in seven species of mammals from a medically curated postmortem database published by Boddy et al., (2020) in addition to 15 species from various sources in the literature reviewed by Albuquerque et al., (2018) and shown in Table   S12. In all studies, the number of individuals with abnormal growth including both benign and malignant tumours (neoplasia) out of the total number of necropsies were used to define the neoplasia prevalence rate (%). For the large treeshrew (Tupaia tana) neoplasia rate is unknown, so we used the estimate reported for the closely related common treeshrew (Tupaia glis). Due to small samples sizes, neoplasia data were combined in Boddy et al., (2020) for the Asian and African elephant (Elephas maximus and Loxodonta africana), so we used the average TL, mass and lifespan of the two species, which are similar (Table S11), in subsequent analyses. Boddy et al., (2020) also provided estimates of lifetime prevalence of malignant neoplasms (invading surrounding normal tissue) of 29 mammal species, allowing us to test if neoplasia prevalence is generally a good predictor of malignancy prevalence (i.e., cancer rate [%]). Most neoplasia data are from wild animals kept in captivity with no age or sex information, and not all tissues were investigated during necropsies, so the data only provide crude estimates of the natural rate of spontaneous cancers.
Species mean TL measured using the TRF method allows comparisons of absolute TLs and were collected from both Gomes et al., (2011) andSeluanov et al., (2007) to increase sample size (Table S12).
Six overlapping species in these two studies show that the TL measurements are highly correlated (R 2 =.9431), but consistently higher in Seluanov et al., (2007). We therefore used an ordinary linear regression of TLs of the overlapping species to adjust measurements of the four species (Table S12) from Seluanov et al., (2007) to match Gomes et al., (2011) using the formula: 0.6617 × TL − 2.4454 kb. TA was measured using different methodologies, so we can only report whether TA was present (1) or absent (0) in the somatic cells of the species. We then test phylogenetic bivariate associations between neoplasia rate and TL, telomerase expression, adult body mass and maximum lifespan (obtained from the AnAge database, Tacutu et al., 2018), and whether TL and/or telomerase expression predict neoplasia rates (response variable) using phylogenetic multiple regression as described above.

| Effect of domestication on telomere length
We assigned a binary variable (0/1) indicating the nine species in Gomes et al., (2011) that have a long history of human domestication and captive-breeding (Scherf, 2000;Wilkins et al., 2014; see Table   S11): cow (Bos taurus), sheep (Ovis aries), pig (Sus scrofa), dromedary camel (Camelus dromedarius), horse (Equus caballus), dog (Canis lupus familiaris), European white rabbit (Oryctolagus cuniculus), laboratory mouse (Mus musculus) and laboratory rat (Rattus norvegicus). We first tested whether TL (response variable) is different in domesticated species using a phylogenetic bivariate regression. Since TL is also hypothesized to be influenced by TA, body mass and lifespan, we also tested the effect of domestication on TL (response variable) while accounting for TA, body mass and lifespan (predictor variables) using phylogenetic multiple regression, as described above.

| Sensitivity analyses
Due to the relatively small sample sizes of the data sets (57 and 22 species, respectively), which may produce uncertainty in the parameter estimates, we performed sensitivity analyses to species (re-)sampling using a jackknifing method. We randomly removed a proportion of species (around 10%, 20%, 30%, 40% and 50%) and then refitted the regression models described above 500 times (for each reduced data set) to estimate sensitivity of β-values to sample size. The power to reliably detect phylogenetic signal may be low for small sample sizes (Münkemüller et al., 2012) and λ was thus calculated for each simulation. These procedures are implemented in the function "samp_phylm" in the "sensiPhy" package (Paterno et al., 2018). All analyses were performed in r (R Core Team, 2020).

| Replication of analyses in Gomes et al., (2011)
Presented in Table 1 are the results reported by Gomes et al., (2011) and our re-analysis-both using two phylogenetic multiple regression models. Consistent with Gomes et al., (2011), we found that both mass and TL significantly predicted lifespan. However, the remaining results are all substantially discordant. While Gomes et al., (2011) found no association between TA and lifespan, we observed a significantly positive association. In predicting body mass, Gomes et al., (2011) found that TA was a significant predictor, but TL was not. We found the opposite: TA did not predict mass, but longer TL predicted decreased body mass. Using current estimates for mass and lifespan from the AnAge database (Table   S11) did not qualitatively alter the findings described above (Table   S2). Gomes et al., (2011) reported a series of bivariate plots (their figure 4) for key variables in their analyses along with p-values.
These p-values are identical to those presented in their figure 2 (our Table 1), strongly suggesting that these significance values were from phylogenetic multiple regressions, not phylogenetic bivariate regressions. A casual reader is likely to interpret these p-values as bivariate associations. Here, we present phylogenetic linear bivariate regression analyses of TL, body mass, maximum lifespan and TA to examine the sensitivity of the multiple regression results and more clearly define these bivariate associations.
The bivariate results (Figure 1) are qualitatively consistent with our multiple regression results, except lifespan did not show any association with TA ( Figure 1d). Based on these bivariate results, a 1% increase in body mass predicts a 0.08% decrease in TL (Figure 1a), while a 1% increase in lifespan predicts a 0.43% reduction in TL ( Figure 1b). The co-evolution of TL and body mass is illustrated in the phylogeny in Figure 2.
In an effort to further reconstruct the reasons for our discrepant results with those of Gomes et al. and explore the sensitivity of our analyses to modelling strategies, we experimented with alternative transformations of the TA measurements (see Supporting Information). Briefly, we find that in multiple regression models, as TA is first log-transformed (Table S3) and then binary transformed (Table S4), TL becomes less and less predictive of either mass or lifespan and the correlation between TL and TA increases (Table S5). When the bivariate association of mass predicting binary TA (0/1) is modelled in a phylogenetic logistic regression (see Supporting Information), complete telomerase repression, which was found in 39 out of 57 species, is significantly associated with higher body mass with an inflection point at 0.8 kg ( Figure S1). Among the 18 species with nonzero TA (see Table S11 these are our best guess at a reconstruction. More complete results from our re-analyses are given in Table S1.

F I G U R E 1
Associations between telomere length, maximum lifespan and body mass (all log 10 -transformed) and telomerase activity (untransformed) across 57 mammal species. Bivariate phylogenetic regressions of these associations are shown. Scatter plots do not depict the phylogenetic controls. Domesticated species are shown in grey there were no significant associations between (untransformed) TA and mass or lifespan (Table S6). Inspection of the model residuals and quantile plots of the associations between TA, mass and lifespan (Figure 1c,d) indicated significant outliers violating the multivariate normal distribution assumptions of the models. To identify these outliers within the phylogenetic regression, we ran a post hoc leave-one-out deletion analysis that detects influential species on parameter estimates while accounting for phylogeny using the function "influ_phylm" in the r package "sensiPhy" (Paterno et al., 2018). After removing the most influential species (see Supporting Information), we reran the main analyses, which revealed a significant negative bivariate association between TA and mass (53 species, β log(mass) = −0.5190 ± 0.2056 SE, p = .0148) that was, however, not robust to controlling for the effect of TL (Table S7).

| Associations between neoplasia incidence and telomere biology
We found that the prevalence of neoplasia is a good predictor of malignancy prevalence (R 2 = .90, β = 0.7735 ± 0.0489 SE, p < .001, Figure S2, see also Figure S3) across 29 mammal species included in Boddy et al., (2020). Because neoplasia prevalence was available for more species, analyses proceeded using this variable.
As expected, neoplasia rate was positively associated with TL ( Figure 3a) and binary TA (Figure 3d, marginally significant trend) across the 22 mammal species shown in Table S11. Thus, a 1% increase in TL predicts an absolute increase in neoplasia incidence of 0.21% (Figure 3a). Neoplasia rate was negatively associated with lifespan ( Figure 3b). A 1% increase in maximum lifespan predicts an absolute decrease in neoplasia incidence of 0.14% (Figure 3b). Gomes et al., (2011) illustrating the co-evolution of telomere length (in kb, red: shorter telomeres, blue: longer telomeres) and body mass (log 10 -transformed, corresponding to the size of the blue spheres). Domesticated species are marked in grey. The timescale shows divergence estimates in million years ago (Mya) Species with some TA (1) are predicted to show an absolute increase in neoplasia incidence of 20.6% (Figure 3d). There was no association between neoplasia rate and mass (Figure 3c). There was only a marginally significant positive association between TL and neoplasia rate, when accounting for binary TA in a phylogenetic multiple regression (Table 2), probably because TA and TL are correlated and some of the variation in TL is explained by TA (Table S5).

F I G U R E 2 Phylogeny of 57 mammal species included in
Neoplasia rates could be biased towards higher estimates in domesticated species due to the easier observation of large numbers of individuals (Boddy et al., 2020;Ewald & Swain Ewald, 2015;Nagy et al., 2007), which could be associated with increased lifespan and senescence in captivity (Tidière et al., 2016). We therefore tested, in a post hoc analysis, whether the effect of TL on neoplasia rate was robust to controlling for potential domestication effects on neoplasia rates (Table S8) or if the effect of TL varied with domestication status (Table S9). TL remained significantly positively associated with neoplasia rate when controlling for domestication status and there was no effect of domestication status on neoplasia rate (Table   S8). The effect of TL on neoplasia rate did not vary with domestication status (Table S9).

| Effect of domestication on telomere length
The arithmetic mean (±SD) TL was 24. Excluding the nine domesticated species from the replication analyses presented in Table 1 did not qualitatively alter our findings (Table   S10).

| Sensitivity analyses
The statistically significant bivariate associations between TL and mass ( Figure 1a) and TL and lifespan (Figure 1b) were generally robust to sample size effects (mean change in β with 49% of the species included was 12%-17%) and the associations remained significant within most of the reduced data sets (85%-100%, The bivariate associations between neoplasia rate and TL, lifespan or TA (22 species, Figure 3a,b,d) were relatively robust to sample size effects, but the significance of the associations decreased considerably as sample size was reduced, with ≥32% corresponding to a simulated sample size of just 15 species (Figures S10-S12). The association between neoplasia rate and mass remained nonsignificant under most simulations ( Figure S13). The multiple regression of neoplasia rate predicted by TL and binary TA was relatively robust to sample size effects ( Figure S14). The phylogenetic signal was not robust in the sensitivity analyses of neoplasia rate (Figures S10-S14).

| DISCUSS ION
Our re-analysis of the data from Gomes et al., (2011) demonstrates that TL has evolved to be shorter with increasing body mass and lifespan, across 57 species of mammals spanning 15 orders ( Figure 2). These results were corroborated by sensitivity analyses of sample size effects ( Figures S4 and S5). This supports the hypothesis that short telomeres evolved in large, long-lived mammals to enable replicative senescence as a mechanism to ameliorate cancer risk conferred by a larger number of cells and longer time to accumulate oncogenic mutations (Gorbunova et al., 2014;Risques & Promislow, 2018). Alternatively, because TL was measured in adults, larger species may erode their telomeres faster than small species, giving rise to the negative correlation between TL and body mass (Monaghan & Ozanne, 2018). However, larger species have been shown to have lower telomere shortening rates across nine tetrapod species, but this association was driven by inclusion of a domesticated mouse strain (Pepke & Eisenberg, 2020). Because telomere loss rates tend to be slower in longer lived species (Dantzer & Fletcher, 2015;Haussmann et al., 2003), telomeres at conception are likely to be even longer in shorter-lived species than estimated in adulthood.
In contrast to our findings, previous studies have shown no significant associations between TL with mass or lifespan across 15 rodent species (Seluanov et al., 2007) or between TL and lifespan of 19 bird species (Tricola et al., 2018). These differences might be explained by these rodents and birds being relatively small and short-lived compared to many of the species included here. Among eight species of relatively large mammals and birds (all >0.5 kg), longer early-life TL predicted longer lifespan (the opposite of our findings), but only when accounting for body mass and telomere shortening rate (Pepke & Eisenberg, 2020 These analyses also showed large uncertainty in the estimation of the phylogenetic signal, which should be interpreted with caution for these small sample sizes. Our re-analysis of Gomes et al., (2011) showed that TA had only inconsistent associations with lifespan and body mass (Table 1; Tables S1-S4), which call into question past results suggesting that high TA has been selected against in larger and longer lived species. We found some hints that TA was inversely associated with mass when telomerase was either log-transformed or binary transformed (Tables S3 and S4 and Figure S1), or when highly influential species (outliers) were excluded from the analyses (Table S7). Thus, the associations between telomere biology, mass and lifespan are sensitive to modelling strategies. While we found an inverse association of binary TA with mass, the biological meaning of binary coded TA is unclear. Perhaps having any TA increases the level of cancer risk (Shay et al., 2001) variable. Furthermore, small levels of TA in normal stem cells is thought to be insufficient to fully maintain TL Sun et al., 2019).
TL was a strong predictor of neoplasia rates across 22 species spanning 11 orders (Figure 3a). Although this conclusion may be limited by a small sample size, the association was relatively robust to the exclusion of a smaller proportion of the species (Table S10). Neoplasia rates are a good proxy for malignancy rates ( Figure S2). In turn, malignancy rates probably reflect the incidence of development of cancers with fitness consequences (Ewald & Swain Ewald, 2015). Thus, natural selection against longer telomeres may be driven by cancer-related mortality. Our findings are consistent with the hypothesis that shorter telomeres have evolved as a cancer resistance mechanism across mammals. Larger and longer-lived species did not have higher neoplasia incidence rates (Figure 3b,c), which is a demonstration of the well-known Peto's paradox. The negative correlation that we found between TL and body mass may partly underlie Peto's paradox; that is, larger species have evolved shorter telomeres as an anticancer mechanism to compensate for otherwise higher cancer risks (Caulin & Maley, 2011). The optimal TL in mammals may therefore have been shaped by a trade-off between the increased cancer risk of long telomeres and the increased susceptibility to age-dependent degenerative diseases conferred by short telomeres (Aviv et al., 2017;Tian et al., 2018). We also found some evidence that binary TA was associated with higher cancer rates (Figure 3d), although this effect was not significant when controlling for TL (Table 2). High TA is a prerequisite for most cancer cells  and may decrease the threshold for cell immortalization , leading to increased cancer risk in species with somatic TA. However, cancer risk does not depend on TA per se, but on TL maintenance (Borah et al., 2015), which may require different TA levels, or occur through alternative mechanisms, across different species (Tian et al., 2018). Although species with short telomeres have lower cancer incidence rates, some species with long telomeres appear to have lower cancer rates than expected given their TL. This may be explained by the evolution of various telomere-independent cancer defence mechanisms in these typically long-lived species Erten & Kokko, 2020). The most striking outlier in our analysis is the tiger (Panthera tigris), having a large body size (162 kg) despite having very long telomeres (50 kb), high TA (10%), but a moderate neoplasia incidence rate (18%), suggesting the potential for discovery of novel tumour suppressor mechanisms within this species.
We provided the first evidence that telomeres are longer across multiple domesticated species spanning five orders (mainly members of Artiodactyla, but also of Perissodactyla, Carnivora, Rodentia and Lagomorpha; Figure 2). However, excluding domesticated species from our analyses did not qualitatively alter the general associations between telomere biology and life-history traits (Table S10).
Long telomeres may be a part of a "domestication syndrome" in mammals (Darwin, 1868;Wilkins et al., 2014), although the mechanisms generating this trait change are not well known (Eisenberg, 2011). Alternatively, long telomeres were ancestral in domesticated species, but somehow connected to the early stages of domestication (Lord et al., 2020). We found that domestication has resulted in an estimated 38.5% increase in TL (Table 3), which is equivalent to the effect that a 98.4% decrease in body mass would have on TL (Figure 1a). There does not appear to be a single directional change in body sizes of domesticated animals (Lord et al., 2020).
However, domesticated species may have experienced changes related to fertility, reproductive lifespan, growth rate, degree of inbreeding, parasite load and energetic constraints (Diamond, 1997;Lord et al., 2020), which are all effects known to shape intraspecific variation in TL (Bebbington et al., 2016;Giraudeau et al., 2019;Manning et al., 2002;Monaghan & Ozanne, 2018;Sudyka, 2019;Weinstein & Ciszek, 2002). Comparisons of TL and other phenotypes in domesticated species with their wild or feral cousins might provide more insights into the particular causes of TL changes with domestication.
Our re-analyses of the valuable data set generated by Gomes et al., (2011) showed considerably divergent results from those shown by Gomes et al., (2011) and that many results are not robust to varying modelling strategies. Unlike Gomes et al., (2011), we showed that TL decreases with increasing body mass. While Gomes et al.'s results suggest a reduction in TA with increasing body mass, we only found weak evidence for this. Our results suggest that the co-evolution of telomere biology with lifespan and body mass across species needs to be reconsidered. The negative association between TL and mass and lifespan could be explained by the increased cancer risk we found in species with longer telomeres. Finally, we found that domesticated species appear to have evolved longer telomeres-suggesting new avenues to understand TL evolution through closer examination of the domestication process and artificial selection experiments.

ACK N OWLED G EM ENTS
We thank Amy Klegarth for collaboration on earlier exploration of these data and Jerry Shay and Chris Venditti for collegial discussion

DATA AVA I L A B I L I T Y S TAT E M E N T
All data are available from the Dryad Digital Repository (Pepke & Eisenberg, 2021: https://doi.org/10.5061/dryad.1c59z w3v0).