Temperature and land use change are associated with Rana temporaria reproductive success and phenology

Data sources

We obtained data on R. temporaria reproductive success and phenology, as well as ammonium and nitrate ion concentrations, from the Environmental Change Network (ECN) (Rennie et al., 2017). Across the UK, there are 10 ECN sites which collect this data (Fig. 1A). The ECN sites differed in the number of ponds sampled and the duration they were sampled for Table 1.

From 1^st January every year (see Table 1 for site-specific sampling coverage), ponds were sampled weekly until males first congregated at the pond. This date was recorded, and subsequent sampling frequencies increased to daily until spawn first hatched. The dates of spawning and hatching were recorded. After hatching, sampling frequencies decreased to weekly until 16 weeks after spawning or when froglets were seen leaving the ponds. The date of leaving was recorded. The areas of spawn present and the percentage of spawn found dead were recorded when each pond was sampled.

The concentrations of ammonium and nitrate ions were measured from the date of spawning by taking 250 ml of pond water to analyse in a laboratory when each pond was sampled. We obtained data on land cover change from the UK Centre for Ecology & Hydrology (UKCEH) (Rowland et al., 2020a, 2020b). Using satellite data, the dominant land cover type in each 25 × 25 m square of the UK in 1990 and 2015 were classified into six classes: woodland, arable, grassland, water, built-up areas, and other. We obtained data on climate change from 37 historic Met Office stations across the UK (Fig. 1B), which have been collecting climate data for at least 44 years (Met Office, 2022). Climate data consists of monthly mean daily maximum temperature (°C), mean daily minimum temperature (°C), days with air frost, and total rainfall (mm).

Data processing

To quantify the land cover change surrounding each ECN location and the distances between each ECN location and each Met Office station, we used QGIS (QGIS Development Team, 2022). Coordinates for the ECN locations were provided by the ECN, and the coordinates for the Met Office stations were taken from the historic station data webpage (Met Office, 2022), using the coordinate reference system OSGB 1996/British National Grid. To calculate the land cover change surrounding each ECN location we created a buffer zone with a radius of 2.5 km surrounding each ECN location, as some common frogs can disperse over 2 km (Kovar et al., 2009). We then used the land cover data obtained from UKCEH to calculate the area of each land cover class within these buffers in both 1990 and 2015. To calculate the distances between each ECN location and each Met Office station, we used the distance matrix tool.

We used the software R (R Core Team, 2020) to create a dataframe for each reproductive success and phenological response variable measured by the ECN (Table 2). Each row represented a pond at an ECN location in a single year from 1994–2015. The reproductive success data frames contained the largest surface area of spawn and percentage of dead spawn at each pond in each year. The phenology data frames contained the earliest date of either congregation, spawning, hatching, or leaving at each pond in each year, all converted to days since the 1^st January (DOY).

We incorporated the data on maximum ammonium and nitrate ion concentrations recorded in each pond in each year whilst breeding adults were present into the spawning dataset. We added data on mean ammonium and nitrate ion concentrations recorded in each pond in each year whilst spawn was present to the data frames on the surface area of spawn, the percentage of spawn dead, and the date of hatching.

To calculate land cover change, we made assumption that the rate of land cover change was constant between years, allowing us to estimate the area of each land cover class surrounding each ECN location in each year from 1994–2015. We identified the nearest Met Office station to each ECN location and used the climate data from these stations as estimates for the climate at the ECN locations. As congregation and spawning dates depend on winter climate (Carroll et al., 2009; Benard, 2015), we calculated the average of each climatic variable at each ECN location in each year over the winter months (December–February). This yielded a single mean temperature datum for that 3 month period per pond per year. Similarly we calculated the average of each climatic variable at each ECN location in each year during the months whilst spawn was present (January–May), giving a single datum for that 5 month period per pond per year. We added these values to the data frames containing information on the surface area of spawn, the percentage of spawn dead, and the date of hatching.

Data analysis

