Designing a Long-term Occupancy Monitoring Plan for a Cryptic Reptile

Worldwide, many reptile populations are declining. For example, seven species of tortoise have recently gone extinct and 48 (73.8%) of the 65 extant species are currently threatened with extinction (Rhodin et al., 2018). Threats to reptile and tortoise populations are numerous and include habitat degradation and loss, commercial and hunting overexploitation, road mortality, and predation by subsidized predators (Rhodin et al., 2018; Stanford et al., 2020). Mojave Desert Tortoises (Gopherus agassizii; Fig. 1), endemic to the Mojave Desert in the southwestern United States, are no exception to declines in Testudinidae. Mojave Desert Tortoises have been federally listed as a threatened species since 1990 (USFWS, 1990), yet they continue to exhibit long-term declines in abundance across most of their range (Allison and McLuckie, 2018). Significant conservation work to recover Mojave Desert Tortoise populations is ongoing (Averill-Murray et al., 2012), and a key component to ensure successful recovery and to identify the need for additional conservation actions is effective monitoring.

Monitoring programs for Mojave Desert Tortoise populations generally use one of three methods: line-distance sampling, occupancy monitoring, or mark–recapture density estimates. Line-distance sampling is the current method used for range-wide monitoring of trends in the density of Mojave Desert Tortoises (Allison and McLuckie, 2018). Line-distance sampling accommodates imperfect detection as a function of distance or visibility from the survey track and with additional effort can adjust for unavailability for detection along the survey line itself (Allison and McLuckie, 2018). Although successful at large scales, line-distance sampling requires a relatively large number (∼60–80) of detected individuals (Buckland et al., 2001) per strata to reliably estimate the distance decay function. The large number of individuals that are required for the line-distance sampling method limits its utility for Mojave Desert Tortoises at local levels because tortoise densities are frequently in the range of 1–2 adults/km² (Allison and McLuckie, 2018), thus requiring extensive transects to achieve 60–80 detected individuals (Smith et al., 2009; Stober et al., 2017). Occupancy monitoring is a more recent technique, where rather than monitoring density or abundance, the proportion of sampling sites that are occupied are monitored (MacKenzie et al., 2002). Occupancy monitoring accommodates imperfect detection because of crypticity and can accommodate unavailability for detection if confounded with imperfect detection (Harju and Cambrin, 2019). Occupancy monitoring has also been demonstrated to be more cost-efficient with higher statistical power to detect declines in occupancy than line-distance sampling can detect declines in abundance for Sonoran Desert Tortoises (Gopherus morafkai; Zylstra et al., 2010). Mark–recapture density estimates can be robust and scaled down to local areas, but are often cost-prohibitive for Mojave Desert Tortoises because of the intensive field effort required to accumulate sufficient capture and recapture samples (Longshore et al., 2003).

Occupancy monitoring is a promising option for monitoring Mojave Desert Tortoises at smaller scales because of its comparative ease of implementation and reflection of a fundamentally important population metric (i.e., occupied range). However, as is the same for other methods of monitoring Mojave Desert Tortoises, it is subject to sporadic bias and noise induced by seasonal unavailability for detection. Seasonal unavailability for detection may be widespread across a variety of taxa, including reptiles and ephemeral aquatic and terrestrial plants (Alexander et al., 1997; Chen et al., 2013; Couturier et al., 2013; Gray et al., 2013; Harju and Cambrin, 2019). For example, perennial forbs and grasses may not produce aboveground parts in a given year, even though the plant itself is truly present at the site (Alexander et al., 1997). Mojave Desert Tortoises present a tortoise version of the problem of seasonal unavailability for detection. The amount of time that Mojave Desert Tortoises spend above ground, and are thus available for detection, varies dramatically in response to weather and vegetation and may be as low as 11–20% of the spring and summer (Duda et al., 1999). This means that for a given survey year, truly present tortoises may remain undetected on all survey occasions, suggesting absence. Harju and Cambrin (2019) dealt with seasonal unavailability for detection of Mojave Desert Tortoises by relaxing the assumed period of demographic closure to encompass multiple years, thereby allowing for estimation of relationships between landscape features and latent occurrence, given large interannual variation in apparent occupancy. But for long-term monitoring of occupancy trends, there is no clear way to deal with seasonal unavailability for detection; thus, long-term monitoring is currently subject to high statistical noise and low precision in trend estimates (Allison and McLuckie, 2018).

