Hot, wet and rare: modelling the occupancy dynamics of the narrowly distributed Dixie Valley toad

Keywords: Anaxyrus williamsi, aquatic conservation, Bayesian hierarchical model, desert springs, endangered species, geothermal energy, Great Basin, multi-state occupancy model.

The size of a species’ range is negatively correlated to its risk of extinction based on fossil records (Payne and Finnegan 2007) and extant species (Newsome et al. 2020). For anurans in particular, range size is the most important predictor of extinction risk, even after excluding species classified as at risk by the International Union for Conservation of Nature only on the basis of the range-size criteria (Cooper et al. 2008). Narrowly distributed species are also likely to face the perils associated with a small population size (Caughley 1994). A small range and lack of redundancy provided by multiple, widely distributed populations mean that any disturbance to the occupied habitat could threaten the entire species with extinction.

Many narrowly distributed endemic species occur in desert springs, which are inherently isolated patches of aquatic habitat in an otherwise xeric landscape (Shepard 1993). These habitats are fragile, given the inhospitable conditions of the surrounding desert, and desert spring communities are under further threat from climate change, introduced species, and groundwater extraction and surface-water diversion to support human populations (Kodric-Brown and Brown 2007; Unmack and Minckley 2008). Disturbances to desert springs have led to the extinction of endemic fish (Contreras-Balderas and Lozano-Vilano 1996) and springsnails (Hershler et al. 2014). Species inhabiting desert springs can be isolated from the nearest suitable habitat by several kilometres, severing connectivity among populations and leading to genetic divergence (Murphy et al. 2012; Seidel et al. 2009). The entire distribution of an isolated species could decrease rapidly in response to changes in the aquatic habitat at a few springs.

Modelling how species occurrence depends on habitat dynamics is crucial for developing management plans to ensure the long-term viability of narrowly distributed species. Dynamic occupancy models can disentangle the observation and occupancy processes and quantify how key environmental factors influence occupancy, colonisation, and extirpation (MacKenzie et al. 2003). Multi-state occupancy models can provide further insight into the dynamics of narrowly distributed species because they enable researchers to estimate how environmental conditions affect different life stages (Mackenzie et al. 2009; Duarte et al. 2020). Knowing the habitat requirements of different life stages is important for conserving species with complex life histories, such as amphibians with a terrestrial adult phase and aquatic larval phase (Semlitsch 2000).

In this study, we used dynamic occupancy models to quantify how spatio-temporal environmental variation influences occurrence of a narrowly distributed species. We surveyed for the recently described Dixie Valley toad (Anaxyrus williamsi) in the desert spring habitat in north-western Nevada, USA, that makes up this species’ entire range (Gordon et al. 2017). Although adults of many toad species spend much of their time outside the breeding season in terrestrial habitat far from water (Bartelt et al. 2004; Forester et al. 2006; Liang 2013), radio-telemetry during the spring and fall has shown that Dixie Valley toads closely associate with water in both seasons (Halstead et al. 2021). Given the close link between Dixie Valley toad behaviour and water, we modelled how the occurrence of adult toads and their larvae were related to water temperature, availability, and depth and estimated how the area occupied by Dixie Valley toads varied over time. Our results have implications for the conservation of this narrowly distributed desert endemic, and show how dynamic occupancy models can be used to link changes in species occurrence or habitat use to local environmental conditions.

