Habitat is the ensemble of resources and conditions that allow an organism to survive and reproduce in a location (Hall et al. 1997). Habitat selection is the behavioral process whereby animals select the habitats they use (Johnson 1980; Hutto 1985). Habitat selection is a foundational concept in ecology because it provides understanding about the essential nature of an organism (Johnson 1980; Hutto 1985; Morris 2003). Habitat selection results in the disproportionate use of one resource over other available resources (Hall et al. 1997; Krausman 1999) and is evaluated by comparing the use of particular resources relative to their availability (Manly et al. 2002). Central to the concept of habitat selection is that the behavioral process is hierarchical across different spatial and temporal scales (Johnson 1980; Mayor et al. 2009). Therefore, evaluating habitat requirements in a multiscale habitat selection framework can provide insight into important habitat characteristics related to an animal's fitness (Uboni et al. 2017).

In contrast, habitat use describes the way an animal uses its environment. Habitat use is not the same as habitat selection because it does not compare resource use versus availability (Krausman 1999). Therefore, habitat use has no implicit consequence to fitness because used resources cannot be interpreted as habitat preferences or requirements (Krausman 1999). As an example, an animal crossed a road, indicating the road was used, but we do not know if the road is an important component of habitat because we do not know the availability of roads on the landscape. However, an analysis of habitat selection showed that while used, the animal avoided the road relative to its availability, suggesting that the road has a negative influence on fitness.

A review of the mammalogical literature, such as those cited in Mammalian Species accounts, demonstrates that available information about species' habitat typically consists of anecdotal descriptions of places where the species has been observed, or infrequent analyses of habitat use. The main exception is for well-studied large mammals of management interest where studies of habitat selection are more common. McGarigal et al. (2016) reviewed > 800 habitat selection papers and < 3% were for small mammals (rodents and shrews). Problems associated with anecdotal descriptions and analyses of habitat use include: imperfect detection (Neu et al. 1974), focus on accessible or well-studied places (Meyer et al. 2015), observer bias (Schooley and McLaughlin 1992), failure to acknowledge scale (e.g., microhabitat descriptions used to characterize landscape-level habitat), and no link to fitness (Uboni et al. 2017). The litany of anecdotal descriptions and studies of habitat use for small mammals usually are interpreted incorrectly as preference, requirements, or quality (Hall et al. 1997). Consequently, it is possible we have a misunderstanding of habitat requirements for the majority of mammal species, which can limit effective conservation or management actions by misdirected efforts (e.g., Peek et al. 1982; Rettie and Messier 2000; Dussault et al. 2005; Bowyer and Kie 2006). To understand a species' ecology or implement meaningful conservation actions, it is necessary to know the resources preferred by the organism that influence their fitness. This cannot be determined using descriptions of habitat or simple analyses of habitat use.

North American chipmunks (Neotamias and Tamias) are a prime example of the lack of studies on habitat selection for small mammal species. Chipmunks are one of the most speciose group of rodents in North America (Burgin et al. 2018), but most literature on habitat for these species is anecdotal or focuses on habitat use. We are aware of only six published studies that evaluated habitat selection in North American chipmunks, none of which evaluated multiple scales. Three studies evaluated microhabitat selection (i.e., third-order selection) in the eastern chipmunk (T. striatus—Geier and Best 1980), the Colorado chipmunk (N. quadrivittatus—Rivieccio et al. 2003), and Palmer's chipmunk (N. palmeri) and the Panamint chipmunk (N. panamintinus—Lowrey and Longshore 2013). Three studies evaluated habitat selection at the landscape scale (i.e., first-order selection) in N. panamintinus (Lowrey et al. 2016) and a subspecies of the Colorado chipmunk, the Oscura Mountains Colorado chipmunk (N. q. oscuraensis—PerkinsTaylor and Frey 2018, 2020). It is possible that mammalogists have a misunderstanding about the habitat needs of many North American chipmunk species because of the lack of habitat selection studies and complete absence of multiscale habitat selection studies. These misunderstandings can influence studies aimed at other aspects of the ecology or evolution of these chipmunk species, particularly any tied to questions about distribution, abundance, or fitness (Uboni et al. 2017).

Neotamias quadrivittatus occurs from Colorado and Utah south to Arizona and New Mexico (Hall 1981). Most of the literature describes it as generally occupying the montane coniferous forest biotic community, which is dominated by large conifer trees (e.g., ponderosa pine [Pinus ponderosa] and Douglas-fir [Pseudotsuga menziesii]—Dick-Peddie 1993), although anecdotal descriptions of habitat suggest a wider array of land cover associations (Best et al. 1994). At the southern distributional limits in the Chihuahuan Desert of New Mexico, there are two subspecies: N. q. oscuraensis and the Organ Mountains Colorado chipmunk (N. q. australis). The state of New Mexico lists these populations as threatened due to their relictual distribution and threats to habitat (New Mexico Department of Game and Fish 2018). Available land cover types in these mountain ranges differ from the typical land cover associations for N. quadrivittatus. For example, the Oscura Mountains contain piñon (Pinus edulis)–juniper (Juniperus monosperma) woodland yet entirely lack montane coniferous forest, whereas in the Organ Mountains, montane coniferous forest is poorly developed and limited to small, scattered stands of ponderosa pine and Douglas-fir. Thus, due to the absence of reliable information about habitat selection at multiple scales, and lack of transferability of most existing information to these populations, it has not been possible to develop sound conservation plans aimed at protecting or enhancing habitat conditions.

