Four decades of phenology in an alpine amphibian: trends, stasis, and climatic drivers

Strong phenological shifts in response to changes in climatic conditions have been reported for many species, including amphibians, which are expected to breed earlier. Phenological shifts in breeding are observed in a wide number of amphibian populations, but less is known about populations living at high elevations, which are predicted to be more sensitive to climate change than lowland populations. The goal of this study is to assess the main factors determining the timing of breeding in an alpine population of the common toad (Bufo bufo) and to describe the observed shifts in its breeding phenology. We modelled the effect of environmental variables on the start and peak dates of the breeding season using 39 years of individual-based data. In addition, we investigated the effect of the lunar cycle, as well as the individual variation in breeding phenology. Finally, to assess the individual heterogeneity in the timing of breeding, we calculated the repeatability of the timing of arrival at the breeding site. Breeding advanced to earlier dates in the first years of the study but the trend continued only until the mid 1990s, and stabilised afterwards. Overall, toads are now breeding on average around 30 days earlier than at the start of the study period. High temperatures and low snow cover in winter and spring, as well as reduced spring precipitation were all associated with earlier breeding. Additionally, we found evidence of males arriving on average before females at the breeding site but no clear and strong effect of the lunar cycle. We only found weak evidence of among-individual variation in shifts in the breeding phenology, as well as a low repeatability of arrival timing. Our findings show that the observed changes in breeding phenology are strongly associated with the environmental conditions. These results contribute to filling a knowledge gap on the effects ssof climate change on alpine amphibian populations. Moreover, we show that changes in phenology, especially in the mountains, can be hard to predict as local microclimatic conditions do not necessarily reflect the observed global climatic trends.

