Breeding Dynamics of Gopher Frog Metapopulations Over 10 Years

The gopher frog Rana capito is a pond-breeding amphibian and one of many at-risk species currently designated for reviews of current and future population trends (i.e., through Species Status Assessments: U.S. Fish and Wildlife Service 2016; Smith et al. 2018) by the U.S. Fish and Wildlife Service (USFWS) to determine a listing decision under the U.S. Endangered Species Act (ESA 1973, as amended). Regional or range-wide population declines have been observed or suspected over the past century due to extensive habitat loss, degradation, and other stressors (Semlitsch et al. 1995; Enge et al. 2014; Crawford et al. 2020a). Subsequently, the species has been a focus of conservation planning efforts by Federal, State, and other partners to stem ongoing declines and sustain the species across its range. State wildlife agencies, Federal agencies, nongovernment organizations, and academic researchers have increased efforts to survey gopher frog populations, study its life history traits, and estimate key demographic rates (e.g., Greenberg 2001; Roznik et al. 2009; Humphries and Sisson 2012; Enge et al. 2014; Richter et al. 2014). These efforts have revealed an immediate need for research to characterize dynamics of breeding and recruitment to estimate extinction risk for gopher frogs at local and range-wide scales, which will inform forthcoming conservation status and management decisions. Gopher frog monitoring often employs various sampling techniques and intensities, so research estimating breeding frequency and recruitment must also account for imperfect detection rates to reduce potential bias in any population parameter estimates (MacKenzie et al. 2003, 2006; Kéry et al. 2009).

Gopher frogs are one of many amphibian species of conservation concern that breed in isolated and ephemeral or semipermanent wetlands embedded within xeric pine savannas in the southeastern United States, including the once extensive longleaf pine Pinus palustris forest (Dodd 1992; Gibbons et al. 2006; Enge et al. 2014; Erwin et al. 2016). Populations of ephemeral pond-breeding amphibians may be particularly sensitive to dynamic and difficult-to-predict extrinsic factors that influence rates of breeding and recruitment (i.e., of larvae into the terrestrial juvenile stage class). Successful annual breeding and juvenile recruitment from ephemeral wetlands requires local climate conditions that allow wetlands to fill and hydroperiods to be sufficiently long for larval development (e.g., Palis 1998; Greenberg 2001; Jensen et al. 2003; Chandler et al. 2016, 2017). For example, Chandler et al. (2016) estimated that isolated wetlands used by endangered reticulated flatwoods salamanders Ambystoma bishopi only had sufficiently long hydroperiods to support larval development in 61% of years and that frequency was declining over recent decades as a result of climate change. For gopher frogs, Greenberg et al. (2015) forecasted that hydroperiods at most isolated wetlands in a central Florida sandhill would be insufficient, including many multiyear periods, for juvenile gopher frog recruitment. Such low probabilities of juvenile recruitment increase the risk of extirpation. Hydroperiods of isolated wetlands across the Southeast can be further altered by urbanization (McCauley et al. 2013), forest cover, and silvicultural practices (Chandler et al. 2017; Jones et al. 2018; Haggerty et al. 2019), which may interact with climate change to increase extirpation rates of amphibian species such as gopher frogs (Snodgrass et al. 2000; Chandler et al. 2017). In addition, other factors, including the periodic colonization or introduction of fish into isolated wetlands, may increase egg and tadpole mortality and decrease recruitment rates (Gregoire and Gunzburger 2008; Liner et al. 2008). Finally, studies of gopher frog breeding frequency and juvenile recruitment rates are often limited in temporal and spatial extent because of logistical constraints and difficulty predicting relatively brief breeding activity, and they rarely account for imperfect detection when estimating demographic rates (i.e., the probability that breeding occurred but was not observed via sampling methods: MacKenzie et al. 2002; Kéry et al. 2009).

There is clearly both a need and an opportunity to better estimate gopher frog breeding dynamics across broader scales. In this study, we applied an occupancy modeling framework to estimate gopher frog breeding dynamics, which integrated multiple data sets collected across the species' range under different sampling intensities. Our objectives were 1) to estimate the influence of terrestrial and aquatic habitat factors and weather on seasonal breeding probabilities, and 2) to use those estimates to characterize seasonal, annual, and regional breeding dynamics at the metapopulation scale over a 10-y period. This analysis builds on previous collaborative research between the University of Georgia and State, Federal, and other partners to develop status and decision models for at-risk amphibian and reptile species in the longleaf pine system. Previous work on gopher frogs included predicting habitat suitability across the species' range (Crawford et al. 2020a), a multipartner decision-making workshop run by the National Conservation Training Center (USFWS), and additional expert elicitation efforts estimating persistence under current conditions and potential management scenarios.

The study extent encompassed 192 ephemeral to semipermanent wetlands sampled across the southeastern coastal plain in Alabama, Florida, Georgia, North Carolina, and South Carolina in the United States (Figure 1) that were known or expected to have historical or recent gopher frog breeding. Sampled wetlands were typically isolated, ephemeral, and had open canopies (Greenberg 2001; Humphries and Sisson 2012; Enge et al. 2014). The surrounding uplands varied across the species' range but included pine savannas, scrub, flatwoods, and sandhill habitats with well-drained, sandy soil, low canopy cover, high herbaceous ground cover, and adequate density of refugia such as gopher tortoise Gopherus polyphemus burrows or stump holes (Greenberg 2001; Roznik et al. 2009; Humphries and Sisson 2012; Enge et al. 2014; Smith et al. 2020).

