For species with variable phenology, it is often challenging to produce reliable estimates of population dynamics or changes in occupancy. The Arizona Toad (Anaxyrus microscaphus) is a southwestern USA endemic that has been petitioned for legal protection, but status assessments are limited by a lack of information on population trends. Also, timing and consistency of Arizona Toad breeding varies greatly, making it difficult to predict optimal survey times or effort required for detection. To help fill these information gaps, we conducted breeding season call surveys during 2013–2016 and 2019 at 86 historically occupied sites and 59 control sites across the species' range in New Mexico. We estimated variation in mean dates of arrival and departure from breeding sites, changes in occupancy, and site-level extinction since 1959 with recently developed multi-season staggered-entry models, which relax the within-season closure assumption common to most occupancy models. Optimal timing of surveys in our study areas was approximately 5–30 March. Averaged across years, estimated probability of occupancy was 0.58 (SE = 0.09) for historical sites and 0.19 (SE = 0.08) for control sites. Occupancy increased from 2013 through 2019. Notably, even though observer error was trivial, annual detection probabilities varied from 0.23 to 0.75 and declined during the study; this means naïve occupancy values would have been misleading, indicating apparent declines in toad occupancy. Occupancy was lowest during the first year of the study, possibly due to changes in stream flows and conditions in many waterbodies following extended drought and recent wildfires. Although within-season closure was violated by variable calling phenology, simple multi-season models provided nearly identical estimates as staggered-entry models. Surprisingly, extinction probability was unrelated to the number of years since the first or last record at historically occupied sites. Collectively, our results suggest a lack of large, recent declines in occupancy by Arizona Toads in New Mexico, but we still lack population information from most of the species' range.

NATURAL history museums and archived data are valuable resources for determining changes in species distributions or how common they are across landscapes (Sullivan, 1993; Shaffer et al., 1998; Tingley and Beissinger, 2009; Meineke et al., 2019). Reliance on surveys of historically occupied sites based on museum records may often overstate declines, however, because assuming a non-zero patch extinction rate, occupancy of historically occupied sites can only decrease over time (Skelly et al., 2003; Fournier et al., 2019). Sampling a combination of historically occupied sites plus control sites selected without prior knowledge of a species' presence can provide reliable information on changes over time and unbiased estimates of current occupancy (Skelly et al., 2003; Weiser et al., 2020).

The challenge of producing reliable population estimates at local and regional scales is magnified for species with variable phenology, such as the Arizona Toad (Anaxyrus microscaphus). The Arizona Toad primarily breeds in slow-flowing streams of the southwestern United States. Their breeding season can begin from February to April, depending on elevation, and is initiated by warming nighttime temperatures rather than rainfall (Schwaner and Sullivan, 2009). Arizona Toads breed during a four- to eight-week period in early spring when males call from streamside habitats to attract females (Blair, 1955; Sullivan, 1992; Ryan et al., 2015). Timing and consistency of Arizona Toad breeding varies greatly throughout its range (Blair, 1955; Sullivan, 1992; Degenhardt et al., 1996), however, making it difficult to predict optimal survey times and the number of surveys required to be confident that males occupying a patch of breeding habitat will be detected.

The Arizona Toad has been petitioned for protection under the U.S. Endangered Species Act because of threats from hybridization with Woodhouse's Toads (A. woodhousii), changes in hydrological regimes, and drought (Bradford et al., 2005; Schwaner and Sullivan, 2009; USFWS, 2015; Wooten et al., 2019). Suspected declines have been reported from Nevada, Arizona, and New Mexico, and the species is a Species of Greatest Conservation Need in all four states where it occurs, which includes Utah (Bradford et al., 2005; AZGFD, 2015; NVDW, 2015; UDWR, 2015; NMDGF, 2016). Due in part to the unusual natural history of Arizona Toads among amphibians in the region, conservation assessments are confounded by a lack of information on population status and insufficient evidence of presumed threats.

Accounting for variation in phenology and detection, and determining optimal timing and survey effort, is critical for producing reliable population status information for the Arizona Toad. To help fill information gaps, we conducted breeding season call surveys during 2013–2016 and 2019 across the range of the species in southwestern New Mexico. We surveyed 86 sites where museum records indicated presence of toads during opportunistic collections from 1959–2009 (hereafter, historical sites) and 59 nearby control sites where we had no knowledge of previous occurrence by toads. We used recently developed staggered-entry multi-season occupancy models to help account for variation in calling activity (Kendall et al., 2013; Chambert et al., 2015). The staggered-entry design is like other multi-season occupancy models that estimate site occupancy, colonization, and extinction while accounting for imperfect detection, but it relaxes the within-season closure assumption by incorporating additional entry and departure parameters. These entry and departure parameters allow estimation of the true state of patches that can transition to occupied (entry) or unoccupied (departure) between surveys (Kendall et al., 2013; Chambert et al., 2015).

To provide guidance for future surveys of Arizona Toads and other species with similar variation in phenology, we compared results from staggered-entry models to those from standard multi-season occupancy models (MacKenzie and Royle, 2005). With the standard multi-season model, entry and departure from patches or variation in calling activity (i.e., true within-season variation in the ability to detect a species) is considered detection error and may provide a less accurate accounting of patch dynamics than staggered-entry models. But when estimating species' conservation status or population trends is of primary interest rather than characterizing variation in phenology, the simple multi-season model—which estimates fewer parameters than the staggered-entry design—may be more desirable because it can produce more precise estimates (Chambert et al., 2015). We compared three sets of dynamic occupancy models with similar covariates to (1) estimate annual occupancy of Arizona Toads during 2013–2019, especially with regard to site type (historical vs. control sites) and habitat type (lotic vs. lentic); (2) describe temporal variation in calling activity and how it affects detection; (3) compare occupancy estimates from staggered-entry and simple multi-season models; and (4) estimate the rate of change in occurrence at historically occupied sites.

