Data and R script for paper 'Citizen observations reveal 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 5 × 5-degree grid cells across the study area and assigned each observation to a grid cell ID. We then estimated the mean (± SD) day of the year (Julian day) for the earliest 5% of observations per species and grid cell, defined as spring emergence day. 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 amphibians1 and the ReptTraits dataset for reptiles2. 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 periods3. We downloaded the minimum and maximum monthly temperature, and monthly precipitation for the years 2010 to 2019 (i.e., coinciding with citizen observations) at 2.5 arc minute resolution (not available at finer resolution), and the average annual temperature at 30 arc second resolution (for the years 1970 to 2000; not available after 2000) from the ‘WorldClim version 2.1’ dataset4. 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, averaged across years. Moreover, we estimated the temperature range during spring, defined as the mean of monthly maximum minus minimum temperatures. For the earliest 5% of observations per species and grid cell, we then extracted the annual mean temperature, minimum and maximum spring temperature, temperature range during spring, cumulative springtime precipitation, and elevation at 3 arc second resolution from SRTM data (https://www.earthdata.nasa.gov/sensors/srtm). We then analyzed the day of the year (response variable) for the earliest 5% of observations, separately by species and grid cell. To account for phylogenetic relatedness among species, we used Bayesian Markov Chain Monte Carlo phylogenetic generalized linear mixed models of the R package ‘MCMCglmm’5. Using raw observations, rather than the mean emergence day per species and grid cell, allowed us to also quantify within-grid cell variation of spring emergence timing in relation to elevation and climate variables. Phylogenetic data were retrieved using the R package ‘rotl’6, generating a phylogenetic tree for all species. In the resulting trees, branch lengths were estimated using the R package ‘ape’7, 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. The analyses were initially run 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. 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. Moreover, 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 the same fixed effects as above. In addition, we included body size for urodela and anura, body mass (log-transformed) for reptiles, hibernation (on land or in water; anura and testudine analyses), breeding type (explosive or prolonged; anura analysis only), and reproductive mode (oviparous, ovoviviparous, or viviparous; lizard and snake analyses). Development (direct versus larval) was highly correlated with body size 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). 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. Parameters whose 90% credible intervals overlapped zero were considered uninformative. References1 Oliveira, B. F., São-Pedro, V. A., Santos-Barrera, G., Penone, C. & Costa, G. C. AmphiBIO, a global database for amphibian ecological traits. Scientific data 4, 1-7 (2017). 2 Oskyrko, O., Mi, C., Meiri, S. & Du, W. ReptTraits: a comprehensive dataset of ecological traits in reptiles. Scientific Data 11, 243 (2024). 3 Wells, K. D. The social behaviour of anuran amphibians. Anim. Behav. 25, 666-693 (1977). 4 Fick, S. E. & Hijmans, R. J. WorldClim 2: new 1‐km spatial resolution climate surfaces for global land areas. International journal of climatology 37, 4302-4315 (2017). 5 Hadfield, J. D. MCMC methods for multi-response generalized linear mixed models: the MCMCglmm R package. Journal of statistical software 33, 1-22 (2010). 6 Michonneau, F., Brown, J. W. & Winter, D. J. rotl: an R package to interact with the Open Tree of Life data. Methods in Ecology and Evolution 7, 1476-1481 (2016). 7 Paradis, E. et al. Package ‘ape’. Analyses of phylogenetics and evolution, version 2, 47 (2019).