Recognizing that gopher frogs can exist in metapopulations of connected wetlands and uplands (Greenberg 2001; Roznik et al. 2009; Humphries and Sisson 2012), we delineated metapopulation boundaries i by buffering all sampled wetlands j and additional terrestrial gopher frog locations (Crawford et al. 2020a) within 3.5 km of each other corresponding to the longest known movement distance (Humphries and Sisson 2012). We subdivided the polygons by features believed by species experts (see Acknowledgments for list of experts and institutions from which communications were obtained in 2019 and 2022) to represent significant barriers to movement, such as corridors of unsuitable habitat or major roads (Figure 1 inset shows seven metapopulations, where five that occur as a contiguous cluster are nevertheless assumed to be demographically isolated from one another). These steps yielded 67 metapopulations (Figure 1). Although some metapopulation boundaries only contained one sampled wetland (see Results), we assumed other nonsampled wetlands could exist and referred to these as metapopulations. Still, we believed the sampled wetlands to be the primary breeding sites within each metapopulation boundary based on previous ground surveys and aerial imagery. We performed all spatial analyses in ArcGIS version 10.4 (ESRI, Redlands, CA).

This study integrated independent detection–nondetection and covariate data sets (Data S2–S7, Supplemental Material) collected by State wildlife agencies and other researchers in all five States to estimate gopher frog breeding dynamics using an occupancy modeling framework (MacKenzie et al. 2002; Royle and Kéry 2007). Although some sampling design components varied between Florida and non-Florida data sets (described below), some methods were consistent. All researchers collected data of evidence of gopher frog breeding—defined as observing egg masses or tadpoles—obtained from timed dipnet and visual surveys of focal wetlands. We structured observation data collected in sampling occasion k into seasons t (an index for the combination of year and season) using the following season designations: Spring (1 March–31 May), Summer (1 June–31 August), Fall (1 September–30 November), and Winter (1 December–28 February). We defined a year as 1 September–31 August. Although breeding can occur year-round, especially in lower latitude or more coastal locations (Enge et al. 2014), summer through early fall generally represents the period of driest conditions and lowest breeding activity that can serve as a reasonable division between one annual term and the next in our model. Season designations aligned with how researchers often structured their monitoring efforts, such that wetlands were targeted for sampling at least once per season. The pooled data set constituted sampling between Fall 2009 and Summer 2019 (40 seasons, across 10 y in total), but each wetland was not sampled every season or year. Researchers noted if wetlands were dry during any visit.

We obtained more detailed detection–nondetection data from the Florida Fish and Wildlife Conservation Commission (FWC: Data S1, Supplemental Material). Researchers from FWC sampled 96 wetlands (grouped into 48 metapopulations) between Fall 2015 and Fall 2018, and they recorded detection–nondetection data from repeat samples of wetlands within a given season. Sampling was sometimes conducted with multiple observers simultaneously dipnetting random areas of a wetland that we treated as separate samples within the same season for this study. Researchers treated samples taken by different observers on the same day or on different days within a season equally as secondary sampling events. Many wetlands were large (between 1 and 47 ha) and remote, posing logistical constraints for field surveys, so researchers often used spatial replication (as opposed to temporal replication) of sampling random locations within a wetland during one visit to obtain information on detection probability.

We also obtained detection–nondetection data from State agencies and additional researchers in Alabama, Georgia, North Carolina, and South Carolina (Data S2). These groups performed opportunistic sampling of 96 wetlands (grouped into 19 metapopulations) between Fall 2009 and Summer 2019. Non-Florida data sets included season- and wetland-specific observations (i.e., whether evidence of breeding was observed), but researchers did not record how many times they visited each wetland each season (see below).