Surveys took place within desert spring habitat in the Dixie Valley in Churchill County, Nevada, USA (Fig. 1), from 2018 to 2021. The desert spring habitat within the range of the Dixie Valley toad consists of shallow wetlands covering a total area of approximately 1.5 km², surrounded by alkali desert scrub, at an elevation of approximately 1040 m. The wetland habitat is fed by >120 springs and seeps of varying discharge rates and temperatures (U.S. Fish and Wildlife Service 2022a). Shallow marshes dominated by bulrushes (Schoenoplectus spp.) are found near the outflows of the larger springs, and seasonally wet meadows dominated by rushes (Juncaceae), salt grass (Distichlis spicata), and composite forbs (Asteraceae) are found farther from the springs. The wetland habitat is surrounded by desert scrub dominated by sagebrush (Artemisia tridentata), saltbush (Atriplex spp.), and greasewood (Sarcobatus vermiculatus; Bureau of Land Management 2011). The Dixie Valley toad was emergency listed as endangered under the US Endangered Species Act in April 2022 (US Fish and Wildlife Service 2022b). In addition to its narrow range, concern over geothermal energy development in the Dixie Valley hot springs and its potential disturbance to the spring-fed wetlands was cited as one reason the species was listed (US Fish and Wildlife Service 2022b).

Fig. 1.  The location of the study area within Churchill County, Nevada, USA. Squares represent grid cells that were sampled for Dixie Valley toads (Anaxyrus williamsi) and have been enlarged slightly to improve visibility. The red star in the inset panel represents the location of the study area in the main panel.

We divided the desert spring habitat within the range of the Dixie Valley toad (Fig. 1) into 20 m × 20 m grid cells and used the ‘spsurvey’ package (Kincaid et al. 2016) in R (R Core Team 2022) to randomly select 60 cells by using a generalized random-tessellation stratified (GRTS; Stevens and Olsen 2004) design that ensured spatial dispersion of selected cells (Fig. 1). We used a stratified sampling design to select cells proportionally on the basis of the area of the spring habitat polygon, with a minimum of three cells per polygon to ensure coverage of each (Fig. 1). We chose 20 m × 20 m as the cell size because this area was small enough to be feasible for repeated surveys within an approximately 1-week primary period, yet large enough that it was reasonable to assume closure (no change in occupancy state) within the primary period (Halstead et al. 2019). Grid cells were surveyed for Dixie Valley toads during the following six primary periods: April 2018, May 2018, May 2019, September 2019, May 2020, and May 2021. Cells were surveyed from one to five times (secondary occasions) within each primary period; starting in 2019, all cells were surveyed on at least two secondary occasions within each primary period (Supplementary Table S1). We surveyed all cells each primary period, even if the surface was dry. Each survey consisted of one or occasionally two observers walking the area of the cell visually surveying for Dixie Valley toads and using a dipnet to capture amphibian adults and larvae for identification when necessary. Cells were surveyed for 5–21 min (mean = 14.4 min, s.d. = 4.1) total (number of observers × number of minutes surveyed per observer) on each secondary occasion. Most surveys took place during the day (n = 582), but we also surveyed cells at night (n = 140) in each primary period except May 2019.

During each survey, we recorded characteristics we expected to influence detection probability of adult and larval Dixie Valley toads, including the length of time spent surveying (survey time, in person-minutes to account for multiple observers on some occasions), time of day, day of year, % cover of emergent vegetation in the cell, and air temperature. We defined emergent vegetation as plants that are normally submerged at ground level with stems that emerge from the water’s surface. Emergent vegetation was primarily composed of bulrushes (Schoenoplectus spp.), cattails (Typha sp.), grasses, and rushes. We also measured and recorded characteristics of the cell that we expected to influence the probability that a cell would be occupied by Dixie Valley toads, including the percentage of the surface of the cell that was covered with water (% wetted surface; wet), the average water depth in the cell (based on a visual estimate by the surveyor and one measurement of the maximum depth in cm; dep), % cover of emergent vegetation (ev), and water temperature (wt). Because a radio-telemetry study of Dixie Valley toads had found that toads often chose brumation sites near springheads (Halstead et al. 2021) and our field observations suggested that Dixie Valley toad larvae were more often observed farther from springheads, we calculated the distance of the centroid of each surveyed cell from the nearest springhead using ArcGIS ver. 10.8.1 (Environmental Systems Research Institute 2020). We then used the distance to the nearest springhead (spdist) as a covariate on all occupancy parameters.