MATERIALS AND METHODS

Study system.—

In New Mexico, the Arizona Toad most often occurs between 1,800–2,700 m elevation in the upper Gila River Basin, including the Gila and San Francisco Rivers, the upper Mimbres River basin of the eastern Mogollon Rim, and in a disjunct population that occurs in the Rio Grande basin east of the Continental Divide in the San Mateo Mountains (Degenhardt et al., 1996; Jennings et al., 2010; Ryan et al., 2017). Most breeding in New Mexico occurs from early March to late April, when males call from stream margins, side channels, or outflows from reservoirs (Sullivan, 1992; Degenhardt et al., 1996). Arizona Toads will also occasionally breed in lentic stock tanks or similar isolated impoundments (M. Ryan, pers. obs.).

Historical records and site selection.—

To identify historical survey locations, we queried eight regional and two national museum collections (Museum of Southwestern Biology, Western New Mexico State University, New Mexico State University, University of Kansas, University of Arizona, United States National Museum, Los Angeles County Museum, Texas Cooperative Wildlife Collection, University of Texas at El Paso, and the American Museum of Natural History) with the largest Arizona Toad collections (Degenhardt et al., 1996; Ryan et al., 2015, 2017). See Ryan et al. (2017) for details on museum queries and specimens. There were six instances where separate collection events were within 1.5 km of each other. In these cases, we drew a centroid around the collection locations and assigned the center point as the historical survey site. The life stage of museum specimens (egg, larva, or metamorph) confirmed the presence of toad breeding activity at 75 percent of historical sites.

For control sites, we selected suitable waterbodies near historical sites. Even small groups of calling toads can be heard from ∼1.0 km (M. Ryan, pers. obs.). Therefore, all but 15 control sites (5 percent) were at least 1.5 km from historical sites. For the 15 exceptions, landscape or topographic features increased the effective distance and ensured sampling independence between sites. The median distance between historical sites and the nearest control site was 4.79 km (sd = 7.76).

Calling surveys.—

We conducted standardized surveys during 2013–2016 and 2019 at historical and control sites throughout most of the known species' range in southwestern New Mexico (mean elevation of sampled sites: 1,879 m; Fig. 1). Night calling surveys were conducted weekly during the breeding season, beginning in early March and ending late May (only 10 surveys occurred during May). We divided the approximately 21,727 km2 sampling area into two broad units: one encompassed the Mimbres River drainage and Middle, East, and West Forks of the Gila River; the other encompassed the San Francisco River drainage and upper Gila River. Sampling occurred throughout the species' approximately eight-week seasonal breeding period. During each survey, a single observer listened for three minutes and noted the detection or non-detection of calling toads (Ryan et al., 2015). We stopped nightly surveys no later than 0030 hours, even though toads would call as late as 0200 hours.

Fig. 1.

Study area in southwestern New Mexico showing historically occupied (n = 86; red circles) and control (n = 59; blue triangles) waterbodies sampled during the 2013–2016 and 2019 breeding seasons for the Arizona Toads (Anaxyrus microscaphus). Filled shapes represent sites where we detected toads at least once, and open shapes represent sites where we never detected toads. Basemap by Stamen Design, under CC BY 3.0.

Fig. 1.

Study area in southwestern New Mexico showing historically occupied (n = 86; red circles) and control (n = 59; blue triangles) waterbodies sampled during the 2013–2016 and 2019 breeding seasons for the Arizona Toads (Anaxyrus microscaphus). Filled shapes represent sites where we detected toads at least once, and open shapes represent sites where we never detected toads. Basemap by Stamen Design, under CC BY 3.0.

The sampling design used during 2013–2016 precluded distinguishing detection error from real variation in calling activity at calling sites. Therefore, in 2019 we modified the survey design by having the single observer record (with a Sony PCM-10 recorder) their 3-minute sample session using a Wildtronics parabolic microphone directed toward the potential calling site (files are archived at the Museum of Southwestern Biology). An observer with no knowledge of field results later listened to each recording and assigned detection or non-detection of calling toads, providing a second survey of the same period that we used to evaluate detection error.

Analysis.—

We used PRESENCE version 12.39 to perform three analyses that estimated initial occupancy, colonization, and extinction while accounting for imperfect detection. Our primary analysis was a staggered-entry multi-season model based on all surveys from 2013 through 2019. For initial occupancy and colonization, we evaluated the importance of site type (86 historical vs. 59 control) and habitat type (21 lentic vs. 124 lotic; Table 1). Given the low initial occupancy in 2013, we did not have enough information to estimate the effect of site or habitat type on extinction. We also evaluated trends (year treated as a continuous covariate) for colonization and extinction (Table 1). To model detection, we included year of sampling as either a trend or factor variable, plus the effects of cloud cover, moon, air temperature, and wind during call surveys. For staggered-entry models, detection is the probability of detecting toads in a site that was truly occupied. We included the effects of year (continuous) and Julian date (day of year) to describe variation in entry and departure parameters (Table 1).

