Recruitment Patterns and Potential Climate Change Impacts on Three Florida Hylids with Different Life Histories

Abstract

Altered weather patterns associated with climate change are likely to adversely affect amphibian recruitment, especially for species dependent on ephemeral, geographically isolated wetlands for breeding. Future changes in temperature and rainfall patterns could affect hydroregimes (periodicity, depth, duration, and timing of water in wetlands) or adult breeding effort. We used 24 years of continuous amphibian trapping, weather, and hydroregime data to identify breeding-to-metamorphosis periods (BMPs) and environmental factors affecting annual recruitment by three hylid species at eight isolated ephemeral limesink ponds in Florida longleaf-wiregrass sandhills. We used standardized climate metrics (Bioclim variables) to predict future precipitation, temperature and hydroregime variables, then used them to predict future recruitment in 2050 and 2070 under two emissions scenarios. We hypothesized that Hyla gratiosa would be more sensitive to short-term pond drying than H. femoralis or H. squirella due to its lower abundance and more specific habitat requirements. Hyla gratiosa recruitment was not explained by adult breeding effort and was more dependent on higher water levels during BMPs than for H. femoralis or H. squirella, independent of rainfall. In contrast, H. femoralis and H. squirella recruitment depended heavily on rainfall independent of pond depth and was positively associated with adult breeding effort. Models predicted moderate decreases in H. gratiosa and H. squirella recruitment by 2050 but projections were highly uncertain for all three species by 2070. Our findings highlight the importance of maintaining wetlands with diverse hydroregimes to accommodate species with different BMPs and hydroregime requirements. Proactive monitoring and conservation measures such as headstarting and creating artificial ponds may be necessary for these and other amphibian species that may suffer reduced recruitment under future climate change.

1. Introduction

Climate change is among the most widespread causes of modern biodiversity loss and poses a severe risk to the persistence of wetland biodiversity [1]. Geographically isolated ephemeral wetlands represent one biodiverse habitat that may be impacted by climate change [2]. These wetlands are known for high amphibian diversity, partly due to seasonal cycles of drying and filling that exclude predatory fish and promote successful amphibian recruitment [3]. However, amphibian species that require ephemeral wetlands are vulnerable to climate change due to their dependence on reliable seasonal rainfall and pond filling for successful breeding and recruitment [4]. Excluding data-deficient species, 42% (n = 1692/3992) of wetland-associated amphibians are currently listed as “Near Threatened” or worse by the International Union for Conservation of Nature (IUCN Red List, 2021). At least 29% of these species (n = 496) face recognized threats due to direct or indirect effects of climate change. Despite this, precise climate-related drivers of declining amphibian recruitment remain poorly understood for many species [4].

The adaptability of wetland-breeding amphibians to climate change depends on complex interactions between local hydroregime variables (i.e., the periodicity, depth, timing, and duration of water), life history, and amount of temperature and precipitation change [5,6,7,8]. For example, in geographically isolated, groundwater-driven ephemeral wetlands (e.g., isolated sinkhole ponds), aspects of hydroregime are critical determinants of suitable breeding habitats. In these habitats, hydroregimes are driven by short-term precipitation and longer-term groundwater discharge [2]. Species with slow larval development rates to metamorphosis, such as most Lithobates species, require long periods of continuous water in ponds to prevent tadpole desiccation before metamorphosis [9]. In contrast, species with rapid rates of larval development, such as many Hyla species, require minimum hydroperiods of shorter duration throughout the breeding, egg, and larval periods (“breeding-to-metamorphosis period” (BMP)) [10,11]. Long-term studies evaluating linkages between temperature, rainfall, hydroregime, and recruitment by ephemeral pond-dependent species with different life histories are needed to predict their responses to climate change and develop resiliency-oriented management plans [12,13].

Geographically isolated wetlands are often ephemeral and used by a wide diversity of amphibians, making them an excellent system for studying long-term recruitment trends related to climate change [14,15]. Because they often are poor dispersers [16,17] juvenile amphibians can serve as accurate proxies for studying on-site recruitment, less confounded by immigration than in more vagile taxa [18]. For example, hylid treefrogs tend to disperse over small distances (<1 km) and are highly affected by habitat-specific dispersal barriers [17,19]. Moreover, hydroregimes can vary considerably even among proximate ponds due to differences in surface topography and, perhaps more importantly, subsurface drainage [20]. This allows for disentanglement of the role of weather or climate from pond-specific hydroperiod variables in shaping recruitment. Some evidence suggests that certain wetlands, such as those at low elevations and closer to the water table, may serve as more viable climate refugia than others (with exceptions) due to a greater chance of holding water for longer periods [21,22]. However, few studies have focused on how amphibians using isolated ponds may respond to altered hydroregimes and a changing climate (but see Greenberg et al., 2017 [9]).