We collected several spatial and temporal covariates during field sampling and previous work that could influence breeding dynamics. We assigned each metapopulation to one of seven regions (Figure 1), which were used when summarizing breeding patterns. In Florida, we assigned sites to one of four regions representing National Oceanic Atmospheric Administration Climate Divisions (https://www.ncdc.noaa.gov/monitoring-references/maps/us-climate-divisions.php [December 2022]). Climate Divisions represent areas of similar patterns of temperature, precipitation, drought, and other climatic factors, based on long-term data. Using Climate Divisions as regions outside of Florida would have resulted in some regions containing only one or two wetlands. Therefore, we assigned sites outside of Florida into one of three larger regions representing grouped EPA IV ecoregions (https://www.epa.gov/eco-research/ecoregions [December 2022]) that had similar ecological and climatic characteristics. We captured site quality, which may influence breeding rates, by calculating the mean habitat suitability index (HSI: 0 = least suitable; 1 = most suitable habitat) within a 1-km buffer of a wetland using predictions from a habitat suitability model (Crawford et al. 2020a). We chose this buffer size to capture habitat quality in the wetland and adjacent uplands while minimizing overlap between buffers of nearby wetlands. The HSI for gopher frog habitat was positively influenced by mean soil drainage, compatible land cover (evergreen forest, mixed forest, shrub–scrub, grassland–herbaceous, and barren land), and fire frequency in the surrounding neighborhood, which ranged from 500 to 2,000 m, depending on the predictor (Crawford et al. 2020a).

We collected additional covariates at each wetland that could vary across seasons and sampling occasions that may influence breeding dynamics. At each sampling occasion, FWC researchers recorded the presence of fish, which can increase tadpole predation risk and decrease the probability of tadpoles surviving to metamorphosis. Fish were originally identified to species, but we pooled data to indicate detection of any species, which included predatory fish (e.g., warmouth Lepomis gulosus) as well as gape-limited fish that may only depredate eggs and small larva (e.g., eastern mosquitofish Gambusia holbrooki). Data indicated high detection probability for fish (i.e., fish were often detected in all or no samples within a season), so we fixed fish presence across all sampling occasions within a season as 0 when fish were not detected and 1 otherwise. For non-Florida sites that did not collect fish presence data, we imputed fish presence each season using a Bernoulli trial with a probability of 0.5. We used the Climate Engine Portal (https://app.climateengine.com/login?ref=https://app.climateengine.com/climateEngine [December 2022]) to download raster data sets summarizing remotely sensed climatic conditions for each season within the study period for all wetlands. We obtained raster data for mean Standardized Precipitation Index (SPI), which measures the deviation of amount of rainfall in a specific period of a year (e.g., March–May 2010) relative to the distribution of rainfall in that same period across many years (McKee et al. 1993). We obtained rasters that calculated SPI for each season relative to long-term climate data from 1980 to 2019. Gopher frogs have a relatively long larval period of approximately 90 to 120 d (Semlitsch et al. 1995). Climatic conditions favorable for breeding in the months preceding and during the season when sampling was conducted may influence the probability of breeding and presence of egg masses and tadpoles; therefore, we developed a final covariate data set where the seasonal SPI represented climate conditions in the sampling season and the one prior (i.e., a 6-mo window). Climate rasters were available at a 4-km resolution, and we extracted seasonal climate values for each wetland.

We developed an integrated model using an occupancy framework to estimate seasonal, annual, and regional responses in gopher frog breeding as a function of wetland- and season-specific covariates. We used a hierarchical state-space formulation fitted in a Bayesian framework (Royle and Kéry 2007; Kéry and Schaub 2012) that separately accounted for the ecological processes of breeding and the observation process of detection, thereby reducing bias and increasing precision of estimates (MacKenzie et al. 2003; Clark and Bjørnstad 2004; Kéry et al. 2009). Specifically, we estimated the breeding probability for wetland j during season t (θ_j,t), and detection probability p_j,t,k during sampling occasion k, given that breeding occurred at a wetland in season t. We assumed closure (i.e., the breeding state of each wetland did not change) within seasons but not between seasons following the “robust design” (Pollock 1982; Kendall et al. 1995) implementation of state-space dynamic occupancy models (Royle and Kéry 2007).

We constrained detection probability (p) to be constant across all Florida wetlands and sampling occasions, and we used p in a cumulative seasonal detection probability for non-Florida wetlands, where K_j,t is a random variable representing the number of sampling occasions (1, 2, 3, or 4) within a season at a wetland modeled as a categorical distribution with probabilities 0.4, 0.3, 0.2, and 0.1, respectively. We chose these probabilities in consultation with data collectors to approximately represent sampling effort.

The model framework required a few key assumptions. Although gopher frogs typically exist in metapopulations with multiple, connected wetlands, we estimated breeding dynamics at the wetland level without explicitly capturing processes in metapopulation models (e.g., patch-specific movement probabilities; Royle and Dorazio 2008). Most metapopulations in our study only included one sampled wetland (see Results), so we included a random metapopulation effect to implicitly account for variation in dynamics that may influence the overall probability of breeding at wetlands within a metapopulation. Our model also assumed metapopulations were persisting across the study period (i.e., breeding was possible each season, when wetlands were not dry) and lacked mechanisms for site survival (1 − extinction) and recolonization typical of dynamic occupancy models (MacKenzie et al. 2003; Bailey et al. 2004; Royle and Kéry 2007; Kéry et al. 2010). Preliminary models that included site survival and recolonization parameters showed poor convergence and produced unreliable estimates that predicted entire metapopulations were becoming extirpated and recolonized multiple times within the 10-y study period. We know that a failure to observe breeding at focal wetlands could be the result of wetlands not filling sufficiently, individuals migrating to proximate alternate breeding sites (Petranka et al. 2004), or most individuals skipping breeding (Bailey et al. 2004; Church et al. 2007; Muths et al. 2010, 2013; Cayuela et al. 2014). Gopher frogs can likely live 6–10 y, similar to longevity estimates of congeners (Richter and Seigel 2002). For these reasons and because metapopulations were separated by long distance (> 7 km) or physical barriers, it was more reasonable to assume gopher frogs were persisting within a metapopulation boundary throughout the study period but exhibiting high year-to-year variation in breeding rather than experiencing multiple extirpation and recolonization events within a 10-y period.

We fit the integrated model using a Bayesian framework with Markov chain Monte Carlo methods in JAGS (Plummer 2003) called from R (R Core Team 2019) via the R2jags package (Su and Yajima 2012). We assigned vague prior distributions for all parameters and hyperparameters governing random effects: breeding probability intercept (expressed on the probability scale) and detection probability—Uniform(0,1); breeding probability covariate parameters (logit scale)—Normal(mean=0, precision=0.37); hyperparameter (standard deviation) governing the random metapopulation effect—Uniform(0,5; see model code in Data S7 for details). We generated three Markov chain Monte Carlo chains using 150,000 iterations where we retained the last 10,000 iterations, yielding a final set of 30,000 samples from posterior distributions of the parameters. We assessed convergence for all models by visually inspecting chain mixing in Markov chain Monte Carlo trace plots, confirming unimodality in posterior distribution plots, and verifying Brooks–Gelman–Rubin statistics < 1.1 for all parameters. We based parameter inferences on posterior means and 95% Bayesian credible intervals (2.5th–97.5th percentile of the distribution).

Of 67 metapopulations, 35 contained one sampled wetland and 32 contained multiple sampled wetlands (range = 2–28). The number of wetlands and metapopulations, respectively, in each region was 26 and 14 (Northwest Florida), 27 and 14 (North Florida), 26 and 12 (North Central Florida), 17 and 8 (South Central Florida), 18 and 3 (Upper Coastal Plain [Carolinas]), 59 and 11 (Atlantic Coastal Plain [Carolinas]), and 19 and 5 (Upper Coastal Plain [Alabama and Georgia]). Researchers of FWC sampled the 96 Florida wetlands between 0 and 11 times within a season between Fall 2015 and Summer 2018, and the mean total number of samples per wetland over the 4-y period was 21.9 (range = 9–37). For the 96 non-Florida wetlands, the mean number of seasons sampled was 10.5 (range = 1–33, out of 40). Over all of the metapopulations studied, gopher frogs were detected on average during 4.16 samples (range = 0–25) at a wetland within the study period. Mean HSI of wetlands was 0.67 (range = 0.02–1.00). Mean seasonal SPI of wetlands was 0.16 (range = −3.07–3.10).

Estimated breeding probability varied considerably by wetland and season and was negatively correlated with fish presence and positively correlated with SPI and prior season breeding state (Table 1). The effect of HSI on breeding was not as strong as other predictors, but we saw some support for an interaction between HSI and SPI (Table 1; Figure 2), indicating that strength of relationship with SPI depended on HSI. The direction of the interaction effect resulted in wetlands with high HSI predicted to have higher breeding probabilities, relative to wetlands with low HSI, during drought conditions (i.e., lower SPI), but this difference became negligible as SPI increased. The estimated mean detection rate from a sample was relatively high (0.69).

Estimated breeding dynamics varied seasonally, annually, and regionally (Figure 3). Mean estimated breeding frequency within the 10-y study period of sampled wetlands was 0.49 (range = 0.11–0.93); mean estimated breeding frequency of metapopulations was 0.68 (range = 0.00–1.00). Mean predicted breeding frequencies within the 10-y study period were within the ranges 0–0.19, 0.2–0.49, 0.5–0.79, and 0.8–1.0 for 6, 13, 16, and 32 metapopulations, respectively. The skewed distribution of metapopulations in the highest breeding frequency category occurred because we calculated breeding frequency as the number of years when breeding was predicted at any wetland within a metapopulation, and 32 of 67 metapopulations contained multiple sampled wetlands. The percentage of wetlands predicted to have breeding varied by year and season but was often highest in Winter and Spring across the species' range (Figure 3). Within a single year cycle (Fall to Summer), mean predicted breeding was highest most often in Winter or Spring in northern (non-Florida) regions. Conversely, some Florida regions had noticeable predicted peaks of breeding in Fall in some years, while breeding was more dispersed across all seasons in other years for North, North Central, and South Central Florida regions (e.g., breeding was highest in Summer, relative to other seasons in a year cycle, for 4 out of 10 y in the study period for these regions: Figure 3). Peaks in predicted breeding often occurred during or in the season following high rainfall events when SPI was high. For example, the increased percentage of wetlands with breeding between Fall 2017 and Spring 2018 in multiple Florida regions corresponded with Hurricane Irma impacting that area in Fall 2017.

Using the estimated detection probability when conducting a sample (timed dipnet sampling and visual surveying), we found that the minimum number of samples needed to achieve a < 5% probability of breeding going undetected at a wetland depended on prior breeding probability of a site (Figure 4). Two samples were necessary to achieve a < 5% probability of failing to detect breeding when prior breeding probability was 0.2; three samples were needed when prior breeding probability was 0.5; and four samples were needed when prior breeding probability was 0.8. In other words, these findings reflect the number of samples needed to be 95% confident that breeding has not occurred when gopher frogs are not observed.

Our work contributes an approach for leveraging multiple available data sets collected by different partners across a species' range, even when sampling timing and frequency vary, to improve estimation of population parameters of interest. It is our experience that State and Federal agencies and other partners often collect data sets for pond-breeding amphibians and other species of concern within their separate jurisdictions that may be compatible for unified analyses similar to those undertaken in this study. Recent efforts have successfully coordinated partners and developed appropriate modeling frameworks for these contexts (e.g., Ahrestani et al. 2017; Crawford et al. 2020b; McGowan et al. 2020). Our study also advanced previous work on gopher frogs by estimating the influence of climatic and wetland-specific factors on breeding dynamics over a 10-y period and across the species' range, which can directly inform site-level management and future models predicting persistence of gopher frog metapopulations. We found that 1) breeding probability at a wetland was positively influenced by seasonal precipitation and negatively influenced by fish presence; 2) there was some evidence that wetlands surrounded by higher quality terrestrial habitat may have increased breeding probability during drought conditions; 3) the percentage of sampled wetlands predicted to have breeding varied seasonally, annually, and regionally across the study, but peaks of breeding in most regions occurred in Winter and Spring and often following high rainfall events (e.g., hurricanes)—with more dispersed breeding across seasons in regions in the Florida Peninsula; and 4) sampling wetlands two to four times within a short period where closure can be assumed (e.g., a few days), depending on prior probability of breeding at a site, can sufficiently reduce the probability of failing to detect breeding events.