Table 1.

Covariates evaluated for each parameter in the staggered-entry multi-season occupancy models for Arizona Toads (Anaxyrus microscaphus) from 141 sites in New Mexico. Julian date was included in the detection parameter for the simple multi-season model, which does not estimate the entry and departure parameters.

Covariates evaluated for each parameter in the staggered-entry multi-season occupancy models for Arizona Toads (Anaxyrus microscaphus) from 141 sites in New Mexico. Julian date was included in the detection parameter for the simple multi-season model, which does not estimate the entry and departure parameters.
Covariates evaluated for each parameter in the staggered-entry multi-season occupancy models for Arizona Toads (Anaxyrus microscaphus) from 141 sites in New Mexico. Julian date was included in the detection parameter for the simple multi-season model, which does not estimate the entry and departure parameters.

We evaluated models in a step-down fashion, using AIC values to select the most parsimonious model (Doherty et al., 2012). We used coefficient effect sizes and confidence intervals to identify uninformative parameters (Arnold, 2010) and assessed the stability of covariate coefficients throughout the model table to identify potential confounding variables. For supported covariates, we used estimated coefficients to predict probabilities for model parameters across the range of collected data (MacKenzie et al., 2006). Confidence intervals for each predicted estimate were estimated based on the delta method (Seber, 1982), using packages msm (Jackson, 2011) and emdbook (Bolker, 2020) in program R v3.6.3 (R Core Team, 2020).

After identifying the best-supported staggered-entry model, we fit two subsequent models with the same covariates. First, given the relative complexity of implementing and interpreting staggered-entry models, we compared our best-supported staggered-entry model to the more common simple multi-season occupancy model. With the standard multi-season design, entry and departure from patches (i.e., true within-season variation in the ability to detect a species) is a violation of model assumptions and is absorbed into detection (Chambert et al., 2015). Therefore, we included Julian date on the detection parameter for the simple multi-season models but otherwise held the model structure the same as the top-ranked staggered-entry model. Second, to provide insight into long-term rates of patch-level extinction for Arizona Toads, we re-fit the top-ranked model using only the 62 historical sites for which we had information on years of historical observations. We evaluated two added covariates for initial occupancy derived from museum records and archived field notes: years since first historical observation for a given site (11–57, mean = 36.5, sd = 11.8) and years since the most recent historical observation for a given site (10–56, mean = 30.2, sd = 11.1).

RESULTS

Calling surveys.—

We conducted 1,986 calling surveys for Arizona Toads at 145 sites during 2013–2016 and 2019. We averaged 13.7 (sd = 8.33) visits per site across all years, although the number of sites and visits varied between years (Table 2). Comparison of detections from in-person surveys and recordings of 278, 3-minute surveys during 2019 revealed almost perfect concordance. There were only two cases that differed in reported toad detection, and both discrepancies were occasions where field personnel heard toads during the site visit but the independent observer that listened to the recording did not. This result indicates observer error in the field was near zero, and thus the most likely source of false negatives (assigning a patch as vacant when it was truly occupied) was from variation in toad calling activity. Based on these low error rates, we did not use data from the recordings in subsequent analyses.

Table 2.

Number of sites, number of visits, and naïve occupancy (proportion of sites where toads were heard calling) for 86 historically occupied and 59 control sites surveyed for Arizona Toads in New Mexico.

Number of sites, number of visits, and naïve occupancy (proportion of sites where toads were heard calling) for 86 historically occupied and 59 control sites surveyed for Arizona Toads in New Mexico.
Number of sites, number of visits, and naïve occupancy (proportion of sites where toads were heard calling) for 86 historically occupied and 59 control sites surveyed for Arizona Toads in New Mexico.

Staggered-entry multi-season model.—

The top model from our staggered-entry analysis indicated that detection varied by year (factor) and percent of moon visible (Table 3). Estimates from this model showed toad detection was highest in the first year of the study (2013) and decreased in all but one of the subsequent years (Table 4). While the covariate for moon visibility was in the top-ranked and several other highly ranked detection models, the estimated effect was small and uninformative (b = –0.0001, 95% CI = –0.0003–0.0001; Arnold, 2010). Based on entry and departure estimates from the best-supported staggered-entry model, the probability that toads were present and calling was highest from approximately 5 March to 30 March and peaked at 0.80 on 14 March (Fig. 2). This top-ranked model included the historical vs. control site covariate to describe variation in occupancy, but not habitat. Estimated probabilities of occupancy in 2013 were similar for historical sites (0.158, 95% CI = 0.080–0.288) and control sites (0.067, 95% CI = 0.016–0.241), but because of higher colonization rate at historical sites, occupancy diverged thereafter (Table 5, Fig. 3A).

Table 3.

Comparison of several staggered-entry occupancy models and the covariates used in each for occupancy (ψ), colonization (γ), extinction (), entry (E), departure (D), and detection (P) parameters for the Arizona Toad, as well as the model structures for the simple multi-season (SMS) model and models based only on historical sites. Covariates included in the table are historical vs. control sites (Hist/Con), lentic vs. lotic sites (Len/Lot), Julian date (Jul), percent moon cover (M), wind (Wi), air temperature (C, Celsius), year as a factor (YrF), and year as a trend (YrT). Staggered-entry models are ordered by delta AIC scores (ΔAIC), model weights (W), and –2 log likelihood values (LogLik).