In this paper, we used 24 years (March 1994–December 2017) of continuous, concurrent amphibian trapping data from eight ephemeral limesink (sinkhole) ponds to examine the influence of weather, hydroperiod variables, and apparent adult breeding effort on recruitment by Hyla femoralis, H. gratiosa, and H. squirella, which vary modestly in their reproductive and larval life histories. Hyla femoralis and H. squirella have less specific breeding habitat preferences and generally higher abundances than H. gratiosa across all life stages [23]. In addition, larval periods are shorter for H. femoralis (50–75 days) and H. squirella (40–50 days) than H. gratiosa (41–160 days) (

Table 1. Size cutoffs (SVLs), breeding-to-metamorphosis periods (BMPs), larval duration [10,24], and sample sizes (first captures) for each study species, Ocala National Forest, Marion County, Florida.

					Metamorphs		Juveniles		Adults
Species	BMP	Breeding Period	Metam. Period *	Larval Duration (Days)	Size (mm)	n	Size (mm)	n	Size (mm)	n	n
H. femoralis	1 May–26 October	1 May–31 August	22 June–26 October	50–75	<13	131	<25	958	≥25	2271	3229
H. gratiosa	15 April–9 October	15 April–28 August	28 May–9 October	41–160	<25	44	<44	285	≥44	60	345
H. squirella	1 April–16 October	1 April–31 August	10 May–16 October	40–50	<15	100	<23	190	≥23	466	656
Total						275		1433		2797	4230

* Metamorphosis periods capture ~90% of observed metamorphs and assume that eggs hatch within 48 h of being deposited [24].

) [10,24]. These species vary somewhat in start and end dates for breeding and metamorphosis, but the BMP broadly spans mid-spring to mid-fall for all three (Table 1). Interspecific differences in some traits (abundance, habitat preferences, larval duration) but similarities in others (spring/summer breeding seasons, shared genus, use of semi-permanent ponds for reproduction) make these hylids valuable models for identifying roles that their non-shared traits play in juvenile recruitment under a changing climate.

We related observed hydroperiod and weekly weather data to climate variables (calculated from daily weather station data aggregated over annual periods) and used these relationships to predict future hylid recruitment under two different climate change scenarios (Figure 1). Our three study species depend heavily on isolated ponds embedded within Southeastern Coastal Plain longleaf pine-wiregrass savannahs [24] and face uncertain futures under climate change. For each species, we examined: (1) if and how recruitment is influenced by adult breeding effort, weather, and hydroregime variables, both immediately before and during their respective BMPs; (2) if future recruitment-related weather and hydroregime variables can be predicted using annual temperature and precipitation-related (Bioclim) climate variables; (3) how future hylid recruitment might be altered under climate change given two different carbon emissions scenarios. We hypothesized that the sensitivity of hylid recruitment to short-term wetland drying would be higher for H. gratiosa than H. femoralis and H. squirella, due to its relatively slower larval development rates.

2. Materials and Methods

2.1. Study Area

We monitored eight small (0.1–0.37 ha) geographically isolated ephemeral limesink wetlands (“ponds” hereafter) fed by groundwater from the Floridan Aquifer System, in Ocala National Forest, Marion County, Florida. Study ponds were surrounded by savanna-like sandhills with wiregrass-forb ground cover and widely-spaced longleaf pine trees with patches of hardwood trees and sand pine (Pinus clausa). Six of our eight study ponds were clustered just north of Salt Springs, FL and located within ~1 km of each other; the remaining two were located ~9.5 km farther south (Figure 2). Average weekly temperatures (1 February 1994–31 December 2017) ranged from 13.5 °C in January to 28.6 °C in August. Average annual precipitation (1 January 1994–31 December 2017) in our study area was 140.7 cm, with more than half occurring during late spring and summer. Tropical systems and thunderstorms occurring in summer and fall, and wet autumn, winter, or spring frontal systems provided most groundwater recharge [26]. Ponds were generally deepest in winter and shallowest in summer [2], due to combined seasonal variation in rainfall, groundwater recharge, and evapotranspiration rates (Knowles Jr. et al., 2002). Common soils were well-to-excessively drained Entisols with <5% silt plus clay in the upper profile, classified in the hyperthermic, uncoated families of Spodic (Paola series) and Typic (Astatula series) Quartzipsamments (Aydelott et al., 1975). Study pond elevations ranged from 4–26 m.

2.2. Field Methods

