DOI: 10.5061/dryad.5hqbzkhbk
title: Habitat distribution models for pygmy rabbits in Idaho

description: 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.
description: 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 &gt;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 (<em>Artemesia</em>
 <em>tridentata</em>) and green rabbitbrush
 (<em>Chrysothamnus</em> <em>viscidiflorus</em>)
 are common shrub species at lower elevations, especially along drainages
 and alluvial fans. On the lower elevation flats, isolated areas of low
 sagebrush (<em>A. arbuscula</em>), bitterbrush
 (<em>Purshia</em> <em>tridentata</em>), shadscale
 (<em>Atriplex</em> <em>confertifolia</em>), and
 spiny greasewood (<em>Glossopetalon</em>
 <em>spinescens</em>) also are common. The dominant tree
 species within the sagebrush-steppe ecosystem are western juniper
 (<em>Juniperus</em> <em>occidentalis</em>) and
 Utah juniper (<em>J. osteosperma</em>), but there are also
 populations of pinyon pine (<em>Pinus</em>
 <em>monophylla</em>) scattered throughout. Many of the
 sagebrush-steppe habitats in southern Idaho include an understory of mixed
 native grasses (e.g., bluebunch
 wheatgrass, <em>Pseudoroegneria</em>
 <em>spicata</em>; Sandberg’s bluegrass,
 <em>Poa</em> <em>secunda</em>; bottlebrush
 squirreltail, <em>Elymus</em> <em>elymoides</em>;
 Indian ricegrass, <em>Achnatherum</em>
 <em>hymenoides</em>), seeded perennial grasses (e.g., crested
 wheatgrass, <em>Agropyron</em>
 <em>cristatum</em>), and annual invasive species (e.g.,
 cheatgrass, <em>Bromus</em>
 <em>tectorum</em>).    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).  <em>Land cover.</em> – 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).  <em>Topography</em>. – 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.  <em>Soils</em>. – 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).
 <em>Climate</em>. – 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 <em>dismo</em> (Hijmans et al. 2020) and a group of packages developed by Adam Smith (<em>enmSdm</em> v0.3.4.6, <em>statisfactory</em> v0.1.7, <em>omnibus</em> v0.3.0.3, and <em>legendary</em> 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 <em>enmSdm</em> 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 (&lt;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 (&gt;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<sup>2</sup> 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<sup>2</sup>, 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 <em>dismo</em> (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. 
description: The data products are raster files that can be opened with
 GIS-capable software (e.g., Esri ArcGIS, QGIS, R).