We focus here on developing a long-term monitoring program for Mojave Desert Tortoise occupancy, specifically because of its applicability to monitoring local populations. Clark County, Nevada, USA, is responsible for managing a 35,000-ha conservation easement (Boulder City Conservation Easement [BCCE]) southeast of Las Vegas, primarily for the conservation of Mojave Desert Tortoises and other species. An occupancy monitoring program (sensu MacKenzie et al., 2002) was implemented from 2013 to 2017. We used the observed results and high interannual swings in apparent occupancy because of seasonal unavailability for detection to simulate power to detect specified declines in Mojave Desert Tortoise occupancy rates under different sampling strategies. We sought to determine the outside bounds of sampling parameter combinations (e.g., number of sample sites and number of secondary period sampling visits per year) and determine the power that different sampling designs would have to detect a range of true changes in occupancy over a 10-yr period. Based on these results, we present recommendations for implementing a long-term occupancy monitoring program with sufficient power to reasonably detect true changes in Mojave Desert Tortoise occupancy rates.

From 2013 to 2014, the Clark County Desert Conservation Program conducted occupancy monitoring at 75 4-ha sites within the BCCE, with three secondary period surveys within each year (March–May), and from 2015 to 2017 at a subset of 60 sites, with seven secondary period surveys within each year (March–June). Dropped sites were excluded at random to reallocate within-year survey effort to the remaining 60 sites. The sampling design is referred to as the standard design, wherein S sites are visited each of K times within each of T years (MacKenzie et al., 2002).

Conducting occupancy sampling for Mojave Desert Tortoises, however, is noisier than simply recording presence and absence. To incorporate potential noise into the simulations, we added two types of variability: sampling noise and process noise. Sampling noise was the error observed as a function of imperfect detection, where some occupied sites do not have any observed individuals within the secondary or primary sampling period even though they truly are present (detection probability ∼ uniform(0.185, 0.734); Harju and Cambrin, 2019). Sampling noise was included via the statistical distributions that we used to generate the simulated datasets. Process noise was that arising from Mojave Desert Tortoise ecology, whereby apparent occupancy within a given primary sampling period was both a function of true occupancy and availability for detection. To model process noise, we used a parametric bootstrap based on the variation in apparent annual occupancy rates observed in a previous Mojave Desert Tortoise occupancy study (Harju and Cambrin, 2019). We used these data to fit a truncated normal distribution on the logit scale, with mean (μ = 0) and standard deviation (σ = 1.063) derived from the data on apparent annual occupancy rate deviations and the bounds of the distribution set to the minimum and maximum observed deviations (a = −1.53, b = 1.647; logit scale). The process ‘noise' parameter was then added to the true across-site occupancy rate, the summed value of which was used as the probability parameter in the sampling distribution used to generate simulated datasets with sampling noise (i.e., the process noise parameter modified the probability of a site appearing occupied in a given year, based on the underlying true occupancy rate in that year; see Supplementary Data for all R code used in the simulations). It is likely that this noise is a weather-related phenomenon, because Mojave Desert Tortoises restrict aboveground movement in drought years (Duda et al., 1999), but for these simulations the noise was treated as a random feature of surveys for Mojave Desert Tortoise occupancy. Furthermore, Harju and Cambrin (2019) assumed constant occupancy during their analysis period (e.g., no trend in occupancy), an assumption that becomes less likely over longer timescales. Here, this means that variability in apparent occupancy is a conservative estimate, in that true variability may be lower and realized power to detect a change in occupancy may be higher in practical settings.

We generated 150 simulated 10-yr occupancy datasets for each of the 270 possible combinations of 10 site sample sizes, three visitation schedules, and nine annual rates of change. For each of those 10-yr occupancy datasets we modeled colonization and extinction rates for individual sample sites as constant over the 10 yr. Colonization is the rate that previously unoccupied sites become occupied (e.g., increase in occupancy rate), and extinction is the rate that previously occupied sites become unoccupied (e.g., decline in occupancy rate; MacKenzie et al., 2003). We extracted colonization rate estimates (from models with positive annual rates of change) and extinction rate estimates (from models with negative annual rates of change) from simulated trajectories and then calculated 85% confidence intervals (CIs) on colonization and extinction rate estimates. We chose 85% CIs to reflect challenges with setting thresholds for “significance” too high given the sampling and process noise, while still seeking sampling combinations that would yield a high degree of power. We defined the power to detect true change, given sampling and process noise, as the rate that 85% CIs in estimated colonization and extinction rates excluded zero. In other words, we defined power as the frequency that the null hypothesis of “constant colonization or extinctions rates” was rejected with α = 0.15. The result compared and contrasted the power of different combinations of sample sites and of secondary visits for detecting a true change in occupancy rates for different magnitudes of true change. Because the true occupancy rate of change is unknown, we also averaged power to detect changes across all modeled rates of change to identify power equivalencies to guide effective allocation of occupancy monitoring budgets.