Our goal therefore was to evaluate habitat selection by N. q. australis at multiple scales to determine whether previous literature and findings were supported in a multiscale habitat selection context. We accomplished this by using radiotelemetry to identify locations used versus randomly sampled available locations. We hypothesized that previous descriptions of habitat in the literature would be upheld in a multiscale habitat selection context. Our objectives were to: 1) estimate home range size; 2) analyze habitat selection at three spatial scales (second order, third order, and microhabitat); 3) interpret results relative to known habitat information about N. quadrivittatus; and 4) provide management recommendations based on our findings.

Materials and Methods

Study area.—We studied N. q. australis (hereafter, chipmunks) in the Aguirre Springs Recreation Area (32°21′21.60″N, 106°33′43.20″W) of the Organ Mountains–Desert Peaks National Monument, Doña Ana County, New Mexico from October 2018 to July 2019 (Fig. 1). The Organ Mountains are an isolated, steep mountain range characterized by abrupt rock outcrops that rise to 2,708 m (United States Geological Survey 2017). Elevations in the study area ranged from 1,689 to 2,333 m. Over a 30-year period (1981 – 2010), mean daily summer temperature was 24.7°C, mean daily winter temperature was 6.6°C, and mean annual precipitation was 404 mm (PRISM Climate Group 2019). Summers were lengthy and hot, with monsoonal rains beginning in July, and winters had brief periods below freezing. Dominant shrub and tree species at lower elevations included mountain mahogany (Cercocarpus montanus), gray oak (Quercus grisea), and alligator juniper (Juniperus deppeana). At higher elevations, dominant trees included Gambel oak (Quercus gambelii) mixed with occasional ponderosa pines (P. ponderosa). Riparian corridors existed within steep-sided drainages (i.e., arroyos) that intermittently contained water, and the dominant vegetation included velvet ash (Fraxinus velutina), desert willow (Chilopsis linearis), and, at lower elevations, Fremont's cottonwood (Populus fremontii). The Bureau of Land Management managed the land and prohibited cattle from grazing within the study area. Hiking trails frequently used by people and their dogs were located in the study area.

Fig. 1.

Location of the study area in the Organ Mountains, Doña Ana County, New Mexico, United States. Coordinates are in UTM Zone 13S. Distributions of the four land cover types where available locations were drawn for the second-order and third-order scale analyses and locations used for the microhabitat scale analysis are illustrated.

Overview of study design.—We used a use versus availability design to investigate habitat selection at three spatial scales: second order, third order, and microhabitat (Tables 1 and 2). In our study, the second order represented selection of the combined home ranges of animals in the population, and we analyzed selection using spatial variables (i.e., study design II of Manly et al. 2002). We defined used locations as telemetry locations for all chipmunks. We defined available locations as random points drawn from a buffered cumulative minimum convex polygon (MCP—Mohr 1947) home range, which was estimated with telemetry locations buffered by the maximum movement distance (175 m) of chipmunks in our study and was created in ArcMap 10.7.1 (ESRI 2019; Redlands, California). We spatially rarified locations used by 10 m to reduce spatial autocorrelation (n = 350). We randomly generated 700 available locations within the buffered cumulative MCP home range. We used logistic regression with binomial distribution and logit link functions to test 120 hypothesis-driven a priori conceptual models ( Supplementary Data SD2 (gyab071_suppl_supplementary_data_sd2.pdf)). We did not use mixed-effects models because preliminary analyses indicated no substantial variation in use among individuals in the population.

Third order represented selection of areas within the home ranges of individual known chipmunks, and we analyzed selection using spatial variables (i.e., study design III of Manly et al. 2002). At the third-order scale, we defined use from telemetry locations of individual chipmunks and defined available locations as randomly generated points within the 100% MCP home range for each chipmunk. The ratio of used:available points was 1:2. We tested 75 hypothesis-driven a priori models using mixed-effects logistic regression with binomial distribution and logit link function and included a random intercept for each chipmunk ( Supplementary Data SD2 (gyab071_suppl_supplementary_data_sd2.pdf); Breslow and Clayton 1993; Gillies et al. 2006). We included a random effect at this scale to allow for the explicit identification of individual variability within the home range (Neter et al. 1996).

Table 1.

Variables included in the second- (home range) and third- (within home range) order analyses of habitat selection by Neotamias quadrivittatus australis in the Organ Mountains, New Mexico, based on radiotelemetry of captured individuals (n = 20), October 2018 to July 2019.

Table 2.

Microhabitat variables collected in the field at locations where Neotamias quadrivittatus australis was observed (n = 56) and paired random locations (n = 56) in the Organ Mountains, New Mexico, June 2018 to July 2019.

In our study, microhabitat represented selection of areas within the home ranges of unknown chipmunks and was analyzed with field-collected and 10-m resolution spatial variables (study design III of Manly et al 2002). At the microhabitat scale, we defined used locations as the area within the immediate vicinity (15 m) of an opportunistically observed chipmunk's location and defined available locations as random locations selected within a 93 m radius of the chipmunk's location, which represented the radius of a 2.7-ha home range for N. quadrivittatus (Bergstrom 1988). We opportunistically observed chipmunks at 56 locations and the ratio of used:available points was 1:1. At the microhabitat scale, we used logistic regression with binomial distribution and logit link function to test 32 hypothesis-driven a priori models ( Supplementary Data SD2 (gyab071_suppl_supplementary_data_sd2.pdf)).