We installed 7.6 m sheet metal drift fences around the perimeter of each pond (n = 9–16 per pond) at the approximate high-water line. Fences were spaced 7.6 m apart and thus enclosed half the perimeter of each pond. Pitfall traps (19-L buckets) were positioned inside and outside of each end of each fence (four per fence), and a double- or single-ended funnel trap (one of each per fence) was positioned at the midpoint of the fence on both sides, to detect directional movement by amphibians to and from wetlands. One vertically-oriented PVC pipe (5 cm width, 1.4 m height) was placed between each fence to attract hylids. We placed sponges in all pitfall and funnel traps and moistened them as needed during surveys to limit amphibian desiccation (see Greenberg et al., 2017 [28] for greater detail). The location of our fences at the upland-wetland interface meant that they captured hylids year-round, though most frequently during hylid BMPs between spring and fall.

We checked drift fences three times weekly in the period 1 February 1994–31 December 2017, or sometimes less frequently (1–2 times weekly) during cold months (November–early March) when capture rates were low. We identified species, measured using calipers (1 mm accuracy), and toe-clipped (to indicate pond and year of capture) upon the first capture, then released individuals on the opposite side of the fence (if captured in a pitfall or funnel trap) or placed them near or in the PVC pipe where captured. We measured pond depths, precipitation, and maximum and minimum air temperatures within our study area weekly. Missing temperature data during periods of equipment malfunction (<10% of the dataset) were imputed using the MissForest package in R [29]. Missing precipitation data were estimated using the pond hydrology model developed by (Greenberg et al., 2015 [2]).

2.3. Determining Breeding-to-Metamorphosis Periods (BMPs)

We focused on patterns and potential drivers of recruitment by three hylid species: Hyla femoralis, H. gratiosa, and H. squirella. These species all breed in spring–summer, but vary somewhat in reproductive phenology and larval periods (Table 1), as well as abundance, habitat specificity, and non-breeding habitat use [10,23,24,30,31]. All three species reside in uplands surrounding wetlands during most of their lives but depend on standing water for breeding and larval development to metamorphosis during their BMPs. We used snout-vent length (SVL) cutoffs based on the literature and professional judgement to define adults, juveniles, and metamorphs, respectively, for each species: H. femoralis (≥25 mm, <25 mm, <13 mm), H. gratiosa (≥44 mm, <44 mm, <25 mm); H. squirella (≥23 mm, <23 mm, <15 mm). Metamorphs were used to estimate the end dates of BMPs based on their capture dates and were treated as a subcategory of juvenile; thus the SVL ranges of the juvenile and metamorph categories overlapped (Table 1). Differences in relative weekly capture rates (peaks) of adults and metamorphs in our study (Figure 3) and the literature [24,30,31] indicate that BMPs span approximately 1 May–26 October for H. femoralis, 15 April–9 October for H. gratiosa, and 1 April–16 October for H. squirella.

2.4. Data Processing

We identified BMPs for each species based on spikes in capture during spring through fall. We used captures of individual juveniles (first captures only) during annual BMPs as an index for estimating annual recruitment at each pond for each species. We assumed that juveniles captured during BMPs through the end of the year represented same-year recruits, due to relatively rapid maturation after metamorphosis in our study species (Conant and Collins, 1998).

Ponds occasionally flooded, causing temporary trap closures. We estimated weekly lost captures due to flooding per pond by combining information about weekly and annual capture rates. For example, suppose that a given pond was “open” (no trap flooding) during week 1 of April during 20 of 24 possible years and there were 1000 total “open” weeks at the pond (all months combined) over the entire study period. Thus, the first week of April accounted for 2% (20/1000) of total trapping effort across all years. If there were 900 total first captures during “open” weeks over the study period with 5% (45/900) occurring during week 1 April, then individuals were captured at a 2.5× faster rate (5% of captures/2% of trapping effort) for that week compared to overall. Thus, if traps were flooded during week 1 April, 2012 and there were 2 first captures per “open” week in 2012 overall (e.g., 80 captures over 40 “open” weeks), we estimated 5 lost captures (2 × 2.5) during week 1 April 2012, given that week’s 2.5× higher-than-average capture rate across all years and an overall 2012 capture rate of 2 captures/open week.

After accounting for weekly flooding, we summed estimates of first captures across all weeks within each year, then accounted for small changes in the number of PVC pipes in place over time to better estimate annual recruitment. To adjust for unequal PVC trap nights over the course of this study (e.g., lost due to prescribed burns or theft until replaced), we scaled annual hylid captures (per pond and year, after adjusting for flooding) to the number of PVC pipes in place during the majority of each given year, as a percentage of the pond’s maximum PVC pipe count. For example, if a pond contained a maximum of 16 PVC pipes at any point during this study, but only 11 were in place for the majority of 2012, then annual hylid captures at that pond were multiplied by 1.45 (16/11) for 2012.

