Digging into Detectability: Uncovering How Temperature Influences Detection Probability of the Fossorial Temblor Legless Lizard

Conserving biodiversity requires accurate information on the distributions of species. For rare species, one or a few surveys might be insufficient to detect its presence at a site. Occupancy models were developed to account for such “false-negative” results by estimating detection probability based on replicated surveys (MacKenzie et al. 2002; Tyre et al. 2003). The nondetection problem is exacerbated if species are often unavailable for detection, such as highly mobile animals or fossorial animals that spend the majority of their time underground. Multiscale occupancy models disentangle species’ availability for detection from the probability of detecting the species given it occurs at the site and is available to be detected (Nichols et al. 2008; Mordecai et al. 2011). For example, many terrestrial salamanders spend most of the year underground and might only be available for detection at the surface when temperatures and humidity are suitable (Halstead et al. 2022). If availability is not explicitly modeled, estimates of occupancy from traditional single-scale models can be biased low (DiRenzo et al. 2022). When availability for detection is much lower than 1, occupancy surveys for rare species can benefit from a hierarchical design that allows for the separation of detection probability and availability.

Designing an occupancy study for rare or understudied species is challenging. The optimal number of sites and replicate surveys depends on expected detection and occupancy probabilities (MacKenzie and Royle 2005). Choosing the sample of sites presents another dilemma. Inference into environmental covariates of species’ occupancy is strongest with a random sampling design, but a random sample could lead to selection of few occupied sites and poor estimates of occupancy (Pacifici et al. 2016). By contrast, only surveying sites expected a priori to have a high occupancy probability can bias estimates of occupancy and detection probability (MacKenzie and Royle 2005) and reduce the likelihood of documenting new localities. For species of conservation concern, acquiring information on occupancy can be valuable to inform time-sensitive species status reviews. A promising method for rare species is the two-phase occupancy design developed by Pacifici et al. (2012). In a two-phase occupancy study, an initial sample of sites is selected at random, an occupancy model is fit to those data, occupancy probability is predicted at a wider pool of sites, and a second sample is selected with the probability of inclusion based on that predicted occupancy probability (Pacifici et al. 2012). By immediately incorporating information learned during the first phase of surveys, the chances of documenting new localities are increased, providing timely data to inform conservation decisions.

Occupancy sampling could provide valuable information on the ecology of North American legless lizards (genus Anniella). In 2013, the California legless lizard Anniella pulchra was split into five species based on genetic and morphological evidence (Papenfuss and Parham 2013). Little is known about the distribution, ecology, and activity of the recently described Temblor legless lizard A. alexanderae (Figure 1) beyond its occurrence at a few localities in the western San Joaquin Valley of California (Parham et al. 2019). Because of its small putative range and potential impacts from large-scale agriculture and fossil fuel extraction, the species status is currently under review by the California Department of Fish and Wildlife (CDFW; California Fish and Game Commission 2022) and the U.S. Fish and Wildlife Service (U.S. Fish and Wildlife Service 2021). Ascertaining the range of the Temblor legless lizard, its prevalence, and the probability of detecting them at occupied sites are all critical for determining whether it warrants conservation protections. Legless lizards are fossorial predators of invertebrates that spend most of their time underground in loose, sandy soils (Miller 1944; Stebbins 2003; Papenfuss and Parham 2013) and therefore could be unavailable for detection if surveys are conducted during conditions that are unsuitable for activity near the surface. Given the paucity of information on the distribution of the Temblor legless lizard and its fossorial habits, combining a two-phase sampling design with a multiscale occupancy model is a promising approach to learn about its activity, distribution, and ecology.

In this study, we conducted occupancy surveys for Temblor legless lizards in the western San Joaquin Valley following a two-phase design. We then used a multiscale occupancy model to quantify availability near the soil surface and detection probability of this fossorial species. We addressed three objectives: 1) quantify how environmental conditions influence the probability of detecting Temblor legless lizards, 2) document new localities for the species, and 3) obtain preliminary estimates of factors influencing the species’ occurrence. Our results shed light on the ecology of this recently described species and point to future directions for research to answer conservation and management questions.

