The fast–slow continuum of longevity among yellow-bellied toad populations (Bombina variegata): intrinsic and extrinsic drivers of variation

Study area

Life-history traits of yellow-bellied toads Bombina variegata were studied in nine localities corresponding to populations 1–8 in Germany and to population 9 in Austria (Fig. 1). Details on geographic location, local climate, survey years and sampling frequency are given in Table 1. Climate data were obtained from the nearby weather stations Montabaur (50.44°N, 7.81°E, 265 m asl, Germany), Marienberg (50.66°N, 7.96°E, 547 m asl, Germany) of the Deutscher Wetterdienst, and Mariabrunn (48.21°N, 16.23°E, 225 m asl, Austria) of the Zentralanstalt für Meteorologie und Geodynamik.

Populations 1–3 inhabited the former military training area Schmidtenhöhe (area: ca. 700 ha) and population 4 a neighbouring clay pit near Koblenz (Rhineland-Palatinate, Germany; Hantzschmann & Sinsch, 2019), now forming part of the NATURA 2000 network (LNatSchG, 2015). The localities occupied by populations 1–4 were about 0.5–2.4 km distant from each other. Since we did not detect any exchange of toad individuals among the four areas, we considered the breeding assemblage at each locality as an isolated and independent population. The annual migratory range of B. variegata is 20–244 m in Bulgaria (Beshkov & Jameson, 1980), 6–600 m in Romania (Hartel, 2008), 203–732 m in Germany (Jacob, Scheel & Buschmann, 2009), and 300 m in Switzerland (Jordan, 2012) supporting the assumed isolation of populations 1+2, 3 and 4. Populations 5-8 were sampled at working and former clay pits 18–48 km distant from populations 1–4 to supplement longevity estimates based on skeletochronology. Population 9 was located in the north-western part of Lainzer Tiergarten, a nature reserve in the west of Vienna (Austria), which is covered by deciduous forest interspersed with meadows. The area (Hochwiese, ca. 1 km²) surveyed for this study consisted of a large wet meadow where small puddles are formed by wallowing wild boars (Sus scofra), and three streams in its immediate vicinity.