Our results offer insights into the effects of climate, habitat, and other wetland-specific factors on gopher frog breeding dynamics that support previous findings of several studies conducted at smaller temporal or spatial scales. Estimated positive relationships between seasonal precipitation (captured with Standardized Precipitation Index in this study) and breeding success agree with previous findings for gopher frogs (Palis 1998; Jensen et al. 2003) and other ephemeral pond-breeding amphibians (e.g., yellow-bellied toad Bombina variegata; Cayuela et al. 2014). Future precipitation patterns will likely influence or restrict successful breeding opportunities for gopher frogs and other pond-breeding amphibians, but these patterns are uncertain in the Southeast. Greenberg et al. (2015) predicted that future precipitation and other climatic factors would shorten hydroperiods for ephemeral wetlands in Florida, whereas Keellings and Engström (2019) predicted increased rainfall and decreased number of consecutive dry days throughout most of the gopher frog's range, but we note the latter study did not include other climatic factors (e.g., temperature) that may influence hydroperiod. Additionally, increased rainfall and flooding events could increase the frequency with which isolated wetlands become temporarily connected to other waterbodies and colonized by fish. Invasive (e.g., African jewelfish Hemichromis letourneuxi) and other predatory fish have already been detected in gopher frog wetlands following floods in southern Florida (A.H.G., personal observation).