Comparison of several staggered-entry occupancy models and the covariates used in each for occupancy (ψ), colonization (γ), extinction (), entry (E), departure (D), and detection (P) parameters for the Arizona Toad, as well as the model structures for the simple multi-season (SMS) model and models based only on historical sites. Covariates included in the table are historical vs. control sites (Hist/Con), lentic vs. lotic sites (Len/Lot), Julian date (Jul), percent moon cover (M), wind (Wi), air temperature (C, Celsius), year as a factor (YrF), and year as a trend (YrT). Staggered-entry models are ordered by delta AIC scores (ΔAIC), model weights (W), and –2 log likelihood values (LogLik).
Comparison of several staggered-entry occupancy models and the covariates used in each for occupancy (ψ), colonization (γ), extinction (), entry (E), departure (D), and detection (P) parameters for the Arizona Toad, as well as the model structures for the simple multi-season (SMS) model and models based only on historical sites. Covariates included in the table are historical vs. control sites (Hist/Con), lentic vs. lotic sites (Len/Lot), Julian date (Jul), percent moon cover (M), wind (Wi), air temperature (C, Celsius), year as a factor (YrF), and year as a trend (YrT). Staggered-entry models are ordered by delta AIC scores (ΔAIC), model weights (W), and –2 log likelihood values (LogLik).
Table 4.

Estimated detection probabilities (SE) for the Arizona Toad from the staggered-entry multi-season occupancy model and the simple multi-season occupancy model.

Estimated detection probabilities (SE) for the Arizona Toad from the staggered-entry multi-season occupancy model and the simple multi-season occupancy model.
Estimated detection probabilities (SE) for the Arizona Toad from the staggered-entry multi-season occupancy model and the simple multi-season occupancy model.
Fig. 2.

Estimated probabilities of Arizona Toads entering sites between subsequent surveys (orange), given that they have not previously entered, based on the estimated relationship between entry and Julian date. Estimated probabilities of Arizona Toad departing sites between subsequent surveys (blue), given that they are already present, based on the estimated relationship between departure and Julian date. Yellow points show the probabilities that toads were available to be detected for each day. Julian values range from 1 March to 30 May.

Fig. 2.

Estimated probabilities of Arizona Toads entering sites between subsequent surveys (orange), given that they have not previously entered, based on the estimated relationship between entry and Julian date. Estimated probabilities of Arizona Toad departing sites between subsequent surveys (blue), given that they are already present, based on the estimated relationship between departure and Julian date. Yellow points show the probabilities that toads were available to be detected for each day. Julian values range from 1 March to 30 May.

Table 5.

Mean probabilities of colonization, extinction, and occupancy (SE) from the staggered-entry multi-season occupancy model and the simple multi-season occupancy model. As per our best-supported model structures, we estimated separate rates of colonization for historical and control sites, but constant rates of extinction for all site types.

Mean probabilities of colonization, extinction, and occupancy (SE) from the staggered-entry multi-season occupancy model and the simple multi-season occupancy model. As per our best-supported model structures, we estimated separate rates of colonization for historical and control sites, but constant rates of extinction for all site types.
Mean probabilities of colonization, extinction, and occupancy (SE) from the staggered-entry multi-season occupancy model and the simple multi-season occupancy model. As per our best-supported model structures, we estimated separate rates of colonization for historical and control sites, but constant rates of extinction for all site types.
Fig. 3.

Estimated probabilities of yearly occupancy (y-axis; ± 95% CI) from the historically occupied (orange) and control (blue) sites for the (A) staggered-entry model (SE) and the (B) simple multi-season (SMS) model. Dotted error bars represent years when surveys were not conducted.

Fig. 3.

Estimated probabilities of yearly occupancy (y-axis; ± 95% CI) from the historically occupied (orange) and control (blue) sites for the (A) staggered-entry model (SE) and the (B) simple multi-season (SMS) model. Dotted error bars represent years when surveys were not conducted.

Simple multi-season model.—

When we re-fit the calling data to a simple multi-season occupancy model, yearly estimates of detection were generally lower than when estimated with the staggered-entry model, as expected (Table 4). Estimates from the simple multi-season model indicated the probability of detection declined through the season (b = –0.006, 95% CI = –0.014–0.0001; Fig. 4). All other parameter estimates from the simple multi-season model were similar to those from the staggered-entry model, resulting in nearly identical occupancy estimates over time (Table 5, Fig. 3B).

Fig. 4.

Estimated detection probability from the simple multi-season model fitted with calling data from historical and control sites, based on the relationship between detection and Julian date from the simple multi-season occupancy model. Julian values represent only the values present in our data set and range from 1 March to 30 May (Julian days 60–150).

Fig. 4.

Estimated detection probability from the simple multi-season model fitted with calling data from historical and control sites, based on the relationship between detection and Julian date from the simple multi-season occupancy model. Julian values represent only the values present in our data set and range from 1 March to 30 May (Julian days 60–150).

Historical site models.—

Re-fitting the top staggered-entry model with only data from 62 historical sites and information on years since the first record for a given site (11–57 yrs) or years since the most recent record for a given site (10–56 yrs) showed that both models provided a similar fit to the data (ΔAIC = 0.05; Table 3). The occupancy coefficient did not differ from 0 based on years since the first observation (b = 0.0006, 95% CI = –0.0045–0.0057) nor years since the most recent observation (b = 0.0006, 95% CI = –0.0045–0.0059; Supplemental Fig. 1; see Data Accessibility). Both historical site models estimated ∼21% of historically occupied sites were still occupied during 2013–2016 and 2019.