We used weekly temperature and precipitation data collected within our study area, adult first capture data (number of individuals), and water depth data for each study pond to compute predictor variables for annual hylid recruitment. We used ten predictor variables to model annual hylid recruitment: pond ID, a distance-based spatial eigenvector to account for spatially autocorrelated recruitment [32], minimum pond depth in the month prior to the BMP, and seven variables measured each year during BMPs (apparent adult breeding effort, mean temperature, total precipitation, minimum pond depth, median pond depth, maximum pond depth, and consecutive days with water in pond). We treated first captures of adults each year during the defined breeding period (Table 1) as a proxy for apparent adult breeding effort. Adults captured outside of the breeding period were considered non-breeding and were not counted towards adult breeding effort. Given the tendency for individuals to regenerate toe-clipped digits over several months (pers. obs.), recaptured adults were likely counted as “first captures” during each new breeding season, but not multiple times within a breeding season.

2.5. Data Analysis

We evaluated potential climate change impacts on hylids in three phases (Phases I, II, and III hereafter). In Phase I, we modeled annual hylid recruitment as a function of adult breeding effort and weekly field-measured weather and hydroregime data. In Phase II, we used finer-scale (daily) climate data from an Ocala, FL weather station to predict field-measured variables to inform recruitment in Phase I. In Phase III, we combined Phase I and Phase II models with climate change projections to predict future recruitment for each study species. We define “climate” hereinafter as the aggregation of weather data over an annual timescale, and “weather” as temperature and precipitation variables over any sub-annual timescale (e.g., daily, weekly, or monthly).

2.5.1. Phase I: Evaluating Patterns of Annual Recruitment

We used boosted regression trees (BRTs) [33] to evaluate the influence of pond conditions and adult breeding effort on annual hylid recruitment. These models use iterative decision trees to identify relationships between predictor variables and a response variable. Each subsequent decision tree works to explain residual deviance remaining from the previous one [33]. To allow for a straightforward interpretation of our results and to maximize model consistency, we fit all BRTs with tree complexity = 2, learning rate = 0.0001, and bag fraction = 0.65–0.75. We assessed model fit using k-fold cross-validation with 10 folds [34]. Because our data were highly zero-inflated (i.e., no recruits captured at some ponds during many years), we constructed BRTs in two steps for each species. First, we constructed BRTs with a binomial error distribution to predict the likelihood of any recruitment occurring at a given pond during a given year (i.e., at least one juvenile captured). Next, we constructed BRTs with a Poisson error distribution to predict the number of juveniles captured in years when recruitment occurred. We then multiplied predictions from these two models together to calculate the total relative recruitment (i.e., scaled juvenile captures) expected for a given pond and year.

2.5.2. Phase II: Evaluating the Relationship between Bioclim and Recruitment-Related Variables

Daily temperature and precipitation data (as opposed to weekly, as measured manually within our study area) were needed to calculate Bioclim variables, which describe the climate of a given location based on several characteristics of its temperature and precipitation regimes. We used daily temperature and precipitation data (1994–2017) from the closest weather station to our study area, located in Ocala, FL [35] to calculate annual values for each of 19 Bioclim variables [25] in the period 1994–2017. After calculating Bioclim variables, we used least absolute shrinkage and selection operator (LASSO) models [36] to assess the relationship between these variables and recruitment-related predictor variables from Phase I. Thus, predictor variables from Phase I were response variables in Phase II. Because we treated year as our unit of replication in Phase II (total n = 24 years), we averaged annual Phase I predictor variables across all eight ponds where necessary (e.g., for maximum pond depth) for Phase II. We did not discard any Bioclim variables before conducting these analyses because (1) our LASSO models were robust to multicollinearity, as they consistently shrunk slope (beta) values of all but 1–3 relatively uncorrelated predictor (Bioclim) variables to zero, (2) we lacked strong a priori hypotheses about specific relationships between Bioclim variables and field-collected variables, and (3) we were more interested in maximizing the accuracy and consistency of Bioclim-derived estimates of Phase I variables than in disentangling precise linkages between the two. We identified an optimal shrinkage penalty parameter (λ) for each LASSO model by comparing the mean-squared error of LASSO models built with different possible λ values using leave-one-out cross-validation [37]. We did not attempt to predict apparent adult breeding effort using Bioclim variables, due to its complicated relationship with recruitment during previous years in addition to unmodeled abiotic factors. Due to their immutability over time, we also excluded pond ID and spatial eigenvectors from all Phase II analyses.