We found some evidence that the effects of seasonal droughts could be mitigated by increasing habitat suitability within and around wetlands. Specifically, high HSI of the surrounding landscape—characterized by well-drained sandy soils, open-canopied land cover types (evergreen or mixed forest, shrub–scrub, grasslands), and frequent fires (Crawford et al. 2020a)—tended to increase the breeding probability during drier periods (when SPI was negative: Figure 2). Burrow and Maerz (2021) demonstrated experimentally that shading and hardwood leaf litter, which occur with succession within wetland basins when fire is excluded and/or suppressed, interact to extend the larval period of gopher frogs, making them more susceptible to wetland drying. Greenberg (2001) observed that metamorph recruitment was significantly lower in wetlands where hardwoods had invaded wetland basins. We also expect that upland habitat suitability, including the maintenance of high densities of refugia (e.g., tortoise burrows, mammal burrows, and stump holes) and ground cover, has a strong effect on gopher frog survival in terrestrial stages (Roznik and Johnson 2009; Roznik et al. 2009; Humphries and Sisson 2012; Burrow 2021; Burrow et al. 2021). The terrestrial stages of many pond-breeding amphibians are key “storage” phases for populations to withstand episodic droughts that limit breeding success (Taylor et al. 2006); therefore, habitat conditions that increase terrestrial fitness of gopher frogs may be more likely to sustain some breeding success during drier years and lead to resilient populations. Field research measuring habitat effects on demographic rates could complement future models that separate population persistence and breeding dynamics to better characterize the effects of climate and terrestrial and aquatic habitat conditions on those dynamics.

Similar to our study, Cayuela et al. (2014) found breeding success in one season was positively influenced by the breeding state in the previous season. This may be, in part, an artifact of the length of gopher frog larval periods (approximately 90–120 d; Semlitsch et al. 1995) relative to the typical length of a “season” in our model (89–91 d). It is reasonable that successful breeding in one season would provide tadpoles to be detected in the subsequent season. This does not preclude the possibility that breeding rates in consecutive seasons and years are not independent as a result of natural cycles of adult abundance, density-dependent mortality of tadpoles (Greenberg 2001), high mortality of first-time breeding adult cohorts (V. Terrell, University of Georgia, unpublished data, 2022), and the propensity for adults to skip breeding seasons, which has been observed among many other amphibians (Bailey et al. 2004; Church et al. 2007; Muths et al. 2010, 2013; Cayuela et al. 2014). Future work could estimate factors that determine the effects of past breeding (e.g., during any season in the previous year) on subsequent breeding while also estimating other demographic processes to better capture longer term population cycles.