We analysed our survey data by using a multi-state dynamic occupancy model adapted from earlier models for amphibian occupancy (Halstead et al. 2019; Duarte et al. 2020). The three possible occupancy states in the model were as follows: (1) occupied by adult or subadult (post-metamorphosis; hereinafter adult) Dixie Valley toads, (2) occupied by larval Dixie Valley toads and adults (larval occupancy is conditional on adult occupancy), and (3) not occupied by Dixie Valley toads. Each grid cell was treated as a ‘site,’ and the probability a cell was occupied during the first survey period is ψ₁ (initial adult occupancy). Cells can transition between occupancy states on the basis of the probability of colonisation (γ), site survival (ϕ), and larval occupancy (r; a measure of reproduction). Colonisation is the transition from the unoccupied state (3) to occupied by adults (with or without larvae), whereas site survival is the probability of remaining occupied by adults, given a cell was occupied during the previous primary period. The transition probabilities between each occupancy state at Time t and each occupancy state at Time t + 1 are defined for each cell (i) by the state transition matrix below.

For the occupancy component of the model, we tested the importance of five environmental covariates that we hypothesised could influence the probability of larval occupancy (r), initial occupancy (ψ₁), colonisation (γ), and site survival (ϕ) on the basis of our experience with Dixie Valley toads in the field (Eqn 1). In Eqn 1, y_i,t represents any of the four occupancy parameters for Cell i in primary Period t.

The three possible detection states in the model were (1) adult/subadult Dixie Valley toad detected, (2) larval Dixie Valley toad detected (with or without adults/subadults), and (3) no Dixie Valley toad detected. There were two distinct parameters for detection probability of adult/subadult Dixie Valley toad (p_A) and larval Dixie Valley toad (p_L). We assumed that no false positives occurred during our surveys (e.g. if the true state of a cell was unoccupied, it was not possible to observe State 1 or State 2). The probability of a survey resulting in each of the three detection states, given the true occupancy state of a cell (i) in Primary period t, is defined by the matrix below.

We included the following survey conditions as covariates on detection probability for adults and larvae: day of year (d; linear and quadratic), survey time (stime), time of day (tod; circular, in radians), emergent vegetation cover (ev), and air temperature (at) at the time of the survey (Eqn 2, p_A shown).

We first fit a full model that included all covariates on detection and occupancy parameters, with no variable selection to estimate the relationships between detection and occupancy probability and each covariate, and to obtain posterior distributions for use as pseudo-prior distributions for covariate effects in our variable selection procedure. Next, we used Gibbs variable selection (Carlin and Chib 1995; Dellaportas et al. 2002) to let the model and data dictate whether covariates were included in the model or not. The coefficient for each Covariate j (β_j) was multiplied by an indicator variable () that allowed each covariate to be ‘turned on’ ( = 1) or ‘turned off’ ( = 0) by the model (Eqn 2). The prior probability that each variable was included was 0.5 for linear terms, and 0.25 for quadratic terms (linear effects had to be included for quadratic effects to be included), and we calculated the posterior probability of inclusion on the basis of the proportion of Markov-chain Monte Carlo (MCMC) samples in which the indicator variable was set to 1 for that covariate (Tables 1, 2). We used a joint prior distribution for the β coefficients, with separate priors for when a covariate is turned on or turned off (Eqn 3), where N(μ,σ) represents a normal distribution with mean = μ and standard deviation = σ. In Eqn 3, μ_tune,j and σ_tune,j are the mean and standard deviation respectively, of the posterior distribution for β_j from the full model, and σ₁ = 1.64 is the default standard deviation to create a diffuse prior when the covariate is turned on.

Table 1.  Summary of Dixie Valley toad (Anaxyrus williamsi) detections by primary period.

Table 2.  The posterior probability of inclusion of each detection (survey) covariate on adult (p_A) and larval (p_L) detection probability.