We conducted our study in the western San Joaquin Valley, in Fresno and Kern counties, California (Figure 2), which is characterized by hot, dry summers and mild winters during which most precipitation occurs. Dominant land uses in our study region include ranchland for grazing cattle, agriculture, and fossil fuel extraction. Vegetation communities include annual grassland (primarily nonnative grasses in grazed areas) and alkali desert scrub (California Department of Forestry and Fire Protection 2015). We defined our study area using a 10-km buffer around the hypothesized range of the Temblor legless lizard (Parham et al. 2019) to expand sampling beyond known localities for the species without expending effort in areas separated by clear geographic barriers. We focused our surveys on public lands owned and managed by the U.S. Bureau of Land Management (BLM) or the CDFW because most privately owned land in this region is challenging to access. Sites ranged in elevation from 98 to 589 m and in latitude from 35.17°N to 36.26°N.

We used a two-phase occupancy design (Pacifici et al. 2012) to select sites to survey. In the first phase, we used a generalized random tessellation stratified (GRTS) sampling design (Stevens and Olsen 2004) to pick a spatially balanced random sample of sites on accessible public lands. We overlaid a 200 m × 200 m square grid onto BLM and CDFW lands within the study area and numbered each grid cell. We removed from our sampling pool parcels of public lands that were inholdings in private land that the management agency could not access. We calculated the percentage of sand in the upper 50 cm of soil for each grid cell using the Gridded Soil Survey Geographic data for California (Soil Survey Staff 2018). Because Anniella are associated with loose, sandy soil suitable for burrowing (Miller 1944; Kuhnz et al. 2005), we dropped cells with <30% sand. We also used vegetation data from the CALFIRE Fire and Resource Assessment Program (California Department of Forestry and Fire Protection 2015) to mask out cells with clearly unsuitable land cover, such as agriculture or urban areas; only cells with >50% combined cover of annual grassland or alkali desert scrub were included in the pool of potential sites based on previously documented association with Temblor legless lizard occurrences (Parham et al. 2019). We selected a sample of 58 random cells using a GRTS design and stratified the sample by landowner and parcel to ensure that nine cells fell within CDFW Ecological Reserves (ERs; four at the Pleasant Valley ER in Fresno County and five at the Lokern ER in Kern County) and 49 within BLM lands, which collectively made up most of the total area available to survey. Two additional cells were manually selected for sampling in 2022 because they were previously occupied by Temblor legless lizards, and these “reference sites” were likely to provide information on detection probability and availability (Halstead et al. 2022). One reference site was in the Pleasant Valley ER and the other was a private ranch in Kern County (Parham et al. 2019). As part of the GRTS selection, we also selected 60 “over-sample” sites as backups in case any of the first-choice sites were inaccessible. If a first-choice site was not accessible, we chose the next accessible over-sample site in the sequence to maintain spatial balance in the GRTS design (Stevens and Olsen 2004). Study sites were concentrated at the northern and southern ends of the hypothesized range of the Temblor legless lizard because few accessible parcels of public land exist in the center of the range (Figure 2).