All data and code to reproduce these analyses is provided at https://github.com/xavharrison/FrogSpawn2023. First, for all variables (spawning, hatching, and percentage dead) we plotted temporal trends from 1995 to 2015 at the site level for the seven sites with sufficient time series. We used Pearson correlations to identify significant temporal trends in these traits after correcting for multiple testing. Sample sizes for the temporal analysis were: spawning (n = 175); hatching (n = 145); percentage dead spawn (n = 178). Then, we used bivariate response mixed effects models on the same dataset to test for correlations between spring temperature (spawning) or spawning temperature (hatching, percentage dead) and phenological variables. Bivariate mixed effect models can estimate the posterior correlation between temperature and phenological outcome variables whilst controlling for site effects, temporal trends and temporal autocorrelation (Houslay & Wilson, 2017; Harrison, 2021). This approach is similar to ‘detrending’ residuals to identify causal effects (e.g., see Votier et al., 2008). By measuring the correlation among the residuals of the two responses (e.g., mean maximum spawning temperature and hatch date), we can identify shifts in timing or changes to mortality over and above any long term trends in advancement of breeding. Here we predict that higher than average maximum temperatures (a positive residual) would be associated with earlier than average hatch dates (a negative residual). This manifests as a negative correlation in the residuals of the model, assesses as significant/important depending on whether the credible intervals cross zero.

Finally, we used generalised linear models (GLMs) and general linear mixed effects models (GLMM) to investigate drivers in variation in date of spawning (n = 23), and hatching (n = 25), as well as the percentage of dead spawn (n = 33) at the landscape scale. Variation in the size of the datasets is a result of ensuring all rows have complete information on all traits, including land use change and pollution. Landscape-scale variables such as land cover, temperature and chemical use are often correlated, so we used a principal components analysis (PCA) on the anthropogenic stressor variables for each of the reproductive success and phenology data frames using the factoextra and FactoMineR packages (Sebastien Le & Husson, 2008; Mundt, 2020). We extracted Principal Component 1 (PC1) and Principal Component 2 (PC2) for each reproductive success and phenological variable to use as predictors in our modelling. We used rotation plots to identify the anthropogenic stressors explained by each principal component (Wei & Simko, 2021), and so crucially the biological interpretation of the importance of each PC changes for each model. For each outcome variable we used Spearman’s rank correlation to detect significant associations with the two primary axes of variation in the PCA (PC1/PC2).

Spawn date

Across seven sites, only T01 showed a significant negative relationship with time, indicative of advancing timing of spawning (cor = −0.7, p.adj = 0.01 Fig. 2A). Here, spawning is advancing at approximately 14.9 days per decade. Sites T02 and T11 also exhibited moderate negative correlations (−0.42 and −0.3 respectively) but were not significant after correction for multiple testing. A bivariate response model controlling for temporal autocorrelation recovered a site-level effect of timing of spawning, where sites with higher winter average temperatures tend to initiate spawning earlier (posterior correlation −0.55%, 95% credible interval −0.82, −0.1; Fig. 2B; Table 3A). However we found no evidence of a residual correlation between winter temperature and spawning date at the datum level (Table 3A), suggesting limited plasticity in spawn date i.e., warmer than average winters do not consistently result in earlier than average spawning.

Univariate correlations with corrections for multiple testing revealed negative relationships between spawn date and both arable area and mean maximum winter temperature (higher values mean earlier spawning), and positive relationships with air frost days and grassland area (Figs. 2C–2F. Table S1A; higher values mean later spawning). Principal components analysis identified that the majority of these variables were collinear (Table S2; Fig. S1A), and recovered a positive correlation between PC1 and spawn date (cor = −0.45, p = 0.03).

Hatch date