We assessed selection of land cover associations at second- and third-order scales using Manly selection ratios (Manly et al. 2002) with the R package adehabitatHS (Calenge 2006). At the second-order scale, we compared the proportion of telemetry locations of all chipmunks in each land cover type to the proportion of each land cover type within the buffered cumulative home range. At the third-order scale, we compared the proportion of telemetry locations in each land cover type to the proportion of land cover type available within the home range of each individual chipmunk.

Chipmunk capture and radiotelemetry.—We captured chipmunks with Sherman live traps (model LFATDG; H. B. Sherman, Tallahassee, Florida) baited with a commercial horse feed mixture of grains and molasses. Upon capture, we transferred the chipmunk to a zippered, mesh handling bag for processing and radiocollar attachment. Total handling time of captured chipmunks was < 5 min to minimize stress.

During processing, we recorded body mass by weighing the chipmunk in the bag with a hanging spring scale (Pesola spring scales, PESOLA Präzisionswaagen AG, Schindellegi, Switzerland). We determined sex by visual observation of genitals and assessed reproductive condition as scrotal or nonscrotal for males and open, closed, pregnant (i.e., swollen abdomen), or lactating for females (Schulte-Hostedde et al. 2002). We classified age as juvenile or adult; we considered chipmunks with pelage appearing fluffier, with a body mass ≤ 50 g, and not displaying signs of reproductive activity to be juveniles. We fitted a 1.8-g radiotransmitter collar (model BD-2C; Holohil Systems, Carp, Ontario, Canada) to adult chipmunks that did not appear stressed and appeared to be in healthy physical condition. Adult chipmunk body mass was 62.5 g ± 8.43 SD (range = 50 – 84 g), so the collar never exceeded the recommended 5% of the chipmunk's body weight (Sikes et al. 2016). Prior to release, we held radiocollared chipmunks in Sherman traps for ∼30 min and then examined them to ensure we properly fitted the collar and each chipmunk had made a full recovery. We recaptured nine chipmunks midway through the study period, at which time we replaced the radiocollar. Capture and handling methods followed recommendations of the American Society of Mammalogists (Sikes et al. 2016) and were approved by the New Mexico State University Animal Care and Use Committee (IACUC 2018-006) and the New Mexico Department of Game and Fish (scientific collecting permit #2868).

We located radiocollared chipmunks with handheld telemetry receivers (model R-1000; Communications Specialists, Orange, California) attached to a 3-element Yagi antenna (Wildlife Materials International, Inc., Murphysboro, Illinois). We recorded observer location with handheld GPS units (model GPSMAP 64st; Garmin, Olathe, Kansas). While recording location data, we made an effort not to disturb the chipmunks or their habitat. We located each chipmunk in random order via triangulation at least twice a day and 5 days/week for the duration of radiotransmitter battery life (∼12 weeks). We ascertained compass bearings for a chipmunk from at least two observer locations based on the direction of the strongest telemetry signal. Observers recorded their location and compass bearing within 1 min to reduce potential error due to moving chipmunks. At times when an observer saw a chipmunk without disturbing it, the observer recorded the exact location of that chipmunk with a GPS unit after the chipmunk had left the vicinity. To maintain temporal independence for each chipmunk, we waited ≥ 1 h before obtaining a subsequent location (Brzeziński et al. 2019).

We estimated chipmunk locations using the software Location of a Signal (Ecological Software Solutions LLC 2019), which estimates the point of intersection for two or more locations with bearings. We evaluated the extent and sources of error associated with radiotelemetry to exclude excessively imprecise locations (Withey et al. 2001). We measured the radiotelemetry error using radiotransmitters placed at known locations and developed an equation to predict linear error (i.e., distance between the true position of the radiotransmitter and the triangulated location) following methods of Withey et al. (2001;  Supplementary Data SD1 (gyab071_suppl_supplementary_data_sd1.pdf)). We used the predicted linear error of each telemetry location to define valid telemetry locations. Valid locations included visual observations of identified individuals and triangulated points with a linear error ≤ 30 m ( Supplementary Data SD1 (gyab071_suppl_supplementary_data_sd1.pdf)). Invalid locations were triangulated points with a linear error > 30 m, nonvisual chipmunk locations with only one point, visual observation of unidentified individuals, or triangulation locations that failed to intersect. We only included valid telemetry locations in statistical analyses at the second- and third-order scales.

Home range analysis.—We estimated home range size based on valid telemetry locations using the MCP (Mohr 1947) and 95% fixed kernel density estimator method (KDE—Worton 1989; White and Garrott 1990). For 95% KDE, we used the ad hoc smoothing parameter (h) because reference bandwidth and least squares cross-validation can result in a Type I error by overestimating home range sizes based on small samples (Kie 2013). We used the R package adehabitatHR for home range calculations (Calenge 2006). We determined the number of telemetry points needed to estimate home range using a home range area curve where the asymptote of the curve was considered to represent the number of locations needed to accurately estimate home range size (Odom and Kuenzler 1955; Haines et al. 2006). We evaluated differences in home range size by sex using an equal variance t-test. We tested for normality by sex using the Shapiro–Wilk normality test and tested for equal variance between sexes using an F-test.