DISCUSSION

Surveys for Arizona Toads during 2013–2016 and 2019 covered approximately 85% of documented historical localities for the species in New Mexico and nearly its entire known range in the state. Based on surveys of 86 historical and 59 control sites, historical sites were approximately three times more likely to be occupied, and occupancy did not vary between lentic and lotic habitats. Analysis of the data using a staggered-entry framework also revealed variation in timing of when toads started (patch entry) or stopped calling (patch departure) and allowed us to identify optimal timing for future surveys. Our results also highlight that conclusions based on naïve results unadjusted for variation in detection probability would have greatly underestimated occupancy by the end of the study, especially for historically occupied sites.

Our main staggered-entry analysis suggested that estimated occupancy of historically occupied sites increased sharply from a low of ∼16% in 2013 to ∼76% in 2019, while occupancy of control sites increased from ∼6% to ∼30% during the same time (Fig. 2). Although this increase is surprising, the result was consistent across all sub-models. Based on modeled estimates and naïve results (Table 2), 2013 was the year that toads were rarest. Year 2013 was also the year with the most intensive per-site survey effort and highest probability of detection (staggered-entry model), lending credence to the low occupancy estimates. Importantly, there was also a trend of lower detection during the study (Table 4), which caused an increasing mismatch between naïve and estimated occupancy over time.

We suspect low initial occupancy in 2013 resulted from extended drought that continued through the winter of 2015, and two wildfires that burned portions of the study area in 2011 (Miller Fire, ∼359 km2) and 2012 (Whitewater-Baldy Fire, ∼1200 km2; Ryan and Latella, 2013; Whitney et al., 2016; Gido et al., 2019). Approximately 50% of surveyed historical toad sites were dry during spring 2013 due to extended drought that reduced stream flows in parts of the Mimbres River drainage and small streams in the western part of the study area (Ryan and Latella, 2013; Gido et al., 2019). Breeding dynamics of the closely related Arroyo Toad (A. californicus) in southern California, which has similar habitat requirements as the Arizona Toad, are also strongly linked with drought-related variation in water availability in small streams (Sweet and Sullivan, 2005; Miller et al., 2012).

In 2013, many of the inundated sites in the West and Middle Forks of the Gila River also had heavy post-fire flooding and channel scouring that eliminated the clear, slow-flowing pools needed for toad breeding and tadpole development (Ryan and Latella, 2013; Gido et al., 2019). Stream and debris flows commonly increase after large wildfires in the southwestern USA (DeGraff et al., 2015; Sedell et al., 2015), sometimes causing loss of important habitat and localized extirpations of amphibians and other aquatic fauna (Hossack and Pilliod, 2011; Whitney et al., 2016; Zylstra et al., 2019). The combined effects of drought and habitat changes from post-fire flooding and sedimentation could have caused extinctions in patches in 2013 that were later recolonized, similar to rapid recovery of other amphibians and stream fishes that have evolved in stochastic environments (Miller et al., 2012; Whitney et al., 2016; Hossack et al., 2017; Gido et al., 2019; Zylstra et al., 2019). This evidence for localized extinction resulting from environmental variability highlights the importance of long-term studies, especially when resampling historically occupied sites in stochastic environments. Had our study concluded in 2013 or 2014, patch dynamics may have confounded long-term trends in occupancy, and low rates of occupancy in those years would have indicated a decline in Arizona Toad occupancy.

The simple multi-season model did not incorporate the entry and departure parameters from the staggered-entry model, and temporal variation in calling activity was considered detection error. As a result, estimated probability of detection averaged 23% lower across years for the simple multi-season vs. staggered-entry models, and detection from the simple multi-season model varied greatly across years even based on the same model (e.g., 0.75 in 2013 vs. 0.23 in 2016; Table 4). However, the nearly perfect concordance of detections from live surveys and recordings of the same 3-minute surveys in 2019 suggests variation in detection during the study was from toad calling activity, not observer error. When planning surveys, the lower end of detection estimates provides a safer margin for ensuring reliable information compared with using the mean, especially if few repeat surveys are conducted within a season or if a study will not extend for several years. Alternatively, removal designs that stop surveys of a site once a target species has been detected might provide a compromise between over-sampling (years with high detection) and under-sampling (years with low detection; MacKenzie and Royle, 2005). Removal designs also allow allocation of effort to a larger number of sites, potentially increasing precision and better reflecting environmental and spatial variation.

Despite clear violation of the within-season closure assumption that is common to simple multi-season models—the assumption that staggered-entry models relax—the simpler model framework provided nearly identical occupancy and vital rate (colonization, extinction) estimates (Table 5). This similarity suggests if scientists are not explicitly interested in describing variation in phenology, the simpler model will suffice. Still, the staggered-entry approach provided useful insight into the variable calling phenology of Arizona Toads and optimal timing to detect them. Based on surveys during 2013–2016 and 2019, the derived probability of patch availability was highest during approximately 5 March to 20 March, which corresponded to the early portion of sampling efforts during most years (Fig. 2). Sampling during this optimal date range would have maximized the likelihood of encountering toads and decreased the likelihood of recording false absences.