We also sought to determine whether there was nonlinearity in the power to detect change as a function of the number of sample sites to determine whether, and where, there was a way to find the most efficient balance between allocating funding to additional sample sizes vs. reduced returns on investment, in terms of power to detect changes in occupancy. We built a multiple linear regression model with the number of sites, number of visits, and quadratic rate of change as predictor variables. We then used segmented regression to identify whether, and where, there was a significant difference in the gains in power to detect change as a function of the number of sites surveyed. All data simulations, result processing, and statistical analyses were performed in program R 4.0 (R Core Team, 2021). The segmented regression was performed using the ‘segmented' package in R (Muggeo, 2003). See Supplementary Data for all statistical R code to create, summarize, and analyze simulated datasets and for an example input file defining combinations of parameters to simulate.

Processing time for the 40,500 simulated datasets, each with 10 yr of simulated data following the MacKenzie et al. (2003) dynamic occupancy model, was considerable, requiring ∼150 h of processing time. Incorporating both sampling and process noise resulted in large swings in apparent occupancy over time for any given simulation, matching the swings in apparent occupancy reported by Harju and Cambrin (2019). Nonetheless, on average, the trends were close to the true calculated trend in across-site occupancy (Fig. 3).

The power analysis showed several general patterns that were expected. Larger annual rates of change in true occupancy resulted in more power to detect true change in occupancy. In other words, as the strength of the decline or increase in the occupancy rate increased in magnitude, more of the 85% CIs on colonization or extinction rates in simulated datasets did not overlap zero (Fig. 4). Also, as the number of sites increased, the power to detect the true change in occupancy increased. By contrast, the number of secondary sampling visits in the simulated dataset did not appear to have a strong effect on increasing power to detect declines in occupancy (Fig. 4a), but was associated with power to detect increases in occupancy (Fig. 4b).

Larger annual rates of population change, both increases and declines, were more likely to exhibit plateaus in the power to detect change (Fig. 4). The highest rate of change analyzed, a 4% per annum decline or increase, generally exhibited equivalent power between 60 and 100 sites. Gains in power with additional sites were linear with lower rates of annual change in occupancy.

Power equivalency calculations, essentially a reorganizing and recombination of the data in Figure 4, provide guidance on selecting between the trade-offs of different sampling designs that yield equivalent power to detect true change (Fig. 5). Averaging across all potential rates of change in occupancy, we found that additional intraseason surveys could offset additional numbers of sites. For example, sampling 60 sites five times per year yielded an equivalent average power to detect change as sampling 80 sites three times per year (Fig. 5). The results could be useful if there are cost differences to prioritizing number of sites vs. number of visits. For example, assuming an average baseline personnel cost of $100 per site visit, sampling 60 sites with five visits each would cost $30,000 [$100 × 60 × 5 = $30,000], whereas sampling 80 sites with three visits each would cost $24,000 [$100 × 80 × 3 = $24,000], yet the power to detect a change in occupancy would be equivalent (52.2% power for 80 sites, three visits vs. 52.8% power for 60 sites, five visits). For reference, the 2017 sampling design on the BCCE of 60 sites with seven visits had an average power of 54.0% to detect a 1–4% change (increase or decline) in occupancy, and 60 sites with five visits would have an average 52.8% power.

The segmented regression found that the rate of increase in power with increasing number of sites declined by ∼60.0% when adding more than ∼47 sites (estimated breakpoint = 46.7 sites, 95% CI = 29.0–64.3; Fig. 6). For every additional 10 sites sampled, there was an increase in percentage points of power to detect change of 0.070 (95% CI = 0.046–0.094; e.g., on average, power was 0.29 with 20 sites and was 0.36 with 30 sites) up to 44 sites. Above 44 sites, an additional 10 sites yielded a power increase of only 0.028 percentage points (95% CI = 0.001–0.056).

Designing an effective long-term monitoring program requires achieving a balance between survey design and cost. We used existing occupancy survey data for Mojave Desert Tortoises to simulate and analyze 150 occupancy datasets over a 10-yr period for each of 270 possible combinations of 10 to 100 sample sites; three, five, or seven secondary visits per year; and under nine different scenarios of true change in the underlying across-site occupancy rate. We found that, as expected, higher numbers of sample sites and larger magnitudes of change in underlying occupancy rates were associated with higher power to detect the true change, given the presence of sampling and process noise (Guillera-Arroita and Lahoz-Monfort, 2012). Moderate declines in occupancy were more reliably detected (e.g., power to detect change >50%), especially with at least 30 sample sites, whereas increases in occupancy were not reliably detected until growth rates were 3–4% per annum with at least 40 sites. Segmented regression found that, across all magnitudes of change in underlying occupancy and number of secondary visits, the power to detect change in occupancy increased strongly with the number of sample sites until 47 sites and then increased much more weakly with additional sites.