The percentage of sampled wetlands predicted to have breeding varied seasonally, annually, and regionally across the study. Peaks of breeding in most regions occurred in Winter and Spring corresponding with high seasonal rainfall, which is consistent with previous studies of gopher frog breeding activity (Palis 1998; Enge et al. 2014). The 10-y period of our study also highlighted the association between peak gopher frog breeding events to peak rainfall events in certain years and regions, such as the peak of breeding in multiple regions in Florida in Fall 2017 to Spring 2018 that corresponded with Hurricane Irma (Figure 3). Breeding tended to be more evenly dispersed across seasons for North, North Central, and South Central Florida regions: these regions had noticeable predicted peaks of breeding in Fall in some years, and breeding was highest in Summer, relative to other seasons in a year cycle, for 4 of 10 y in the study period for these regions (Figure 3). These breeding patterns in the Florida Peninsula corroborate previous findings documenting frequent breeding across all seasons in this area (Enge et al. 2014). We also note that no distinct peaks of breeding were predicted for the Atlantic Coastal Plains (Carolinas) region. This region had substantially more wetlands than any other region, so the trend of lower season-to-season variation may reflect a degree of population and regional resilience where breeding is likely to occur in at least some wetlands each season.

Across the study period, mean predicted annual breeding frequency of metapopulations was high, on average (0.68), and > 0.8 for 32 of 67 metapopulations. We note that surveyed metapopulations tended to be “strongholds” for the species where semiconsistent breeding activity had been observed previously on sites located on public or private lands managed for wildlife conservation. Other multiyear studies of gopher frogs have observed lower breeding frequencies (0.06–0.67) among years at wetlands, but they did not account for detection probability. Regardless, they have suggested that, although metapopulations could be stable over long periods with infrequent breeding and recruitment, these conditions likely increase the risk of local extirpation (Semlitsch et al. 1995; Greenberg 2001).

Recent work has demonstrated the value of developing integrated models that leverage multiple available data sets for data-limited species (Crawford et al. 2018, 2020b; Robinson et al. 2018; Zipkin and Saunders 2018). Integrating multiple data sets in unified analyses can increase precision of estimates, reduce effects of potential bias of individual data sets, and allow for estimation of additional parameters of interest (Schaub and Abadi 2011; Zipkin and Saunders 2018). However, developing integrated analyses requires that the underlying demographic processes and observation protocols be shared across data sets, or that ancillary information is available to analytically resolve differences among data sets. A key assumption of our model was that sampling protocols (and the probability of detecting egg masses and tadpoles with a sample) were equivalent across sites that were sampled by different partners. We confirmed that partners used standard dipnet and visual surveys at all sampled sites, but we were not able to incorporate additional details related to sampling (e.g., sampling time) into the analysis. We also had to include uncertainty around the number of samples per season at all wetlands outside of Florida because partners did not record this information systematically. Our analysis overcame these data gaps in non-Florida samples by leveraging information from Florida samples that aided in estimation of detection probability. However, more accurate estimates of breeding dynamics and other population parameters can be achieved if partners conducting sampling, including Federal and State agencies, systematically document and store detection–nondetection data and additional information related to sampling each time a wetland is visited. We hope this study illustrates the importance of keeping records of sampling efforts and findings for future needs.

Larger scale monitoring programs often collect detection–nondetection data of species of interest using spatial replication (i.e., observers sample multiple areas at one time within a site) rather than temporal replication (i.e., observers repeatedly sample one or more areas within a site) in cases where limited monitoring resources prevent repeated visits to large, remote sites. Estimates of detection and occupancy are prone to bias when spatial replication is used in place of temporal replication; however, bias is negligible under certain sampling and ecological contexts (MacKenzie et al. 2006; Kendall and White 2009; Hines et al. 2010). MacKenzie et al. (2006) discussed that the use of spatial replication at a single time (e.g., on a single day) leads to negligible bias in model estimates if multiple sample locations within a site are selected at random with replacement (i.e., during selection of random sample locations at a site of interest, any location may be selected by chance for sampling multiple times at the same visit). The FWC data set satisfies these criteria because it represents a combination of 1) spatial replication when multiple observers randomly surveyed different areas of a wetland on a single day, and 2) for a subset of wetlands, temporal replication when one or more observers returned to survey a wetland (including resampling the same areas) later within the closed season. Bias from substituting spatial for temporal replication is also minimal when the modeled species is highly mobile or distributed randomly across sampling locations within a site, which allows some probability of observing all detection histories within a season (Kendall and White 2009; Hines et al. 2010). Although localized patchiness of gopher frog egg masses and tadpoles can occur during breeding events with a small number of breeders, our partnering species experts indicated that egg masses can be distributed widely across a wetland (see Acknowledgments for list of experts). Given the contexts of data collection and species' ecology in our study, we expect bias of model estimates to be inconsequential. For monitoring data from other contexts that used spatial replication without replacement or surveyed less mobile species, our model would likely overestimate occupancy (i.e., breeding) probability. Other model structures could be used in these contexts to account for potential biases (see Hines et al. 2010). For example, estimating occupancy with time-to-detection methods—where detection probability is estimated using the time between the start of a survey and first observation of the species—provides an efficient alternative with similar performance to methods using temporal replication for a variety of ecological contexts, although these methods performed best for widespread species with high detection probabilities similar to our study (Halstead et al. 2021). Future monitoring efforts could collect data using spatial replication with replacement, which 1) efficiently uses resources even under logistical challenges of field surveys, and 2) minimizes bias of population parameters that can inform conservation decisions (Kendall and White 2009).