Broad-scale surveys focused on historically occupied sites have provided some of the most compelling evidence of severe amphibian declines, including in the western USA (Clarkson and Rorabaugh, 1989; Corn et al., 1989; Drost and Fellers, 1996; Fisher and Shaffer, 1996). For example, during 1983–1987, the (now) federally threatened Chiricahua Leopard Frog (Rana chiricahuensis) was found at only two of 36 sites in Arizona where it was documented in the 1960s and 1970s (Clarkson and Rorabaugh, 1989). We expected re-fitting the top-ranked staggered-entry model with information on years since the first observation (11–57 yrs) and years since the most recent observation (10–56 yrs) of Arizona Toads at historically occupied sites would result in a steep extinction curve, with older sites the least likely to be occupied. However, occupancy was unrelated to time since the first or last record, indicating a site where toads were last documented 50 years ago was just as likely to be occupied as a site where toads were last documented 15 years ago (Supplemental Fig. 1; see Data Accessibility). Imprecise location information for some museum specimens, and a lack of information on whether some historical sites hosted breeding activity historically or simply had a vagrant adult, adds uncertainty to estimates for historically occupied sites. Nevertheless, the near-zero relationship between time and occupancy for these sites suggests a lack of widespread, recent population declines in the system. The lack of a relationship between age of record and extinction is notable because surveys of historically occupied sites are biased towards documenting extinction. Whether by chance or from changes in habitat or environmental conditions, populations are prone to extinction (Skelly et al., 2003; Fournier et al., 2019).

Overall, our results indicate calling activity (and thus detectability) of Arizona Toads in southwestern New Mexico is highly variable within and among years, but our results do not indicate large declines in patch occupancy in recent decades. These occupancy trends were reinforced by multiple model types and substantiated by a flat time-occupancy curve at historically occupied sites. Further, large increases in occupancy of historical and control sites since 2013, the year following an extended drought and post-fire habitat changes that likely made many breeding patches temporarily unsuitable, suggest that local breeding sites can be re-occupied by adults. Yet the large differences in detection probabilities among years also suggest populations may be small, and we lack information on whether most breeding sites reliably produce tadpoles that survive to metamorphosis. Importantly, we also still lack information on the status and trends of populations over most of the species' range.

DATA ACCESSIBILITY

Supplemental material is available at https://www.ichthyologyandherpetology.org/h2020133. Unless an alternative copyright or statement noting that a figure is reprinted from a previous source is noted in a figure caption, the published images and illustrations in this article are licensed by the American Society of Ichthyologists and Herpetologists for use if the use includes a citation to the original source (American Society of Ichthyologists and Herpetologists, the DOI of the Ichthyology & Herpetology article, and any individual image credits listed in the figure caption) in accordance with the Creative Commons Attribution CC BY License.

ACKNOWLEDGMENTS

We thank T. Chambert for comments that improved the manuscript. This work was conducted under New Mexico Department of Game and Fish Scientific Collection Authorization #3329 and University of New Mexico Institutional Animal Care and Use Committee Protocol 13-100983-MC. This research was funded by the New Mexico Department of Game and Fish State Wildlife Program Grants T-32-4 and T-32-P-3, and the U.S. Fish and Wildlife Service State Wildlife Grant T-32-P3, 18, and the U.S. Geological Survey Amphibian Research and Monitoring Initiative (ARMI). This manuscript is ARMI product no. 784. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

LITERATURE CITED