Anniella are commonly sampled using artificial cover objects (i.e., cover boards) made from plywood or cardboard (Papenfuss and Parham 2013). In September and October 2021, we set 25 cover boards at each site. We made 13 cover boards from 23/32-in-thick (1.8-cm-thick) plywood sheathing measuring 2 ft × 2 ft (0.6 m × 0.6 m) and made 12 cover boards by folding flat a double-wall corrugated cardboard box (S-4731; ULINE, Pleasant Prairie, WI) lengthwise such that it measured 18 in × 36 in (45.7 cm × 91.4 cm), gluing the two halves together, and applying wood sealant to the top surface to retain moisture under the cardboard. On BLM lands where grazing occurred, we covered cardboard cover boards with fiberglass mineral-surfaced roofing material to discourage disturbance by cattle. Within each selected cell, we set cover boards in the habitat that appeared most likely to be used by Temblor legless lizards based on the presence of vegetation and loose soil for burrowing, and cover board arrays varied in shape depending on the habitat. If the habitat was homogeneous, we set boards in a square grid covering approximately 20 m × 20 m. If a wash or patchily distributed shrubs occurred in a cell, we positioned cover boards around those features as we expected Temblor legless lizards to prefer these features (Miller 1944; Kuhnz et al. 2005). We recorded spatial coordinates of the center and corners of each cover board array using a handheld Global Positioning System (Etrex 10; Garmin Ltd., Olathe, KS) in Universal Transverse Mercator North American Datum 1983. We placed one iButton Hygrochron (DS-1923-F5#; Maxim Integrated, San Jose, CA) under the central cover board at each site to record temperature and relative humidity hourly. iButtons measured temperature to the nearest 0.5°C and relative humidity on a percentage basis (0–100) to the nearest 0.64%. We then surveyed these cover boards for Temblor legless lizards in winter and spring 2022 (see Field data collection below).

After preliminary occupancy modeling based on survey results from 2022, we dropped the 20 sites with the lowest predicted probability of occupancy from our model. The predicted probability of occupancy was based on both ecological covariates and the number of surveys conducted at the site (more surveys with nondetection results in lower predicted probability of occupancy). We selected 20 new sites for 2023 using the GRTS algorithm with the probability of selection weighted by the predicted occupancy probability from the output of the phase 1 model (Pacifici et al. 2012). Predicted occupancy probability ranged from 0.11 to 0.29 based on distance to the nearest ephemeral stream (i.e., wash) in the National Hydrography Dataset Plus version 2.0 (U.S. Geological Survey 2019) and soil composition from the Gridded Soil Survey Geographic data for California (Soil Survey Staff 2018). We then manually selected 10 additional nonrandom sites to sample based on an expected high probability of Temblor legless lizard occupancy from local biologists’ expertise and review of stream courses and satellite imagery. In September 2022, we set 25 cover boards (13 plywood and 12 cardboard) at these 30 new sites, following the same procedure described above.

In 2022, we surveyed cover boards for Temblor legless lizards across 60 sites from 15 February to 15 April. In 2023, we surveyed 68 sites from 8 February to 12 April, along with 4 additional surveys on 23 May 2023. We surveyed sites a variable number of times by design (see Occupancy modeling below). We conducted all surveys between 0830 and 2030 hours because we expected Temblor legless lizards to be diurnally active during the winter and early spring. We recorded substrate temperature (to the nearest degree Celsius) at the time of surveying using a handheld thermometer (Model 9842; Taylor Precision Products, Oak Brook, IL) that we placed on the soil surface under the central cover board. We also collected data on the microhabitat surrounding each cover board (grass, bare ground, tree, shrub, or litter). We categorized soil texture (i.e., particle size) by hand into 1 of 12 categories following a protocol for field characterization of soil texture by feel (Thien 1979), which, when used by trained technicians in the field, can produce classifications comparable to laboratory analysis (Salley et al. 2018). Each category is a qualitative measure of the percent composition of sand, silt, and clay that is defined by the U.S. Department of Agriculture Soil Textural Triangle (Soil Survey Staff 1951).

We used fingers or a 3-tine hand cultivator to search the underlying substrate for legless lizards. We recorded the Universal Transverse Mercator coordinates, cover board type, and microhabitat where each Temblor legless lizard was found. We measured snout-to-vent length and total length of captured individuals to the nearest millimeter using a meter stick. We measured the mass of individuals to the nearest gram using micro-line spring scales (PESOLA, Switzerland). We identified individuals as Temblor legless lizards based on the distinctive light gray ventral coloration unique to this species (Papenfuss and Parham 2013). We took photographs of dorsal and ventral patterns to share with T. Papenfuss (author of the species description for Temblor legless lizards) for confirmation of species identity; we confirmed that all photographed Anniella were Temblor legless lizards. The species description states that genetic analysis is not required to distinguish Temblor legless lizards from other Anniella; morphological characters are sufficient (Papenfuss and Parham 2013; Parham et al. 2019). After sampling, we released all Temblor legless lizards under the board under which we had found them.

To assist CDFW in its status review of Temblor legless lizards, we collected tail tips from up to five Temblor legless lizards per site in 2023 for genetic analysis by the CDFW Wildlife Health Lab (Rancho Cordova, CA). We obtained tissue by cutting approximately 5–6 mm of the tail tip with sterilized surgical scissors. CDFW geneticists sequenced two previously defined genetic markers for differentiating species in the genus Anniella (Parham and Papenfuss 2009): a mitochondrial DNA sequence (NADH dehydrogenase subunit 2 and five adjacent tRNAs) and a nuclear gene (brain-derived neurotrophic factor precursor). All 16 genetic samples were confirmed as belonging to A. alexanderae based on neighbor joining methods (M. Buchalski, unpublished data); these sequences grouped with clade “B” in Parham and Papenfuss (2009), which was later described as A. alexanderae in Papenfuss and Parham (2013).

We collected iButtons on the last survey date for each site and downloaded data covering the time period that the iButton was deployed in the field. Of 60 sites with iButtons deployed in 2022, we recovered data from 57. Of 68 iButtons deployed in 2023, we recovered data from 54, whereas the remaining iButtons were lost or rendered inoperable due to flooding. For sites with missing iButton data, we used Bayesian imputation (Bonner and Schwarz 2006) to impute missing temperature and relative humidity data. For temperature, we defined priors for missing data using a mean equal to the temperature collected under a cover board with a handheld thermometer at the time of checking and a standard deviation calculated using iButton temperatures collected at nearby sites at the same date and time. For relative humidity, we defined priors for missing data using the mean and standard deviation of relative humidity recorded from nearby sites at the same date and time.

Although we selected the majority of sites using GRTS, our inclusion of nonrandom sites (n = 12) in our sample (n = 89) could bias estimates of ψ. We fit the following two occupancy models to evaluate the importance of distinguishing random and nonrandom (manually selected) sites: 1) a model with a single parameter for ψ and 2) a model with two parameters, one for random sites and one for manually selected sites. Given the nonrandom selection of some sites, we focus on finite sample occupancy (the number of occupied sites within our sample) and do not generalize to occupancy probability throughout the study region.

