Biogeographic history moulds population differentiation in ageing of oxidative status in an amphibian Free

Ageing is a ubiquitous process of decline in physiological function, survival and fecundity with advancing age (de Magalhães and Passos, 2018). In recent times there has been growing interest in elucidating the physiological mechanisms of ageing in free-living animals or non-model laboratory organisms (Gaillard and Lemaître, 2020). The oxidative stress theory of ageing postulates that age-associated reductions in physiological functions are caused by slow steady changes of oxidative status with age (e.g. accumulation of molecular oxidative damage), that are associated with life expectancy of organisms (Finkel and Holbrook, 2000; Partridge and Gems, 2002). However, we still know relatively little about age-related changes in oxidative status markers within individuals in wild animals. This research gap is particularly unfortunate if we consider the bourgeoning interest in understanding the functional and fitness consequences of within- and among-individual variation in oxidative status (Costantini, 2014). A pillar of this research is that the expression of oxidative status traits might integrate the individual history or genetic background and underlies life-history trade-offs (Costantini, 2014). For example, variation of oxidative status markers may be linked to immune function (Sorci and Faivre, 2009), the expression of secondary sexual traits (Alonso-Alvarez and Galvan, 2011), personality traits (Herborn et al., 2011), chronic stress exposure (Hau et al., 2015) or demographic traits (Costantini et al., 2015). This emerging view of oxidative status biology as an integrative research ground suggests potential implications of physiological ageing in eco-evolutionary processes that mould the integrative phenotypes, such as dispersal.

Dispersal is a key context-dependent process that contributes to shape spatial patterns of biodiversity at all levels of organization, from genes to ecosystems (Hanski and Gilpin, 1991; Hansson, 1991; Bohonak, 1999; Ronce, 2007; Clobert et al., 2012; Saastamoinen et al., 2018). Recent theoretical and empirical evidence, however, suggests that dispersal could also depend on individual state or phenotypic quality, which can either influence or be influenced by the evolution of virtually all the features of an organism, including behaviour, morphology, physiology and genetics (Bonte et al., 2012; Clobert et al., 2012; Canestrelli et al., 2016a,b).

Importantly, dispersal out of the range-edge of an expanding population can impose ecological pressures on dispersers, resulting from novel biotic interactions (Sakai et al., 2001; Lee and Klasing, 2004; Brown et al., 2013) or abiotic conditions (e.g. temperature, humidity) at which they need to physiologically adapt (Wilson et al., 2000; Johnston and Temple, 2002; Steinhausen et al., 2008). Several lines of evidence have shown that metabolic activity (Myles-Gonzalez et al., 2015; Louppe et al., 2018) or immune function (White and Perkins, 2012; Cornet et al., 2016) may differ between edge and core populations. This is particularly relevant because both metabolic rate and immune function may affect certain traits of oxidative status (Costantini, 2014), strengthening the idea that oxidative status and its role in physiological ageing might vary among populations across an expansion gradient.

In this study, we hypothesized that, if individual variation in physiological ageing reflects variation in the capacity to cope with the new demographic and/or environmental conditions encountered during a range expansion, a spatial structure in oxidative status might emerge along the range of a historically expanded population. We explored this hypothesis by means of a longitudinal common garden experiment lasting 1 year, using the Tyrrhenian tree frog Hyla sarda as the study species. This tree frog underwent a major northward range expansion during the Late Pleistocene, whose time line, source area and expansion route have been well characterized by previous studies, using multiple lines of evidence (Bisconti et al., 2011a,b; Spadavecchia et al., 2020 preprint). Specifically, we studied the spatial differentiation along the inferred route of the species' past expansion in four blood-based markers of oxidative status, including damage level and both non-enzymatic and enzymatic antioxidants.

