Environmental relationships can differ across the geographic range of species, especially for widespread generalists.  Because habitat specialists are more vulnerable to environmental changes, incorrect assumptions about consistent habitat associations could hinder strategic conservation efforts.  We used species distribution models (SDMs) to evaluate intraspecific variation in habitat associations for a habitat specialist of conservation concern, the pygmy rabbit (Brachylagus idahoensis), which is endemic to the sagebrush biome of the western USA.  Our goal was to model habitat associations for pygmy rabbits across a portion of their range to evaluate regional variation and contrast predictions with results from a model developed at the rangewide extent.  We created inductive SDMs using maximum entropy methods within five ecological regions that encompassed about 20% of the species rangewide distribution and spanned diverse environmental gradients.  We included a suite of environmental predictor variables representing topography, vegetation, climate, and soil characteristics.  Results of the regional models identified substantial variation in habitat associations across the five regions, with each retaining a unique set of environmental predictors.  Bioclimatic variables were the most influential environmental parameters in all five regions, but the specific variables differed.  The models developed at regional extents predicted smaller areas of habitat (an average of 15% less for suitable habitat and 80% less for primary habitat) than predictions generated from a model developed at the rangewide extent.  Because bioclimatic variables were effective in discriminating areas used by pygmy rabbits, they also provided an opportunity to assess potential changes in habitat distribution by incorporating future climate projections.  Distributions modeled under two mid-century emission scenarios projected substantial reductions in suitable habitat for pygmy rabbits across most regions and pronounced variation among regions in the magnitude and direction of the climate effects.  Collectively, results of this work underscore the need to incorporate regional variation in habitat associations into planning for current and future conservation and management strategies.

Study area

Our study focused on the distribution of pygmy rabbits in Idaho, USA, in the central portion of the species’ geographic range. We defined five ecologically based regions occupied by pygmy rabbits across southern Idaho for development of region-specific SDMs (Fig. 1). To delineate the regional extents, we first examined the Bureau of Land Management (BLM) fine-scale habitat management polygons for the greater sage-grouse, which reflect hydrologic units, topographic features, vegetation communities, barriers to sage-grouse movements based on telemetry data, and expert knowledge (Stiver et al. 2015). We overlaid known occurrences of pygmy rabbits and then combined adjacent fine-scale polygons to create five regions that encompassed separate areas occupied by pygmy rabbits (Regions: 1 - Owyhee Desert, 2 - Southern Desert, 3 - Central Desert, 4 - Bear Lake, and 5 - Lemhi-Salmon). The boundary between adjacent Regions 1 and 2 was defined by the Bruneau River Canyon, which is a prominent landscape feature >60km in length and up to 370m deep. The division between Regions 3 and 5 in part reflects differences in habitat between mountain valleys and watersheds associated with the Northern Rockies Ecoregion from the lower elevation areas that are more closely associated with the Snake River Plains Ecoregion (Bailey 2016). Three of the five regions overlap with the neighboring states of Oregon, Nevada, Utah, and Wyoming (Fig. 1). 

The climate in the study area is characterized by cold winters and hot, dry summers with most precipitation occurring during the cool season months (November-May; Runkle et al. 2017). The low-elevation areas of southern Idaho are shielded by mountains to the east and west resulting in generally limited precipitation (Runkle et al. 2017), however, variation in precipitation and temperature patterns exists among the five modeling regions (Table 1; WWRC 2018). In general, the southwestern portion of the study area has the lowest historical precipitation values and highest mean annual temperatures, and the eastern portion has the highest precipitation and lowest mean annual temperatures (RAWS data; Western Regional Climate Center, 2018; Table 1). 

