Data and R script for paper 'The role of climate and species traits on the timing of spring emergence in amphibians and reptiles'

We downloaded citizen observations of amphibians and reptiles from the GBIF database (gbif.org). We included human observations (i.e., excluding machine observations, material citations and samples, and preserved specimens) from the year 2010 to 2024. Only records with valid coordinates, coordinate uncertainty of < 5 km, and located at a latitude > 30° N were retained. Moreover, we removed duplicate observations, defined as species records with identical latitude and longitude on the same date. We created 100 × 100 km grid cells across the study area and assigned each observation to a grid cell ID. After excluding year-round active species, we estimated the mean day of the year for the 5th percentile of observations per species and grid cell, defined as spring emergence day, using the R package ‘Phenesse’ (Belitz et al. 2020). To validate the use of citizen observations for estimating spring emergence timing, we compared the obtained mean spring emergence day to published literature. To evaluate the correlation between spring emergence based on literature and citizen observations, we used a linear mixed effect model of the r package ‘lme4’ with emergence day based on citizen observations as response variable, emergence date from literature as predictor, and species as random effect. We then assigned ecological traits to each species using the AmphiBIO database for amphibians (Oliveira et al. 2017) and the ReptTraits dataset for reptiles (Oskyrko et al. 2024). For amphibians, traits included maximum body size and offspring development mode (larval versus direct development). For reptiles, the traits were maximum body mass and reproductive strategy (oviparous, ovoviviparous, viviparous). Additionally, for anurans, urodeles, and testudines, we compiled information on whether species generally (i) hibernate on land or in water, and (ii) are explosive breeders, defined by a short breeding period and scramble competition, or prolonged breeders, defined by longer breeding periods (Wells 1977). We downloaded the minimum and maximum monthly temperature (in °C), and monthly precipitation (in mm) for the years 2010 to 2019 (i.e., coinciding with citizen observations) at 2.5 arc minute resolution (277.5 km × 196 km at 45° latitude; not available at finer resolution) from the ‘WorldClim version 2.1’ dataset (Fick and Hijmans 2017), and the average annual temperature at 30 arc second resolution (0.93 km × 0.65 km at 45° latitude) from the ‘CHELSA-SWB’ dataset (Karger et al. 2023). From the monthly raster layers, we calculated the minimum and maximum temperature and cumulative precipitation for the months of February to April (hereafter defined as spring), i.e., when most species emerge from brumation (Table S1), averaged across years. Moreover, we estimated the temperature range during spring, defined as the mean of monthly maximum minus minimum temperatures. Elevation was downloaded at 3 arc second resolution from SRTM data (https://www.earthdata.nasa.gov/sensors/srtm). We then extracted the elevation and climate variables for each observation and averaged them per species and grid cell to be included in the analyses. We analyzed the factors affecting variation in spring emergence (day of the year; response variable) on grid cell level. To account for phylogenetic relatedness among species, we used Bayesian Markov Chain Monte Carlo phylogenetic generalized linear mixed models of the R package ‘MCMCglmm’ (Hadfield 2010). Phylogenetic data were retrieved using the R package ‘rotl’ (Michonneau et al. 2016), generating a phylogenetic tree for all species. In the resulting trees, branch lengths were estimated using the R package ‘ape’ (Paradis et al. 2019), which was also used to create an inverse variance-covariance matrix for subsequent analyses. To account for spatial autocorrelation of the data, we computed spatial distances among the grid cell centroids (latitude and longitude) and performed a principal coordinate analysis (PCoA) on the resulting distance matrix. The PCoA axes were then included as random effects in the MCMCglmm model to account for spatial autocorrelation. Moreover, species was included as a random effect to account for multiple observations. We initially ran an analysis for all species combined, including taxonomic group (see below), elevation (log-transformed), mean annual temperature, spring temperature range, and cumulative precipitation in spring as fixed effects, as well as the interaction of taxonomic group and mean annual temperature (the most important variable; see results) to quantify if the different groups respond differently to geographic variation in temperature. Latitude was highly correlated with mean annual temperature and spring temperature range (Pearson correlation coefficient r: >-0.7) and thus not included as fixed effect. Similarly, mean annual temperature was highly correlated with minimum and maximum spring temperature (Pearson’s r > 0.9), which were consequently excluded from the analysis. We then conducted separate analyses for 5 taxonomic groups to evaluate if their spring emergence timing was affected differently depending on differences in physiology and reproductive ecology: (1) urodela, (2) anura, (3) testudines, (4) lizards (all squamata except snakes), and (5) snakes, including elevation and the same climate variables as above as fixed effects. In addition, we included body size for anura, body mass (log-transformed) for reptiles, hibernation (on land or in water; anura and testudine analyses), breeding type (explosive or prolonged; anura analysis), offspring development (direct versus larval; urodela analysis) and reproductive mode (oviparous, ovoviviparous, or viviparous; lizard and snake analyses). Finally, we included the above-mentioned species traits as interaction term with mean annual temperature to quantify if these traits affect spring emergence timing conditional on geographic variation in temperature. Body size was highly correlated with development (with larger species typically having direct development) and thus excluded from the urodela analysis. Moreover, hibernation type was not included in the urodele and squamata analyses, because there was too little variation in the data (most species hibernating on land and having a prolonged breeding season). None of the other variables were highly correlated (Pearson’s |r| < 0.6 and variance inflation factors < 3). Numeric fixed effects were scaled and centered to increase convergence and obtain comparable estimates. We used weakly informative priors for fixed and random effects (variance = 1, nu = 0.02) for the Bayesian analyses, which were carried out with 1,000,000 iterations, with burn-in of 10,000 and thinning interval of 500. For additional information, please refer to our manuscript. References Belitz, M. W., E. A. Larsen, L. Ries, and R. P. Guralnick. 2020. The accuracy of phenology estimators for use with sparsely sampled presence‐only observations. Methods in Ecology and Evolution 11:1273-1285. Fick, S. E., and R. J. Hijmans. 2017. WorldClim 2: new 1‐km spatial resolution climate surfaces for global land areas. International journal of climatology 37:4302-4315. Hadfield, J. D. 2010. MCMC methods for multi-response generalized linear mixed models: the MCMCglmm R package. Journal of statistical software 33:1-22. Karger, D. N., S. Lange, C. Hari, C. P. Reyer, O. Conrad, N. E. Zimmermann, and K. Frieler. 2023. CHELSA-W5E5: Daily 1 km meteorological forcing data for climate impact studies. Earth System Science Data 15:2445-2464.Michonneau, F., J. W. Brown, and D. J. Winter. 2016. rotl: an R package to interact with the Open Tree of Life data. Methods in Ecology and Evolution 7:1476-1481. Oliveira, B. F., V. A. São-Pedro, G. Santos-Barrera, C. Penone, and G. C. Costa. 2017. AmphiBIO, a global database for amphibian ecological traits. Scientific data 4:1-7. Oskyrko, O., C. Mi, S. Meiri, and W. Du. 2024. ReptTraits: a comprehensive dataset of ecological traits in reptiles. Scientific Data 11:243. Paradis, E., S. Blomberg, B. Bolker, J. Brown, J. Claude, H. S. Cuong, R. Desper, and G. Didier. 2019. Package ‘ape’. Analyses of phylogenetics and evolution, version 2:47. Wells, K. D. 1977. The social behaviour of anuran amphibians. Animal Behaviour 25:666-693.