The Tyrrhenian tree frog Hyla sarda (De Betta 1853) is a small, cryptically coloured amphibian endemic to the Tyrrhenian islands of Sardinia, Corsica, the Tuscan archipelago (Capraia and Elba), and some neighbouring islets (western Mediterranean Sea). Among the European tree frogs, it is the most closely linked to breeding ponds, spending most of the year in the adjacent vegetation, although it can tolerate even prolonged periods of aridity (reviewed by Lanza et al., 2007). Hyla sarda shows a lek mating pattern, with most breeding sites observed in lowland and mid-altitude areas, although it can be found from the sea level up to more than 1500 m.a.s.l. (Lescure and de Massary, 2012). Breeding choruses can be heard from late March to late July (Lescure and de Massary, 2012; Lanza et al., 2007). Adult tree frogs mostly feed on arthropods, and are predated by various species of birds and snakes (Lanza et al., 2007). A study of the sexual size dimorphism in H. sarda (Cadeddu et al., 2012) showed that, while a positive relationship between age and size was observed for both sexes, males and females differ both in reproductive life span (1–4 years in males and 2–5 years in females) and body size (on average up to 32.8±2.4 and 38.2±3.2 mm in males and females, respectively). Age classes of breeding adults also differ between sexes (Cadeddu et al., 2012), with most males (80%) being 2–3 years old, and most females (71%) being 3–4 years old.

Previous phylogeographic (Bisconti et al., 2011a; Spadavecchia et al., 2020 preprint) and population genetic studies (Bisconti et al., 2011b) showed that H. sarda colonized Corsica and the Tuscan archipelago from an ancestral area located along the central-eastern coast of Sardinia. This range expansion event probably occurred in two steps (Spadavecchia et al., 2020 preprint). The first step allowed the species to colonize northern Sardinia at the end of the last glaciation, taking advantage of a wide – albeit transient – land bridge connecting Sardinia to Corsica throughout the glacial epoch. In the second step, H. sarda colonized the whole of Corsica and the other islands, following the post-glacial amelioration of climatic conditions in these northern islands. Remarkably, early in the post-glacial epoch, the transient land bridge connecting Sardinia to Corsica was lost, preventing subsequent gene flow between the two islands (Bisconti et al., 2011a).

To minimize potentially confounding effects, our sampling design was carefully planned along the following strategies: (i) we limited collection to adult males only (snout–vent length, mean±s.d.: 34.8±2.8 mm), thus controlling for sex, life stage and, at least in part, age (see Cadeddu et al., 2012); (ii) despite the fact that this species is common in freshwater environments at almost any altitude, we restricted our sampling sites exclusively to coastal pools, thus avoiding any effect due to environmental factors associated with altitude; (iii) sampling was carried out in two geographic areas per island (Fig. 1; Table S1), in two breeding sites per geographic area, and in a minimum of two distinct pools within each site, thus controlling for both local habitat effects and demographic aspects such as kinship or inbreeding; (iv) the sampling area was limited to a coastal strip at the eastern side of the species range, thus minimizing bioclimatic differences between geographic areas and across time epochs (current versus last glacial maximum) – in fact, in an earlier species distribution model, this side obtained, by comparison, higher bioclimatic suitability scores under both current and past bioclimatic conditions (Bisconti et al., 2011a); and (v) we intentionally excluded from sampling the northern islands of Capraia and Elba, where the species is also present, because they were most likely colonized through jump dispersal rather than through a spatial diffusion process, and thereby via founder events and consequent genetic drift (Bisconti et al., 2011a,b).

A total of 55 males of H. sarda were included in this study. Tree frogs were collected with hand nets following mating calls at night in late May to early June, transported into our laboratory and housed in individual cages (25 cm×25 cm×25 cm) in a climatically controlled room at a temperature of 24–25°C (monitored using a Hobo MX2201 temperature data-logger), relative humidity of 60–80%, and natural photoperiod, following suggestions from Staniszewski (1995). Cages were provided with a dechlorinated water tank, a plant and oak wood as shelter, and were placed randomly with respect to population origin. Twice a week, cages were cleaned, water was renewed, and tree frogs were fed with crickets (Acheta domestica; see Udomsil et al., 2019, for information on cricket nutritional content).

Sampling procedures were performed under the approval of the Institute for Environmental Protection and Research ‘ISPRA’ (protocol no. 5944), Ministry of Environment ‘MATTM’ (protocol no. 8275), Regione Sardegna (no. 12144) and Corsica (nos 2A20180206002 and 2B20180206001). Permission to temporarily house amphibians was granted by the Local Health and Veterinary Centre with licence code 050VT427. All handling procedures outlined in the present study were approved by the Ethical Committee of the University of Tuscia for the use of live animals. After transportation and housing of the frogs in the experimental facilities, we ascertained that day–night activity rhythms were normalized (higher activity at night time) by hourly visual scan sampling for 2 weeks before any experiment started; observations also included reactivity to food and calling activity at night. Moreover, we measured the body condition index at the onset of the experiment (following Baker, 1992; mean±s.d.: 0.003±0.012 and 0.002±0.013 for Corsican and Sardinian frogs, respectively) and monthly mass variation. No adverse effects on the overall health of tree frogs were observed during the procedures. The animals were released in the original sampling locations at the end of the experiment.