Collection of spatial variables.—We used spatial variables created in a geographic information system (GIS) that were selected based on their hypothesized importance to habitat selection, including land cover type, normalized difference vegetation index (NDVI), distance to drainage, aspect, hill shade, slope, vector ruggedness measure (VRM), topographic position index (TPI), and cliff (Table 1). Land cover types that included conifer trees have been associated with increased relative likelihood of occurrence for N. q. australis (Frey and Kopp 2013). We included NDVI as a variable because land cover types with conifer trees have been shown to have the greatest NDVI values in the Organ Mountains (Frey and Kopp 2013) and using NDVI would help distinguish among different cover types (i.e., conifer, nonconifer, and nonvegetated areas). Drainages sometimes hold water and water provides more opportunity for drinking and might increase the production of food-producing plants (Moir and Ludwig 1979; Block and Finch 1997). Previous studies have suggested that N. q. australis is adapted to cool climates and associated vegetation types (Patterson 1980; Frey and Kopp 2013); we therefore hypothesized aspect and hill shade would influence habitat selection. Areas with a southwest-facing aspect and decreased hill shade are expected to have more sunlight across the year, hence will not be as likely to support conifers and other cryomesic vegetation (Moir and Ludwig 1979; Block and Finch 1997). Topographic position can influence ecological characteristics of a site, and N. q. australis is thought to be associated with steep slopes (Rivieccio et al. 2003); we therefore hypothesized that slope, VRM, TPI, and cliff would influence habitat selection.

We created all spatial variables in ArcMap 10.7.1 (ESRI 2019). For each variable (except land cover type and distance to drainage), we calculated average values within a 30-m buffer of the location to account for uncertainty associated with the telemetry locations. We selected 30 m as the buffer based on the mean linear error calculated for the study area ( Supplementary Data SD1 (gyab071_suppl_supplementary_data_sd1.pdf)). We created a map of four land cover types—arid, montane scrub, riparian, woodland—based on 1-m resolution satellite imagery available from the National Agriculture Imagery Program (United States Department of Agriculture, Farm Service Agency 2015) and ground-truthed ground cover data ( Supplementary Data SD3 (gyab071_suppl_supplementary_data_sd3.pdf)). The woodland land cover in our study was a ponderosa pine and Gambel oak association and considered a montane coniferous forest biotic community type (Dick-Peddie 1993). We calculated the proportion of each land cover type within the 30-m buffer of the location. We rescaled land cover proportions from 10 m to 30 m and 90 m resolutions. We used NDVI data available at 30 m resolution from Landsat 8 surface reflectance data (United States Geological Survey 2019). We calculated average NDVI over the entire year using Landsat 8 surface reflectance data available monthly for 2019 (Vermote et al. 2016). We rescaled NDVI from 30 to 90 m. We defined drainages, which represented all perennial and intermittent streams, using the National Hydrography Dataset (United States Geological Survey 2019), and calculated the distance in meters from locations to nearest drainage.

To create topographic variables, we used a 10-m Digital Elevation Model raster available from The National Map (United States Geological Survey 2017). We created each topographic variable at three spatial resolutions by using the aggregate tool to rescale variables from 10 m to 30 m and 90 m resolution or from 1-cell to 3-cell and 9-cell moving window (Table 1). We calculated aspect using the aspect tool. We folded aspect to the northeast–southwest line following the equation in McCune and Keon (2002). To calculate hill shade, we obtained data on the sun's altitude and azimuth from the United States Navy Astronomical Applications Department for 0900, 1200, and 1500 h on the 2019 solstices and equinoxes and calculated the average of these values. We calculated the degree of slope using the slope tool and VRM using the Benthic Terrain Modeler toolbox add-on (Walbridge et al. 2018). We calculated TPI as a continuous variable using the Topography Tools toolbox 10.3 add-on (Dilts 2015). We defined cliff as a pixel that had a TPI of –2 to 2 and a slope > 15°.

Collection of field variables for the microhabitat scale.—For the microhabitat scale analysis, we collected data on microhabitat characteristics that could be important for foraging or reducing predation risk. We included all of the spatial variables (except folded aspect, cliff, and NDVI) at a 10-m spatial scale and a suite of field-collected variables. Previous microhabitat use studies found N. q. australis to be associated with increasing litter cover and decreasing shrub and grass cover (Rivieccio et al. 2003), rocks (Sullivan 1996; Rivieccio et al. 2003), and burned areas (Johnson et al. 1998). We measured microhabitat variables along four perpendicular 15-m transects radiating from each location. Variables included tree and shrub species composition and abundance, distance to nearest shrub, tree, and boulder, presence of burned vegetation, and ground cover (Table 2). We defined trees, shrubs, and forbs, using descriptions from the United States Department of Agriculture PLANTS database (United States Department of Agriculture, Natural Resources Conservation Service 2019). Trees were considered woody plants that, at maturity, had a trunk. Shrubs were considered multistemmed woody plants. Forbs were considered vascular plants that lacked significant secondary woody growth. We used the United States Department of Agriculture PLANTS database for plant species names (United States Department of Agriculture, Natural Resources Conservation Service 2019).