2.5.3. Phase III: Projecting Future Recruitment under Different Climate Change Scenarios

We used a CMIP5 general circulation model [38] under two realistic emissions scenarios (RCP 4.5 and RCP 6.0; Meinshausen et al., 2011 [39]), along with our Phase II models, to predict values for variables related to mean annual hylid recruitment (i.e., predictor variables from Phase I) for 2050 and 2070. Predicted variables were then used as inputs in our Phase I models to predict future annual mean recruitment for each species. Predictor variables from Phase I that could not be explained using Bioclim variables in Phase II were held constant at their mean values in all Phase III simulations. To assess model consistency and generate distributions of predictions for mean future recruitment, we replicated the above approach 300 times for each species using bootstrapping [34]. We also used 10,000 bootstraps to construct distributions for mean annual recruitment in the period 1994–2017, given our observed data for each study species.

2.5.4. Decadal Capture Trends

To enable visual interpretation of long-term capture trends, we calculated 10 year running averages of annual new captures for each study species in the period 2003–2017. For each year and species, we calculated an average of the number of first captures (all age classes combined) across the previous 10 years. For example, the 10 year running average for new captures of H. femoralis in 2015 was calculated by summing all new captures from 2006–2015 and then dividing by 10.

3. Results

We captured 4230 total individual hylids, including 3229 (10,647 recaptures) H. femoralis, 345 (31 recaptures) H. gratiosa, and 656 (3186 recaptures) H. squirella. Hylid species varied widely in recruitment from year to year, among ponds, and in the proportions of adult versus juvenile captures (Figure 4). Overall, juveniles comprised 29.8% of first captures for H. femoralis, 82.6% for H. gratiosa, and 30.1% for H. squirella (Table 1). Most H. femoralis (65%) and H. squirella (79%), but few H. gratiosa (5%) individuals, were captured in PVC pipes (as opposed to pitfall or funnel traps). Among PVC pipe captures, recaptures (excluded from our analyses) comprised the majority of total captures (83% for H. femoralis, 61% for H. gratiosa, and 86% for H. squirella). Adult and juvenile captures peaked in similar weeks for each species, but juveniles lagged behind adults slightly. A smaller, secondary peak in juvenile captures was observed months later in the year for H. femoralis and H. squirella (Figure 4).

3.1. Phase I: Influence of Pond-Level Conditions on Annual Recruitment

Recruitment was highly dependent on environmental or pond conditions during the BMP for all three species, but the relative importance of each varied among species. Recruitment was highest during rainy years for all three species, but heavy rains during the BMP benefitted H. femoralis and H. squirella far more than H. gratiosa (Figure 5, Figure 6 and Figure 7). Consistently-filled ponds during the month before breeding were associated with increased recruitment for H. femoralis and H. gratiosa but did not explain H. squirella recruitment. Recruitment increased concomitant with apparent adult breeding effort for H. femoralis and H. squirella, but not H. gratiosa. Conversely, minimum and maximum pond depth during BMPs were poor predictors of recruitment for H. femoralis and H. squirella, but were positively correlated with H. gratiosa recruitment, with its maximum recruitment occurring above thresholds of 70 cm (for minimum pond depth) and 150 cm (for maximum pond depth).

3.2. Phase II: Predicting Pond-Level Variables Using Bioclimatic Data

Three of seven environmental predictor variables in Phase I (total precipitation, maximum pond depth, and consecutive days that ponds held water during the BMP) were estimable using Bioclim variables for all study species (

Table 2. Correlations (r) between observed annual values of each of our Phase I predictor variables (1994–2017) and values predicted by Bioclim variables in Phase II using LASSO models, Ocala National Forest, Marion County, Florida.

Field-Derived Variable *	Correlation between Past Annual Values (1994–2017) and Predictions from Bioclim Variables in LASSO Model **
	H. femoralis	H. gratiosa	H. squirella
Average weekly mean temperature	—	—	—
Total precipitation	0.95	0.96	0.95
Median pond depth	0.45	—	—
Minimum pond depth	—	—	—
Maximum pond depth	0.62	0.68	0.64
Minimum pond depth in month prior to breeding season	—	0.56	0.57
Maximum consecutive days with water in pond	0.72	0.69	0.69

* Field-derived variables were calculated during BMP periods unless stated otherwise. ** Dashes indicate Phase I variables that could not be reliably predicted using Bioclim variables without resulting in model overfitting.

). In addition, Bioclim variables were able to predict median pond depth during the BMP for H. femoralis (r = 0.45), and minimum pond depth in the month prior to breeding for H. gratiosa (r = 0.56) and H. squirella (r = 0.57). The remaining two Phase I environmental variables (mean temperature and minimum pond depth) were not accurately predicted by Bioclim variables (i.e., all beta values were shrunk to zero in the LASSO model), and were thus held at their mean values when simulating the potential influence of climate change on hylid recruitment.