Site T04 (Upper Teesdale, North Pennines) exhibited a negative trend over time in hatch date, indicative of earlier hatching every year (p = −0.4, p.adj = 0.045; Fig. 3A). Sites T01, T02 and T011 also showed negative trends (correlations −0.23 to −0.44), though these were not significant after correction for multiple testing. A bivariate response model controlling for temporal autocorrelation recovered a negative correlation between hatch date and mean maximum spawning temperature at the site level (i.e., sites with higher spawning temperatures tend to hatch earlier, though the credible intervals crossed zero (posterior correlation −0.46%, 95% credible interval −0.83, 0.03; Fig. 2B; Table 3B). We did find evidence of a residual correlation between spawning temperature and hatch date at the datum level (Table 3B), suggesting earlier hatching in warmer years consistent with faster development of eggs.

Univariate correlations with corrections for multiple testing revealed negative relationships between hatch date and both arable area and mean maximum spawning temperature (higher values mean earlier hatching), a positive, potentially non-linear relationship with air frost days, and a positive relationship with grassland area (Figs. 3C–3F. Table S1A; higher values mean later hatching). Principal components analysis identified that the majority of these variables were collinear (Table S3; Fig. S1B). Increaases in PC1 (representing fewer air frost days, higher temperatures and greater amounts of arable area) resulted in earlier hatching (cor = −0.7, p < 0.001).

Reproductive success

Though several sites showed increases over time in the percentage of dead spawn, only site T11 (Yr Wyddfa, Wales) had a significant trend when correcting for multiple testing (cor = 0.54, p.adj = 0.02; Fig. 4A). Bivariate response modelling identified a negative posterior correlation between temperature during spawning and spawn death (cor = −0.8%, 95% credible interval −0.96, −0.45; Table 3C) suggesting that sites with lower maximum spawning temperatures tend to show greater degrees of spawn mortality, on average. We did not detect correlation between the residuals for spawn death and spawn temperature (posterior mean correlation 0.04%, 95% credible intervals −0.12, 0.21), meaning that within sites there was no evidence that years with colder than average spawning temperatures suffered higher than average mortality. Univariate analyses revealed a negative relationship between arable area and percentage of dead spawn (more arable area reflects lower mortality), and a positive relationship with grassland area (Figs. 4C–4E. Table S1C). Principal components analysis identified strong collinearity among these traits (Table S4; Fig. S1C). Increases in PC1 (representing fewer air frost days, higher mean temperatures and more arable area. Table S4) resulted in lower proportions of dead spawn (cor = 0.49, p = 0.007).

Phenology, land use change and temperature

This work is consistent with previous studies linking temperature shifts to alteration in amphibian breeding phenology (e.g., Carroll et al., 2009; Phillimore et al., 2010; Blaustein et al., 2010; Benard, 2015). Consistent with our predictions, at the landscape scale we found that variation in the onset of spawning and hatching was linked to warmer average winter and spawning temperatures, respectively. Similarly, within populations we found that the magnitude of the slope of phenological advancement was population specific, as did a previous study on R temporaria (Phillimore et al., 2010). Spatial heterogeneity in the slope of phenological trends may be common (e.g., Primack et al., 2009), and not solely associated with latitudinal variation. For example, Sparks et al. (2007) found considerable variation in the slope of shifts in spawning initiation dates with mean Jan-March temperature between populations of R temporaria in the UK and Poland, even though they were found at similar latitudes. Both studies (Sparks et al., 2007; Phillimore et al., 2010) attribute the observed among-population variation in slope to local adaptation to climate. Spatial variation in the strength of phenological (Primack et al., 2009; Phillimore et al., 2010; this study) or phenotypic (Sheridan et al., 2018) adjustments to climate change highlight that ‘space-for-time’ studies of species traits under climate change should be used with caution (Phillimore et al., 2010; Sheridan et al., 2018).

We still lack a comprehensive understanding of the costs of a population not ‘keeping up’ with rates of phenological advancement, nor how this may vary in concert with the mean and variance of local climatic conditions. For example, warmer winters have been shown to advance phenology in congeners like wood frogs (Rana sylvatica), but also delay larval development time when followed by cooler post-breeding temperatures (Benard, 2015). Thus, even if environmental trait means shift in a consistent direction (winters get warmer, earlier), increased variation around that mean may erode selection on changes in phenological traits. In common toads Bufo bufo in the United Kingdom, warmer than average temperatures correlated with decreased body condition and survival (Reading, 2007). Recent work on wood frogs showed that phenological shifts can expose individuals to colder temperatures and resulted in lower tolerance of offspring to pollutants like NaCl (Buss, Swierk & Hua, 2021). Collectively these data suggest that in some species, warming can induce cryptic cost of breeding plasticity in multiple life stages that are not immediately apparent if looking at phenological variables alone (see Blaustein et al., 2010).

Using bivariate models revealed that within sites, warmer than average winter temperatures for a particular year were not associated with earlier than average spawning, indicating limited short-term (rather than inter-annual) plasticity in the timing of spawning. However, we did uncover such a pattern for hatching, consistent with more rapid development times under warmer spawning temperatures. Previous work on natterjack toads has shown that devleopment time is reduced under warming temperatures (Sanuy, Oromí & Galofré, 2008), similar to the plasiticty we have observed here. The lack of signal of short-term plasticity in spawning could arise if mean maximum winter temperature is too coarse a measurement to capture the climatic cues used within years to determine exact spawning date. Conversely, such plasticity might be the property of individuals and not populations. For example, previous work on two terrestrial newt species found marked among-individual variation in exploratory behaviour linked to temperature shifts (Hubáček & Gvoždík, 2023). This suggests individual frogs may be capable of adjusting their reproductive behaviour in response to fine-scale among-year variation in climate (plasticity), but against a background of strong among-individual heterogeneity, making the pattern difficult to detect at the population level. Addressing this possibility requires that we track individual frogs across years, for example by using unique dorsal markings (Petrovan et al., 2024), to measure individual behavioural adjustments over time.

We also uncovered associations between land use change and breeding phenology. Increased arable and built-up land cover is associated with earlier hatching, whilst grassland is associated with later hatching. Arable and built-up areas tend to be warmer than the surrounding natural habitats due to unvegetated ground, whilst grassland areas, with higher vegetation levels, are cooler (Lembrechts, Nijs & Lenoir, 2019; Schmidt, Lischeid & Nendel, 2019). The warming effects of arable and built-up areas can even increase the temperature of the wider landscape, leading to earlier breeding phenology in the areas surrounding human-dominated land (Tian et al., 2020). These effects could be reduced by increasing vegetated areas in these land cover types, thus providing cooler microclimates (Greenwood et al., 2016). Previous work has shown that anthropogenic habitat alteration can alter breeding phenology in multiple Australian amphibian species, specifically that increased urbanisation is linked to earlier breeding (Liu et al., 2022). However, the nature of our data prevents us from disentangling the relative influence of land use change and temperature regime shifts. Future work could investigate microclimatic variation at breeding sites to explore how urbanisation changes climatic envelope experienced by amphibians, and compare these trends to sites with similar mean temperatures in (more southerly) rural areas.

Spawn mortality, land use change and temperature

We found no clear association between mean winter maximum or minimum temperatures and spawn mortality, when temperature was used as the sole predictor in models. Instead, using Principal Component Analysis we found that the percentage of dead spawn was linked to a composite measure of frost days, rainfall, land use change, and temperature. Increases in the proportion of arable and urban areas, nitrate levels, and increased mean maximum temperature, were associated with higher spawn survival. Conversely, increased numbers of frost days, rainfall and grassland areas were associated with decreasing reproductive success. A key finding was that some populations are exhibiting consistent increases in the proportion of dead spawn over time, which may indicate that long term viability of these locations may be compromised even if no further land use or environmental changes occur.

The correlations between land use, pollution and spawn survival did not align with our predictions. We expected increased farmland (and associated pollution such as nitrates), and increased urbanisation to be detrimental rather than associated with higher spawn survival. However, though human-driven land use change is often associated with lower reproductive success, some studies have shown that they can support declining populations. For example, the average reproductive success of multiple amphibian species in America breeding in arable ponds was no different to in natural wetlands (Knutson et al., 2004). Though some species did respond negatively to arable land use (Knutson et al., 2004), R. temporaria have been found to use arable ponds more than other amphibians (Hartel, Băncilă & Cogălniceanu, 2011). Therefore, arable land cover may be able to support robust populations of R. temporaria, leading to higher spawn survival in the ECN sampling locations. Arable ditches may confer landscape connectivity for amphibians, reptiles, and mammals; Maes, Musters & De Snoo (2008) found that ditches in agricultural environmental schemes supported similar Rana esculenta abundances as nature reserves, demonstrating how these features of arable land cover can serve as important avenues for dispersal between breeding sites. This is crucial in preventing population isolation, perhaps explaining why we found a negative relationship between the arable area and the maximum percentage of spawn that died (Allentoft & O’Brien, 2010).

Built-up areas have also been shown to increase the fitness of some wildlife populations, with some threatened birds (Kettel et al., 2019) and amphibians (Iglesias-Carrasco, Martín & Cabido, 2017) having greater reproductive success ad body condition respectively in built-up areas than in the countryside. There are multiple theories for why this could be Saenz, Hall & Kwiatkowski (2015) suggest that the occurrence of chytridiomycosis, an amphibian disease that can reduce fitness, could be lower in urban areas. However, chytridiomycosis is not common in British amphibians (Garner et al., 2005) so this is unlikely to explain our findings. Hall & Warner (2017) suggest that the high densities of prey insects in urban areas or lower predation pressures that allow adults to spend more time hunting, could increase fitness. Further evidence for this is from Germany, where R. temporaria adults were found to be bigger in urban greenspace than in the surrounding countryside (Niemeier et al., 2020). Larger adult body size is likely to increase offspring survival (Hall & Warner, 2017), providing an explanation as to why built-up areas are associated with lower R. temporaria spawn mortality.

Many of the studies that demonstrate the value of human-dominated land cover types also acknowledge the need for management. In built-up areas, it is important to reduce barriers to movement by connecting urban greenspace (Mazgajska & Mazgajski, 2020; Niemeier et al., 2020). These suggestions highlight that local management for R. temporaria could be more important that the broad-scale land cover type in determining this amphibian’s reproductive success. Landscape management practices could also explain why grassland is negatively associated with R. temporaria reproductive success (i.e., positively associated with mortality). In the UK, only 2% of grassland is classed as diverse (Bullock et al., 2011), due to widespread “improvement” (Vickery et al., 1999) and livestock grazing (Fuller, 1987; Bullock et al., 2011). However, grazing can be detrimental to amphibians. Livestock can cause high levels of wetland bank erosion (Trimble, 1994), leading to increased sediment deposition. When investigating the impact of cattle on Bufo achalensis, a toad species endemic to Argentina, Jofré, Reading & di Tada (2007) found that increased sediment levels reduced algal growth, a key food source for larval amphibians, leading to higher mortality in the B. achalensis larvae. This could result in increased isolation between R. temporaria populations, potentially leading to the increased percentage of spawn death observed in ECN locations surrounded by grassland. High densities of livestock can also directly reduce amphibian spawn reproductive success through disruption and trampling (Knutson et al., 2004).

Our modelling could not conclusively implicate temperature regime changes as a potential driver of difference in amphibian reproductive success, though maximum temperature and number of frost days were part of the composite measure (PC1) associated with spawn survival. Although climate change is leading to warmer temperatures, it is also causing extreme weather conditions (Huber & Gulledge, 2011). One example is the occurrence of spring cold-snaps, a phenomenon that has been shown to have detrimental effects on a wide range of organisms (Augspurger, 2013; Benard, 2015; Turner & Maclean, 2022). Benard (2015) observed an increase in Rana sylvatica larvae being exposed to cold-snaps from 2006 to 2012, leading to altered development. Freezing is known to kill R. temporaria (Pasanen & Karhapää, 1997), making it likely that cold-snaps could lead to high spawn mortality. Warmer springs and fewer frost days likely explain the lower proportions of dead spawn observed in this dataset under these conditions.