Our analysis had limitations that should be considered when interpreting the results. Our data were insufficient to construct a model that could reliably separate breeding dynamics from metapopulation persistence or other mechanisms typical of metapopulation models (e.g., wetland-specific immigration and extinction rates). We assumed that populations were persisting over the 10-y period and breeding could potentially occur. Gopher frogs have an expected lifespan of 6–10 y and were detected at least once during the study in 56 of 67 metapopulations, which supported the assumption of persistence at most sites. The 11 metapopulations where gopher frogs were not detected had only one or two sampled wetlands, no other known breeding sites, and estimated breeding frequencies of < 0.3 (i.e., < 3 of 10 y with breeding). Although it is possible that gopher frogs used nonsampled, nearby wetlands for breeding during this period, it is likely that some of these metapopulations are extirpated, given that partners have observed degraded wetland and upland habitat conditions and no detections of gopher frogs during additional surveys at some of these sites. Including extirpated metapopulations in our analysis could have potentially underestimated detection and overestimated breeding probabilities; however, we believe any bias introduced by their inclusion was limited because only 166 of 2,910 surveys (5.7%) came from wetlands in the 11 potentially extirpated metapopulations.

Although sampled wetlands represented what partners believed were the primary breeding sites within metapopulation boundaries based on previous monitoring and habitat characteristics, the true number of potential breeding wetlands was unknown. Therefore, we included a random effect in the model to account for metapopulation-to-metapopulation variation in number of breeding wetlands and other factors that may influence breeding dynamics but were not explicitly captured. We also note that egg and tadpole observation data used in this study is informative about breeding dynamics and adult occupancy, but additional information about direct rates of recruitment into terrestrial life stages would help determine if individual breeding wetlands are serving as population sources or sinks, which could better inform management decisions. Otherwise, the presence of declines might be masked by unsuccessful breeding attempts of long-lived adults. Future efforts could use an integrated modeling framework that operates over multiple scales of effort to estimate persistence and other metapopulation dynamics; for example, detection–nondetection data from dipnet sampling of many wetlands and metapopulations could be leveraged with capture–mark–recapture data from a subset of metapopulations (Schaub and Abadi 2011; Zipkin and Saunders 2018).

Our results have direct applications for site-level monitoring and management aimed at increasing successful breeding opportunities of gopher frogs and other associated pond-breeding amphibians. First, partners could tailor monitoring of breeding activity to align with predicted seasonal peaks and conduct two to four samples within a season (ideally within a few days when population closure can be assumed), depending on prior knowledge or estimates of breeding frequency at a site, to minimize the probability of missing breeding events. Gopher frog wetlands are often dispersed across a state or region making it logistically difficult to monitor with repeated visits; therefore, conducting two to four samples could be accomplished using multiple observers during a single visit, if sample locations can be selected randomly with replacement.

Second, partners can consider efforts to monitor and prevent fish colonization or remove fish when detected in known gopher frog breeding wetlands. Limiting fish presence in ephemeral breeding wetlands has become an interest of amphibian conservation efforts by our partnering species experts and their management partners (see Acknowledgments for list of experts), because predatory and invasive fish are commonly detected in wetlands (Enge et al. 2014) and predation can decrease egg and larval survival (Gregoire and Gunzburger 2008). In some cases, this may require hydrological restoration of uplands or wetlands, if ditching or other alterations have increased fish colonization or lengthened the wetland hydroperiod. Removal of invasive and predator fish with piscicides (e.g., rotenone), applied carefully when rare amphibians are absent from the wetland to avoid nontarget mortality of amphibians (e.g., Shulse and Semlitsch 2014), is a strategy under consideration in some conservation efforts.

Third, our analysis suggests that improving upland habitat conditions surrounding breeding wetlands by providing a relatively open canopy and midstory with mechanical removal or prescribed fire may increase wetland hydroperiods in periods of drought and provide more opportunities for successful breeding (Chandler et al. 2017; Simpson et al. 2021). This may be particularly important for gopher frogs as increasing temperatures and changes in the hydrologic cycle fueled by climate change are predicted to increase the duration and intensity of droughts (USGCRP 2018). Furthermore, a recent habitat suitability model classified only 8.6% of habitat with the gopher frog's range as either moderate or highly suitable for the species, suggesting a need for habitat restoration throughout the species' range (Crawford et al. 2020a).

This work precedes research that will 1) incorporate breeding frequency into a population viability analysis for the gopher frog that could be adapted for other pond-breeding amphibians, and 2) use the population viability analysis to evaluate the effects of stressors and management alternatives on breeding frequency and population persistence to inform conservation decisions at local to range-wide scales. Although we expect the population viability analysis to initially require expert opinion and thus have high uncertainty around certain parameters, effective monitoring of breeding activity using information from this study could allow future predictions of population persistence to be iteratively updated and improved. Rigorous predictions of gopher frog population persistence, which account for breeding dynamics, can inform forthcoming assessments of extinction risk and designation of the species' conservation status by the U.S. Fish and Wildlife Service.

Please note: The Journal of Fish and Wildlife Management is not responsible for the content or functionality of any supplemental material. Queries should be directed to the corresponding author for the article.