Using this joint prior ensures that the value for β_j stays within the range of reasonable values when the covariate is turned off thanks to the informative prior, whereas a broad, uninformative prior is used when the covariate is turned on, and the informative prior does not influence the posterior for β_j (Hooten and Hobbs 2015).

Including the indicator variables () for each covariate, the linear models for occupancy parameters (Eqn 4) and the probability of detection of adults in Cell i during Survey j in Year t then become (Eqn 5):

The same covariates were included in the linear predictor for the probability of detecting adult () and larval Dixie Valley toads in Cell i during Survey j in Primary period t (). We fixed p_L and r to zero for the September 2019 survey period, because on the basis of the biology of the species, we did not expect larvae to be present in the late summer and fall. No larvae were observed during surveys in September 2019.

For four environmental covariates (emergent vegetation cover, water temperature, percentage wetted surface, and mean water depth), we included a submodel such that our estimates from each individual survey are measured with error around some ‘true’ unobserved mean value of that characteristic for that cell during that primary survey period (e.g. May 2018). This helps account for the fact that some cells have noisy estimates of environmental covariates for multiple surveys within the same month. The unobserved ‘true’ value for each covariate is then included as a covariate on the detection and occupancy parameters. We did not record the cover of emergent vegetation in surveyed cells in the September 2019 and May 2020 primary periods. To impute missing emergent vegetation data for these two primary periods, we set an informative prior for the true emergent vegetation cover in each primary period with a mean equal to the mean of all emergent vegetation measurements for that cell, and a standard deviation equal to the standard deviation of all emergent vegetation measurements for that cell (Eqn 6).

For water temperature, % wetted surface, and mean water depth, the prior for the true value in each cell in each primary period was based on the mean and standard deviation of field measurements of that covariate from that primary period (Eqn 7, only water temperature shown).

On the basis of the results of the Gibbs variable selection, we then fit a ‘final’ model that included only parameters for which the posterior probability of inclusion was greater the prior probability of inclusion (0.5 for linear effects, 0.25 for quadratic effects), following Duarte et al. (2020). For the final model, the priors for covariate effects on detection and occupancy were N(0,1.64). We based our inference about the number of cells occupied by each life stage in each primary period, detection probabilities, and occupancy parameters on this final model.