We identified to species all trees and shrubs that had trunks or stems located within a 1-m belt on either side of the transects. We recorded the number of each species of shrub by 1-m height classes (< 1, 1 – 2, 2 – 3, 3 – 4, > 4) and the number of each species of tree by 10-cm-diameter size classes (< 10, 10 – 20, 20 – 30, 30 – 40, 40 – 50, > 50) along each transect. We determined diameter at breast height of trees with a measuring tape. We measured distance from the location to the nearest shrub and to the nearest tree and identified them to species. We calculated woody plant diversity with the Shannon–Weaver index (Shannon and Weaver 1948; Margalef 1957). We measured distance from the location to the nearest boulder, which we defined as a rock ≥ 1 m in diameter. We recorded the degree woody plants had been burnt by wildfire: no burn (no evidence of fire), partial burn (< 75% burn coverage), or complete burn (≥ 75% burn coverage) for shrubs, trees, and dead standing trees.

We collected ground cover data every 2 m along the transects. We used standard 20 × 50 cm Daubenmire plot and classing categories (1: 0 – 5%, 2: 6 – 25%, 3: 26 – 50%, 4: 51 – 75%, 5: 76 – 95%, 6: 96 – 100%) to estimate ground cover by type at 1 m above the ground (Bonham et al. 2004). Ground cover types included forbs, grasses, woody plants, leaf litter, fine woody debris, coarse woody debris, bare, and rock by type. We defined dead leaves as leaf litter. We defined fine woody debris as logs or branches with a diameter < 10 cm and coarse woody debris as logs or branches with a diameter ≥ 10 cm (Harmon and Sexton 1996). We classified rock as exposed bedrock, boulders, or small rocks (< 1 m in diameter).

Statisticalanalysis.—We carried out statistical analyses in Program R (R Foundation for Statistical Computing, Vienna, Austria). For the three habitat selection scales, we created scatterplots of each variable to examine potential outliers and data transformations. We performed a univariate logistic regression on each variable and used the R Package MuMIn (Bartoń 2019) to calculate Akaike's information criterion adjusted for small sample sizes (AIC_c) to assess model fit to the chipmunk location data (Burnham and Anderson 2002). For spatial variables calculated at three resolutions (10, 30, 90 m), we retained the resolution with the lowest ΔAIC_c value (Burnham and Anderson 2002). We then removed any variable that had either 85% CI overlapping zero or were less biologically relevant when the Pearson correlation coefficient between two variables was > 0.70 (Arnold 2010). We used 85% CI to ensure we maintained variables that could be included in best-approximating models (Arnold 2010). Based on the scatterplots, we tested for quadratic effects on slope and hill shade at the second-order scale and used AIC_c to assess model fit. We retained the quadratic term if the quadratic model had a lower ΔAIC_c than the model with the linear term. We hypothesized interactions among NDVI, aspect, hill shade, slope, and VRM at the second- and third-order scales. We retained only those interactions that had 85% CI that did not overlap zero. We standardized the remaining continuous variables around the mean and standard deviation to compare the relative influence of resources on habitat selection (Zar 1986).

We calculated Manly selection ratios and 95% Bonferroni-adjusted confidence intervals at the second- and third-order scalestodetermineifproportionofalandcovertypeuseddiffered from proportion available (Johnson 1980; Manly et al. 2002). We considered selection ratios > 1 to indicate selection for a specific land cover type and < 1 to indicate avoidance (Manly et al. 2002). At the third order, we evaluated differences in Manly selection ratios between sexes using equal variance t-tests. We tested for normality by sex using the Shapiro–Wilk normality test and tested for equal variance between sexes using an F-test.

Table 3.

Ranking of logistic regression models for second-order selection (home range) by Neotamias quadrivittatus australis (n = 20) in the Organ Mountains, New Mexico, October 2018 to July 2019. Model variables, number of parameters in the model (K), difference in Akaike's information criterion corrected for small sample sizes (ΔAIC_c), Akaike weights (w_i = estimated probability of model i being the best model given data and model set), and model deviance for candidate models developed to explain differences in selection between used and available radiotelemetry locations. See  Supplementary Data SD3 (gyab071_suppl_supplementary_data_sd3.pdf) for the full model set; only models that cumulatively made up 95% AIC_c weight are included in this table.

For each habitat selection scale, we considered all models with a ΔAIC_c value of < 2 as competitive and selected the top model based on ΔAIC_c (Burnham and Anderson 2002). We calculated variance inflation factors (VIF) using the R package car (Fox and Weisberg 2019) for variables in the top model to assess multicollinearity (James et al. 2014). At the second- and third-order scales, we carried out k-fold cross-validation partitioned into 10 folds to test the predictive capabilities of our top models (Boyce et al. 2002). At the microhabitat selection scale, we undertook k-fold cross-validation partitioned into five folds due to the small sample size. For the second-order scale, we created a map of predicted use based on unstandardized coefficients of variables in our top model (Long et al. 2009).

Results

Home range.—We recorded 1,256 locations from 20 adult chipmunks (n = 10 M, n = 10 F), of which 637 locations were valid (31.85 locations ± 17.38 SD for each animal;  Supplementary Data SD4 (gyab071_suppl_supplementary_data_sd4.pdf)). The home range area curve (100% MCP) appeared to stabilize for most chipmunks at 30 locations ( Supplementary Data SD5 (gyab071_suppl_supplementary_data_sd5.pdf)). We estimated MCP and KDE home range for 10 individuals (n = 5 M, n = 5 F) that had ≥ 30 locations. Home ranges of each sex were normally distributed, and variances were equal between sexes (95% MCP: F_4,4 = 2.81, P = 0.34; 100% MCP: F_4,4 = 0.98, P = 0.98; KDE: F_4,4 = 2.20, P = 0.46). Mean home range did not differ between males and females (95% MCP: t₈ = 0.72, P = 0.49; 100% MCP: t₈ = 0.57, P = 0.58; KDE: t₈ = 1.87, P = 0.10). Estimated mean home ranges for the study population were 95% MCP of 2.55 ha ± 1.55 SD (range = 1.09 – 4.00 ha); 100% MCP of 3.25 ha ± 1.32 SD (1.59–5.26 ha); and KDE of 2.09 ha ± 1.21 SD (0.39–3.96 ha).