Tree frogs were allowed to acclimate to captivity conditions for 1 month before the beginning of experimental procedures. Blood sampling was carried out in the morning (09.00–12.00 h), 24 h after the last meal, by cardiac puncture with a 0.3 ml syringe. Blood samples were collected from 55 individuals (32 from Sardinia and 23 from Corsica). One year later, another sample of blood was collected from 46 individuals that were still alive (31 from Sardinia, 15 from Corsica). Before the collection of blood, animals were anaesthetized by immersion for 2 min in a solution of MS-222 (tricaine methanesulfonate, 0.5% m/v). Blood samples were immediately centrifuged at 10,000 rpm for 3 min in order to separate the serum from the erythrocytes into two different tubes that were immediately stored at −80°C until laboratory analysis. After blood collection, each individual was placed in a wet box; we returned frogs to their cages when they were fully awake.

Four markers of oxidative status were measured according to established protocols for vertebrates (Costantini et al., 2011a,b). Briefly, we used the d-ROMs test (Diacron International, Grosseto, Italy) to quantify oxidative damage in serum. The ROMs (reactive oxygen metabolites in serum) concentration is expressed as mmol l⁻¹ of H₂O₂ equivalents. The OXY Adsorbent test (Diacron International) was used to quantify the ability of serum non-enzymatic antioxidant compounds (vitamins, carotenoids, uric acid, thiol proteins) to react in vitro with a known concentration of hypochlorous acid (HOCl). OXY (serum non-enzymatic antioxidant capacity) values were expressed as mmol l⁻¹ of HOCl neutralized mg⁻¹ protein. The Ransel assay (Randox Laboratories, Crumlin, UK) was used to measure the activity of the antioxidant enzyme glutathione peroxidase (GPx) in haemolysates; the activity of GPx was expressed as U mg⁻¹ protein. The Ransod assay (Randox Laboratories) was used to quantify the activity of the antioxidant enzyme superoxide dismutase (SOD) in haemolysates; the activity of SOD was expressed as U mg⁻¹ protein. The concentration of proteins in haemolysates was quantified using the Bradford test. All analyses were run in duplicate, and the mean coefficient of variation ranged between 6.88% and 8.82%.