Vegetation characteristics vary across the study area (Appendix S1: Fig. S1). The five ecological regions we defined coincide with five divisions of Idaho’s floristic provinces (Ertter & Moseley 1992). Within the study area, big sagebrush (Artemesia tridentata) and green rabbitbrush (Chrysothamnus viscidiflorus) are common shrub species at lower elevations, especially along drainages and alluvial fans. On the lower elevation flats, isolated areas of low sagebrush (A. arbuscula), bitterbrush (Purshia tridentata), shadscale (Atriplex confertifolia), and spiny greasewood (Glossopetalon spinescens) also are common. The dominant tree species within the sagebrush-steppe ecosystem are western juniper (Juniperus occidentalis) and Utah juniper (J. osteosperma), but there are also populations of pinyon pine (Pinus monophylla) scattered throughout. Many of the sagebrush-steppe habitats in southern Idaho include an understory of mixed native grasses (e.g., bluebunch wheatgrass, Pseudoroegneria spicata; Sandberg’s bluegrass, Poa secunda; bottlebrush squirreltail, Elymus elymoides; Indian ricegrass, Achnatherum hymenoides), seeded perennial grasses (e.g., crested wheatgrass, Agropyron cristatum), and annual invasive species (e.g., cheatgrass, Bromus tectorum).   

The study area is characterized by a diversity of topographic features. In general, the eastern portion of the study area is characterized by higher elevations and more complex topography (Table 2). The Lemhi-Salmon Region (Region 5) in east-central Idaho has broad valleys bounded by high mountain ranges including the Continental Divide of the Rocky Mountains, which runs along the eastern border. The Owyhee Desert Region (Region 1) is characterized by a combination of mountains and deep river canyons.  In addition, the distributions of land ownership and land use differ markedly across the five regions. In general, populations of pygmy rabbits in Idaho occur on multiple-use public rangelands managed by the BLM and, at higher elevations, on open lands managed by the US Forest Service (FS). Public rangelands are often embedded in a mosaic of private lands managed as pastures or agricultural land with native habitat adjacent or nearby. Across public and private lands, livestock grazing is one of the most common land uses.  Like much of the sagebrush ecosystem in the western USA, fires have increasingly impacted sagebrush communities in our study area (Mata-González et al. 2018; Miller & Heyerdahl 2018).

Environmental data

Environmental data used for the SDMs included bioclimatic variables, soil characteristics, vegetation properties and topography (Appendix S1: Table S1). All layers were projected to the same coordinate system (Albers equal-area conic projection, North American Datum (NAD) 1983, Contiguous United States) with the same geographic bounds and cell size (30m resolution), and converted to ASCII file format for input into R 3.5.2 (R Core Team 2018) and Maxent 3.4.0 (Phillips et al. 2017b). 

Land cover. – Land cover variables used for modeling habitat for pygmy rabbits included vegetation types known to be used (e.g., sagebrush), or alternately, avoided by the species (e.g., trees). Tree canopy cover has been negatively associated with pygmy rabbits, presumably because trees provide perches for avian predators, and increasing tree cover is associated with a reduction of the sagebrush or shrub understory (Larrucea & Brussard 2008b; Woods et al. 2013; Edgel et al. 2014). Agricultural lands also are expected to be avoided by pygmy rabbits, although some rabbits may use sagebrush habitats adjacent to agricultural lands. During this project, national land cover datasets were being updated, and consequently, we incorporated components from both newly updated and previous versions. Primarily, we used the updated Rangeland Fractional Components dataset from the National Land Cover Database (USGS 2020; Rigge et al. 2019). In addition, we included variables from the 2016 land cover database (i.e., agricultural land cover; USGS 2016a), the 2016 NLCD USFS Tree Canopy Cover (i.e., tree canopy cover; USGS 2016b), and the Provisional Remote Sensing Shrub/Grass NLCD Base Products for the Western US (i.e., sagebrush canopy cover; USGS 2016c).  