Second-order selection.—The only competitive model for second-order selection had proportion of montane scrub land cover at the 90-m scale (MS-90), NDVI at the 90-m scale (NDVI-90), distance to drainage, folded aspect at the 90-m scale (aspect-90), slope² at the 90-m scale (slope-90), and an interaction between NDVI-90 and aspect-90 (Table 3;  Supplementary SD6 (gyab071_suppl_supplementary_data_sd6.pdf)). Probability of selection by chipmunks increased with increasing MS-90, NDVI-90, and aspect-90 (Table 4; Fig. 2). Probability of selection by chipmunks decreased as distance to drainage increased (Table 4; Fig. 2). Probability of selection by chipmunks increased until the slope of an area reached ≥ 20° and then began to decrease (Table 4; Fig. 2). The interaction between NDVI-90 and aspect-90 suggested that chipmunks selected for northeasterly facing slopes but also selected more southwesterly facing aspects if they were greener (Table 4; Fig. 2). Ten k-fold cross-validation indicated that the second-order selection model had strong predictive power (ρ = 0.81, P < 0.0001). The second order predictive map indicated there were very few areas on the landscape with > 25% probability of selection (Fig. 3).

Table 4.

Variables, standardized beta coefficient estimates, SE, 85% CI, and variance inflation factors (VIF) for the best-fitting model for second-order selection (home range) by Neotamias quadrivittatus australis (n = 20) in the Organ Mountains, New Mexico, October 2018 to July 2019.

Fig. 2.

Predicted probability of Neotamias quadrivittatus australis second-order (home range) selection for distance to drainage (m), proportion of montane scrub land cover type, folded aspect (radians; zero represents NE-facing and 3.14 represents SW-facing), and slope (degrees) in the Organ Mountains, New Mexico, October 2018 to July 2019; shading represents 85% CI.

Fig. 3.

Projected second-order (home range) habitat selection function for Neotamias quadrivittatus australis in the Organ Mountains, New Mexico, October 2018 to July 2019 (90 m resolution; n = 20). Coordinates are in UTM Zone 13S. Dark areas represent increase in probability of selection; inset map is the density of valid telemetry locations.

At the second-order scale, the highest Manly selection ratio was for riparian land cover (selection ratio = 2.64). Chipmunks also selected for montane scrub land cover (selection ratio = 1.96) and avoided arid and woodland land cover types (selection ratio = 0.16 for both; Fig. 4A).

Third-order selection.—The only competitive model for third-order selection included NDVI at 90-m scale (NDVI-90), distance to drainage, hillshade at 10-m scale (hillshade-10), slope at 90-m scale (slope-90), and TPI at 10-m scale (TPI-10; Table 5;  Supplementary Data SD6 (gyab071_suppl_supplementary_data_sd6.pdf)). Probability of selection by chipmunks increased with increasing NDVI-90 and slope-90 (Table 6; Fig. 5). Probability of selection by chipmunks increased for strongly negative TPI-10, which we considered as locations within home ranges that were near the base of cliffs in drainages (Table 6; Fig. 5). Probability of selection by chipmunks decreased with increasing distance to drainage and hillshade-10 (Table 6; Fig. 5). The k-fold cross-validation indicated that the top third-order selection model had poor predictive power (ρ = 0.40, P < 0.0001).

Fig. 4.

Selection ratios and Bonferroni-adjusted 95% CI for ten female (F) and ten male (M) radiotracked Neotamias quadrivittatus australis in the Organ Mountains, New Mexico, March–July 2019. Dashed line indicates no selection (i.e., ratios = 1). A) Second-order (home range) selection ratios (n = 20) and B) third-order (within home range) selection ratios for chipmunks with ≥ 30 telemetry locations (n = 10).

Table 5.

Ranking of mixed-effects logistic regression models for third-order selection (within home range) by Neotamias quadrivittatus australis (n = 10) in the Organ Mountains, New Mexico, March–July 2019. Model variables, number of parameters in the model (K), difference in Akaike's information criterion corrected for small sample sizes (ΔAIC_c), Akaike weights (w_i = estimated probability of model i being the best model given data and model set), and model deviance for candidate models developed to explain differences in selection between used and available radiotelemetry locations. See  Supplementary Data SD3 (gyab071_suppl_supplementary_data_sd3.pdf) for full model set; only models that cumulatively made up 95% model weight are included in this table.

Table 6.

Variables, standardized beta coefficient estimates, SE, 85% CI, and variance inflation factors (VIF) for the best-fitting model for third-order selection (within home range) by Neotamias quadrivittatus australis (n = 10) in the Organ Mountains, New Mexico, March–July 2019.

At the third-order scale, there were no differences in Manly selection ratios by sex (t₃₅ = –0.14, P = 0.89). All chipmunks selected for montane scrub land cover (selection ratio range = 1.05–1.40; Fig. 4B), and all chipmunks avoided arid land cover (selection ratio range = 0–0.60; Fig. 4B). Most (8) chipmunks avoided riparian land cover, although one chipmunk selected it and one chipmunk did not either select or avoid it (selection ratio range = 0.00–1.06; Fig. 4B). All seven chipmunks with woodland land cover available avoided it (selection ratio range = 0–0.18; Fig. 4B); three chipmunks did not have woodland land cover type within their home range.