3.3. Phase III: Predicting Changes in Recruitment under Different Climate Change Scenarios

Climate change-related effects on recruitment appeared possible for all three hylid species, but with varying effect sizes and levels of confidence (Figure 8). Moreover, potential impacts on recruitment varied only weakly between RCP 4.5 and RCP 6.0 emissions scenarios. Mean annual recruitment appeared more likely to be negatively impacted by climate change for H. gratiosa and H. squirella than for H. femoralis under both emissions scenarios, but the distribution of possible predicted mean recruitment values in 2050 was shifted at least slightly downward for all three species compared to in the period 1994–2017 (except for H. femoralis under the RCP 4.5 emissions scenario). The distributions of predicted mean annual recruitment in 2070 were broader than 2050 for all three species, rendering the 2070 projections relatively uninformative.

3.4. Decadal Capture Trends

The three study species varies widely in long-term capture rates over the course of our study (Figure 9). For H. femoralis, the (likely) most abundant species, 10 year capture rates vacillated over the course of this study, but were roughly the same at its start in 1994 and conclusion in 2017. In contrast, H. gratiosa and H. squirrela had much lower capture rates throughout this study. However, these two species had opposite trajectories; decadal H. gratiosa captures declined around the mid-to-late 2000s and did not recover, while decadal H. squirella captures increased around the early 2010s (surpsassing H. gratiosa captures) and did not level off by the end of this study.

4. Discussion

We investigated potential climate change impacts on hylid recruitment by first assessing the influence of field-measured weather and hydroregime variables on recruitment (Phase I), then assessing the relationship between these field-measured variables and annual climate variables (Phase II), and finally synthesizing the two above models with future climate projections (Phase III). Precipitation predicted hylid recruitment most accurately when considered alongside pond-specific hydroperiod variables (minimum, median, and maximum pond depth, and consecutive days that ponds held water). For Phase II, climate variables explained >62% of the variation in maximum pond depth and >69% of the variation in consecutive days that ponds held water, but poorly predicted minimum and median pond depths (Table 2), suggesting that local factors affecting hydroperiod such as groundwater discharge (Boswell and Olyphant, 2007) and subsurface drainage (Greenberg et al., 2015)—thus hylid recruitment—differ among ponds. Our Phase III results are likely conservative estimates of how climate change may impact recruitment, as we assumed that hydroregime variables related to long-term groundwater storage (unmeasured) would remain constant under climate change scenarios despite changes in annual climate variables (Table 2). Nonetheless, our results suggest, at a minimum, potential slight-to-moderate future decreases in overall hylid recruitment by 2050 under both climate scenarios (Figure 8). While our models were designed to maximize predictive capabilities rather than causal inference about individual predictor variables, predictions of reduced recruitment by 2050 may reflect increased rates of evaporation and (or) evapotranspiration rates (thus less water in ponds) as temperatures continue to increase [9].

Adult capture rates (Table 1) indicated that annual breeding effort was relatively low for H. gratiosa compared to H. femoralis and H. squirella. Additionally, the relationship between adult breeding effort and juvenile recruitment was weaker for H. gratiosa than H. femoralis or H. squirella. Given differences in recruitment patterns among our study species, non-shared traits among these species (e.g., larval duration, abundance, and breeding condition preferences) were likely equally or more important for shaping their recruitment than shared traits (e.g., breeding season and genus). Our study species and other pond-breeding anurans within our study area also show high variability between adult breeding effort and successful juvenile recruitment, likely due in part to unseen (underwater) factors affecting egg or larval survival, such as predation or disease [9].

Potential differences in detectability among study species could have biased interspecific comparisons of adult breeding effort. Most adult and juvenile H. femoralis (65%) and H. squirella (79%) individuals were captured when voluntarily using PVC pipes adjacent to ponds. A high recapture rate in PVC pipes (83% and 86%, respectively) further suggests their use of PVC pipes as daytime refugia. In contrast, 95% of H. gratiosa individuals were captured in pitfall and funnel traps, and only 27 of the 44 total captures in PVC pipes were recaptures. This trend was mainly driven by juveniles; only 5 of 20 juvenile H. gratiosa captures in PVC pipes were recaptures. Although all three hylid species were able to climb out of pitfall traps (pers. obs.), much greater voluntary use of PVC pipes by H. femoralis and H. squirella suggests that detectability of H. gratiosa was much lower. This, combined with an apparently lower abundance of H. gratiosa than H. femoralis or H. squirella [23], could at least partially explain low adult capture rates.