Arnold,
T. W.
2010
.
Uninformative parameters and model selection using Akaike's Information Criterion
.
The Journal of Wildlife Management
74
:
1175
1178
.
AZGFD (Arizona Game and Fish Department).
2015
.
Species of greatest conservation need in Arizona. State Wildlife Action Plan
.
Arizona Game and Fish Department
,
Phoenix, Arizona
.
Blair,
A. P.
1955
.
Distribution, variation, and hybridization in a relict toad (Bufo microscaphus) in southwestern Utah
.
American Museum Novitates
1722
:
1
38
.
Bolker,
B.
2020
.
emdbook: ecological models and data in R. R package version 1.3.12.
Bradford,
D. F.,
Jaeger,
J. R.
and
Shanahan.
S. A.
2005
.
Distributional changes and population status of amphibians in the eastern Mojave Desert
.
Western North American Naturalist
65
:
462
472
.
Chambert,
T.,
Kendall,
W. L.
Hines,
J. E.
Nichols,
J. D.
Pedrini,
P.
Waddle,
J. H.
Tavecchia,
G.
Walls,
S. C.
and
Tenan.
S.
2015
.
Testing hypotheses on distribution shifts and changes in phenology of imperfectly detectable species
.
Methods in Ecology and Evolution
6
:
638
647
.
Clarkson,
R. W.,
and
Rorabaugh.
J. C.
1989
.
Status of leopard frogs (Rana pipiens complex: Ranidae) in Arizona and southeastern California
.
Southwestern Naturalist
34
:
531
538
.
Corn,
P. S.,
Muths,
E.
Kissel,
A. M.
and
Scherer.
R. D.
2011
.
Breeding chorus indices are weakly related to estimated abundance of boreal chorus frogs
.
Copeia
2011
:
365
371
.
Corn,
P. S.,
Stolzenburg,
W.
and
Bury.
R. B.
1989
.
Acid precipitation studies in Colorado and Wyoming: interim report of surveys of montane amphibians and water chemistry
.
U.S. Fish and Wildlife Service Biological Report 80 (40.26).
Degenhardt,
W. G.,
Painter,
C. W.
and
Price.
A. H.
1996
.
Amphibians and Reptiles of New Mexico
.
University of New Mexico Press
,
Albuquerque, New Mexico
.
DeGraff,
J. V.,
Cannon,
S.
and
Gartner.
J. E.
2015
.
Timing of susceptibility to post-fire debris flows in the western USA
.
Environmental & Engineering Geoscience
250
:
173
186
.
Doherty,
P.,
White,
G.
and
Burnham.
K.
2012
.
Comparison of model building and selection strategies
.
Journal of Ornithology
152
:
317
323
.
Drost,
C. A.,
and
Fellers.
G. M.
1996
.
Collapse of a regional frog fauna in the Yosemite area of the California Sierra Nevada, USA
.
Conservation Biology
10
:
414
425
.
Fisher,
R. N.,
and
Shaffer.
H. B.
1996
.
The decline of amphibians in California's great Central Valley
.
Conservation Biology
10
:
1387
1397
.
Fournier,
A. M. V.,
White,
E. R.
and
Heard.
S. B.
2019
.
Site-selection bias and apparent population declines in long-term studies
.
Conservation Biology
33
:
1370
1379
.
Gido,
K. B.,
Propst,
D. L.
Whitney,
J. E.
Hedden,
S. C.
Turner,
T. F.
and
Pilger.
T. J.
2019
.
Pockets of resistance: response of arid-land fish communities to climate, hydrology, and wildfire
.
Freshwater Biology
64
:
761
777
.
Hossack,
B. R.,
Honeycutt,
R. K.
Sigafus,
B. H.
Muths,
E.
Crawford,
C. L.
Jones,
T. R.
Sorensen,
J. A.
Rorabaugh,
J. C.
and
Chambert.
T.
2017
.
Informing recovery in a human-transformed landscape: drought-mediated coexistence alters population trends of an imperiled salamander and invasive predators
.
Biological Conservation
209
:
377
394
.
Hossack,
B. R.,
and
Pilliod.
D. S.
2011
.
Amphibian responses to wildfire in the West: emerging patterns from short-term studies
.
Fire Ecology
7
:
129
144
.
Jackson,
C. H.
2011
.
Multi-state models for panel data: the msm package for R
.
Journal of Statistical Software
38
:
1
29
.
Jennings,
R. D.,
Christman,
B. L.
Painter,
C. W.
and
Frye.
J. L.
2010
.
Near-stream herpetofauna of the upper Gila River of southwestern New Mexico
.
The New Mexican Botanist (Special Issue)
2
:
88
115
.
Johnson,
D. H.,
and
Batie.
R. D.
2001
.
Surveys of calling amphibians in North Dakota
.
Prairie Naturalist
33
:
227
247
.
Kendall,
W. L.,
Hines,
J. E.
Nichols,
J. D.
and
Grant.
E. H. C.
2013
.
Relaxing the closure assumption in occupancy models: staggered arrival and departure times
.
Ecology
94
:
610
617
.
MacKenzie,
D. I.,
Nichols,
J. D.
Royle,
J. A.
Pollock,
K. H.
Bailey,
L. L.
and
Hines.
J. E.
2006
.
Occupancy Estimation and Modeling: Inferring Patterns and Dynamics of Species Occurrence
.
Academic Press
,
San Francisco
.
MacKenzie,
D. I.,
and
Royle.
J. A.
2005
.
Designing occupancy studies: general advice and allocating survey effort
.
Journal of Applied Ecology
42
:
1105
1114
.
Meineke,
E. K.,
Davies,
T. J.
Daru,
B. H.
and
Davis.
C. C.
2019
.
Biological collections for understanding biodiversity in the Anthropocene
.
Philosophical Transactions of the Royal Society B: Biological Sciences
374
:
20170386
.
Miller,
D. A. W.,
Brehme,
C. S.
Hines,
J. E.
Nichols,
J. D.
and
Fisher.
R. N.
2012
.
Joint estimation of habitat dynamics and species interactions: disturbance reduces co-occurrence of non-native predators with an endangered toad
.
Journal of Animal Ecology
81
:
1288
1297
.
NMDGF (New Mexico Department of Game and Fish).
2016
.
Species of greatest conservation need in New Mexico. State Wildlife Action Plan
.
New Mexico Department of Game and Fish
,
Santa Fe, New Mexico
.
NVDW (Nevada Department of Wildlife).
2015
.
Species of greatest conservation need in Nevada. State Wildlife Action Plan
.
Nevada Department of Wildlife
,
Reno, Nevada
.
R Core Team.
2020
.
R: a language and environment for statistical computing
.
R Foundation for Statistical Computing
,
Vienna, Austria
.
Ryan,
M. J.,
Giermakowski,
J. T.
Latella,
I. M.
and
Snell.
H. L.
2017
.
No evidence of hybridization between the Arizona Toad (Anaxyrus microscaphus) and Woodhouse's Toad (A. woodhousii) in New Mexico, USA
.
Herpetological Conservation and Biology
12
:
565
575
.
Ryan,
M. J.,
and
Latella.
I. M.
2013
.
Current status of the Arizona Toad (Anaxyrus microscaphus) in New Mexico: identification and evaluation of potential threats o its persistence
.
New Mexico Department of Game and Fish
.
Ryan,
M. J.,
Latella,
I. M.
Giermakowski,
J. T.
and
Snell.
H. L.
2015
.
Final Report: Status of the Arizona Toad (Anaxyrus microscaphus) in New Mexico
.
New Mexico Department of Game and Fish
.
Schwaner,
T. D.,
and
Sullivan.
B. K.
2009
.
Fifty years of hybridization: introgression between the Arizona Toad (Bufo microscaphus) and Woodhouse's Toad (B. woodhousii) along Beaver Dam Wash in Utah
.
Herpetological Conservation and Biology
4
:
198
206
.
Seber,
G. A. F.
1982
.
The Estimation of Animal Abundance and Related Parameters
.
MacMillan
,
New York
.
Sedell,
E. R.,
Gresswell,
R. E.
and
McMahon.
T. E.
2015
.
Predicting spatial distribution of postfire debris flows and potential consequences for native trout in headwater streams
.
Freshwater Science
34
:
1558
1570
.
Shaffer,
H. B.,
Fisher,
R. N.
and
Davidson.
C.
1998
.
The role of natural history collections in documenting species declines
.
Trends in Ecology & Evolution
13
:
27
30
.
Skelly,
D. K.,
Yurewicz,
K. L.
Werner,
E. E.
and
Relyea.
R. A.
2003
.
Estimating decline and distributional change in amphibians
.
Conservation Biology
17
:
744
751
.
Sullivan,
B. K.
1992
.
Calling behavior of the Southwestern Toad (Bufo microscaphus)
.
Herpetologica
43
:
383
389
.
Sullivan,
B. K.
1993
.
Distribution of the Southwestern Toad (Bufo microscaphus) in Arizona
.
The Great Basin Naturalist
53
:
402
406
.
Sweet,
S. S.,
and
Sullivan.
B. K.
2005
.
Bufo californicus Camp, 1915 Arroyo Toad
,
p.
396
400
.
In:
Amphibian Declines: The Conservation Status of United States Species
.
Lannoo
M. J.
(ed.).
University of California Press
,
Berkeley, California
.
Tingley,
M. W.,
and
Beissinger.
S. R.
2009
.
Detecting range shifts from historical species occurrences: new perspectives on old data
.
Trends in Ecology & Evolution
24
:
625
633
.
UDWR (Utah Division of Wildlife Resources).
2015
.
Species of greatest conservation need in Utah. State Wildlife Action Plan
.
Utah Division of Wildlife Resources
,
Salt Lake City, Utah
.
USFWS (U.S. Fish and Wildlife Service).
2015
.
90-day finding on a petition to list the Arizona Toad under the Endangered Species Act
.
Federal Register
80
(126)
:
37568
37579
.
Walls,
S. C.,
Waddle,
J. H.
and
Dorazio.
R. M.
2011
.
Estimating occupancy dynamics in an anuran assemblage from Louisiana, USA
.
Journal of Wildlife Management
75
:
751
761
.
Weir,
L. A.,
Royle,
J. A.
Gazenski,
K. D.
and
Villena.
O.
2014
.
Northeast regional and state trends in anuran occupancy from calling survey data (2001–2011) from the North Amphibian Monitoring Program
.
Herpetological Conservation and Biology
9
:
223
245
.
Weiser,
E. L.,
Diffendorfer,
J. E.
Lopez-Hoffman,
L.
Semmens,
D.
and
Thogmartin.
W. E.
2020
.
Challenges for leveraging citizen science to support statistically robust monitoring programs
.
Biological Conservation
242
:
108411
.
Whitney,
J. E.,
Gido,
K. B.
Pilger,
T. J.
Propst,
D. L.
and
Turner.
T. F.
2016
.
Metapopulation analysis indicates native and non-native fishes respond differently to effects of wildfire on desert streams
.
Ecology of Freshwater Fish
25
:
376
392
.
Wooten,
J. A.,
Sullivan,
B. K.
Klooster,
M. R.
Schwaner,
T. D.
Sullivan,
K. O.
Brown,
A. D.
Takahashi,
M.
and
Bradford.
P. R.
2019
.
Thirty years of hybridization between toads along the Agua Fria River in Arizona: part II: fine-scale assessment of genetic changes over time using microsatellites
.
Journal of Herpetology
53
:
104
114
.
Zylstra,
E. R.,
Swann,
D. E.
Hossack,
B. R.
Muths,
E.
and
Steidl.
R. J.
2019
.
Drought-mediated extinction of an arid-land amphibian: insights from a spatially explicit dynamic occupancy model
.
Ecological Applications
29
:
e01859
.

Author notes

1

U.S. Geological Survey, Northern Rocky Mountain Science Center, Missoula, Montana 59801; Email: (MJF) mikejforzley@gmail.com; and (BRH) blake_hossack@usgs.gov. Send reprint requests to MJF.

2

Wildlife Biology Program, W.A. Franke College of Forestry and Conservation, University of Montana, Missoula, Montana 59812.

3

Arizona Game and Fish Department, 5000 W Carefree Hwy, Phoenix, Arizona 85086; Email: MRyan@azgfd.gov.

4

Arizona Game and Fish Department, 9140 E 28th St., Yuma, Arizona 85365; Email: ilatella@azgfd.gov.

5

Division of Amphibians and Reptiles, Museum of Southwestern Biology, Department of Biology, MSC03 2020, University of New Mexico, Albuquerque, New Mexico 87131; Email: tomas@unm.edu.

6

U.S. Geological Survey, Fort Collins Science Center, Fort Collins, Colorado 80526; Email: muthse@usgs.gov.

7

U.S. Geological Survey, Southwest Biological Science Center, 520 N. Park Avenue, Tucson, Arizona 85719; Email: bsigafus@usgs.gov.

Associate Editor: M. J. Lannoo.