Microhabitat selection.—The only competitive model for microhabitat selection included woody plant diversity, distance to nearest tree, coarse woody debris ground cover, rock ground cover, and TPI at 10 m resolution (TPI-10; Table 7). Probability of selection by chipmunks increased with increasing woody plant diversity, coarse woody debris ground cover, and rock ground cover (Table 8; Fig. 6). Probability of selection by chipmunks increased as TPI-10 values became more negative, reflecting areas near the base of slopes (Table 8; Fig. 6). Probability of selection by chipmunks decreased as the distance to nearest tree increased (Table 8; Fig. 6). The k-fold cross-validation indicated that the top microhabitat selection model had moderate predictive power (ρ = 0.46, P < 0.0001).

Discussion

We found that habitat selection by N. q. australis was fundamentally different from descriptions of habitat use by the species in the literature. Prior studies described N. quadrivittatus as primarily associated with montane coniferous forests (Patterson 1980; Best et al 1994). In contrast, N. q. australis in our study area avoided pine–oak woodland land cover type and selected a suite of characteristics associated with relatively deep, steep-sided arroyos. At the second order, the study population selected particular areas of arroyos that were relatively green, contained more montane scrub, and had moderate slopes. Similarly, at the third order, chipmunks selected locations within home ranges near the base of arroyo sides that were topographically shadier, steeper, and greener. The second-order predictive map demonstrated that selected habitats only occurred in certain areas of arroyos rather than forested slopes, reinforcing the key difference between our habitat selection results and prior descriptions of habitat use. Because we based our results on a comparison of used versus available locations, we can have more confidence in these results and can interpret them as habitat requirements for N. q. australis in our study area.

Our results demonstrated that failure to consider multiple scales could lead to potential misinterpretation of habitat information. For example, based on the second-order Manly selection ratio results, chipmunks selected riparian land cover nearly twice as strongly as montane scrub land cover, suggesting an association with riparian vegetation. Nevertheless, when we included the context of the third-order and microhabitat scales, we know that chipmunks selected for certain conditions of arroyos—shady, steep sides that were more green—rather than riparian vegetation itself. Without the context of additional scales, managers might misapply conservation efforts for N. q. australis by enhancing riparian vegetation in the Organ Mountains. Finally, we demonstrated that extrapolating microhabitat results to a broader, landscape-level scale can lead to misinformation about habitat. At the microhabitat scale, chipmunks selected diverse woody areas with coarse woody debris and rock cover near the base of slopes. Without the information from the coarser scales, we could misinterpret these results and incorrectly conclude that N. q. australis is associated with coniferous forests, as suggested in the previous literature on N. quadrivittatus habitat use. However, given the context of the second- and third-order scales, we know that these microhabitat conditions of rocky sites with downed logs and diverse, woody vegetation near the base of slopes were selected within arroyos. Rivieccio et al. (2003) provided another example of the importance of context when extrapolating microhabitat information. Their study suggested that N. q. australis was positively associated with logs and negatively associated with grass and shrub cover, which could again lead to the inappropriate conclusion of an association with coniferous forests. Nevertheless, that study lacked the ability to draw conclusions from the microhabitat information to the location of chipmunks on the landscape because it did not consider broader scales.

Fig. 5.

Predicted probability of Neotamias quadrivittatus australis third-order (within home range) selection for topographic position index (TPI), distance to drainage (m), hill shade, slope (degrees), and normalized difference vegetation index (NDVI) in the Organ Mountains, New Mexico, March–July 2019; shading represents 85% CI.

Extrapolating results of a habitat selection study to a different geographic area might lead to misinterpretation of habitat information. Prior studies appropriately used occupancy modeling (Perkins-Taylor and Frey 2018) and species distribution modeling (Perkins-Taylor and Frey 2020) to investigate first-order habitat selection by N. q. oscuraensis (Meyer and Thuiller 2006). N. q. oscuraensis selected areas at high elevation with piñon woodland and escarpments (Perkins-Taylor and Frey 2018, 2020). In the Oscura Mountains, these selected habitats are broadly distributed due to the relatively uniform topography of the range. In contrast, the Organ Mountains are much more topographically rugged than the Oscura Mountains, resulting in complex patterns of environmental variation and patchy vegetation. This, and the selection of particular areas of arroyos by N. q. australis, suggested that the distribution of N. quadrivittatus in the Organ Mountains might be more restricted and fragmented than in the Oscura Mountains. For example, piñon woodland is rare in the Organ Mountains and exists only as small relict stands. Although we can only hypothesize what we might expect to find in a first-order model because we did not address first-order selection in this study, this suggests that applying N. q. oscuraensis selection for piñon woodland to N. q. australis habitat information could be incorrect.

Table 7.

Ranking of logistic regression models for microhabitat selection by Neotamias quadrivittatus australis (n = 56) in the Organ Mountains, New Mexico, July 2018 to July 2019. Model variables, number of parameters in the model (K), difference in Akaike's information criterion corrected for small sample sizes (ΔAIC_c), Akaike weights (w_i = estimated probability of model i being the best model given data and model set), and model deviance for candidate models developed to explain differences in selection between used and available locations. See  Supplementary Data SD3 (gyab071_suppl_supplementary_data_sd3.pdf) for full model set; only models that cumulatively made up 95% model weight are included in this table.