Age-specific detection differences notwithstanding, a higher proportion of juvenile captures (versus adult) suggests that recruitment was less dependent on annual adult breeding effort for H. gratiosa than for H. femoralis or H. squirella. Assuming low year-to-year variation in detectability of breeding adults within species, our “adult breeding effort” metric should accurately reflect yearly breeding effort within, but not necessarily among species. Therefore, we are confident in the above finding of a weaker relationship between breeding effort and recruitment for H. gratiosa than the other two species. Despite far fewer adult captures than other species, juvenile captures of H. gratiosa (n = 285) outnumbered those of H. squirella (n = 190), and the proportion of juvenile:adult first captures was far higher for H. gratiosa (4.7:1; 82.6% juveniles) than for H. femoralis (0.3:1; 29.8% juveniles) or H. squirella (0.3:1; 30.1% juveniles). These proportional differences are potentially explained in part by adult H. gratiosa climbing out of pitfall traps more often than juveniles. Moreover, adult H. gratiosa are larger than adult H. femoralis and H. squirrela and may climb of pitfall traps more easily than these two species. However, it seems unlikely that frequent escapes from pitfall traps by adult H. gratiosa could fully explain a 14x higher juvenile:adult ratio than for H. femoralis or H. squirella. Alternative or additional explanations for a high juvenile:adult ratio for H. gratiosa relative to H. femoralis or H. squirella include potentially larger clutch sizes or higher per capita combined egg and larval survival, despite potentially lower adult breeding effort. Evolution of a single, brief breeding event per year may have strengthened selection for higher offspring survival, and is theoretically an explanation for the evolution a wide variety of amphibians breeding only once per year. If so, adult breeding effort may not be a helpful metric for predicting H. gratiosa (or other similar species that breed once annually) recruitment or population viability.

Hyla femoralis and H. squirella were more dependent on heavy rainfall for successful recruitment than H. gratiosa, despite having relatively shorter larval periods [10,40,41]. Possible reasons are (1) greater reliance on rainfall as a cue to initiate breeding by H. femoralis and H. squirella than H. gratiosa and (or); (2) a greater reliance of H. gratiosa on groundwater-driven high water levels (>70 cm pond depth; Figure 6), independent of recent precipitation, for breeding and/or larval development. The need for rainfall independent of pond depth may be especially true for H. squirella, as only 0.5% of juveniles (n = 1/190) were captured during the relatively dry year of 1998 despite that year’s high water levels. Moreover, water level was not an important predictor of the species’ recruitment (Figure 7).

A delayed association between water levels and annual rainfall can occur in karstic regions [20], including our study area [2], which may impact hylid recruitment. For example, a remarkable 56.5% (n = 161/285) of H. gratiosa juvenile captures across our study occurred during 1998; despite below-average annual precipitation during hylid BMPs, 1998 ranked second among study years for average pond depth, with minimum and maximum pond depths falling well above thresholds (70 cm and 150 cm, respectively, averaged across ponds) predicted to maximize H. gratiosa recruitment. This trend was likely driven by two tropical storms (Jerry and Josephine) in 1995 and 1996 followed by heavy El Niño rainfall in December, 1997 causing high groundwater levels. As a result, modest precipitation in 1998 resulted in pond flooding.

Unlike weekly hydroregime [2], some annual hydroregime variables (median and minimum pond depth) were poorly predicted using climate variables in our Phase II models, highlighting the long-term confounding influence of groundwater-driven lags in water depth following precipitation. This delayed association between heavy rainfall and pond depth was corroborated by annual mean groundwater discharge at nearby Silver Springs, Florida, which was higher in 1998 than any other year of our study [42]. Groundwater-driven filling of ponds may therefore shield some amphibians from short-term drought impacts, provided that the drought is preceded in previous years by heavy rains associated with groundwater recharge. Notably, the frequencies of both droughts and tropical storms are increasing in Florida under climate change [43]. Reducing water withdrawals from the Floridan aquifer, especially during drought years, may represent an essential tool for conserving isolated ephemeral pond-breeding amphibians. This is particularly true for H. gratiosa and other species whose recruitment success is associated with deeper water. Fortunately, our study sites were not located in an area with heavy groundwater withdrawals, which may aid in future hylid persistence.