Data S1. Dataframe of detection–nondetection data for Florida wetlands (N = 96) sampled between 2015 and 2018 in the southeastern United States. Fields (columns) include wetland identification number (ID), state name, and detection histories by season and sample for gopher frogs Rana capito.

Available: https://doi.org/10.3996/JFWM-21-076.S1 (11 KB CSV)

Data S2. Dataframe of detection–nondetection data for non-Florida wetlands (N = 96) sampled between 2009 and 2019 in North Carolina, South Carolina, Georgia, and Alabama in the southeastern United States. Fields (columns) include wetland identification number (ID), state name, and detection histories by season and sample for gopher frogs Rana capito.

Available: https://doi.org/10.3996/JFWM-21-076.S2 (14 KB CSV)

Data S3. Dataframe of wetland-specific information for sampled wetlands in Florida (N = 96) between 2015 and 2018. Fields (columns) include wetland identification number (ID), metapopulation identification number (Pop_ID), state name, mean Habitat Suitability Index (HSImean_1000), standardized precipitation index data for each season of the study (column names begin with “spi”), the site's regional group for the study (Region), and the site's EPA ecoregion group (EPA_LICAP).

Available: https://doi.org/10.3996/JFWM-21-076.S3 (38 KB CSV)

Data S4. Dataframe of wetland-specific information for sampled wetlands in North Carolina, South Carolina, Georgia, and Alabama (N = 96) between 2009 and 2019. Fields (columns) include wetland identification number (ID), metapopulation identification number (Pop_ID), state name, mean Habitat Suitability Index (HSImean_1000), standardized precipitation index data for each season of the study (column names begin with “spi”), the site's regional group for the study (Region), and the site's EPA ecoregion group (EPA_LICAP).

Available: https://doi.org/10.3996/JFWM-21-076.S4 (35 KB CSV)

Data S5. Dataframe summarizing the year, season, and number of samples for the Florida data (sampled for 96 wetlands between 2015 and 2018).

Available: https://doi.org/10.3996/JFWM-21-076.S5 (2 KB CSV)

Data S6. Dataframe of wetland identification number (ID), state name, and detection histories by season and sample (see column names) for any fish observed during dipnet sampling for Florida wetlands (N = 96) sampled between 2015 and 2018.

Available: https://doi.org/10.3996/JFWM-21-076.S6 (11 KB CSV)

Data S7. Program R code to reproduce models and outputs for estimating gopher frog Rana capito breeding dynamics in wetlands (N = 192) between 2009 and 2019 in the southeastern United States. The code uses all dataframes included in Data S1–6.

Available: https://doi.org/10.3996/JFWM-21-076.S7 (17 KB R)

Reference S1. Enge KM, Farmer AL, Mays JD, Castellón TD, Hill EP, Moler PE. 2014. Survey of winter-breeding amphibian species. Final report. Florida Fish and Wildlife Conservation Commission, Fish and Wildlife Research Institute.

Available: https://doi.org/10.3996/JFWM-21-076.S8 (8.993 MB PDF) and https://f50006a.eos-intl.net/ELIBSQL12_F50006A_Documents/14_Enge.pdf

Reference S2.U.S. Fish and Wildlife Service. 2016. Long-term listing transformation and five guiding principles of the unified listing team. Director's memorandum 16 March 2016.

Available: https://doi.org/10.3996/JFWM-21-076.S9 (98 KB PDF) and https://dx.doi.org/10.3996/052017-JFWM-041.S9

This work is funded and supported by the U.S. Fish and Wildlife Service through the former Gulf Coastal Plains and Ozarks Landscape Conservation Cooperative (Cooperative Agreement Award F16AC00985) and the Natural Resources Conservation Service (Award 68-7482-16-545 and Agreement NR183A750001C001). We thank M. Harris and D.T. Jones-Farrand for supporting the development of this work and multiple individuals at the Southeast Climate Adaptation Science Center for technical guidance. We thank all individuals who generously contributed data and expert input on population boundaries (in 2019) from the North Carolina Wildlife Resources Commission (J. Hall, J. Humphries), South Carolina Department of Natural Resources (A. Grosse), University of Georgia Savannah River Ecology Laboratory (S. Lance), Georgia Department of Natural Resources (M. Elliott, D. Sollenberger), Jones Center at Ichauway (L. Smith), Florida Fish and Wildlife Conservation Commission (M. Sisson), and Alabama Department of Conservation and Natural Resources (E. Shelton-Nix, D. Colbert). Additional data were also generously provided by Archbold Biological Station, Auburn University, Department of Defense, U.S. Army (Ft. Stewart), and U.S. Air Force (Eglin Air Force Base. We thank B. Harris, M. Fedler, L. Tirado, P. Hill, and many others who contributed to field sampling in Florida. We thank P. Howell, B. Halstead, and two other reviewers for providing comments that improved the manuscript. The Georgia Cooperative Fish and Wildlife Research Unit is jointly sponsored by U.S. Geological Survey, the University of Georgia, U.S. Fish and Wildlife Service, Georgia Department of Natural Resources, and the Wildlife Management Institute.

Any use of trade, product, website, or firm names in this publication is for descriptive purposes only and does not imply endorsement by the U.S. Government.