Table 8.

Variables, standardized beta coefficient estimates, SE, 85% CI, and variance inflation factors (VIF) for the best-fitting model for microhabitat selection by Neotamias quadrivittatus australis (n = 56) in the Organ Mountains, New Mexico, June 2018 to July 2019.

We can have further confidence that our results are reliable because we tested and incorporated telemetry error, which many radiotelemetry studies ignore (e.g., Harris et al. 1990; Withey et al. 2001; Bartolommei et al. 2012). Telemetry error is study site-specific and can lead to incorrect assignment of locations, thereby influencing habitat selection results (White and Garrott 1990). To control for telemetry error in our study, we excluded locations with excessive error from all analyses, we calculated spatial variables for the second- and third-order selection within a circular area around each telemetry location, and we excluded triangulated telemetry locations from the microhabitat analysis. We also selected scales of study based on biologically meaningful levels, which some habitat selection studies fail to consider (Manning et al. 2004; Bowyer and Kie 2006). We defined our scales based on known home range information for N. quadrivittatus instead of human preconceptions to be certain we did not select arbitrary scales for analysis (Manning et al. 2004; Bowyer and Kie 2006). We also calculated spatial variables at multiple spatial scales and used information theory to determine the most informative resolution for analysis (Thompson and McGarigal 2002). We used statistical analyses that were appropriate for our study design and type of attribute data, which some researchers fail to justify adequately (Thomas and Taylor 2006; Jenkins et al. 2019). We chose logistic regression as our method of statistical analysis because the nature of our response variable was binomial: used or available (Garson 2016). We also met the assumptions of logistic regression because we accounted for outliers by using standardized continuous predictors, for multicollinearity among predictors by using Pearson correlation analysis and VIF, and for independence between observations by spacing observations temporally and rarefying locations spatially.

Conservation implications.—Across most of its distribution, N. quadrivittatus occurs in relatively cool, mesic environments of major mountain ranges (Best et al. 1994). In contrast, the Organ Mountains are a small, isolated mountain range located in the Chihuahuan Desert at the species' southern distributional limits. The Organ Mountains have relatively hot and dry conditions compared with other ranges where the species occurs. We found that N. q. australis selected specific areas of certain arroyos that were relatively deep and steep-sided. These arroyos provide a relatively cool microclimate due to topographic shading. In a post hoc test, we found that locations within the arroyos were cooler than locations outside of the arroyos in all seasons except winter ( Supplemental Data SD7 (gyab071_suppl_supplementary_data_sd7.pdf)). The selected areas of the arroyos also were unique in that they held water more consistently than other areas of the study area, due to occurrence of natural springs and precipitation runoff (Blake et al. 2020). Studies have found that water is important for reproductive success in other species of chipmunks (Heller and Poulson 1970; Hirshfield 1975). Arroyos also might have greater diversity of woody vegetation due to the increased availability of water. These physical and ecological conditions of arroyos allow them to function as a “climate change refugia” by increasing the chances of retaining surface water and milder environmental conditions (Morelli et al. 2016:3). That N. q. australis selected relatively cool, sparsely distributed mesic areas suggests it is existing on the edge of its ecological tolerances in our study area.

Fig. 6

Predicted probability of Neotamias quadrivittatus australis microhabitat selection for coarse woody debris (CWD) ground cover, woody plant diversity, rock ground cover, topographic position index (TPI), and distance to nearest tree (m) in the Organ Mountains, New Mexico, June 2018 to July 2019; shading represents 85% CI.

Despite carrying out our study in the supposed core area for chipmunk distribution in the Organ Mountains (Patterson 1980), we found that N. q. australis selected for a very specific landform that is unique on the landscape. A post hoc analysis of arroyos within the study area found that they made up 22.52% of the landscape ( Supplementary Data SD7 (gyab071_suppl_supplementary_data_sd7.pdf)); however, based on the second-order predictive map, there only were four 90 m pixels with predicted habitat selection > 50% and only three with > 75%. As a result, N. q. australis could have a small and fragmented distribution in the mountain range, thus increasing vulnerability to habitat disturbance and population instability in occupied patches from recreational activities, military activities, and disturbances such as wildfire (Sullivan 1996; Morrison et al. 2006). Given the limited habitat available to this chipmunk, it therefore is imperative to protect pockets of habitat from disturbance and mitigate impacts of climate change.

As a relict taxon that selects refugial habitat, conservation of N. q. australis could focus on protecting arroyos as climate change refugia, following the steps recommended by Morelli et al. (2016). Briefly, beneficial conservation and management plans would be those that focus on maintaining arroyo chipmunk habitat and implementing priority actions such as minimizing disturbance from wildfire and other causes. Future research is needed to determine the climatic buffer that arroyos provide to understand and manage refugia features. Our study did not consider first-order selection; therefore, we are uninformed as to how our habitat selection findings relate to chipmunk distribution across the mountain range. Based on the narrow set of habitat requirements for N. q. australis and the uniqueness of the selected portions of arroyos on the landscape, we hypothesize that first-order selection also will include arroyos and that these conditions will continue to be rare on the landscape. We suggest implementing an occupancy study to examine first-order habitat selection. If the chipmunks are as limited in habitat across their distribution as our study suggests, the conservation stakes for N. q. australis could be grim.