at the breeding site but no clear and strong effect of the lunar cycle. We only 52 found weak evidence of among-individual variation in shifts in the breeding 53 phenology, as well as a low repeatability of arrival timing. 54 55 5. Our findings show that the observed changes in breeding phenology are 56 strongly associated with the environmental conditions. These results 57 contribute to filling a knowledge gap on the effects of climate change on 58 alpine amphibian populations. Moreover, we show that changes in phenology, 59 especially in the mountains, can be hard to predict as local microclimatic 60 conditions do not necessarily reflect the observed global climatic trends. 61 Introduction and the environment. More specifically, our goal is to (i) identify the environmental 140 variables (e.g., temperature, snow cover, moon cycle) that could be driving the 141 observed breeding phenology of this population (both the start and the peak of the 142 breeding season), (ii) analyse if there is significant variation in the phenological shifts to breed at the study pond. We then marked (first by toe-clipping, then starting in 155 1993 by implanting PIT tags), measured, and released them in the same place 156 (Hemelaar, 1988;Grossenbacher, 2002). To make sure we captured both early and 157 late arrivers, we repeated this procedure for on average 5-6 nights, with breaks in-158 between of about 2-4 days (i.e, the data conform to Pollock's (1982) robust design). 159 The length of the fieldwork period usually covers the breeding season duration, 160 which typically lasts about two weeks at our study pond. This design also had the 161 advantage of not overly stressing the toads. In total, for the period 1982-2020, 3053 162 uniquely recognizable individuals have been caught, of which 1852 were males and 163 1201 females. For each individual we have a record of presence for each capture 164 night over the study period. Given the reduced size of the pond and the repeated 165 capture rounds within a capture night, we assumed high capture probabilities 166 (capture probability p ≈ 0.85 per year based on a preliminary analysis of the mark-167 recapture data). At the population level we determined for each year a start, a peak, 168 and an end date of breeding (i.e., first capture night, the capture night when most 169 toads were captured, and last capture night, respectively). These calendar dates 170 were all transformed into days of the year (where January 1st = 1), to facilitate

Population trend 194
Moreover, to check how the standard deviation (SD) of the start or the peak 217 breeding dates changes over time, we calculated for both start and peak breeding 218 the SD of the residuals of each of the 1000 piecewise regressions, using a rolling 219 window approach (with a 10-year window) with the function rollapply of the package 220  In addition to these five climatic variables, the lunar cycle has also been 252 identified to be an important factor for the timing of breeding in amphibians, with in  To quantify both the effects of climate and of the moon cycle on the breeding 261 phenology, we modelled two separate linear regressions on the day of the breeding 262 start and the day of peak breeding over the period 1982-2020. As explanatory 263 To further study the association between the moon cycle and breeding 272 phenology we tested if start and peak breeding tended to happen more frequently 273 under certain moon phases. To do this, we used the rayleigh.test function of the 274 circular R package (Agostinelli & Lund, 2017)  angles. Also in this case we first ran the tests on the originally assigned dates and 281 then we ran them on 1000 simulated datasets of start and peak breeding dates and 282 obtained the 2.5 th and 97.5 th percentile of the p-values. uniquely marked individuals, as many individuals were breeding in multiple years 294 (mean = 2.21 years, SE = 0.02)), using a linear mixed model (package lmerTest; 295 assigned arrival dates, and then, to account for uncertainty in the assignment of the 297 dates of arrival to the pond we simulated 1000 new datasets where every individual 298 arrival date is newly sampled from an uniform distribution and can be as early as the 299 capture night preceding the original arrival date, or if it was the first capture night of 300 the season, up to seven days before. Using these 1000 new datasets we ran 1000 301 models and obtained the 2.5 th and the 97.5 th percentile values for each parameter. 302 As a random effect, applied on both the intercept and the slope of both PC1 and 303 PC2, we included individual identity (ID). This was done not only to observe if 304 individuals react differently to changing environmental conditions, but also to account 305 for the non-independence of the data. Moreover, we also included year as a random 306 effect on the intercept, to account for unexplained year-specific variation in the data. 307 Finally, we included the effect of sex to account for differences between males and 308 female. To properly be able to compare the effects of continuous variables (i.e., the 309 two PCs and the cosine of the lunar angles) with the effect of a categorical variable 310 (i.e., sex), we standardised the three continuous variables by subtracting the mean 311 and dividing by two times the standard deviation (Gelman, 2008). Finally, as a 312 measure of model fit, we calculated the conditional R 2 value using the 313 r.squaredGLMM function from the package MuMIn (Barton, 2019). Finally, we also estimated repeatability (i.e., the upper limit of heritability) of arrival 318 dates at the breeding site. High values of repeatability (r) mean that individuals are 319 consistent in their relative arrival timing (e.g., always among the first ones), and vice 320 versa. To calculate r, we used for each individual the date of first capture for each year that it was captured. This date is a relatively good proxy for the date of arrival at 322 the breeding site, as the data collection usually starts every year approximately when 323 the first toads arrive at the pond. The date was converted to the day of the year 324 (where January 1st = 1), and then standardised by subtracting the year-specific 325 mean and dividing by the year-specific standard deviation. We then used the 326 function rpt from the package rptR to calculate r using individual ID as the group 327 variable (Stoffel et al., 2017), and bootstrapping 1000 times to obtain the 95% CI. 328 As for all the other analyses, to account for the uncertainty in the assignment of the 329 dates, we repeated the calculation of r 1000 times, sampling different arrival dates

Population trend 344
Both the breeding start dates and the dates of peak breeding show very similar 346 trends (Pearson's correlation coefficient = 0.91), with both also showing marked 347 between-year variation over the study period. Nonetheless, a shift towards earlier 348 breeding dates is observable, with breeding happening now on average around 30 349 days earlier compared to the start of the study period ( Figure 1). The piecewise 350 regression on breeding start dates identified a single breakpoint in the temporal trend 351 in the year 1993 with a pre-1993 steep advancement of breeding dates followed by a 352 post-1993 almost flat trend ( Figure 1A; Table 1). The analysis of the robustness of 353 the piecewise regression, done by simulating data and running 1000 piecewise 354 regressions, performed very similarly, with 910 cases out of 1000 where the year 355 1993 was identified as breakpoint and the model coefficients were very close to the 356 piecewise regression conducted on the originally assigned breeding dates (Table  357   S1). The piecewise regression on peak breeding dates also identified 1993 as a 358 breakpoint year ( Figure 1B; Table 1). In this case, the analysis of the robustness 359 showed slightly more variation, with the breakpoint years mostly obtained being 1993 360 and 1996 (274 and 283 out of 1000 respectively) (Table S1). Moreover, we found the 361 standard deviation (SD) of the residuals of the piecewise regressions on both start 362 and peak breeding dates to vary considerably, with higher SDs at the start and the 363 end of the study period ( Figure S2). To further check the pattern in the residuals we 364 split them in four different decades and checked their distribution ( Figure S3).  384 Table S4 in the Appendix shows the summary of these five piecewise regressions.

Determinants of variation in the breeding phenology in the population 386
The first two principal components (PC) of the principal component analysis (PCA) 388 described together more than 70% of the variation in the data, and both had a 389 standard deviation (i.e. the squared root of their eigenvalue) above one ( Figure S1; 390 Table S5). Therefore, applying the Kaiser rule, we kept the scores of these two PCs 391 (PC1 and PC2) as explanatory variables in the following linear regressions on the 392 start of the breeding season and on peak breeding (also including the scaled cosine 393 of lunar angle). PC1 was mostly determined by winter temperature (+0.45 loading) 394 and winter and spring SWE (-0.61 and -0.64, respectively). PC2 was mostly 395 determined by spring weather conditions. Spring temperature had a negative loading 396 (-0.68), while precipitation had a positive loading (+0.68) ( Figure S1; Table S6). 397 Regarding the start of the breeding season, the model (adjusted R 2 = 0.41) 398 indicated a significant negative relationship with PC1 and a significant positive 399 relationship with PC2 ( Table 2). The cosine of the lunar angle had a non-significant 400 effect. Similarly, for the regression on the dates of peak breeding, we found a 401 significant negative relationship with PC1 and a significant positive relationship with 402 PC2, while the cosine of the lunar angle had a small and non-significant effect (Table  403 2). The adjusted R 2 was 0.54. In both cases the outcome is that warmer 404 temperatures in winter and spring, less snow cover, and weaker precipitations are all 405 associated with an earlier start and peak of the breeding season. Both the 1000 406 linear regressions on the simulated dates of the start of the breeding season and the 407 1000 on the simulated dates of peak breeding performed similarly to the two 408 regressions on the originally assigned dates (Table S2), indicating that our analysis 409 is robust to possible imperfect assignment of dates of start and peak breeding. Rayleigh's test on the moon phases on breeding season start and on peak dates. In 416 both cases we obtained a non-significant p-value (0.27 and 0.08 respectively), 417 indicating that we could not confidently reject the null-hypothesis of the data being 418 uniformly distributed in the circular space. In addition, the outcome of the Hermans-419 Rasson test for multivariate deviations indicated that the null hypothesis could not be 420 rejected for both start and peak breeding (p-value = 0.38 and 0.21 respectively). To

Determinants of individual variation in breeding phenology 436 437
To better understand if there are among-individual differences in the phenological 438 response to changing climatic variables, we used a linear mixed model to test for the 439 effect of climatic variables on the individual breeding start dates (i.e., the date on 440 which an individual was first captured). We found only a small difference in the 441 response of breeding phenology to climatic variables among individuals (i.e., low 442 values for the random effect ID, both on intercept and slopes, Table 3). We found a 443 strong significant positive effect of PC2 on the breeding dates (17.51 ± 3.27 SE), 444 meaning that stronger precipitation and lower minimum spring temperatures are 445 associated with a delay in the breeding. We also found a significant and strong 446 negative effect of PC1 (-10.14 ± 2.85 SE), indicating that colder winter temperatures 447 and higher SWE are associated with a delay in the breeding. We also found a 448 significant but weak effect of the cosine of the lunar angle (1.57 ± 0.14 SE), 449 suggesting a possible small role of the lunar cycle. Finally, we observed an effect of 450 sex indicating that males arrived on average earlier than females (-1.45 ± 0.14 SE) 451 (Table 3). The 1000 models on the 1000 simulated datasets, ran to assess the 452 robustness of the analysis to imperfect assignment of arrival dates, showed a similar 453 outcome to the main model (Table S3).

Discussion
Our results show that variation in the breeding phenology is strongly associated with 467 climatic conditions, which vary substantially among years but also show trends 468 across times. We also found low repeatability values and low variability in individual data, as our measure of precipitation included both snow-and rainfall. We found that 486 a higher amount of precipitation in spring (combined with a decrease of spring 487 temperature) was associated with a later breeding date. In fact, at low temperatures, 488 precipitation in the form of snowfall or freezing rain can delay the melting of the snow cover, therefore leading to a delay in the breeding. The observed negative 490 association between snow water equivalent (SWE) and breeding timing is in line with 491 the rest of the findings. In fact, SWE depends considerably on temperatures and 492 precipitation, as well as other aspects such as exposition, and it is a key factor that 493 influences phenology (Corn, 2003). The very similar trend observed for peak activity 494 in breeding indicates that both start and peak breeding are influenced mostly in the 495 same way by the same climatic variables. 496 When looking at the individual timing of arrival we still found an important effect on 498 the breeding phenology of PC1 (TW and SWESp/W) and PC2 (TSp and PrecSp) ( Table  499 3). However, we found only non-significant and small among-individual variation in 500 phenological response to changing climatic conditions (Table 3). As reproduction 501 happens only once a year in explosive breeders living in temperate zones, 502 synchronisation in breeding could be key to maximise reproductive output (Ims, 503 1990). Such an accurate synchronisation can be achieved more easily when all 504 individuals hibernating close to each other express similar responses to external 505 cues triggering their migration to the breeding pond, instead of responding 506 individually in different ways, highlighting once more that the breeding phenology is 507 mainly driven by climatic conditions. 508 Moreover, the low values of r (i.e., the upper limit of heritability) that we found Finally, we found that on average males tend to arrive earlier than females (Table 3), 530 similarly to what has been found in lowland populations of B. bufo (Loman & 531 Madsen, 1986;Höglund & Robertson, 1987; but see Gittins et al., 1980). In 532 these studies, males, especially bigger ones, were observed to arrive on average 533 earlier at the breeding pond. Smaller males, on the other hand, were observed 534 intercepting females on their way to the pond, betting on the fact that the females 535 would lay the eggs as soon as they arrived at the pond, avoiding competition from 536 the other bigger males. A more detailed future analysis of body size and its effects 537 on the timing of migration to the breeding site could confirm this theory also for our 538 study population.
Climate change is leading to on-average increasing temperatures both globally but 540 also at smaller scales such as in the European Alps (Vitasse et al., 2021)  Despite the observed stabilisation of the trend of the breeding dates (Figure 1), the 583 study population appears to experience increased variation in the dates of the start 584 of the breeding season ( Figure S2 and Figure S3). This increased variation could be 585 explained by extreme weather events whose occurrence is expected to increase