Although pygmy rabbits are associated with relatively high densities of sagebrush and other shrub species, many studies also have noted the importance of an understory of grasses and forbs (e.g., Heady & Laundré 2005; Larrucea & Brussard 2008a). The Rangeland Fractional Components dataset included a layer depicting percent herbaceous cover, along with a nested subset layer for cover of annual herbaceous plants. Because grasses and forbs comprise roughly half of the summer diet of pygmy rabbits (Shipley et al. 2006; Schmalz et al. 2014), we also included a measure of vegetation productivity, Normalized Difference Vegetation Index (NDVI), estimated for the late summer period. We used cloud-free eMODIS 7-day composite NDVI images (USGS 2016d) from July, August, and September to calculate the mean monthly maximum NDVI values, and the intraseasonal average and standard deviation from 2000 to 2016. This variable was used in previous work to model distributions for both sage-grouse and pygmy rabbits in Idaho (Smith et al. 2021). Lastly, because studies have linked presence of pygmy rabbit burrows to areas of reduced ground litter and microbiotic crust (e.g., Weiss & Verts 1984; Himes & Drohan 2007; Edgel et al. 2014) or reduced herbaceous cover (Gabler et al. 2001), we also included percent litter and percent bare ground variables to reflect these ground covers from the updated Rangeland Fractional Component dataset (Rigge et al. 2019). 

Topography. – We included topographic features in the regional models because they also can influence animal habitat, usually by altering soil deposition, vegetation composition and growth, thermal environments, and properties of precipitation (Appendix S1: Table S1). We resampled a 1/3 arc-second, 10-m resolution national elevation dataset (NED; USGS 2019) to create 30-m resolution raster data layers for elevation, slope, and curvature using Topography Tools (ArcGIS 2010; Dilts 2015) in ArcMap 10.7.1. We created a data layer for a topographic position index (TPI) by calculating the normalized difference between elevation at a central pixel and the surrounding average elevation, and we selected 500m focal radius to represent slope position and general landforms (Weiss 2001). We also estimated an index of terrain roughness at a 200-m scale by calculating the standard deviation of elevation to represent terrain ruggedness associated with the approximate size of an individual’s home range. 

Soils. – We included soil properties because pygmy rabbits are obligate burrowers, and relatively deep, loamy soils that retain integrity of burrows are associated with their presence (Green & Flinders 1980). Soil depth and texture also serve as indicators for ecological conditions that influence potential vegetation, and relatively deep, well-drained soils in areas used by rabbits also support greater shrub growth (Weiss & Verts 1984). We included soil parameters in the models from the POLARIS soils database (Chaney et al. 2016).  We calculated mean values for two soil depth bins (0–30 and 30–100 cm) for eight soil characteristics (bulk density, pH, organic matter percentage, clay percentage, sand percentage, silt percentage, saturated water content (Theta-S), and depth to restrictive layer; Appendix S1: Table S1).

Climate. – Pygmy rabbits are small-bodied endotherms that can be affected by climate directly because seasonal thermal extremes are often outside of their thermal neutral zone (Katzner & Parker 1997; Milling et al. 2018). Pygmy rabbits do not migrate or shift space use seasonally, but instead use burrows and above-ground micro-sites as thermal refuges during both winter and summer seasons (Milling et al. 2018). Climate also can influence the distribution of pygmy rabbits indirectly through effects on vegetation and soil development.  We modeled the influence of temperature and precipitation by including 19 bioclimatic variables (Hijmans et al. 2005) from long-term datasets describing average conditions during 1981–2010 (PRISM Climate Group, 2012; Appendix S1: Table S1). These variables are commonly used to model the bioclimatic conditions shaping species distributions (e.g., Elith et al. 2006; Anderson & Gonzalez 2011; Stanton et al. 2012). 

Location data