Sites were spatially clustered in some areas because of limited accessibility to public lands. We tested if spatial clustering of sites influenced occupancy probability by modeling spatial covariance with distance between sites using a Gaussian process (Johnson et al. 2013). We found no support for covariance in occupancy varying as a function of distance; the posterior distribution for the effect of distance was nearly equal to the prior distribution. Given this lack of evidence of spatial covariance, we used a nonspatial occupancy model for inference.

In 2022, we did not strictly follow a hierarchical sampling design for a multiscale occupancy model as we refined our survey protocol. We sampled 60 sites on a minimum of three secondary occasions (i.e., distinct weeks) each. We performed three tertiary surveys (i.e., cover board checks) at 42 sites and four or more tertiary surveys at 18 sites (Table S1, Supplemental Material). We did not perform more than one tertiary survey per secondary occasion at most sites in 2022, but at our reference site in the Pleasant Valley ER and selected sites at which Temblor legless lizards were detected, we performed more than one tertiary survey per secondary occasion to obtain additional information on detection probability and availability (Table S1, Supplemental Material).

In 2023, we sampled all sites following a multiscale occupancy design. We sampled each site on at least two secondary occasions with two tertiary surveys within each secondary occasion to model availability and detection probability (Table S1, Supplemental Material). We separated secondary occasions by 3 to 4 weeks (mean = 25 d, range = 11–39 d), and we separated tertiary surveys by 1–3 d (mean = 1.75 d). Two sites that we sampled in 2022 could not be sampled in 2023 because winter flooding washed out cover boards, resulting in a total of 68 sites sampled in 2023. Of 68 sites, we sampled 35 on two secondary occasions (four total surveys), we sampled 33 sites on a third secondary occasion, with either one (n = 17 sites) or two (n = 16) tertiary surveys during the third secondary occasion. Finally, we surveyed three sites once more during a fourth secondary occasion in late May 2023 during a visit to retrieve iButtons.