Statistical analyses were run using R (version 3.5.3; http://www.R-project.org/). Generalized linear mixed-effect models (GLMMs; lme4 package) were used to compare frogs from the two islands for each blood-based oxidative status marker. In each model, we included island (Sardinia and Corsica), sampling period (beginning and end of the experiment) and their interaction as fixed factors; individual was included as a random factor to control for non-independence of multiple measurements. As dependent variables, we entered each marker singly; we also ran a model using the SOD/GPx ratio as a dependent variable because prior work on humans suggested that an unbalanced ratio of the two enzymes (i.e. higher values) might be associated with cellular ageing (De Haan et al., 1996). This is because an overabundance of SOD may convert more superoxide free radical to hydrogen peroxide, which in turn is quenched by GPx. However, other enzymes (e.g. catalase) also contribute to protect against oxidative damage, and thus other types of ratio can also prove valuable to estimate physiological senescence (De Haan et al., 1995). A Gamma distribution and a log link function were used for ROMs, OXY and SOD/GPx, while a Gamma distribution and an identity link function were used for GPx and SOD. We relied on the Akaike information criterion (AIC) to identify the distribution type and the link function that improved model fitting.

Preliminary models showed that entering sampling location within each island as a random factor did not improve the fit of the model (i.e. the AIC value was not reduced beyond 2); thus, it was not subsequently considered. Preliminary models also showed that model fitting based on the whole dataset was better than that based on the subset of individuals for which we had data for both the pre- and post-experimental period (data not shown).

GLMM-based repeatability models (glmer in package ‘lme4’) were then used to quantify the within-individual repeatability of all oxidative status markers from the first to the second sampling period. Gamma distribution with a log link function was used. The observation-level variance was obtained by using the trigamma function (Nakagawa et al., 2017). We entered each marker singly as a dependent variable and individual as a random factor. Preliminary models showed that the adjusted repeatability calculated by entering location as a fixed factor did not improve the fit of the models.

We compared the mortality percentage observed from the first to the second sampling of blood between islands using [P=(p_s×n_s+p_c×n_c)/(n_s+n_c)] in Statistica version 10 (StatSoft, Tulsa, OK, USA), where P is the proportion of dead frogs, n is the sample size, s is Sardinia and c is Corsica. We also tested whether there was any selective mortality between the first and second sampling period of individuals from a given island having different values of each marker at the beginning of the experiment as compared with the overall distribution. We ran a generalized linear model (glm in package lme4) with a binomial error distribution and a logit link function that included a given marker as main factor and survival probability (0=dead, 1=alive) until the end of the experiment as the response variable. These analyses were restricted to Corsica because only one frog from Sardinia died before the second sampling of blood.

Finally, we analysed bioclimatic features of the collection sites, and their potential correlation with the observed pattern of variation of oxidative status markers. We selected four bioclimatic variables deemed important for H. sarda, based on the available literature (Lanza et al., 2007): BIO 9 (mean temperature of the driest quarter), BIO 11 (mean temperature of the coldest quarter), BIO 14 (precipitation of the driest month) and BIO 19 (precipitation of the coldest quarter). We extracted values of these variables for each sampling site from Chelsa Climate, a database of climatic variables averaged over about 35 years (1979–2013) and provided at ∼1 km resolution (Karger et al., 2017). Geographic variation of these bioclimatic variables was reduced by extracting the first principal component (explaining 78% of the variation), in a principal component analysis (loadings of variables onto PC1, based on correlations: BIO 9=0.80, BIO 11=0.94, BIO 14=−0.86, BIO 19=−0.93; all P<0.05) and used this variable as a bioclimatic index.

We found a significant interaction between island and sampling period for ROMs, OXY, GPx and SOD/GPx ratio (Table 1). Post hoc analyses by Tukey test showed that the islands had similar values of all tested markers at the beginning of the experiment (Fig. 2). Conversely, at the end of the experiment, Corsican tree frogs had significantly higher ROMs (P=0.013) and SOD/GPx ratio values (P=0.022), and lower OXY (P=0.005) and GPx (P=0.038) than those from Sardinia (Fig. 2). From the beginning to the end of the experiment, ROMs decreased in Sardinian frogs (P=0.003) and increased in Corsican ones (P=0.028); GPx increased in Sardinian frogs (P=0.007) and decreased in Corsican ones (P=0.015); and OXY did not change in Sardinian frogs and decreased in Corsican ones (P<0.001). SOD increased significantly in both Sardinian (P=0.006) and Corsican frogs (P=0.038) over the course of the experiment; the SOD/GPx ratio increased significantly in Corsican frogs (P=0.001) (Fig. 2). All markers of oxidative status were significantly repeatable over the course of the experiment; estimates of repeatability coefficients (and the associated observation level variance) were as follows: ROMs: 0.40 (0.017); OXY: 0.25 (0.061); GPx: 0.48 (0.139); SOD: 0.26 (0.160); SOD/GPx: 0.36 (0.216).

We recorded a higher mortality of frogs from Corsica (8 out of 23, 34.8%) than for those from Sardinia (1 out of 32, 3.1%; P=0.002) at the end of the experiment (i.e. sampling time 2). We do not know the cause of death; however, variation in oxidative status markers was not associated with survival. Values of each oxidative status marker at the first blood sampling time did not predict the mortality probability until the second blood sampling time in frogs from Corsica (P≥0.41), indicating that results were robust for any potential bias introduced by selective mortality of individuals with either higher or lower values of a given marker at the first sampling period.

Values of all markers at the first sampling period were not significantly correlated with the bioclimatic index (ROMs, r=−0.05; OXY, r=−0.02; GPx, r=−0.05; SOD, r=0.23; SOD/GPx, r=0.23), indicating a negligible effect of local climate in explaining among-site variation in oxidative status.

Using a longitudinal common garden experiment, we show a significant differentiation in age-related changes of oxidative status markers between populations of tree frogs along a historical range expansion route from Sardinia (core population) to Corsica (edge population). At the beginning of the experiment, we found that Sardinian and Corsican frogs had similar concentrations of all oxidative status markers and that this variation was not explained by bioclimatic characteristics of the sampling sites. One year later (i.e. at the end of the experiment), Sardinian and Corsican frogs differed significantly in most markers analysed. Specifically, frogs from Corsica had significantly higher ROMs (a marker of oxidative damage) and lower antioxidant molecules (OXY and GPx) than those from Sardinia. In contrast, we found that expression of the antioxidant enzyme SOD was increased at the end of the experiment in frogs from both islands. Also, Corsican frogs had higher SOD/GPx ratio values than Sardinian frogs. In humans, higher values of the SOD/GPx ratio indicate faster cellular ageing (De Haan et al., 1996); however, SOD and GPx were not correlated (r=0.17, P=0.13), indicating that a direct functional interaction between these two enzymes might have been small or that additional enzymes (e.g. catalase; De Haan et al., 1995) would have been necessary to measure in order to have a more robust ratio of antioxidant enzymatic activity. We also found that frogs from Corsica suffered higher mortality than frogs from Sardinia over the course of the experiment. Our results were robust for any selective mortality of frogs having particular values of markers (very high or very low) compared with the overall distribution because oxidative status markers were not associated with mortality probability. Also, frogs were maintained under the same common garden conditions, thus eliminating any noise due to environmental variation. Overall, these results suggest a differentiation in the pace of physiological ageing between conspecific populations that have been geographically and genetically isolated since the end of the Pleistocene epoch (Bisconti et al., 2011a).

Mitochondrial activity is one important source of the free radical superoxide anion that, in turn, is a precursor of other pro-oxidants, such as hydrogen peroxide (Halliwell and Gutteridge, 2015). SOD, an enzyme that protects against the superoxide anion, was upregulated over the course of the experiment in frogs from both islands. Thus, we might exclude a strong role of mitochondrial superoxide anion in explaining our results if the level of SOD expressed by frogs from the two islands was similarly effective in controlling pro-oxidant activity of superoxide anion. The upregulation of SOD with time is consistent with prior longitudinal and cross-sectional studies on captive zebra finches (Marasco et al., 2017) and free-ranging cheetahs (Costantini et al., 2017), indicating that exposure to certain free radicals might increase with age as shown by research on laboratory strains (Costantini, 2014).

Alternatively, Sardinian and Corsican frogs might have differed in the activity of the hypothalamic–pituitary–interrenal axis that releases corticosterone (Sapolsky et al., 2000; Romero, 2004). Corticosterone promotes a suite of metabolic processes that increase gluconeogenesis, catabolism of energy reserves and inhibition of energy storage, and increase in locomotory and foraging activity (Sapolsky et al., 2000). Corticosterone can also induce changes in oxidative status, especially when animals are being exposed to increased levels of corticosterone for long periods of time (Costantini et al., 2011a,b). In the house sparrow (Passer domesticus), individuals from the range edge produced more corticosterone, possibly in response to their exposure to novel environmental conditions (Liebl and Martin, 2013; Martin and Liebl, 2014). Corticosterone has also long been considered a promotor of dispersal by stimulating activity in birds and, recently, in movement-related survival in cane toads (Belthoff and Dufty, 1998; Jessop et al., 2018). Further work will be needed to elucidate the role of corticosterone as a pacemaker of age-related changes of oxidative status in tree frogs.

The results of our work raise the important question of why these population differences in oxidative status have evolved. One explanation might lie with the legacy of a past positive selection on those life-history strategies that were more successful during the range expansion. Rollins et al. (2015) found a substantial upregulation of genes involved in metabolism and cellular repair, including mechanisms of response to oxidative stress, in cane toads at the edge of an invasion front. These results suggested that dispersing cane toads might have prioritized investment in dispersal capacity and rapid demographic expansion rather than in self-maintenance mechanisms (Rollins et al., 2015). Although the cane toad system is very different from ours, the strategy of giving priority to protection of dispersal capacity would be successful whether the costs of reducing protection against oxidative damage are low as compared with those incurred to sustain active defence mechanisms. This explanation appears to be congruent with the disposable soma hypothesis of ageing, which states that optimization of trade-offs in the allocation of limited resources among self-maintenance and other activities inevitably comes at a cost for the soma in the long term (Kirkwood, 1977; Kirkwood and Rose, 1991). Similarly, the antagonistic pleiotropy hypothesis of ageing posits that natural selection would favour the evolution of life histories in which alleles have beneficial effects on fitness early in life even when they are detrimental in old individuals (Williams, 1957; Nussey et al., 2013). The cost of ageing would be low because it would emerge at a phase of life (old age) when the force of natural selection is expected to be low. Our longitudinal data covered a non-negligible portion of the expected lifespan of a tree frog. The selection pressures produced by extrinsic mortality factors are thought to play a major role in the evolution of ageing rates (Monaghan et al., 2008). The risk of death from extrinsic causes, such as disease or predation, during dispersal is particularly important (Clobert et al., 2012). Investing resources in a robust and long-lasting body might be maladaptive if extrinsic causes of mortality are high. In these circumstances, a much better strategy, in evolutionary terms, would be to invest more resources into traits that favour dispersal. Thus, natural selection might have favoured frogs that invested more into mechanisms that favour dispersal and colonization, and this strategy came at a cost for the soma (higher oxidative stress and mortality) of Corsican frogs in the long term. While there is no reason to think that predators would differ between the source and the newly colonized areas, predation risk could have indeed played a role in population differentiation, as, being density dependent, it is unlikely to have been constant during the range expansion phase. Unfortunately, we have no information about past or current predation pressure in our sampling sites. Also, we lack robust data about among-population variation in life-history traits; thus, we cannot infer about the adaptive meaning of these population differences in oxidative status and in mortality rate.

Variation in climatic conditions among sites might also be very relevant because the oxidative status of ectotherms is sensitive to variation in abiotic factors, like temperature or humidity (e.g. Costantini, 2014; Dupoué et al., 2020). Selection owing to different climatic regimes might have significantly affected the oxidative status of frogs after the dispersal came to an end. Also, this climate-induced selection could have affected the response of frog populations to our laboratory climatic conditions. However, this scenario is not well supported by our data because values of oxidative status markers after 1 month of acclimation to our housing conditions were similar among populations and were not correlated with the bioclimatic conditions recorded in the study sites over the past 35 years.

Finally, we found a moderate individual repeatability of all measured markers (0.25–0.48) from the first to the second sampling period. Estimates of within-individual repeatability of oxidative status markers are rarely reported in the literature and the few available are in line with our estimates (e.g. Saino et al., 2011; Récapet et al., 2019). The mechanistic reasons for and biological meaning of these intrinsic individual differences in oxidative status markers have not been elucidated yet. We argue that, together with the clear spatial structure we found, this significant individual consistency in oxidative status suggests that the pace of physiological ageing might be heritable and that selection on oxidative status might contribute to drive the evolution of individual reaction norms.

In conclusion, our experiment provides evidence for a role of a historical range expansion in explaining the different pace of physiological ageing between populations of tree frogs. These results also indicate that, while age-related changes in SOD might be a shared (or conserved) mechanism (sensu Partridge and Gems, 2002) of ageing, those in other oxidative status markers (ROMs, serum non-enzymatic antioxidants and GPx) might be private (i.e. peculiar to particular evolutionary lineages) because they were expressed only in Corsican frogs. The possible conserved SOD expression across the two islands is also supported by the lack of genetic variation between Sardinian and Corsican frogs in the gene SOD-1 (Bisconti et al., 2011b). Further studies will be needed to elucidate the endogenous mechanisms responsible for the observed differences in oxidative status, and the correlated functional and fitness consequences, and whether such differences are owing to genetic variants or developmental plasticity. In particular, it will be important to characterize the life histories of the different populations in order to infer about the adaptive meaning of the differences in oxidative status we observed.

We warmly thank all the collaborators that assisted in housing and sampling: Andrea Chiocchio, Mauro Zampiglia, Armando Macali, Alessandro Carlini, Cinzia Mastrogiovanni, Ambra Pazzani and Lorenzo Latini. We also thank two anonymous reviewers for providing valuable comments that helped to improve the presentation and interpretation of our results.