We acquired occurrence records for pygmy rabbits collected during 2000–2019 within the five regional extents across southern Idaho and portions of Oregon, Nevada, Utah, and Wyoming from multiple sources (Fig. 1). Most of the data were obtained from the Idaho Species Diversity Database (ISDD: https://idfg.idaho.gov/species/). We used additional locations from recent field surveys conducted by BLM biologists in southwestern Idaho (Region 1). Location data from neighboring states were provided by individual state agencies and used in a rangewide modeling effort for this species (Smith et al. 2019b). We vetted all data to ensure reliability and retained only records for which species identification was confirmed via photographs, sightings by biologists, collection of field specimens, or confirmed pellets at burrow sites. The verified occurrences were randomly subsampled with a minimum distance of 250m using SDMToolbox (Brown et al. 2017), which is the approximate diameter of a female home range (Estes-Zumpf & Rachlow 2009),. Lastly, we removed occurrences within the perimeter of any fire that burned during 2000-2019. In total, 1,584 presence locations were included within the five regional extents (n = 421, 103, 167, 159, 734 in Regions 1-5, respectively). We randomly generated 10,000 background points, or pseudo absences, with a minimum separation distance of 250m within each region to model potential habitat (Phillips and M. Dudík 2008; Barbet-Massin et al. 2012). Individual background points were eliminated if land cover was classified as urban, agriculture, open water, or masked data; or if they occurred within 250m of presence locations. The final number of background points for each region ranged from 9,053 to 9,904.

Model development and testing

We modeled predicted habitat for pygmy rabbits within the five ecological regions using maximum entropy modeling (Maxent 3.4.0, Phillips et al. 2006). This inductive approach uses species presence data, randomly generated background locations, and environmental information in the form of spatial raster data to estimate the relative probability of occurrence across the modeled region based on similarity with habitat conditions at known locations. 

Prior to running the Maxent models, we optimized model parameters (feature type and regularization multiplier) using the R package dismo (Hijmans et al. 2020) and a group of packages developed by Adam Smith (enmSdm v0.3.4.6, statisfactory v0.1.7, omnibus v0.3.0.3, and legendary v0.0.0.9; Smith 2019) to promote model fit and reduce model complexity. The default set of feature types, which are mathematical transformations of the variables used in the modeling process that constrain the variable response curves, in Maxent are linear, quadratic, product, threshold, and hinge (Elith et al. 2011; Merow et al. 2013). We used the occurrence data to determine which combination of feature types was optimized to fit the relationship between occurrences and environmental characteristics. Regularization multipliers penalize the inclusion of parameters that result in little or no “gain” to the model (Merow et al. 2013). We tested a range of multipliers (values 0.5 to 5 in increments of 0.5, and values 6 to 20 in increments of 1) with different feature types to identify the best-performing combination based on Akaike Information Criteria corrected for small sample size (AICc; Warren & Seifert 2011; Wright et al. 2015).

We modeled the distribution and suitable habitat for pygmy rabbits independently for each region. Using the optimized enmSdm parameters, we created full models starting with 53 environmental variables (Appendix S1: Table S1). We selected the Cloglog output, which scales the model results (i.e., probability of species presence) to values between 0 and 1 (Phillips & Dudik 2008). We assessed model performance using a 10-fold cross-validation, in which a random subset (20%) of locations was withheld for model testing. We used the area under the receiver operating curve (AUC) to evaluate threshold-independent model fit using the withheld samples to assess model success in discriminating between presence locations and background points (Elith et al. 2011; Merow et al. 2013).  

Once the full models were completed, we analyzed the contributions of each environmental variable to the model for each region and iteratively removed highly correlated variables and those that did not contribute to the model. Permutation importance (PI) is a relative value that measures the contribution of a single variable to the full strength of the model. When the PI for a variable was low (<2%), we removed the variable from the model for that region; we chose this threshold to be consistent with the model-building process used to create the rangewide model for pygmy rabbits (Smith et al. 2019b). To accomplish a reduction of correlated variables, we constructed a pairwise correlation matrix to identify pairs that were highly correlated (>0.8) and eliminated the variable with the lower PI. We chose a relatively high threshold because Maxent is robust to collinearity (Elith et al. 2011; Feng et al. 2019). We iteratively reduced variables, reoptimized model parameters, and reran models until correlation and contribution criteria were reached. This interactive process resulted in a final “reduced” model for each region that best fit occurrence data and habitat conditions within the regional extent. Twenty-two variables were not retained in any of the final regional models (Appendix S1: Table S1).

We used variable response curves to interpret how environmental variables influenced the predicted distribution of pygmy rabbits and their habitat. Two sets of response curve plots are produced by Maxent: one set displays the predicted habitat suitability as a function of one environmental variable in the model, while setting all other environmental variables in the model to the average sample values. The second set of curves represents a model that only uses a single environmental variable, showing the marginal effect of changing exactly one variable, whereas the plot including the average values of all other variables may take advantage of sets of variables changing together. Both plots reflect the dependence of predicted presence (or habitat suitability) both on the selected variable and on dependencies induced by correlations between the selected variable and other variables (Phillips et al. 2016).  

Mapping predicted habitats

Thresholds are needed to create and map classes of habitat suitability, and we used two types of thresholds to categorize habitat into categories of nonhabitat, suitable habitat, and primary habitat. To distinguish nonhabitat from suitable habitat, we chose the value that maximized the sum of testing sensitivity plus specificity (MTSS), which is often used for presence-only SDMs because it is relatively unaffected by changes in the ratio of occurrences to background points (Liu et al. 2005). For the higher threshold separating suitable from primary habitat, we used the average predicted value of known occurrence locations. Consequently, all primary habitat was predicted to be the same or higher quality than the average value in places where the species was confirmed to be present. This threshold was chosen for consistency in comparing our regional models to a rangewide model and because literature supports this “conservative” approach (Liu et al. 2016). We evaluated correlations between the size of the regions and percent of the region predicted as suitable or primary habitat using Spearman rank correlation tests. Although we did not have an independent data set to test final model accuracy, we used the threshold-dependent True Skill Statistic (TSS; Allouche et al. 2006) to compare the ability of each model to distinguish suitable and unsuitable habitat.

Projections under future climate scenarios

To explore how potential changes in climate might influence the distribution of pygmy rabbits and their habitats, we projected the final model for each region using climate variables from mid-century (2050s) climate projections. This approach assumes that both species-environment relationships and relationships among environmental variables are consistent through time. We used projections from an ensemble of 15 General Circulation Models (GCMs) of the Coupled Model Intercomparison Project Phase 5 (CMIP5) given two greenhouse gas emission scenarios (Representative Concentration Pathways [RCPs] 4.5 and 8.5) (AdaptWest Project 2015, Wang et al. 2016). A stabilization scenario, RCP 4.5, assumes that radiative forcing level conditions will stabilize at 4.5 Watt/m² before 2100 (Taylor et al. 2012). In contrast, RCP 8.5 represents a more extreme emissions scenario in which the radiative forcing level reaches 8.5 Watt/m², which is typical for projections in the literature resulting in high greenhouse gas concentration levels (Taylor et al. 2012). The dataset was downloaded from the AdaptWest Project (2015) in Lambert Conformal Conic projection, at 1-km resolution, covering North America. The original dataset consisted of 48 monthly temperature and precipitation variables; we used the R package dismo (Hijmans et al. 2020) and the function “biovars” to calculate 19 projected bioclimatic variables (Appendix S1: Tables S1, S2). Both sets of future climate projections (RCP 4.5 and RCP 8.5) were reprojected to Albers equal-area conic projection, NAD 1983 Contiguous United States, resampled to 30-m resolution, and clipped to the five regional extents. 

We compared results of our regional models with predictions generated by a rangewide model for pygmy rabbits (Smith et al. 2019b). Using the same methods to identify thresholds, we quantified geographic overlap for both suitable and primary habitats predicted by each model (i.e., rangewide and regional) within the five regions. In addition, we mapped areas of overlap and non-overlap between the regional and rangewide models to examine differences in spatial distributions of model results. 

We used a similar approach to compare predictions of the current and future regional models. Under future climate scenarios, we expected a reduction in the amount of suitable habitat and a shift in the distribution of habitat within regions based on projected warmer summer temperatures and reduced winter precipitation across Idaho (Abatzoglou et al. 2021). Given uncertainty in climate projections, we examined overlap between current and potential future habitat by combining the suitable and primary habitat categories. We used the same MTSS threshold to separate nonhabitat and habitat, then identified areas of habitat contraction, persistence, and expansion under the two future climate scenarios. 

The data products are raster files that can be opened with GIS-capable software (e.g., Esri ArcGIS, QGIS, R).