We fit a Bayesian implementation of a multiscale occupancy model with Markov chain Monte Carlo sampling using JAGS version 4.3.0 (Plummer 2003) accessed through R version 4.2.3 (R Core Team 2023) using the “runjags” package (Denwood 2016). We based our model on the multiscale model used by Halstead et al. (2022) and adapted it to the survey design used for Temblor legless lizards. We fit models on four chains for 250,000 sampling iterations after discarding the initial 10,000 iterations as burn-in. We thinned the resulting chains by a factor of 10, resulting in a final posterior sample of 100,000 iterations for inference. We evaluated model convergence and mixing of chains by visually inspecting trace plots and calculating the potential scale reduction factor (Brooks and Gelman 1998); all parameters had potential scale reduction factor values of <1.01, indicating convergence. We ran a posterior predictive check based on comparing the expected number of detections to observed data and replicate data generated by the model (Kéry and Schaub 2012). The Bayesian P value was 0.615, indicating that there was not a lack of fit to the observed data. Unless otherwise noted, we summarized parameter estimates with the mean and 95% equal-tailed credible interval of the posterior distribution. Archived data are available on ScienceBase (Rose et al. 2024; https://doi.org/10.5066/P9R7UIV7), and model code to reproduce analyses are available on GitLab (Rose and Halstead 2024; https://doi.org/10.5066/P95O3WW2; Data S1, Supplemental Material).

We conducted a total of 239 cover board surveys in 2022 and 324 surveys in 2023 (Table 1). We observed Temblor legless lizards at 9 sites out of 60 sampled in 2022, with a total of 31 detections (unique site and date combination). Of these nine sites, the first detection occurred during the first tertiary survey at five sites, and we did not detect Temblor legless lizards until the third tertiary survey at four sites. We observed Temblor legless lizards at 10 sites out of 68 sampled in 2023, with a total of 23 detections. Of the 10 sites with observations in 2023, 2 were first sampled in 2023, 1 was sampled on three tertiary surveys in 2022 without detecting Temblor legless lizards, and the remaining 7 also had Temblor legless lizard detections in 2022. We did not detect Temblor legless lizards in 2023 at two sites where the species was observed in 2022. One site in the Pleasant Valley ER was flooded by Los Gatos Creek during heavy rains, and the habitat was highly disturbed. The other site in the Kreyenhagen Hills was not disturbed between years. Of the 10 known occupied sites in 2023, the first detection occurred during the first tertiary survey at four sites, during the second tertiary survey at one site, during the third tertiary survey at four sites, and during the fourth tertiary survey at one site. Notably, at a new locality for the species southwest of Coalinga, California, in both 2022 and 2023, we only detected a single Temblor legless lizard during the first of six tertiary surveys each year.

Combining the results over both years, we documented Temblor legless lizards at 12 out of 89 sites. Of these 12 occupied sites, we manually selected 4 for inclusion in the study based on the apparent suitability of the habitat, 1 site in the Pleasant Valley ER was previously known to be occupied by Temblor legless lizards, and 7 were random sites selected by the GRTS algorithm. We observed a total of 74 Temblor legless lizards under cover boards, with 41 observations under cardboard cover boards and 33 under plywood. We found most Temblor legless lizards at sites with soil characterized as loamy sand (n = 43), followed by sand (n = 16), loam (n = 7), silty clay (n = 5), and sandy loam (n = 3). Temblor legless lizards occurred at elevations ranging from 127 to 427 m, covering most of the range of elevation at surveyed sites (98 to 589 m). All sites at which we observed Temblor legless lizards had loose or friable soils. We observed most Temblor legless lizards under cover boards placed under shrubs or trees (n = 57), followed by grass (n = 13) and no vegetation (n = 4). We observed Temblor legless lizards under cover boards placed beneath screwbean mesquite Prosopis pubescens, California juniper Juniperus californica, saltcedar Tamarix ramosissima, and saltbush Atriplex sp. Daily mean temperature under cover boards fluctuated by day and week, with an overall positive trend from early February to mid-April in both years (Figure S1, Supplemental Material). Daily mean relative humidity under cover boards declined overall during the sampling period in 2022, with temporary increases following rainfall (Figure S1, Supplemental Material). In 2023, the relative humidity was higher on average and remained consistently high from early February to early April, only showing a clear decline at two sites where iButtons remained in the field until late May (Figure S1, Supplemental Material).

The occupancy model with one ψ parameter produced similar posterior estimates as the model with separate ψ parameters for random and manually selected sites (Table 2). Therefore, we focus on results from the model with one ψ parameter for all sites but do not extrapolate occupancy beyond our finite sample of sites. Detection probability was positively related to substrate temperature at the time the cover board was checked (Pr[β_temp,p > 0] = 0.997). The relationship between detection probability and substrate temperature was more uncertain for temperatures below 10°C and above 20°C, because fewer data points were available outside the range of 10–20°C (Figure 4). There was weak support for an effect of relative humidity under the cover board at the time it was checked on p (Pr[β_rh,p > 0] = 0.796). There was high uncertainty in the relationship between p and humidity for relative humidity of <60%, with most observations occurring at higher humidity values (Figure 4). The probability of detecting a Temblor legless lizard during a single survey at average temperature and humidity, given that they occupied the site and were available for detection at the time of the survey, was 0.514 (95% credible interval = 0.344–0.696). The mean availability of Temblor legless lizards for detection was high (mean = 0.771, 0.569–0.958), with little variation in availability among secondary periods (Figure S2, Supplemental Material). The product of θ and p (comparable to p in a single-scale occupancy model) for a single survey under average conditions was 0.392 (0.274–0.520).

There was a negative relationship between ψ and the distance to the nearest ephemeral stream (−0.872; −1.903 to −0.048); sites closer to ephemeral streams were more likely to be occupied (Pr[β_st.dist,ψ < 0] = 0.982; Figure 5). Sites without high clay content in the soil were more likely to be occupied than sites classified to have high clay content (Pr[β_clay,ψ < 0] = 0.995; Figure 5). Only 1 of 12 known occupied sites had high clay content in the soil, with texture classified as silty clay. Of the 89 sites surveyed in 2022 and 2023, the model estimated that a median of 13 sites (12–17) were occupied by Temblor legless lizards.

The number of sequential negative surveys (n*) required to have high confidence that a site was truly unoccupied by Temblor legless lizards depended on the substrate temperature at the time of the survey, soil texture, and the distance that the site was located from an ephemeral stream (Figure 6). Sites closer to ephemeral streams required more surveys because their baseline occupancy probability was higher. Likewise, sites with low clay content soils required more surveys with nondetection to have high confidence that the site was unoccupied. More surveys were required to have high confidence that a site was unoccupied if the substrate temperature at the time of the survey was low (<15°C) because p was lower at colder temperatures (Figure 6).

Our study provided the first quantitative estimates of detection probability for any Anniella species and identified substrate temperature as a strong predictor of the probability of observing Temblor legless lizards where they occur. Support for relative humidity under cover boards as a predictor of detection probability was weak, but humidity was consistently high during our sampling period from mid-February to mid-April, particularly in 2023, which had high winter rainfall. We expect that detection probability decreases at lower levels of humidity under cover boards (as a proxy of soil moisture) later in the spring and in summer, given the importance of soil moisture for Anniella (Miller 1944). Likewise, although we found a positive relationship between detection probability and temperature during the winter and early spring, our survey period ended before temperatures increased to the extreme heat that predominates in the region from June to September. It is likely that detection probability exhibits a unimodal relationship with temperature, with p decreasing at high temperatures later in the year based on the temperature preferences of Anniella (Bury and Balgooyen 1976; Miller 1944). Continuing surveys into the late spring and summer could help elucidate how high temperatures affect detection of Temblor legless lizards and when the species is not available for detection near the surface. It is likely that the species is available for detection later than mid-April when we ended most surveys; we observed a Temblor legless lizard on 23 May 2023 when daytime high temperatures were above 30°C and humidity under the cover board was approximately 50%. The detectability of Temblor legless lizards earlier in the winter also remains unknown, and nighttime low temperatures could influence activity near the surface (Miller 1944).

Based on the influence of substrate temperature on detection probability, the relationship between occupancy and distance to ephemeral streams, and estimated availability from our multiscale occupancy model, we calculated the number of sequential negative surveys needed to have high confidence that a site is unoccupied by Temblor legless lizards. It is clear from these calculations and the raw detection data that a few negative surveys are not sufficient to conclude that this species is absent from a site. In addition to prior belief about detection probability and occupancy, the estimated number of negative surveys presented above is also highly dependent on the availability of Temblor legless lizards for detection near the surface, which was generally high during our late-winter to early-spring survey period. If the probability that Temblor legless lizards were available for detection was lower (as expected in the hot, dry summer), then a larger number of negative surveys would be required to have high confidence that the species was absent from a site. We note that our estimates of detection probability are potentially dependent on the number of cover boards that we used at each site. Increasing the number of cover boards used to sample each site would likely increase detection probability, which could be valuable for determining occupancy of Temblor legless lizards. Our results also show that surveys over multiple years might be needed to document the presence of Temblor legless lizards, particularly during droughts when conditions at the surface are unsuitable.

We sought to strike a balance between primarily sampling sites selected at random while not wasting effort at locations that were clearly unsuitable for Temblor legless lizards. The inclusion of reference sites in 2022 and manually selected sites that appeared suitable in 2023 provided benefits in terms of data on detection probability and two new localities for the species. The cost for including these reference sites in our model is that our estimate of occupancy probability was potentially biased, although relationships between ψ and covariates were unchanged when random and manually selected sites were modeled with separate ψ parameters. Still, given the potential bias in ψ, we refrained from making predictions about occupancy probability throughout the study area. We used a two-phase sampling design to increase our chances of finding Temblor legless lizards at new locations in 2023, but the three sites with new observations in 2023 were either selected in the first phase and sampled in 2022 without detecting the species (one site) or manually selected to sample in 2023 (two sites).

Our ability to make strong inferences about ecological factors influencing occupancy of Temblor legless lizards was constrained by the low number of known occupied sites (12 out of 89), but we found some evidence for soil texture and distance to ephemeral streams affecting occupancy probability. The proximity of Temblor legless lizard occurrences to washes fits with knowledge of the ecology of other Anniella species (Miller 1944; Papenfuss and Parham 2013). Alluvial fans and washes are associated with loose, friable soils suitable for burrowing by Anniella, and the soils of these ephemeral drainages can retain greater moisture than surrounding habitats (Balding and Cunningham 1974; Bull 1977). Given the species’ fossorial nature and its method for burrowing, it is unsurprising that we found support for lower occupancy in sites with soils that had high clay content. Previous studies of Anniella emphasized the affinity for sandy soils (Miller 1944; Kuhnz et al. 2005; Papenfuss and Parham 2013). Therefore, we screened out areas with low sand content in the soil from the pool of potential sites, and within our pool of sites, most were classified as a type of sand or loam. Still, within a given site, there was microgeographic variation in soil texture and composition. A future study characterizing soil texture, composition, moisture, and temperature at a finer scale with multiple samples per site could lend further insight into habitat selection and activity of Temblor legless lizards.

Surrounding land use could influence the occurrence of Temblor legless lizards, although we observed the species at sites actively grazed by cattle, within 100 m of fossil fuel extraction, and within 500 m of agricultural fields. The scale at which land use affects Anniella occupancy is unknown, but given the small home ranges of individuals (mean = 71 m² for A. pulchra; Kuhnz 2000), populations might persist in undisturbed habitat adjacent to human land uses as long as the patch size is sufficient to support a viable population. Increased access to privately owned lands could facilitate sampling a larger portion of the species’ putative range, particularly the central two-thirds in Kings County and northern Kern County. Further refinement of our sampling design could be necessary to increase the likelihood of documenting the species at new localities and acquire a better understanding of the abiotic and biotic conditions affecting its occurrence at macro- and microgeographic scales.

Much remains to be learned about the ecology and conservation status of Temblor legless lizards. Nondetection of Temblor legless lizards during one or a few surveys cannot be conclusively interpreted to indicate the absence of this species from a site. At many sites, we did not observe Temblor legless lizards until we had conducted three or more surveys, and at one site we only observed a single individual at the first survey out of six. Furthermore, we demonstrated that Temblor legless lizards are far less likely to be detected when substrate temperatures are low (<15°C), and our results provide guidance for optimizing survey protocols during the winter and spring. Our findings also expand the species’ range beyond that hypothesized by Parham et al. (2019) into the foothills of the Diablo Range in southwestern Fresno County. It is likely that more localities outside the hypothesized range could be found by sampling habitat further outside this range polygon, which could help elucidate boundaries between the range of Temblor legless lizards and other Anniella species. One question for managers considering the status of the Temblor legless lizard is whether it is truly rare or whether the low occupancy rate in this study reflects the difficulty in identifying suitable habitat for this small fossorial species from landscape-scale data available as geographic information system layers. Our results indicate that answering this question could be facilitated by increased spatial and temporal coverage of sampling and better characterization of habitat features that influence the occurrence of Temblor legless lizards.

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

Reference S1. Parham JF, Koo MS, Simison WB, Perkins A, Papenfuss TJ, Tennant EN. 2019. Conservation assessment of the California legless lizard (Anniella). Sacramento, CA: Prepared for the California Department of Fish and Wildlife.

Data S1. Site data and survey data for Temblor legless lizard Anniella alexanderae surveys in 2022 and 2023 in western San Joaquin Valley, California, are available at the U.S. Geological Survey ScienceBase (Rose et al. 2024; https://doi.org/10.5066/P9R7UIV7). File “TLL_site_data.csv” contains data for each site sampled, with one row per site. File “TLL_survey_data.csv” contains data for each individual survey for Temblor legless lizards in 2022 and 2023. The code to reproduce the analyses is available on the U.S. Geological Survey GitLab (Rose and Halstead; https://doi.org/10.5066/P95O3WW2).

Table S1.The number of surveys per secondary occasion (Sec) at each site that we sampled for Temblor legless lizards Anniella alexanderae in 2022 and 2023 in western San Joaquin Valley, California. A value of “NA” indicates that the site was not sampled in that year.

Figure S1. Daily mean temperature (°C; A and C) and relative humidity (%; B and D) under cover boards at each site that we sampled for Temblor legless lizards Anniella alexanderae in 2022 (top row) and 2023 (bottom row) in western San Joaquin Valley, California. Thick black lines represent the daily mean averaged across all sites, and thin lines in color represent daily means for individual sites.

Figure S2. Posterior estimates of Temblor legless lizard Anniella alexanderae availability (θ) by week of the field season for 2022 (black) and 2023 (red) in western San Joaquin Valley, California, based on occupancy models fit to survey data. Points represent means of the posterior distribution, thick lines represent 50% equal-tailed credible intervals, and thin lines represent 95% equal-tailed credible intervals. For 2022, week 1 began on 15 February 2022 and week 9 began on 13 April 2022. In 2023, week 1 began on 7 February 2023, week 10 began on 10 April 2023, and week 11 was 23 May 2023.

We thank C. Irons and K. Smith for assisting with cover board deployment, L. Parker and L. Stewart for administrative support, and E. Tennant, J. Battistoni, M. Westphal, S. Bullock, D. Meade, and J. Dart for coordinating access to sites. We thank the CDFW, U.S. BLM, and B. Grant for allowing our research on their lands. We thank M. Koo for sharing the range polygon for the Temblor legless lizard and T. Papenfuss and J. Parham for providing advice on sampling for Temblor legless lizards. We thank M. Martinez for reviewing an earlier draft of this manuscript and three anonymous reviewers and the Associate Editor for comments that improved this paper. This work was supported by the U.S. Fish and Wildlife Service and the U.S. Geological Survey Ecosystems Mission Area. We performed all work under CDFW Scientific Collecting Permit SC-10779.

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