These findings provide some guidance for developing a long-term monitoring program for Mojave Desert Tortoises in similar populations (Allison and McLuckie, 2018). Given the observed variability in apparent occupancy and detection rates reported in Harju and Cambrin (2019), no study design had strong power to reliably and consistently detect small trends of 1–2% per annum declines or increases in across-site occupancy. This agrees with occupancy power analyses for Sonoran Desert Tortoises and Gopher Tortoises (Gopherus polyphemus), where power to detect 1% declines in occupancy was unachievable (Zylstra et al., 2010; Chandler et al., 2020). For slightly larger per annum trends of 2% decline in Mojave Desert Tortoise occupancy, at least 50 sites would be needed to detect declines with a moderate level of power (i.e., 50–60% power). Therefore, it should be considered that if one of these assessed standard study designs is used, managers should not expect high power to detect small changes in occupancy.

For moderate-to-large trends in occupancy, such as 3 or 4% per annum increases or declines in occupancy rates, higher levels of power to detect the change (i.e., 70–80% power) were seen with 50 to 60 sample sites. By comparison, segmented regression showed a marked diminished return on field effort with more than 47 sites, across all numbers of visits and magnitude of occupancy rate change. Equivalency calculations also showed that across the modeled rates of change, sampling 80 sites for only three visits each season provided the same power level as sampling 60 sites with five visits each.

The largest single constraint inherent in these results is the application of the standard design. Other occupancy sampling simulation studies have found that rather than applying the standard sampling design, it is more efficient to apply a removal sampling design (MacKenzie and Royle, 2005) or a conditional sampling design (Specht et al., 2017). We argue not against these conclusions, but rather that such a flexible budget and effort allocation is not feasible for many organizations, where use-or-lose budgets are allocated before sampling or where shuffling of field technician responsibilities cannot accommodate high uncertainty in the scope of follow-up efforts (MacKenzie and Royle, 2005; Specht et al., 2017). It may be easier to accommodate these more efficient sampling designs for in-house studies, but it becomes much more complicated when contracting surveys outside of an organization because of the need for cost certainties and potential contracting constraints that may exist within an organization. Thus, although we have conducted these simulations within the framework of the standard design, we acknowledge that more flexible sampling paradigms may be available, albeit with more complicated budgeting, planning, and analysis requirements than the standard design.

We treated process error in availability for detection as a random process, but in reality such noise is likely a function of weather and climate. Mojave Desert Tortoises are known to restrict aboveground movements during drought years when annual forbs are not available as forage (Duda et al., 1999). Future work to evaluate the relationship between apparent occupancy and weather or climate could improve occupancy estimation and via inclusion as a covariate could reduce the apparent noise in occupancy estimates (Harju and Cambrin, 2019). We also note that dynamic occupancy models are potentially a function of multiple simultaneous drivers (i.e., colonization and extinction rates), and our simulations specified a “true” rate of occupancy increase or decline and then estimated colonization and extinction rates rather than fixing these rates themselves. In reality, changes in colonization or extinction are the driving patterns in observed annual occupancy, and multiple combinations of the two could result in the same observed occupancy pattern (MacKenzie et al., 2003). For example, either a reduction in colonization of 0.1 or an increase in the extinction rate of 0.1 would have the same impact on change in apparent occupancy. Including covariates behind the mechanism of occupancy rate change can help improve understanding and precision of the apparent change in occupancy rate.

These results suggest that for designing a long-term occupancy monitoring program for Mojave Desert Tortoises using the standard design, at a minimum, 47 sites should be surveyed to maximize the gain in power for each additional site. Furthermore, to maximize the total power equivalency and efficiency, we suggest that the number of sample sites for Mojave Desert Tortoise occupancy at the BCCE be increased from the current 60 to 80, but that the number of secondary sample visits be reduced to 3 to maintain power while reducing cost. These changes are expected to yield an equivalent power to the current sampling scheme, but at a lower cost. If a higher level of power is desired, local managers should even further increase the number of sample sites.

We thank M. Henson and other members of the Clark County Desert Conservation Program for support of field data collection for parameterization in this project. Numerous field technicians assisted with data collection and we extend our thanks. This work was funded through grants from the Southern Nevada Public Land Management Act and through funds collected as part of implementing the Multiple Species Habitat Conservation Plan. Project funders played no role in study design, data collection, statistical analysis, or manuscript preparation. All initial fieldwork was completed under the MSHCP 10(a)1(B) permit TE 034927-0 from the U.S. Fish and Wildlife Service.

Supplementary data associated with this article can be found online at http://dx.doi.org/10.1670/ 21-087.S1.