Because of their generally higher abundance and relatively rapid rates of larval development (Krysko et al., 2019), H. femoralis and H. squirella likely capitalize more on transient rainfall or periods of continuous water in ponds for breeding opportunities than H. gratiosa. This likely explains our finding that decadal capture rates were generally stable or increased throughout our study period for H. femoralis and H. squirella but decreased for H. gratiosa (Figure 9). This finding also agrees with that of a recent study suggesting that H. gratiosa may be more sensitive to environmental stressors than H. femoralis or H. squirella [23]. However, observed lower sensitivity of H. femoralis and H. squirella to past environmental stressors (1994–2017) may not directly translate into future climate resiliency. Adult breeding effort—an important predictor of recruitment for these species—was highly variable, and could be affected by novel climatic changes not yet experienced during our study [44]. For example, variability in annual precipitation decreased over the course of our study but is expected to increase dramatically in future decades [43].

Interestingly, we found only a weak a relationship (|r| < 0.25) between annual adult breeding effort and Phase I climate variables for all three study species. Detection of a strong climate-mediated breeding response could be confounded by wide population fluctuations among ponds and years [45]. Because we were unable to model the effects of climate change on adult breeding effort, our Phase III results should be interpreted cautiously. Adult breeding effort should be carefully monitored as an indicator of population persistence, and active management approaches such as headstarting should be considered if warranted based on future reductions in populations. Climate change is likely to alter breeding effort in amphibian species reliant on weather cues to initiate breeding, but impacts are likely to vary widely among species [44,46]. Moreover, our results showing that precipitation was highly predictive of H. femoralis and H. squirella recruitment—even when pond depths were high—suggests that increasingly erratic precipitation patterns under climate change [43] could impact future populations. Additionally, this correlation suggests both the importance of rainfall in filling breeding ponds and potentially eliciting behaviors such as adult breeding activity or larval partitioning of suitable microhabitats that favor successful recruitment [47].

Spatial isolation of ponds such as those in our study is likely to affect the climate change vulnerability of hylids in multiple ways. Isolated, spatially distant breeding sites in conjunction with poor amphibian dispersal capabilities [15,16] may limit colonization after local extinctions, potentially limiting the likelihood of population persistence. Thus, the importance of frequent on-site recruitment for ensuring population persistence is amplified for species that depend on isolated ponds for breeding and larval development. Our finding that hylid recruitment may decrease under future climate change scenarios highlights the importance of maintaining multiple isolated ponds at varying distances apart so that some populations can be “rescued” by immigration. Equally important is the role of groundwater aquifers in driving pond surface water availability and thus suitable breeding sites for amphibians [48,49]. Aquifer levels will eventually drop after repeated years of low recharge, but can nonetheless provide a critical buffer against short- and medium-term droughts that are expected to increase in Florida under climate change [43]. We recommend monitoring water levels for ponds in karstic regions to identify and prioritize those with more reliable groundwater discharge during amphibian BMPs, recognizing that occasional pond drying is important to preclude fish and reduce aquatic insect predators [50,51]. Regardless, poor recruitment during drought years will likely remain problematic for H. femoralis, H. squirella, and other pond-breeding species whose annual recruitment cycles depend on high same-year rainfall, unless greenhouse gas emissions are reduced rapidly enough to mitigate future climate change.

5. Conclusions

Our results highlight several important takeaways for managing pond-breeding amphibians under climate change and point to critical future research needs. Groundwater conservation and reduced greenhouse gas emissions to slow climate change will increase the likelihood of population persistence for our study species and other pond-breeding amphibians in groundwater-driven wetlands. Our study species differed in recruitment patterns despite their similar breeding seasons and taxonomic position, indicating that other traits (e.g., larval duration, breeding condition preferences, or abundance) were potentially more important for shaping hylid recruitment. Creating artificial ponds tailored to required breeding conditions may aid in future recruitment [52,53], especially for species with specific hydroregime requirements for successful recruitment, such as pond depth (e.g., >70 cm minimum BMP pond depth for H. gratiosa). During long-term droughts when water cannot be maintained in ponds, and for species not resilient to long-term rainfall deficits (e.g., H. femoralis and H. squirella), labor-intensive approaches (e.g., headstarting) may be needed to enhance the likelihood of population persistence. We also found mixed support for a relationship between pre-breeding season rainfall or pond depth and recruitment. We suggest that future studies evaluate impacts of early season conditions on factors impacting subsequent amphibian breeding and egg/larval survival and development, such as food availability for larvae, predator and competitor regulation, and disease transmission. Moreover, additional research is needed to disentangle complex relationships among weather, climate, adult breeding effort, amphibian dispersal distance, and the surrounding upland matrix on amphibian population connectivity among geographically isolated ponds.