We fit occupancy models using JAGS software ver. 4.3.0 (Plummer 2003) implemented through R ver. 4.1.3 (R Core Team 2022), by using the ‘runjags’ package ver. 2.2.0-3 (Denwood 2016). We centred and standardised all covariates (except time of day) to have a mean of 0 and standard deviation of 1 to improve the performance of the MCMC sampler. We ran models on four independent chains, with the first 10 000 iterations discarded as burn-in, followed by 25 000 sampling iterations that were thinned by a factor of 10, resulting in a final posterior sample of 10 000. We inspected trace plots for evidence of mixing and convergence of chains. We also calculated the potential scale reduction factor () (Brooks and Gelman 1998) as a measure of the convergence of model parameters; all parameters had values of ≤1.01 and showed no evidence of failure to converge. When summarising posterior distributions, we present medians and 95% equal-tailed credible intervals (CI). Data are available at Rose et al. (2022a; https://doi.org/10.5066/P9QCIC87) and model code to reproduce occupancy analysis can be accessed at Rose et al. (2022b; https://doi.org/10.5066/P97DSXJM).

The number of cells in which only adult Dixie Valley toads were detected varied among primary periods from a minimum of 12 in April 2018 to a maximum of 31 in May 2018 (Table 1). No larval Dixie Valley toads were observed in September 2019, but larvae were detected in at least two surveyed cells in all five spring primary periods (Table 1). Of 60 cells surveyed, Dixie Valley toads of any life stage were detected in 48 cells over the entire study period. Of these 48 cells, adult toads were detected within 45 cells, and larvae were detected within 12 cells at some point in the study (Table S2, Supplementary Figs S1–S6).

The detection probability of larval Dixie Valley toads (p_L) was more uncertain and more variable among primary periods than for adults (p_A). The highest estimate of p_L was in May 2020 (median = 0.88, 95% CI = 0.41–0.99) and the lowest was in May 2018 (0.18, 0.02–0.75; Fig. S7). The estimated detection probability for adult Dixie Valley toads (p_A) was more consistent among primary periods, with a median varying from 0.32 to 0.57 (Fig. S7). The variable selection results showed support for the inclusion of effects of survey time, time of day, air temperature, and emergent vegetation cover on p_A (Table 2). Adult detection probability (p_A) was positively related to survey time, air temperature, and emergent vegetation cover, and exhibited a u-shaped relationship to time of day; adult detection probability was highest in the morning and evening (Fig. 2). The variable selection results showed support for including an effect of day of year, survey time, and time of day on p_L. As for adults, larval detection probability (p_L) was also positively related to survey time (Fig. 3). Larval detection probability was negatively related to day of year; larvae had a higher probability of being detected in surveys conducted earlier in the year. Larval detection probability exhibited a quadratic relationship with time of day; detection probability was high in the morning and early afternoon, before decreasing in the late afternoon and being lowest at night.

Fig. 2.  Relationships between detection probability for adult Dixie Valley toads (p_A) and survey covariates, with posterior probability of inclusion >prior probability of inclusion. Solid lines are median predictions from the posterior distribution, dashed lines are 95% equal-tailed credible intervals, and grey shading represents the posterior probability density. Vertical tick marks above the x-axis represent observed values for survey covariates, jittered to improve visibility of the density of observed values.

Fig. 3.  Relationships between detection probability for larval Dixie Valley toads (p_L) and survey covariates, with posterior probability of inclusion >prior probability of inclusion. Solid lines are median predictions from the posterior distribution, dashed lines are 95% equal-tailed credible intervals, and grey shading represents the posterior probability density. Vertical tick marks above the x-axis represent observed values for survey covariates, jittered to improve visibility of the density of observed values.

The number of surveyed cells estimated to be occupied by Dixie Valley toads of any life stage was lowest in September 2019 (median = 30, 95% CI = 21–42) and highest in May 2018 (45, 41–51) and May 2020 (39, 33–47; Fig. S8). The number of surveyed cells occupied by larval Dixie Valley toads was similar in April and May 2018, but dropped in May 2019, May 2020, and May 2021 (Fig. S8).

Estimates of ψ were largely stable during our study period, with the highest estimate for May 2018 and lowest for April 2018 and September 2019 (Fig. S9). Estimates of r in an average cell were low during each primary period, whereas the probability of an occupied cell remaining occupied (ϕ) was generally high throughout the study (Fig. S9). The probability of a cell being colonised (γ) varied greatly between primary periods; γ was high for May 2018 and May 2020, and low for May 2019, September 2019, and May 2021 (Fig. S9). The variable selection indicated that ψ₁ was influenced by % wetted surface and mean water depth. There was a positive relationship between ψ₁ and the mean water depth in the cell; the probability of initial occupancy increased rapidly with mean water depth for cells up to approximately 20 cm deep (Fig. 4). There was support for a positive relationship between ψ₁ and the area of surface water in a cell (Tables 3, 4, S3; Fig. 4).

Fig. 4.  Relationships between probability of occupancy during initial primary period (ψ₁) for Dixie Valley toads (Anaxyrus williamsi) and environmental covariates, with posterior probability of inclusion >prior probability of inclusion. Solid lines are median predictions from the posterior distribution, dashed lines are 95% equal-tailed credible intervals, and grey shading represents the posterior probability density. Vertical tick marks above the x-axis represent observed values for environmental covariates, jittered to improve visibility of the density of observed values.

Table 3.  The posterior probability of inclusion of each occupancy (cell) covariate on initial adult occupancy (ψ₁), site survival (ϕ), site colonisation (γ), and larval occupancy (r).

Table 4.  Prior distributions and summary statistics for posterior distributions of parameters in final dynamic occupancy model.

The variable selection showed support for relationships between larval occupancy (r) and % wetted surface, emergent vegetation cover, water temperature, and distance to nearest spring (Fig. 5, Tables 3, S3). Larval occupancy was positively related to % wetted surface of the cell, and negatively related to % emergent vegetation cover. Larval occupancy exhibited a quadratic relationship with water temperature, with the highest probability of larval occupancy at water temperatures near 20–28°C (Fig. 5). As distance from the nearest spring increased, r also increased, indicating larvae were more likely to occur in cells farther from spring heads (Fig. 5).

Fig. 5.  Relationships between probability of larval occupancy (r) for Dixie Valley toads (Anaxyrus williamsi) and environmental covariates, with posterior probability of inclusion >prior probability of inclusion. Solid lines are median predictions from the posterior distribution, dashed lines are 95% equal-tailed credible intervals, and grey shading represents the posterior probability density. Vertical tick marks above the x-axis represent observed values for environmental covariates, jittered to improve visibility of the density of observed values.

Two covariates influenced colonisation of previously unoccupied cells by Dixie Valley toads on the basis of the variable selection. The probability of a cell being colonised (γ) was positively related to % wetted surface and emergent vegetation cover (Tables 3, 4, S3; Fig. 6). The variable selection indicated support for water temperature influencing the probability of an occupied cell remaining occupied (ϕ) (Table 3). As water temperature increased from 10°C to ~25°C, ϕ increased; the relationship between ϕ and water temperature became much more uncertain for temperatures above 25°C and below 10°C. There was equivocal support for the inclusion of an effect of emergent vegetation cover on ϕ and the positive relationship was highly uncertain (Table 3, Fig. 7).

Fig. 6.  Relationships between probability of colonisation (γ) for Dixie Valley toads (Anaxyrus williamsi) and environmental covariates with posterior probability of inclusion >prior probability of inclusion. Solid lines are median predictions from the posterior distribution, dashed lines are 95% equal-tailed credible intervals, and grey shading represents the posterior probability density. Vertical tick marks above the x-axis represent observed values for environmental covariates, jittered to improve visibility of the density of observed values.

Fig. 7.  Relationships between probability of site survival (ϕ, i.e. the site remains occupied) for Dixie Valley toads (Anaxyrus williamsi) and environmental covariates, with posterior probability of inclusion >prior probability of inclusion. Solid lines are median predictions from the posterior distribution, dashed lines are 95% equal-tailed credible intervals, and grey shading represents the posterior probability density. Vertical tick marks above the x-axis represent observed values for environmental covariates, jittered to improve visibility of the density of observed values.

The distribution of Dixie Valley toads was closely tied to the availability of water in their desert spring-fed wetland habitat. Dixie Valley toad larvae were most likely to be found in less-vegetated areas with high surface water coverage that were farther from spring heads and were warm, but not too hot (20–28°C). The cells in which Dixie Valley toad larvae occurred were likely to be warmed by solar radiation rather than hot spring water, on the basis of their distance from spring heads. Adult Dixie Valley toads were more likely to occur in areas with more surface water, water depth of at least 10 cm, and a greater emergent vegetation cover. The importance of thermally stable surface water to both pre- and post-metamorphic life stages of this species is not surprising, given the xeric terrestrial habitat surrounding the spring-fed wetlands. In a radio-telemetry study, adult Dixie Valley toads were most often located in water, and rarely moved more than 14 m from aquatic habitat (Halstead et al. 2021). The Dixie Valley experiences cool spring and fall temperatures, and radio-telemetry data have shown that adult toads select warmer water and substrate temperatures during these seasons (Halstead et al. 2021). Other Anaxyrus species behaviourally thermoregulate to maintain more stable body temperatures than their environment (Rausch et al. 2008; Forget-Klein and Green 2021). Development of aquatic amphibian larvae is closely tied to water temperatures, with larvae commonly reaching metamorphosis faster when reared in warmer water (Smith-Gill and Berven 1979). Given the dependence of adult and larval life stages on thermally stable surface water, disturbances to the hydrology of the springs in Dixie Valley could affect the area of occupancy for this endemic toad species.

Emergent vegetation cover exhibited contrasting relationships with the occupancy of each life stage. Dixie Valley toad larvae were more likely to occur in less vegetated cells, and cells with a greater vegetation cover were more likely to be colonised by adults. The increased probability of occurrence for Dixie Valley toad larvae in cells with less vegetative cover could reflect a preference for warmer, more sunlit areas of the wetland during the spring to facilitate faster growth. In lab experiments, Anaxyrus [Bufo] boreas tadpoles have been shown to select water temperatures at or above 28°C (Beiswenger 1978; Bancroft et al. 2008). Alternatively, larvae could prefer areas of open water to reduce predation risk from insect predators (Gunzburger and Travis 2005), which can reach greater densities among aquatic vegetation in lentic habitats (Tolonen et al. 2003) and are likely to depredate Dixie Valley toads (US Fish and Wildlife Service 2022a). Changes in vegetation management practices have been shown to affect the density of emergent vegetation in wetland habitat for a congeneric desert spring endemic, the black toad (Anaxyrus exsul) (Murphy et al. 2003). Overgrowth of emergent vegetation has the potential to negatively affect amphibian reproduction (Greenberg and Green 2013; Bansal et al. 2019), as overgrowth has affected fish inhabiting desert springs by reducing open water (Kodric-Brown and Brown 2007). Because of the conflicting relationships adult and larval Dixie Valley toad occupancy had with emergent vegetation cover, maintaining a mix of vegetated and open-water habitats might be key to meet the needs of all life stages of the Dixie Valley toad.

There was substantial turnover in which cells were occupied by Dixie Valley toads between April 2018 and May 2018 and between September 2019 and May 2020, with higher colonisation rates in the latter period of each pair. In contrast, turnover was low from May 2018 to May 2019 and from May 2020 to May 2021. The high turnover from early to late spring and between fall and spring is likely to reflect differences in habitat use by adult toads between seasons, because individuals emerge from brumation to mate and forage in spring and prepare to enter brumation in the fall (Halstead et al. 2021). It is possible that seasonal changes in activity meant that some toads were not available for detection during the early spring and fall primary periods, biasing estimates of occupancy during these periods. Still, toads were more likely to be found in cells with surface water in all periods. Despite the turnover among occupied cells between seasons, the site-survival rate and overall occupancy rate were fairly stable throughout our study, with between 55% and 70% of cells estimated to be occupied by Dixie Valley toads in any given period. Dixie Valley toads have a narrow distribution restricted to a single desert spring-fed wetland complex, but utilise much of the available habitat. Within the study area, adult toads were detected within 75% of all surveyed cells at some point during our study, whereas larvae were detected only in 20% of cells. The more restricted distribution of Dixie Valley Toad larvae underlines the importance of warm, sparsely vegetated aquatic habitat to facilitate recruitment and persistence of the species.

Our use of occupancy models differs from that in studies of other amphibians, in which occupancy is classified at the scale of discrete wetlands separated by non-breeding terrestrial habitat. Dixie Valley toads inhabit a small number of wetlands within a few square kilometres, which makes modelling occupancy at the wetland scale largely irrelevant to understanding what influences occurrence of this narrowly distributed species. Modelling occupancy at sites within continuous habitat changes the interpretation of the estimated occupancy parameters (Efford and Dawson 2012). Because we divided largely continuous wetland habitat into 20 m × 20 m cells and modelled occupancy at the cell scale, the estimates of ψ presented here might be best interpreted as the proportion of the wetland area used by Dixie Valley toads in that primary period. The cell size in this study (400 m²) is of similar magnitude to the core activity areas of Dixie Valley toads estimated from radio-telemetry in the fall (mean = 173 m², 95% equal-tailed interval = 8–823 m²), but core activity areas in the spring months are unknown (Halstead et al. 2021). If core activity areas differed between spring and fall, it could lead to different estimates of ψ during each season. Also, it is possible that some toads moved between cells within a single primary period, which could inflate our estimates of ψ in the classical interpretation of area occupied by a species, although ψ could still represent the area used by toads during that period.

The effects of survey conditions on the detection probability of adult and larval Dixie Valley toads have implications for planning future monitoring for this species. For both life stages, surveys of 400 m² cells should last at least 15–20 min to have a satisfactory detection probability. Surveys near midday in the spring were more likely to detect larval Dixie Valley toads, whereas surveys targeting adults during the morning or evening and at higher ambient air temperatures are likely to have a greater success. Our realised detection probabilities for larval Dixie Valley toads were highly variable among primary periods. This temporal variation could reflect differences in the timing of surveys and interannual variation in the timing of breeding; smaller, younger larvae could be harder to detect than are larger, older larvae. Larval abundance could also vary inter-annually, because of females skipping breeding in some years (documented in A. boreas, Muths et al. 2010), and variation in abundance can result in variable detection probabilities (Royle and Nichols 2003). If documenting reproduction is important for monitoring the viability of this endangered species, future surveys could be targeted towards conditions with a higher detection probability for larval Dixie Valley toads.

Our results quantified relationships between occupancy rates of Dixie Valley toads and environmental conditions within wetlands. Desert springs can exhibit inter-annual changes in discharge in response to precipitation trends (Weissinger et al. 2016) and seasonal fluctuations in surface water in spring-fed wetlands can affect the habitat use of resident aquatic vertebrates (Grover 2019). The suitability of wetland habitat for adult and larval Dixie Valley toads will likely vary temporally in response to natural and anthropogenic changes in water availability and temperature. A key future step for modelling Dixie Valley toad occupancy will be to model seasonal and inter-annual variation in water availability and temperature. A dynamic habitat model could be linked to the occupancy model presented in this study to project changes in Dixie Valley toad occupancy under future scenarios of water availability and thermal regimes. Integrated habitat suitability–occupancy models can show relationships between metapopulation processes (e.g. colonisation, extinction) and habitat dynamics (Mackenzie et al. 2011).

Globally, desert spring habitats are under threat from groundwater extraction, surface-water diversion, climate change, and invasive species (Shepard 1993; Unmack and Minckley 2008; Davis et al. 2017; Rossini et al. 2018). Such disturbances have led to the extirpation and extinction of endemic fish and invertebrates from desert springs (Contreras-Balderas and Lozano-Vilano 1996; Minckley 1999; Kodric-Brown and Brown 2007; Hershler et al. 2014). In addition to these stressors, the development of geothermal energy in the Great Basin of the western USA could pose an acute threat to spring-fed wetlands by causing springs to run dry, as has occurred in nearby Jersey Valley, Nevada, USA (Bureau of Land Management 2019). For species dependent on desert springs, such as the Dixie Valley toad, loss and degradation of wetland habitat could be damaging to their long-term viability. Conserving narrowly distributed desert spring endemics requires an understanding of the habitat requirements of all life stages of a species. Our study has demonstrated that multi-state dynamic occupancy models have promise for informing conservation and management actions for narrowly distributed endemic species.

Supplementary material is available online.

Data are available on ScienceBase (Rose et al. 2022a; https://doi.org/10.5066/P9QCIC87) and model code to reproduce analyses are available on GitLab (Rose et al. 2022b; https://doi.org/10.5066/P97DSXJM).

The authors declare no conflicts of interest.

This project was funded by the US Bureau of Land Management, US Fish and Wildlife Service, US Navy, and the US Geological Survey Ecosystems Mission Area.