The Ethics committee of FB3, Department of Biology, University of Koblenz–Landau (#Bv 01/V/2004), approved research. Research and collecting permits referring to the directive #92/43/EWG (FFH-directive) and to the German Federal Law of Nature Conservation (Bundesnaturschutzgesetz, BNatSchG) were issued by the Struktur- und Genehmigungsdirektion Nord in Koblenz (425-104.111.0501; 0603; 0702; 1602; 1711), and by the municipality of Vienna (MA22–2752/96, MA22–1451/06 and MA 22–984143-2015-5) in Austria. All applicable institutional and national guidelines for the care and use of animals were followed (Kreisverwaltung Mayen-Koblenz, Az 39183-04).

Field procedures and individual toad identification

At each site surveyed, B. variegata specimens were hand-caught at the local ponds and their terrestrial vicinity on a weekly to three-weekly base during the activity period. Capture sessions consisted in an area-constrained, exhaustive search for post-metamorphic individuals without time limit (duration: 0.5 to 5 h). Sex, snout-vent length (SVL, measured with callipers, rounded to the nearest mm), body mass (measured with an electronic balance, rounded to the nearest 100 mg), and ventral spot pattern of each specimen were recorded for individual recognition. Criterion for the distinction among immature, male or female individuals was the minimum SVL of males possessing nuptial pads. Smaller individuals were considered immature, same- or larger-sized individuals without nuptial pads females. In populations 1–8, the photographing procedure involved placing the toads in a petri dish with a lid containing a foam insert (method adapted from Vörös, Szalay & Barabás, 2007; Gollmann & Gollmann, 2011; Cruickshank & Schmidt, 2017). Inverting the petri dish revealed the ventral pattern, which was photographed using different types of digital cameras. In population 9, toads were held in the hand for photography. Additionally, we collected phalange bones from all adult individuals captured in populations 1–8. For age determination, the 3rd or 4th digit of a forelimb was toe-clipped using ethanol-sterilised scissors. We did not sample recaptured adults (i.e., those lacking a digit) again. The wound was cleaned with 70% ethanol and the individual was released in situ. The digit was stored in 100% ethanol or in 4% formaldehyde solution at room temperature until skeletochronological laboratory processing.

Handling times were less than one minute per animal, and we released individuals at their original capture location once photographs had been taken. Matching of images to identify recaptured individuals was done by a combination of searching a photographic database with all records using the Wild-ID software (Bolger et al., 2012) and subsequent confirmation of identity of the best 20 potential matches by eye (Cruickshank & Schmidt, 2017). In population 9, identification of individuals was exclusively done by eye, partly in the field using a printed catalogue of ventral patterns (Gollmann & Gollmann, 2012).

Demographic life-history traits examined

Statistical analyses

We compared the normally distributed local RRR and CRR estimates following the first year after initial capture using an ANOVA with multiple group comparison by the post-hoc Tukey-HSD test. Since cohort recapture rates decreased in time during the following years at each locality, we tested, if the differences among populations persisted using an ANCOVA with the covariate time (year of recapture − year of initial capture). The effects of toad age (senescence) and study year on CRR were tested for in a 2-factor ANOVA. Calculations were run using the program package Statgraphics Centurion version 18.1.01. The significance level was set at alpha = 0.05.

Longevity in the surveyed populations

Longevity was assessed in all populations. Estimates were obtained using the CMR data and skeletochronological age determination for each adult cohort separately. CMR data served as a conservative estimate of the lifespan by counting the number of consecutive years in which an adult was recaptured and adding two more years for the post-metamorphic period before sexual maturity (Gollmann & Gollmann, 2012). To obtain an exact estimate, this method requires that the period covered by recapture events is longer than local longevity. This was true for some cohorts of populations 1 and 9. Additionally, we calculated a probabilistic estimate of the minimum lifespan of the eldest 1% adults of a cohort for population 9 in which the recapture period was often shorter than local longevity. We fitted a square-root(x) regression model to the corrected recapture rates [P = a + b*sqrt(x)], where P is the probability that an individual is recaptured x years or more after the first capture. Then, the longevity estimate is recapture year x at P = 0.01 plus 2 years for the juvenile period.

Skeletochronological age estimation is based on quantifying the number of dense hematoxylinophilic narrow growth marks (Lines of Arrested Growth, LAGs, each representing a hibernation) separated by faintly stained broad growth zones in the round bones of amphibians (Smirina, 1972). Individuals of up to eight years of age are correctly aged by LAG-counting, whereas the lifespan of older individuals is systematically underestimated, due to the increasing rapprochement of LAGs at the periphery of bones (review: Sinsch, 2015). Standard laboratory protocol included embedding of the phalange samples in Historesin™ (JUNG), cross-sectioning of diaphysis (12 µm) using a JUNG RM2055 rotation microtome, and staining of cross sections with 0.5% cresylviolet. Cross sections were examined under a light microscope for the presence of growth marks at magnifications of 400x using an OLYMPUS BX 50. The first and the senior author assessed the number of LAGs independently. Age is the maximum number of LAGs detected in an individual plus one for the year of capture, i.e., the activity period following the last hibernation. Then, longevity at a given locality is the maximum age that was recorded throughout the whole study period.

Statistical analyses

We compared the log10-normalized longevity estimates using an ANCOVA with the fixed factor locality and the continuous covariate year. To test for group differences, we used multiple group comparison by the post-hoc Tukey-HSD test (Statgraphics Centurion version 18.1.01, significance level alpha = 0.05).

Potential impact factors on survival rates and longevity

Weather (local scale)

Local weather conditions vary between years with cold wet winters and hot dry summers and may affect annual survival rates (Weinbach et al., 2018). To quantify the weather/survival relationship, we retrieved data on air temperature and precipitation from weather stations near to populations 1–4 and 9. We calculated two sets of influence parameters: (1) Annual means of average (T_X), minimum (T_Min) and maximum air temperature (T_Max), and annual sum of precipitation (R_Sum). (2) Average and maximum summer air temperature (June, July, August; ST_X and ST_Max) and the corresponding sum of precipitation (SR_Sum) as a surrogate for summer heat and drought suffered during the activity period in the year of first capture. Analogously, average and minimum winter air temperature (December, January, February; WT_X and WT_Min) and the corresponding sum of precipitation (WR_Sum) served as surrogates of winter severity.

Statistical analyses

We applied a multiple regression analysis to quantify the effect of each weather data set complemented with the condition index (see below) on the CRR estimates. To obtain the smallest significant model, we used the procedure: backward selection at F < 4.0 to include those weather variables into the model, which explained a significant amount of variance observed in the CRR of the age-mixed cohorts (Statgraphics Centurion version 18.1.01, significance level alpha = 0.05).

Climate (continental scale)

To test for a general pattern underlying the fast–slow continuum of longevity throughout the geographical range, we identified published information on longevity estimates for B. variegata basing on CMR or on skeletochronology. CMR studies qualified, if the study period covered at least two successive years and raw data on recapture were included. These criteria were met by Beshkov & Jameson (1980), Gollmann (1981), Kapfberger (1984), Seidel (1988), Seidel (1992), Seidel (1993), Plytycz & Bigaj (1993), Marchand (1993), Plytycz et al. (1996), Sy (1998), Sy & Grosse (1998), Di Cerbo (2001), Gollmann (2005), Dino, Milesi & Di Cerbo (2010), Di Cerbo et al. (2011), Gollmann & Gollmann (2012), Primus (2013), and Cayuela et al. (2019a). Age estimates by skeletochronology are limited to Seidel (1988), Plytycz & Bigaj (1993), and Bülbül et al. (2018). Longevity data were obtained for 21 localities including our own study sites (Tables 1 and 2). As surrogate measures of local climate, we used the combinations of (1) latitude, longitude and altitude (coarse-grained approach), and of (2) average annual air temperature and annual sum of precipitation (fine-grained approach). To obtain empirically measured data on temperature and rainfall, we retrieved data from nearby weather stations (within a radius of 20 km) at a similar altitude (within a range of ±100 m asl) as that of the study locality during the study periods at each site.

Statistical analyses

We used a multiple regression analysis to assess the correlation between log10-normalized longevity and log10-transformed data sets 1 and 2. Power analyses showed that 21 observations allowed for the detection of significant multiple correlation at R² = 0.51 (data set 1) and R² = 0.47 (data set 2), respectively. As study periods reported were short in several case studies, we cannot exclude that published longevity does not represent precisely actual maximum life expectancy. Therefore, we used a supplemental discriminant analysis to distinguish the geographical features of localities inhabited by short-lived (reported longevity less than 10 years) or long-lived populations (reported longevity at least 20 years), i.e., the extremes of the fast–slow continuum. Multiple regression and discriminant analyses were performed using the program package Statgraphics Centurion version 18.1.01, power analyses using G*Power version 3.1.9.2. (Faul et al., 2009). The significance level was set at alpha = 0.05.

Energetic investment to somatic growth and to reserves for winter survival

The allocation of ingested energy to somatic growth and subsequently, the growth pattern (age-size relationship) depend on food availability in the habitat and has consequences for longevity (Kozłowski, 1992; Sztatecsny & Schabetsberger, 2005; Hecht et al., 2019). We quantified the seasonal variation of the condition index as a surrogate measure of the nutritional state at capture for populations 1–4 and 9. The growth patterns of populations 1–9 reflect the long-term investment of food to somatic growth. The condition index of an individual is calculated as the studentized residual of the SVL-mass relationship using a multiplicative model ln(mass) = a + b × ln(SVL), with a = intercept and b = slope (residual index: Băncilă et al., 2010; Scheele et al., 2014; Mikoláš, 2016). Growth following metamorphosis was estimated using the Von Bertalanffy equation (1938): SVL_t = SVL_max − (SVL_max − SVL_met) * e^−k∗t, where SVL_t = average body length at age t; SVL_max = asymptotic body length; SVL_met = body length at metamorphosis; t = number of growing seasons experienced (n LAGs), and k = growth coefficient (i.e., shape of the growth curve). SVL_met varied between 11 mm and 19 mm (Schäfer et al., 2018) and was set to an average value of 15 mm (Di Cerbo & Biancardi, 2010). We analysed age/size data for populations 1–8 based on skeletochronological estimates, and for population 9 based on CMR age data.

Statistical analyses

To detect differences in the nutritional state of toads at the study sites, we analysed the variation of the condition index in a 2-Factor ANCOVA with sex and locality as categorical factors, and activity season (March 26 to October 10, i.e., day 85 to 283) as continuous co-variable. The von Bertalanffy growth model was fitted to the average growth curve using the least square procedure (nonlinear regression). Estimates of SVL_max and k are given with the corresponding 95% confidence interval. If confidence intervals did not overlap in pairwise comparisons, estimates differed significantly at P < 0.05. Statistical procedures were performed using the program package Statgraphics Centurion version 18.1.01.

Pathogen prevalence in short-lived populations

Disease- and parasite-mediated reduction of longevity may affect B. variegata populations (Spitzen-van der Sluijs et al., 2017; Campbell et al., 2018). Therefore, we looked for dead specimens and living individuals showing external signs of disease such as skin lesions, ectoparasites, and abnormal behaviour in populations 1–8. Dead individuals, which were not run over by tanks or cars or with marks of predator attack, were dissected to assess the presence of endoparasites (Sinsch, Kaschek & Wiebe, 2018). The skin and nostrils of Bombina were examined for the presence of eggs and larvae of the blowfly Lucilia bufonivora Moniez, 1876, infecting common toads Bufo bufo at the Schmidtenhöhe. In 2017, we sent phalange samples of 49 individuals collected in populations 1–4 to the Zoological Society of London (Institute of Zoology) to detect infections with the chytrid fungus Batrachochytrium dendrobatidis and Ranavirus.

Prediction 1: Annual survival rates vary locally among years, but the average magnitude differs among populations

We compared recapture rates (X ± SE) of age-mixed adult populations one year following the first capture among populations 1–4 (Schmidtenhöhe, Germany) and 9 (Lainzer Tiergarten, Austria) to quantify the variation of the corresponding survival rates. Recapture rates differed significantly among populations 1–4 (ANOVA: RRR, F_4,31 = 7.25, P = 0.0004; CRR, F_4,25 = 10.49, P = 0.0001). A closer inspection showed that population 9 differed from all other populations (Tukey-HSD test, P < 0.05). Average RRR were 0.32 ± 0.04 (populations 1–4) and 0.52 ± 0.03 (population 9), corresponding average CRR were 0.41 ± 0.06 and 0.71 ±0.03.

When using the number of years after first capture as a covariate of recapture rates (ANCOVA: RRR, F_1,250 = 397.4, P < 0.0001; CRR, F_1,237 = 614.0, P < 0.0001), the difference persisted (ANCOVA: RRR, F_1,250 = 178.6, P < 0.0001; CRR, F_1,237 = 177.4, P < 0.0001) (Fig. 3). Comparing RRR and CRR at each locality, again using the number of years after first capture as a covariate, CRR was significantly greater than RRR in population 9 (ANCOVA: F_1,398 = 52.0, P < 0.0001), but not in populations 1–4 (ANCOVA: F_1,67 = 1.3, P = 0.2611). In summary, annual survival rates varied within each population, but differed significantly in magnitude between the four short-lived populations and the long-lived one.

Prediction 2: Annual survival rates do not decrease with age within a population

To test for senescence in populations 1–4 (Schmidtenhöhe, Germany), we analysed recapture rates of known-age toads, i.e., skeletochronologically studied individuals. Corrected recapture rate did not differ among the age classes (1–5 LAGs; 2-factor ANOVA: F_4,49 = 1.4, P = 0.2687; Fig. 4A), but varied among the survey years (2-factor ANOVA: F_4,49 = 3.2, P = 0.0223). Survival to the next year was significantly higher in 2016 (CRR=0.58) than in 2005 (0.26) and 2017 (0.33; Tukey-HSD test, P < 0.05). Average recapture rate (X ± SE) of known-age adults was 0.38 ± 0.04. Survival rates (SR) estimated from the static life tables did not differ significantly among the main adult age classes (Fig. 4B). Average survival rate was 0.46 ± 0.03. Estimated annual survival rates tended to be slightly higher on average than corresponding CRR, but the 95% confidence intervals overlapped widely: 0.40–0.52 (SR) vs. 0.31- 0.52 (CRR). In summary, we did not detect senescence in the short-lived populations.

Prediction 3: The variation of annual survival rate correlates with the among-years variability of local weather and condition index

In populations 1–4 (Schmidtenhöhe, Germany), the most explanatory regression model based on yearly averages of temperature descriptors and of precipitation regime explained 63.4% (R² corrected for df) of the among-years variation of CRR (Multiple Regression Analysis: F_2,11 = 10.51, P = 0.0044). The model included only the average annual air temperature T_X and average annual maximum of daily air temperature T_Max, whereas precipitation and condition index did not contribute significantly to CRR variation. Using descriptors of winter and summer conditions instead of yearly averages, local weather variation accounted almost entirely for the corresponding variation of CRR (Multiple Regression Analysis: R²(corrected for df)=90.7%, F_5,11 = 22.56, P = 0.0008). The refined model included five weather variables following the equation CRR = -0.165 - 0.209*WT_Min + 0.009*WR_Sum –0.605*ST_X + 0.216*ST_Max –0.002*SR_Sum. Again, condition index did not account significantly for CRR variation. Cold winters and warm summers affected CRR positive, whereas precipitation was of minor importance (Fig. S2). In contrast, variation of weather at population 9 (Lainzer Tiergarten, Austria) did not account for a significant amount of among-years variation in CRR as none of the weather variables correlated significantly (Regression analyses, P > 0.05). Overall, weather regime at the Lainzer Tiergarten (Austria) was more continental than that at the Schmidtenhöhe (Germany), i.e., minima during winter were lower and average summer temperature was higher (Fig. S2).

Prediction 4: Longevity variation among populations does not reflect differential energetic investment to somatic growth

Longevity of age-mixed cohorts varied in all populations among years (Table 3). In the short-lived populations 1–8, variation was from 5 to 8 years (skeletochronological and CMR data), in the long-lived population 9 (CMR data, probabilistic estimate for the cohorts 1998-2008) from 13 to 23 years. In May 2019, the recapture of a female of the 1998 cohort (SVL₁₉₉₈ = 47 mm) in population 9 provided the first empiric evidence for a realized 23-years longevity. Log10-transformed longevity estimates differed significantly among populations (ANCOVA, F_8,35 = 10.2, P < 0.0001), which formed one homogeneous group including populations 1–8 and another represented by population 9 (Tukey HSD-test, P < 0.05). The covariate study year was marginally significant (ANCOVA, F_1,35 = 4.34, P = 0.0471) showing a slight trend of decreasing longevity during the study period.

The condition index (descriptor of nutritional state) did not differ between males and females (2-Factor ANCOVA, F_1,1956 = 0.48, P = 0.4863) and its seasonal fluctuations between spring and autumn did not co-vary with time (2-Factor ANCOVA, F_1,1956 = 1.55, P = 0.2126). In contrast, there were significant differences among the populations (2-Factor ANCOVA, F_4,1956 = 41.35, P < 0.0001; Table 3). The distinction between short- and long-lived populations by annual survival rate and longevity did not correspond to the complex pattern of condition index variation (Table 3).

The age-SVL relationship (descriptor of long-term investment to somatic growth) did not differ significantly in SVL_max or k among short- and long-lived populations (Table 3; Fig. S3). The asymptotic maximum SVL ranged from 46.6 mm to 50.6 mm, the growth coefficient k was similar in all populations. The recalculated growth pattern (model parameters: SVL_max = 52.2 mm [47.5–56.9 mm], k = 0.327 [0.235–0.419]) of the short-lived population at Enez, Turkey (Bülbül et al., 2018), did not differ significantly from the populations studied by us. In summary, the invested energy to somatic growth was similar from Germany to Turkey, and in short- and long-lived populations.

Prediction 5: Longevity variation among populations is related to climate gradients

Longevity estimates demonstrated a fast–slow continuum of populations throughout the geographical range (Tables 2 and 3). In 13 populations (10 in Germany, 1 in Bulgaria, 1 in Switzerland, 1 in Turkey), B. variegata were short-lived, i.e., longevity estimates ranged from 5 to 9 years. In three populations (1 in Austria, 2 in Italy), toads were long-lived with a realized longevity of 20 to 23 years. In-between these extremes of the continuum, we detected five populations (1 in Germany, 1 in Switzerland, 2 in Austria, 1 in Poland) with a longevity ranging from 11 to 15 years, partially overlapping with the within-population range of longevity variation in population 9 (Lainzer Tiergarten, Austria). The log10-transformed maximum longevity estimate per locality (n = 21) did not correlate significantly with local geographical features (multiple regression analysis: R² = 0.137, F_3,20 = 0.9, P = 0.4630), or with local temperature and rainfall regime (multiple regression analysis: R² = 0.199, F_2,20 = 2.25, P = 0.1346). When the separate cohort estimates for populations 1–9 were included, log10-transformed longevity correlated significantly with local geographical features (multiple regression analysis: R² = 0.266, F_3,47 = 6.68, P = 0.0008), and also with local temperature and rainfall regime (multiple regression analysis: R² = 0.115, F_2,47 = 4.05, P = 0.0247). The presumptive significance is an effect of pseudo-replication that disappears, if we consider exclusively the number of localities for the degrees of freedom. As the calculated multiple correlations were all below the detection thresholds for significance (power analyses), we cannot exclude that local climate accounts for up to 26% of longevity variation. Geographical features (data set 1) and temperature and rainfall regimes (data set 2) of the localities inhabited by short- or long-lived population, did not distinguish between the two ends of the fast–slow continuum (discriminant functions: data set 1, Eigenwert = 0.229, Wilks-Lambda = 0.813, Chi-squared = 2.6; df = 3, P = 0.4608; data set 2, Eigenwert = 0.378, Wilks-Lambda = 0.726, Chi-squared = 4.2; df = 2, P = 0.1242).

Prediction 6: Low longevity is associated with high pathogen prevalence

In the short-lived populations 1–8 (Westerwald, Germany), dead specimens were rarely encountered in the field and with the exception of one individual, all (n = 5, 2004-2018) were run over by tanks or cars. A dead female of population 3 found at the shore of a pond in 2017 and without external lesions was checked for endoparasites. The digestive system contained a 35 mm Megalobatrachonema terdentatum female (Nematoda). Mass die-offs of B. variegata or of any syntopic amphibians were never observed. Blowfly (Lucilia bufonivora) eggs were found on the skin of a single female of population 3 in 2016. Indication of a former infection was detected in form of dark spots at the nostrils of another female of population 4 in the same year. There was evidence for the presence of Batrachochytrium dendrobatidis in 6 individuals of populations 1 and 2, and by Ranavirus in one individual of population 2. In summary, pathogens in low prevalence were present in most of the short-lived populations studied.

Sources of local variation in survival rates and longevity

All populations studied varied with respect to these life-history traits confirming predictions 1 and 4. Senescence, i.e., age-specific mortality and the timing of its increase, may be a source of variation in survival rates, as postulated by current theories of aging (antagonistic pleiotropy, disposable soma; Williams, 1957; Kirkwood, 1977; Kirkwood & Austad, 2000). While there is growing evidence for these theories in species with determinate growth (birds and mammals), our data demonstrate that in the short-lived populations the annual survival rate is independent of age. In the long-lived population, the annual survival rate remains high during the entire lifetime of a cohort suggesting that there is no decrease of survival with increasing age. In agreement with a study on palearctic salamanders (Cayuela et al., 2019b), we conclude that senescence is negligible in the variation of annual survival rates confirming prediction 2.

Local temperature and precipitation regimes, i.e., weather, show great among-years variation and may affect survival by affecting water balance and metabolism of amphibians (Hillman et al., 2008). In the short-lived populations of B. variegata, weather variables accounted indeed for almost all of the variability in annual survival rates indicating that duration of the activity period (mild winters, Fig. S2A) and temperature-modulated intensity of activity (cold summers, Fig. S2B) play a major role for survival. This is in agreement with prediction 3 and with variation of survival rates in Capensibufo rosei (Becker et al., 2018). Still, survival rates never reached the level of that in the long-lived population indicating the presence of an intrinsic longevity component defining the overall magnitude of survival. In this population, variation of annual survival rates was unrelated to local weather variability suggesting that major mortality factors differ in nature between localities. For example, Cayuela et al. (2017) propose that antipredator and antipathogen skin secretions and modification of reproduction behaviour, i.e., skipping reproduction under harsh conditions, increase adult survival in a French population of B. variegata. With respect to global climate change, the marginally significant trend to decreasing longevity in our study populations during the past 20 years may suggest the presence of an additional source of variation. Yet, database is still too poor to lean hard evidence to this trend. If the trend proves true, prolongation of the activity period and corresponding longer exposure to predators may underlie the putative reduction of longevity.

The nutritional state is related to the availability and quality of food in the toads’ habitat and determines the availability of energy reserves for metabolism during hibernation and subsequently survival (Hillman et al., 2008). Using the condition index as a rough estimate for nutritional state and reserves (mass contribution of fat body), we did not detect a seasonal change during the activity period at any study site. If the condition index does reflect the amount of fat body (Hecht et al., 2019), metabolic reserves were probably not depleted substantially during winter or increased to the end of the activity season to prepare for hibernation. As the condition index did not correlate with annual survival rates as well, we conclude that the nutritional state of toads was probably not a major source of mortality. Alternatively, the momentary condition index may be too rough to mirror correctly nutritional state in B. variegata.

If there is a trade-off between allocation of resources to reproduction in early life and somatic maintenance (disposable soma theory of aging; Kirkwood, 1977; Kirkwood & Austad, 2000; Räsänen et al., 2008), we expect that the long-term investment into somatic growth should lead to distinct growth pattern in short- and long-lived populations. Yet, we do not find any evidence that toads of long-lived population allocate more resources to somatic maintenance/growth than do toads of short-lived population. Contrary to the predictions of life-history theory, we find that age at sexual maturity is about be same at continental scale. Since sexual maturation requires a threshold SVL of about 30–32 mm, the age of first reproduction is usually the third year of life, exceptionally the second (Gollmann & Gollmann, 2012; Bülbül et al., 2018; this study). In conclusion, local growth pattern seems to be invariant throughout the range of distribution confirming prediction 4.

Pathogen prevalence may considerably affect survival and cause local extinction (Adams, Pessier & Briggs, 2017; Finnerty, Shine & Brown, 2018). For example, infections with Batrachochytrium dendrobatidis and Ranavirus decrease longevity of Palearctic amphibians (Spitzen-van der Sluijs et al., 2017; Campbell et al., 2018). Despite the occasional detection of Bd and Ranavirus in few individuals of the short-lived populations, pathogen prevalence seems too low to indicate substantial disease-driven mortality. Data on parasites are too scarce to be conclusive. Consequently, we have to leave open a possible contribution of pathogens to longevity variation leaving prediction 6 unresolved.

Among the sources of local variation in survival and longevity of B. variegata, the only extrinsic factor identified is local weather variability in the short-lived populations, explaining trait variance almost entirely. In contrast, none of factors evaluated turned out to influence variation in the long-lived population. An additional role of climate change in constraining longevity remains unresolved at present.

Sources of continental variation in longevity

At the continental scale, altitudinal and latitudinal clines in demographic life-history traits are an often-observed pattern in amphibians (Morrison & Hero, 2003; Morrison, Hero & Browning, 2004). Palearctic examples are the widespread anurans Rana temporaria (Sinsch, Pelster & Ludwig, 2015) and Epidalea calamita (Leskovar et al., 2006; Oromi, Sanuy & Sinsch, 2012). Yet, our study demonstrates that the contribution of geographic clines and climate features to the variance of the fast–slow longevity continuum in B. variegata populations is less than 26%, if there is any. Altitude seems to play an insignificant role as exemplified in the Italian populations (Albino, Parco dei Colli di Bergamo; Table 2). In the two populations, longevity is the same (20 years; Di Cerbo et al., 2011) despite of inhabiting habitats at 450 and 1,148 m asl, respectively. We conclude that general climate clines related to latitude and altitude of habitats do not contribute substantially to the marked differences in realized longevity in this species, falsifying prediction 5.

Excluding geographical clines, local climate, resource allocation pattern and pathogens as determinants for the nature of the component determining the magnitude of longevity, i.e., the position of a population within the fast–slow continuum of longevity, the question remains, which factors enable some populations to reach more than 20 years of life expectancy? The key to an answer is probably the fact that longevity in the field reaches almost the life expectancy of captive individuals (Mertens, 1964; Mertens, 1970) indicating that predation plays a minor role in long-lived populations. If we assume that there are no dramatic local differences in predator abundance, we expect that improved antipredator protection of the local toads reduces mortality caused by predators. Skin toxins as chemical defences of amphibians are well known to determine unpalatability, to decrease predation, and to increase survival (Darst, Cummings & Cannatella, 2006; Kowalski et al., 2018; Hettyey et al., 2019).

To stimulate future research, we present here a new hypothesis that can explain the broad range of among-population variation in longevity—the palatability hypothesis. Skin secretions of B. variegata include a cocktail of more or less bioactive peptides such as Bombesin, Bombinins and Bradykinin (Simmaco, Kreil & Barra, 2009). Besides anti-bacterial properties, the main effect of these peptides seems to be an interference with the hormonal system of the predator, as insulin-, gastrin- and TRH-releasing properties in mammals and birds have been evidenced (e.g., Marenah et al., 2004; Gonzalez et al., 2008; Simmaco, Kreil & Barra, 2009; Xu & Lai, 2015). Skin peptides of B. variegata are apparently not toxic in a sense of causing death to predators but probably have disruptive metabolic effects combined with irritation of mucosa. If these effects are uncomfortable enough for the predator, it will learn to avoid this prey. The palatability hypothesis proposes that the antipredator protection by bioactive skin secretions varies among B. variegata populations determining their position in the fast–slow continuum of longevity. If so, longevity variation reflects mainly the effect of predation on more or less chemically protected toads.

Investigations on the quantity of bioactive substances in the skin secretions of B. variegata are still lacking. Yet, there is ample evidence that the concentrations of Bradykinin in Bombina pachypus, of Bufogenins in Bufo gargarizans and of Tetrodoxins in Notophthalmus viridescens vary considerably among conspecific populations (Zhang et al., 2005; Yotsu-Yamashita et al., 2012; Ma et al., 2018; Coppari et al., 2019). Multiple, still unknown factors may cause this variation, from environmental variability to differences in the local gene pool (Bókony et al., 2019). Future measurements of skin peptide profiles and their biological activity on potential predators will show, whether there are significant differences among individuals of short- and long-lived populations.