Predicting species distributions is increasingly important in conservation biology and, in the context of contemporary climate change, can be particularly informative for populations at the periphery of the range. Peripheral populations can exhibit unique patterns of habitat use in response to marginal conditions and can provide colonists adapted to novel or extreme environments. We conducted surveys for American pikas (Ochotona princeps) during 2007–2009 in 144 sites in Craters of the Moon National Monument and Preserve, Idaho. This species appears vulnerable to global climate change, and Craters of the Moon is an extensive area of low-elevation lava habitat situated on an interior edge of the species' range. We found pikas to be readily detectable when both direct and indirect sign were used. An estimate of detection probability from a subset of 72 sites that were visited twice was 0.92. We detected pikas in 31% of survey sites overall but only at sites above 1,605 m. We used logistic regression to model the distribution of pikas as a function of elevation, substrate, and vegetation cover. Pikas were most likely to occur on structurally complex pahoehoe lava flows above 1,600 m. The odds of pika occurrence on pahoehoe lava was >10 times that of aa lava flows and increased by 2 times for each SD increase in elevation. Pikas were also more likely to occur on lava flow sites with higher structural complexity and forb cover. An area of pahoehoe lava encompassing >250 km2 in the northern portion of Craters of the Moon contained 91% of pika detections and all predicted site-occurrence probabilities >0.38, an optimal cutoff value determined by examining model receiver operating characteristic curves. Craters of the Moon may provide long-term refugia for the species, given the extent of lava habitat there. However, the importance of elevation in our models suggests that accelerated climate change could erode suitable pika habitat in the park. Most research on pikas has relied on censuses or nonrandom convenience surveys, but we demonstrate an efficient probabilistic sampling approach that has broad application for pika monitoring and research.
Understanding species distributions has been a fundamental pursuit for ecologists, and distribution modeling remains an active area of development (Elith and Leathwick 2009; Guisan and Thuiller 2005; Guisan and Zimmermann 2000). Probabilistic distribution models explicitly embrace uncertainty and provide insights into influences on probabilities of occurrence and persistence over time (Araújo and Williams 2000; Guisan and Zimmermann 2000). Species–environment models of distributions have become fundamental tools for conservation biology and are increasingly important to evaluate the pace and impact of global climate change on ecosystems (Araújo and Williams 2000; Elith and Leathwick 2009; Parmesan 2006; Thomas et al. 2001). Models informed by data from populations at the periphery of a species' range can be particularly useful in this regard (Thomas et al. 2001), and peripheral populations can exhibit unique patterns of habitat use and occurrence in response to marginal conditions that have range-wide conservation implications (Parmesan 2006). Peripheral populations may face higher extinction risks than core populations but can also be important reservoirs for genetic diversity and provide colonists adapted to novel or extreme environments (Channell and Lomolino 2000; Kark et al. 2005; Sexton et al. 2009; Spector 2002).
The talus-dwelling American pika (Ochotona princeps) is distributed in patches in mountainous regions of western North America across broad latitudinal and elevational gradients (Galbreath et al. 2009; Millar and Westfall 2010; Smith and Weston 1990). The species is restricted primarily to habitats associated with permafrost (Hafner 1993, 1994; Millar and Westfall 2010), and evidence of recent extirpation and elevational retreat as a result of warming temperatures and reduced precipitation and snowpack has increased concern for the long-term viability of the species (Beever et al. 2003, 2010; Grayson 2005; McDonald and Brown 1992; Smith et al. 2004). Pikas are thermally sensitive (MacArthur and Wang 1973; Smith 1974), and the species appears to have responded to climate change with rapid range contractions during the Holocene and over the last century (Beever et al. 2003; Galbreath et al. 2009; Grayson 2005; Hafner 1994; Hafner and Sullivan 1995; Moritz et al. 2008). Elevational range contractions in the Great Basin appear to be particularly pronounced (Beever et al. 2003; Grayson 2005; Millar and Westfall 2010). These patterns follow a general trend of relaxation among montane fauna since the last glacial epoch, particularly for species with poor dispersal abilities like the pika (Brown 1971; Lawlor 1998; McDonald and Brown 1992; Smith 1987; Smith and Ivins 1983). The species was recently considered for protection under the United States Endangered Species Act primarily in response to concerns about the impacts of climate change (United States Fish and Wildlife Service 2010).
American pikas are vulnerable to heat stress and have been shown, when caged during experiments, to expire under brief exposure to temperatures as low as 25.5°C (Smith 1974). Pikas modify activity patterns during summer to avoid daytime temperature extremes and spend considerable amounts of time under talus and rock, which provide insulated dens, runways, and areas for food caching (Smith 1974; Smith and Weston 1990). Despite these heat-avoidance behaviors, pika distributions generally are restricted to cool montane environments and are not typically found below 2,500 m in the interior mountain ranges of the western United States (Smith and Weston 1990; but see Beever et al. 2008; Millar and Westfall 2010; Simpson 2009). One exception is the relatively low-elevation lava flows of Craters of the Moon National Monument and Preserve (hereafter, Craters of the Moon) in southern Idaho. These extensive lava flows cover 1,025 km2 and provide some of the most extensive spatially contiguous potential habitat available to the American pika. Craters of the Moon encompasses the entire range of the formerly recognized subspecies O. p. goldmani (Hall 1981), which was recently subsumed into the Northern Rocky Mountains phylogenetic unit, O. p. princeps (Galbreath et al. 2009; Hafner and Smith 2010).
Craters of the Moon is located immediately north of the topographic Great Basin and shares many Great Basin attributes including vegetation and climate. Maximum July air temperatures regularly exceed 28°C (Davey et al. 2006), well above the 25.5°C upper lethal limit for pikas described by Smith (1974). Surface temperatures >60°C have been recorded on the lava (M. Munts, pers. comm.). Precipitation is low, averaging <390 mm annually (Davey et al. 2006), and is lower than what has been reported at many locations where pikas persist (Hafner and Sullivan 1995; Millar and Westfall 2010). Likewise, elevations spanned by Craters of the Moon lava flows trend southward from 1,985 to 1,301 m, well below the average minimum elevations of extant Great Basin populations (Beever et al. 2003; Grayson 2005; Millar and Westfall 2010). These elevations are also >250 m below the lower limits predicted by Hafner (1993) for the geographic coordinates of Craters of the Moon. Craters of the Moon, isolated on the edge of the Snake River Plain, is clearly at an elevational range periphery for the species (Galbreath et al. 2010; Hafner and Sullivan 1995), and the persistence of pikas there remains enigmatic. Continued persistence of the species in Craters of the Moon is of immediate concern given the increasing rate of putative extirpations observed in the Great Basin and the isolation of the area from other suitable habitat.
An understanding of factors influencing pika distribution and habitat use is a prerequisite to evaluating the long-term viability of the species at Craters of the Moon. More broadly, insight into the factors controlling pika distribution can be gained by investigating the species at Craters of the Moon because at this site pikas appear to be persisting outside of the typical climatic envelope and in an unusually large and contiguous complex of potential habitat. Island size has been identified as a primary factor influencing faunal extinction risk in the Great Basin (Brown 1971), and the size and spatial configuration of talus habitat patches influence population dynamics of pikas in montane environments (Franken and Hik 2004; Moilanen et al. 1998; Peacock and Smith 1997; Smith and Gilpin 1997). Habitat attributes related to patch quality, such as substrate type and vegetation cover, are important for the collared pika (Ochotona collaris) and other organisms (Fleishman et al. 2002; Franken and Hik 2004) but generally have not been considered in studies of the American pika. Hafner (1994) and Millar and Westfall (2010) have demonstrated the importance of glacial and periglacial landforms, such as rock–ice features and permafrost, to pika occurrence patterns in montane environments. The United States Fish and Wildlife Service decision not to list the species under the Endangered Species Act was due in part to a recognition that pikas are able to seek out subsurface thermal refugia. Understanding how substrate, vegetation, and other site attributes influence pika distribution and habitat selection is therefore important across a wide range of environments. Most of the available information about the American pika has been collected from classic high-elevation montane environments, but recent reports on low-elevation populations (Beever et al. 2008; Millar and Westfall 2010; Simpson 2009) indicate the need for information from a broader range of environmental conditions.
The objectives of our study were to identify key environmental factors associated with the occupancy patterns of pikas and test a predictive model of the species' distribution across Craters of the Moon. Using presence–absence surveys and an a priori suite of logistic regression models, we evaluated the strength of evidence for competing models that included various combinations of elevation, substrate, and vegetation predictors. These predictors represented resource (vegetation) and indirect (elevation, substrate) ecological gradients (Guisan and Zimmerman 2000), with the former providing food and the latter serving as a proxy for the direct gradients of temperature, precipitation, and snowpack presumed to act on the physiological limitations of pikas. We generated a distribution map of predicted occurrence probabilities and evaluated model performance with a holdout sample. Finally, we summarized the amount and spatial configuration of lava flow habitat according to predicted occurrence probabilities and considered the implications of this summary for long-term persistence of pikas in Craters of the Moon and elsewhere.
Materials and Methods
Study area
We conducted our study within the boundary of Craters of the Moon National Monument and Preserve, located in the Snake River Plain of southern Idaho (42°39′–43°35′N, 113°9–58′ W; Fig. 1). The park consists of extensive basalt lava fields (Fig. 1) that cover approximately 1,025 km2. These fields consist of numerous lava flows that erupted from a series of fissures along a central north–south-trending rift and that range in age from approximately 2,100 to 15,000 years. The lava flows in Craters of the Moon span an elevational gradient of 684 m that runs from 1,985 m at the northern edge down to 1,301 m near the Snake River in the south. The climate is semiarid, and weather patterns predictably reflect the elevational gradient (Davey et al. 2006). Most precipitation occurs during winter, often falling as snow in the northern portion of the area (Davey et al. 2006). Summer droughts extend from June through September. Thirty-year mean annual precipitation reported from the northern end of the study area is 390 mm (weather station 102260), but only 250 mm on the south end of the study area (weather station 105972). Average summer precipitation (June–August), a metric found by Hafner and Sullivan (1995) to be associated with pika genetic diversity, is 66 mm in the north and 38 mm in the south, lower than most extant sites studied by Hafner and Sullivan (1995). Thirty-year average January temperatures are −7°C in the north and −5°C in the south. Thirty-year average July temperatures are 20°C in the north and 20.5°C in the south. July maximum temperatures average 29°C in the north and 30°C in the south. January snowpack averages 508 mm in the north and only 25 mm in the south.
Figure 1
Map of Craters of the Moon National Monument and Preserve, Idaho, including the 2 major lava fields that comprise the National Park Service (NPS) lands and the 4 sampling frames used in the study. The inset shows the location of Craters of the Moon in relation to the Snake River Plain and the state of Idaho.
Vegetation on the lava flows is highly variable, depending on the age and type of flow. Aa lava forms when viscous lava cools and fractures into broken fragments, creating sharp-angled boulder fields (Fig. 2a). Pahoehoe lava, produced by fluid lava flowing under a cooling surface, is generally smooth and ropy, with large coherent features such as deep fissures and collapsed pits (Fig. 2b). Sparsely vegetated flows support a variety of shrubs, including rock spirea (Holodiscus dumosus), mock orange (Philadelphus lewisii), and fernbush (Chamaebatiaria millefolium), sparse forbs, and small clumps of perennial grasses. Aa and pahoehoe lava cover approximately 406 and 619 km2 of Craters of the Moon, respectively (Bell et al. 2009). Fifty-three percent of mapped aa lava was classified as barren (as opposed to sparsely vegetated aa), whereas only 5% of mapped pahoehoe was classified as barren, as determined from a recent land cover map of the study area (Bell et al. 2009).
Figure 2
Photographs of A) aa lava and B) pahoehoe lava in Craters of the Moon National Monument and Preserve, Idaho.
Craters of the Moon appears to be isolated from other suitable pika habitat. Most of the surrounding area is heavily grazed rangeland and irrigated agricultural lands that are devoid of talus. Natural vegetation surrounding and between lava flows is sagebrush steppe. The northern portion of Craters of the Moon abuts the foothills of the Pioneer Mountains, where some talus habitat occurs. Typical dispersal distances recorded for pikas are <300 m (Franken and Hik 2004; Moilanen et al. 1998; Peacock and Smith 1997; Smith 1974), but longer migrations of >1 km also occur under suitable conditions (Hafner and Sullivan 1995; Peacock and Smith 1997). The nearest known recorded extant pika locality is in the Pioneer Mountains, approximately 35 km north of Craters of the Moon. Few patches of talus habitat occur in the intervening lands, and the amount of interchange among these 2 areas is unknown but suspected to be low (D. J. Hafner, pers. comm.).
Sampling design
Pika detectability generally is considered high (>90%), and studies often attempt to census populations within talus patches (Beever et al. 2003, 2008; Moilanen et al. 1998; Morrison and Hik 2008; Peacock and Smith 1997; Smith 1974). Millar and Westfall (2010) recently have proposed a nonrandom judgment or convenience sampling approach. Because of the extent and inaccessibility of Craters of the Moon, a census was infeasible, and we used a probabilistic sampling design involving 4 discrete sampling frames that were accessible from roads (Fig. 1) and captured a range of environmental conditions. Each frame was constructed as a 1-km buffer around a section of selected road, and overlain on a pre-existing parkland cover map (created circa 2001). Road segments were selected that either penetrated the lava flows (i.e., the Monument sampling frame; Fig. 1) or that provided the best access to a range of lava types and elevations. The interior portion of the study area has no roads or trails beyond those in the Monument sampling frame (Fig. 1), and safety and travel time constraints precluded more extensive sampling. The access roads we used included an infrequently traveled secondary road into the Monument area and rugged backcountry jeep trails. In a post hoc exploratory analysis we tested for an association of distance to road with pika detection patterns (reported in “Results”) but did not include this factor in our set of a priori candidate models.
We identified target lava areas and excluded nonlava land cover within the 2001 map and clipped those target areas to sampling frame boundaries. The 4 sampling frames encompassed 2,282, 2,123, 605, and 579 ha of mapped lava suitable for sampling. We used the generalized random tessellation stratified (GRTS) spatially balanced sampling design (Stevens and Olsen 2004) to draw random survey locations in each frame. The GRTS algorithm is an improvement over simple random or systematic sampling designs because it maintains spatial balance even in the event of survey site deletions and additions, such as when map errors or other logistical issues arise. It also is a more efficient design in the presence of spatial pattern or correlation because it enforces regularity in the spacing of survey sites, within the constraints of sampling frame boundaries, somewhat minimizing the effects of spatial autocorrelation (Stevens and Olsen 2004). The average nearest-neighbor distance between sampling sites within frames was 484 m, beyond typical dispersal distances recorded for the species.
We surveyed 72 sites in the Monument sampling frame (Fig. 1) during September 2007. Each site was surveyed twice with independent observers and on different days to permit estimation of detection probabilities. Presence–absence surveys involving species that are imperfectly detected can be negatively biased (MacKenzie et al. 2002; Royle and Dorazio 2008; Tyre et al. 2003), although pikas are considered readily detectable when present, particularly when their characteristic feces and food caches (“haypiles”) are used as indirect signs of site occupancy, in addition to the direct detection possible by sighting individuals or hearing their distinctive calls (Beever et al. 2003, 2008, 2010; Moilanen et al. 1998; Smith 1974). We revisited 68 sites in the Monument sampling frame in June 2008 to estimate annual turnover in site occupancy. We surveyed 23 and 30 sites, respectively, in the Laidlaw Park and Huddle's Hole sampling frames in September 2008 (Fig. 1). We surveyed 19 sites in the Eastside frame during July 2009 (Fig. 1). Ancillary searches were conducted periodically during the survey period in several additional portions of the study area, including the Wapi lava field (Fig. 1) in 2007.
Survey methods and site characterization
We used both direct sign (aural or visual detection of individual pikas) and indirect sign (detection of pika feces and haypiles) for determining site occupancy in our surveys (Beever et al. 2003; Morrison and Hik 2008; Smith 1974). We made an effort to distinguish between old and fresh feces and hay and, for modeling, restricted detections to include only fresh sign. Fresh feces were dark brown in color and contained green plant material when opened and were often stuck together in small mounds by dried urine and mucus, whereas old feces were gray in color, dry, and scattered. Likewise, fresh hay also was green, and recent additions to old haypiles were conspicuous. Each survey lasted 20 min and included both active searching and quiet listening within that period. Survey site dimensions were defined to be circles with 12-m radii (452 m2), and signs outside that distance were not considered a detection for modeling purposes, although incidental detections were noted. Plot dimensions were established temporarily in the field with metric tapes and flagging. The 12-m plot dimension was selected through consideration of typical territory sizes of pikas and of a logistically tractable area to search within the 20-min survey period. In other locations pikas defend territories of approximately 400–700 m2 (Smith and Weston 1990), and the 452-m2 plot was small enough to fall within an average territory and large enough to be used regularly during daily activities, helping ensure that occupied sites had nonzero detection probabilities (MacKenzie 2006).
At each site we measured elevation and categorically described vegetation and geologic attributes considered to have predictive utility of pika occurrence. Elevations for the center of each site were obtained from a 10-m digital elevation model of the entire study area produced by the United States Geological Survey ( http://www.usgs.gov). Elevations of survey sites ranged from 1,434 to 1,846 m. Site elevations were mean-centered and standardized (i.e., Z-scores, subtracted from mean elevation and divided by SD) for use as a predictor in models, requiring model parameter estimates to be interpreted in terms of SD. Use of standardized regression covariates improves numerical performance of model-fitting algorithms when raw values of covariates are on different scales, and improves interpretability and comparability of parameter estimates among studies from different study areas (Gelman and Hill 2007; Schielzeth 2010).
Vegetation attributes included visual estimates of percentage cover within the 12-m-radius plot for graminoids (grasses and sedges), forbs, shrubs, and trees. Cover estimates were categorized using Daubenmire's (1959) scale modified with the addition of 0%, trace (<1%), and 100% categories to discriminate better among sparsely vegetated and barren lava, resulting in 9 cover classes. We included only graminoid and forb cover as predictors in models. Forbs differ from graminoids in physical and chemical features and nutritive properties, and these differences have been used to explain foraging behaviors of pikas (Dearing 1997; Huntly et al. 1986).
Geologic attributes included lava type and structural complexity. Lava types at survey sites were classified in the field as either pahoehoe or aa. We also used a new land cover map produced by the National Park Service (NPS—Bell et al. 2009) to assign pahoehoe or aa to each survey site to support predictions across the entire study area. This map was produced following the National Vegetation Classification System (Comer et al. 2003); classification accuracy of relevant land cover classes was ≥80% (Bell et al. 2009) and closely matched our field-based designations (88%). Lava complexity at survey sites was ranked into 1 of 3 classes. Class 3 sites were the most complex encountered during the study and included pahoehoe with an abundance of large broken surfaces containing fissures, pits, overhangs, and aa boulder fields with large fragments (>0.5 m2) distributed across the landscape in ridges and pits with relief >2.5 m (Fig. 2). Class 2 sites were intermediate in structural complexity, with few broken surface features, and aa fragments <0.5 m2. Class 1 sites were the least complex, with no broken surface features, and with only narrow crevices, smooth pavement-like pahoehoe features, and aa boulder fields with small fragments and no surface relief. For models, we collapsed complexity classes 1 and 2 together to create a single indicator variable identifying the most complex sites. Observers trained together at the beginning of each survey session and calibrated decisions concerning old and new pika signs, vegetation cover estimates, and lava complexity classifications. Representative photographs and cover estimation guides were carried during surveys for reference and calibration.
Data analyses
We estimated detection probabilities from the 2007 survey data as an initial step in model fitting. Detection probabilities were reported as >90% in Lava Beds National Monument (C. Ray, pers. obs.), in a low-elevation study area in Nevada (Beever et al. 2008), and across the hydrographic Great Basin (Beever et al. 2010), and we had expected equally high detection rates in our study. We used the R programming environment (R version 2.8.1, 2008; R Development Core Team, Vienna, Austria) and modified code from Royle and Dorazio (2008) to estimate detection probabilities from the site × survey detection history matrix. Site turnover between 2007 and 2008 in the Monument sampling frame was estimated by computing the proportion of pika detection sites gained and lost in 2008.
After considering detectability and site turnover in the Monument sampling frame (reported in “Results”), we pooled all data from all years (n = 144 independent sites) and randomly selected 94 sites (65%) for model development (training data set), and reserved 50 (35%) for model validation (Guisan and Zimmermann 2000; Van Houwelingen and Le Cessie 1990). Twelve sites (8%) were classified incorrectly by the NPS map (Bell et al. 2009) as occurring on steppe vegetation and were removed from analyses of geographic information system (GIS) models but were retained for use in ground-based models because these sites were classified successfully to lava type in the field. We used logistic regression models, which required assumptions that detection probabilities were ≈1 and that sites were independent and of constant occupancy status throughout the study period.
We constructed 2 sets of models for evaluation. The first set included 17 models representing a series of hypotheses about factors associated with pika distribution, drawn from published literature and from a prior pika survey conducted in Lava Beds National Monument in 2005–2006 (C. Ray, pers. obs.). These models were constructed with site attribute data collected in the field and included combinations of elevation, graminoid and forb cover, and lava type and structural complexity (substrate) predictors thought to be ecologically relevant. The following 4 hypotheses were addressed. First, we used a null model to represent the possibility that pikas in Craters of the Moon are unconstrained by elevation because of their ability to use extensive cool microhabitats within the lava, and that they do not discriminate among different types of lava and vegetation conditions on lava flows. Second, we constructed models with lava type and lava structural complexity parameters to represent the hypothesis that pikas can persist in the harsh climate of Craters of the Moon because of cool microhabitats deep within the lava. Structurally complex lava also might reduce weasel predation, which has been suggested as a potentially important factor in pika population dynamics elsewhere (Clinchy et al. 2002; Smith and Gilpin 1997). Third, we included models with forb and graminoid cover parameters because pikas are selective in the kinds of plants consumed during grazing and haying activities (Gliwicz et al. 2006; Huntly et al. 1986; Johnson 1967), and because forb and graminoid cover were shown to be predictive of pika site occupancy in another lava-dominated habitat (C. Ray, pers. obs.). Finally, given the apparent importance of elevation in current patterns of persistence and extirpation elsewhere in the species range (Beever et al. 2003; Grayson 2005; Moritz et al. 2008), we developed models with an elevation parameter. Elevation itself is an indirect predictive gradient that represents the compound influences of temperature and precipitation on species occurrence probabilities (Guisan and Zimmermann 2000; Körner 2007).
The second (GIS) model set included only a null model and 3 models with covariate combinations of elevation and lava type. Sample size and the post hoc classification of our survey sites using the NPS land cover map required us to aggregate barren and sparsely vegetated lava map classes, somewhat confounding the effects of lava type and vegetation in our GIS models. The GIS models were less precise but allowed us to estimate the probability of pika occurrence across the entire lava flow complex in Craters of the Moon and generate a predictive map of the species' distribution. Because our hypotheses were not mutually exclusive, we considered a global model and reduced models for each set that contained particular combinations of covariates simultaneously evaluating the relative importance of elevation, substrate, and vegetation using information–theoretic model selection (Burnham and Anderson 2002).
Candidate models were fit to the training data set and evaluated using Akaike's information criterion adjusted for small sample size (AICc—Burnham and Anderson 2002). Model parameters were estimated by maximum-likelihood procedures, and covariates were modeled with the logit link (i.e., logistic regression) using R. We used both Kendall's τ and Spearman's ρs rank correlation coefficients to test for correlation between lava type and vegetation cover and to guard against multicollinearity in models. Estimates of correlation were <0.5 for all site covariate pairs. We evaluated model goodness-of-fit by estimating ĉ (residual deviance/degrees of freedom—Burnham and Anderson 2002) for the global model in each set. A value substantially >1 (e.g., 2) would have triggered the use of a quasi-likelihood adjustment (QAICc—Burnham and Anderson 2002). We summed model weights to obtain the 95% confidence set of models and used model averaging (Burnham and Anderson 2002) to account for model selection uncertainty. Model weights were summed for each predictor variable to obtain relative importance values, which can range in value from 0 to 1 (Burnham and Anderson 2002). Models within 2 AICc units of the top model (with smallest AICc value) were considered competitive and averaged (Burnham and Anderson 2002). We computed Nagelkerke's generalized R2 to provide a measure of the proportion of explained variance for each model (Nagelkerke 1991). Final models were validated with the reserved data and performance summarized by computing the area under the curve (AUC) of the receiver operating characteristic (ROC—Fielding and Bell 1997; Guisan and Zimmerman 2000). We tested for residual spatial autocorrelation after accounting for the effect of habitat by computing Moran's I spatial correlation index (Cliff and Ord 1981), although this was expected to be minimal given the low rate and short distances of typical pika dispersal and the absence of other obviously spatially structured habitat variables not included in models. We calculated Moran's I with inverse distance weights within each sampling frame where pikas were detected. We used a 1,000-m threshold above which sites were not considered within the spatial “neighborhood” and given a null weight. Finally, we summarized the amount and configuration of predicted habitat by summing the area of lava within classes of predicted occurrence probabilities and examining a map of the probability surface across Craters of the Moon.
Results
Surveys and ancillary searches
We detected pikas in 45 of 144 sites (31%) surveyed during 2007–2009. Signs of pika occurrence determined to be fresh were found only in the Monument and Huddle's Hole sampling frames (Fig. 1). Pika presence was determined primarily through the discovery of indirect sign during searches. Pikas were seen or heard (within 100 m of a survey site center) on only 6 occasions during formal surveys. At 6 sites fresh sign was encountered nearby (>12 m and <100 m from the site center), but these sites were not treated as occupied for modeling purposes. At one of these sites a pika call also was heard nearby. All of these 6 sites were within the Monument sampling frame and located near other occupied sites, and their classification as nondetection sites (in accordance with the 12-m-radius site boundary) did not strongly influence results (i.e., negligible changes in parameter coefficient estimates, R2N, and AUC values). Neither fresh nor old sign was encountered anywhere in the Eastside sampling frame. Fresh feces were not found in the Laidlaw sampling frame, and old feces were found at only 2 of the northernmost sites in that frame. We found no other evidence of pika occurrence in these 2 frames. We also found no evidence that detections were related to the distance between survey sites and roads (τ = 0.13, n = 94, P = 0.11).
Patterns of detection varied spatially across the study areas along elevation and vegetation gradients and among substrate types. The lowest elevation survey site where pikas were detected was 1,631 m in the Huddle's Hole sampling frame. Old feces were found at 1,605 m in 2 adjacent sites in the northernmost (and highest) portion of the Laidlaw sampling frame. Auxiliary searches outside formal sampling frames failed to detect pika sign below 1,600 m. An incidental aural pika detection was recorded immediately north of the Laidlaw sampling frame at an elevation of approximately 1,600 m. Forty-one of 45 detection sites (91%) were classified on the ground as occurring on pahoehoe lava. Of those, 23 (56%) also were classified as having a complexity level of 3 (highly complex). Most detections occurred on sites with relatively high forb cover. Thirty-eight (82%) detection sites were classified with forb cover ≥1%, and of those, 18 (40%) were classified with >5% forb cover. Conversely, only 12 detection sites were classified with graminoid cover ≥1%. Haypiles we observed contained mostly forb material and some shrub material, but rarely graminoid material.
Occurrence models
Detection probability of pikas in the Monument sampling frame in 2007 was estimated to be 0.92 (95% CI = 0.74–0.98). The increase in estimated site occupancy from the naive estimate after accounting for this small amount of detection bias was <1%. Parameter estimates and SE of a naive logistic regression model β0 + β1*(lava type) were identical, within 2 significant digits, to the same model adjusted to account for detectability.
Site turnover in the Monument sampling frame between September 2007 and June 2008 was 16%. Ten sites determined to be unoccupied in 2007 were occupied in 2008, and 1 previously occupied site in 2007 was found unoccupied in 2008, although an aural detection was made within 100 m of this site during the 2008 survey, suggesting that the turnover there was not actual.
Modeling results showed that pika occurrence probabilities varied spatially across measured environmental gradients. The 95% confidence set of ground-based models included the 5 most complex of 17 candidate models (Table 1A). Elevation and lava type were included in all 5 models, and both variables had relative importance values of 1. Lava complexity and forb cover were included in 3 models, with relative importance values of 0.93 and 0.69. Graminoid cover was included in only 2 of the top models and had a low relative importance value of 0.31. Three ground-based models were well supported using the AICc criterion (ΔAICc < 2) and had similar model weights (wi > 0.25) and R2N values (0.39–0.40; Table 1A). These models were averaged to generate parameter estimates (Table 2). Selection among GIS models provided similar evidence for the importance of elevation and lava type in predicting pika occurrence. The full GIS model containing both elevation and lava type parameters had 79% of model weights (Table 1B).
Table 1
Model selection results for occurrence of pikas in Craters of the Moon, Idaho. Models reported are the 95% confidence set of A) ground-based models and B) geographic information system (GIS) models, ranked using Akaike information criterion (AICc). wi indicates the model weight, k is the number of model parameters, and R2N is Nagelkerke's generalized R2. Ground-based models included substrate and vegetation predictor variables measured at each site, whereas GIS models included only predictors obtained from a digital elevation model (DEM) and land cover map that enabled prediction of pika occurrence probabilities across the study area. Model names reflect the predictor variables for elevation (standardized, interpreted in terms of SD, in m), substrate (lava type and complexity indicators), and vegetation (forb and graminoid cover).
Table 2
Parameter estimates and margins of error (±2 SE) for averaged ground-based models and the highest-ranked geographic information system (GIS) model describing occurrence of pikas in Craters of the Moon, Idaho. Variable names are defined as in Table 1.
Parameter estimates from averaged ground-based models and from the top GIS model indicated that all variables were associated with pika site occurrence, although considerable uncertainty surrounded the model-averaged estimates for elevation, lava complexity, and vegetation cover (Table 2). For each 1-unit increase in SD above the mean elevation it was estimated that the odds of pika occurrence increased at least 2-fold (e.g., e0.69 = 1.99; Table 2). Mean elevation of survey sites was 1,679 m, and SD was 109 m. The estimated coefficient of elevation from the ground-based model was 1.5 times that estimated from the GIS model after accounting for lava complexity and vegetation cover, although the uncertainty in this estimate also increased considerably and the 95% confidence interval for this estimate included zero (Table 2). Lava type was associated most strongly with pika distribution, and coefficient estimates were similar for both the averaged ground-based model and the top-ranked GIS model (Table 2). The odds of pika occurrence at pahoehoe sites were estimated to be more than 10 times those of aa sites. Sites ranked as most complex (i.e., class 3) were estimated to be >4 times more likely to be occupied by pikas than less complex sites. The odds of pika occurrence were estimated to increase by 1.6 times for each 1-unit increase in forb cover class. Conversely, graminoid cover was negatively associated with pika occurrence, but the importance and magnitude of the association was uncertain (Table 2). Occurrence probabilities increased along the elevation gradient at a faster rate for sites containing pahoehoe lava, and detection sites were clustered on the upper end of the elevation gradient above 1,600 m (Fig. 3).
Figure 3
Occurrence probabilities, with 95% confidence interval bands (dashed lines), for pikas by elevation for Craters of the Moon, Idaho. Probabilities and intervals were computed from parameter estimates from the highest ranked (lowest Akaike information criterion [AICc]) GIS-based logistic regression model. Sites (points) are coded according to lava type and observed detection status. Data shown include both training and validation data (n = 132).
Model fit and predictive utility were adequate for both ground-based and GIS models, although ground-based models outperformed GIS models as expected because of the inclusion of additional predictors. Global goodness-of-fit, as measured by ĉ, was ≈1 for both model sets. AUC for the averaged ground-based model from training data was 0.84 and was 0.81 for the highest-ranked GIS model. AUC for ground-based and GIS models using validation data was 0.76. The true positive prediction rate (sensitivity) of the GIS model used for predicting pika occurrence throughout Craters of the Moon was 100% for training data and 94% for validation data. True negative rates of prediction (specificity) were 64% for training data and 57% for validation data, suggesting that our model overestimated pika occurrence somewhat. The overall accuracy for the GIS model using training data was 75%, and error was 24%. Accuracy of this model using validation data declined to 72%, and error increased to 27%. Spatial correlation among GIS model residuals was low, as expected. Moran's I was estimated to be 0.06 and 0.07 for the Monument and Huddle's Hole sampling frames, respectively, although these estimates did not differ significantly from zero (P ≥ 0.50).
Amount and spatial configuration of predicted habitat
We used the highest-ranked GIS model to predict pika occurrence probabilities throughout Craters of the Moon (Fig. 4). Mean predicted probability across the entire study area was 0.2. However, as previously described, considerable variation existed in the spatial distribution of predicted occurrence (Fig. 4). Probabilities declined along the elevation gradient from north to south and differed for pahoehoe and aa lava types (Fig. 3). Most of the study area was predicted to have low occurrence probabilities (Figs. 4 and 5). Sixty-two percent (638 km2) of the study area had below-average (<0.2) probability of occurrence and 40% (418 km2) had <0.1 probability of occurrence (Fig. 5). Only 24% (250 km2) of the study area had occurrence probabilities ≥0.38, the optimal cutoff value identified from the ROC curve “shoulder,” the point along the convex hull of the ROC curve where the correlation (φ) between true positive and true negative rates is largest. All 41 pahoehoe detection sites were included among sites of high occurrence probability when using this cutoff value. Areas containing occurrence probabilities ≥0.38 were restricted to pahoehoe flows on the northern portions of the study area above 1,600 m and included 1 very large area representing 195 km2, 2 areas >5 km2, 4 others >1 km2, and numerous (>250) small areas <500 m2 (Fig. 4).
Figure 4
Distribution map for pikas in Craters of the Moon, Idaho, on the basis of occurrence probabilities estimated from the highest-ranked geographic information system (GIS)-based model. Pikas were most likely to occur on pahoehoe lava flows in the northernmost, highest-elevation portions of the study area.
Discussion
The pika appears to be largely restricted in its distribution to an approximately 250-km2 area of structurally complex, sparsely vegetated pahoehoe lava flows in the northern portion of Craters of the Moon National Monument and Preserve above 1,600 m. This is considerably less area than appears available to the species there superficially but remains a substantial area of contiguous habitat. Assuming that pikas maintain 700 m2 territories on average (the upper end of estimates reported by Smith and Weston 1990), and that at least 38% of territories are occupied within the 250 km2 area of northern pahoehoe lava (Fig. 4), >135,000 occupied pika territories may exist in Craters of the Moon. The next-largest area of lava habitat occupied by the American pika of which we are aware is at Lava Beds National Monument (C. Ray, pers. obs.) in northern California, containing approximately 22 km2 of surface lava. Although Craters of the Moon is a particularly spectacular example, given its size and isolation from surrounding montane habitat, pikas persist in many other locations throughout the intermountain west where extensive lava flows (e.g., Lava Beds National Monument, Newberry National Volcanic Monument) and other low-elevation habitats (Beever et al. 2008; Millar and Westfall 2010; Simpson 2009) are found. Our results provide insight into some of the factors that might influence pika persistence, especially in peripheral habitats.
Lava type had the strongest influence on occurrence probabilities from both our ground-based and GIS-based sets of candidate models. The complexity ranking of lava was influential also in ground-based models. Sites comprised of highly complex pahoehoe lava had higher occurrence probabilities. Aa lava sites, regardless of complexity, were used rarely by pikas. Most previous studies of pika habitat in montane talus environments have not differentiated among substrate attributes and have assumed implicitly that all habitat patches were of equal quality (Moilanen et al. 1998; Smith and Gilpin 1997). However, Franken and Hik (2004) considered an index of habitat quality on population dynamics of collared pikas and found it important to their incidence function results. Variation in the extent and timing of snow cover has been related to pika recruitment and survival (Kreuzer and Huntly 2003; Morrison and Hik 2007). C. Ray (pers. obs.) found crevice depth positively associated with pika occurrence in lava flow sites in Lava Beds National Monument. Most recently, Millar and Westfall (2010) found that rock-ice features in montane environments were associated with patterns of pika occurrence. Without inferring habitat quality directly, our results suggest that substrate does matter, presumably by ameliorating predation rates, the impacts of harsh surface weather conditions, or both. Some of the pahoehoe lava flows in Craters of the Moon are >10 m thick (Kuntz et al. 1988), and we suspect that complex pahoehoe flows offer pikas better protection from predators (weasels) and weather conditions. Additionally, pahoehoe lava flows appear to support more (albeit sparse) vegetation than aa flows and therefore might offer pikas a better foraging environment. It is also possible that pikas inhabiting pahoehoe flows benefit from thermal anomalies in which air flowing through rock can be cooler than expected in the summer and warmer than expected in the winter. This phenomenon has been suggested as a mechanism to explain the close association between pikas and rock-ice features (Millar and Westfall 2010) in montane environments and also might influence differential pika selection among lava types and structural characteristics in Craters of the Moon. Anecdotally, air temperatures in pahoehoe crevices felt noticeably cooler than in the interstices of aa boulders during our summer surveys.
In some situations, such as occurs at Craters of the Moon, substrate features could enable pikas to persist outside their typical bioclimatic envelope and might reduce the effects of accelerated climate change, which in our region is predicted to include a reduced snowpack and elevated summer temperatures (Karl et al. 2009). However, our results showed that elevation was associated with patterns of pika distribution, even after accounting for substrate. Elevation reflects the temperature and precipitation gradients across Craters of the Moon (Davey et al. 2006), and the physiological constraints of pika metabolism and thermoregulation might present the ultimate controls on pika distribution. If so, accelerated climate change could erode the area of lava habitable to pikas in Craters of the Moon and reduce their long-term probability of persistence, despite the mitigating effects of substrate. This is consistent with the findings of Beever et al. (2003, 2010), Galbreath et al. (2009), Grayson (2005), and Hafner (1994) that pikas inhabiting interior regions of the western United States are experiencing range contractions along an elevational gradient.
Elevation and both lava type and complexity presumably influence vegetation cover in Craters of the Moon, despite the low to moderate correlations (<0.50) observed among our data. However, even after accounting for these other factors, our ground-based models suggested a positive, albeit weak, influence of forb cover and a negative influence of graminoid cover on pika site occurrence probabilities. Johnson (1967) analyzed pika diets in Colorado and found that forbs and shrubs, rather than graminoids, comprised most of the food ingested throughout the year. Other studies have reported that pikas tend to graze or eat grasses directly while selecting forbs, possibly on the basis of their nutritional value, for haying (Dearing 1997; Gliwicz et al. 2006; Huntly et al. 1986). Our results were similar to those observed in Lava Beds National Monument (C. Ray, pers. obs.), where a positive relationship between pika occurrence and forb cover and a negative relationship with graminoid cover were documented. Pikas are central-place foragers and can modify local patterns of vegetation (Huntly et al. 1986; Roach et al. 2001). However, given that pikas preferentially select forbs, particularly when haying, this impact should have resulted in an observed negative relationship between site occupancy and forb cover in our study area. Our results therefore suggest that pikas might be selecting for sites with relatively high forb cover at Craters of the Moon. The evidence for a negative relationship between pika occurrence and graminoid cover was less clear, however, and could be a spurious effect of localized pika grazing (Huntly et al. 1986; Roach et al. 2001).
The diagnostics of accuracy and predictive performance for our models showed adequate results. Fielding and Bell (1997) suggested that an AUC score >0.7 indicates acceptable predictive power for presence–absence models, and our models exceeded this level, even with a validation data set not used in model parameter estimation. Also, the small amount of residual spatial autocorrelation suggests that our models included the primary sources of spatial variation in patterns of pika occurrence. However, we note that both training and validation data sets were obtained from only a small proportion of Craters of the Moon. Also, both our ground-based and GIS models had relatively low model specificity and could have overestimated pika occurrence probabilities in some of the study area, particularly in the southern portion where we failed to document pikas. Accordingly, our results might be overly optimistic and should be viewed with an appropriate degree of caution in terms of the amount of area containing high occurrence probabilities. Another concern is the amount of annual site turnover (16%) observed between 2007 and 2008. Site turnover might have introduced additional error into our models by altering occupancy patterns and influencing parameter estimates. Inferring process from patterns of occupancy is problematic (Clinchy et al. 2002; Tyre et al. 2001), and future studies of this kind would be improved by explicitly accounting for local rates of turnover as described by MacKenzie et al. (2003) and Royle and Dorazio (2008) with longitudinal, repeat presence–absence surveys. Repeat surveys from a spatially extensive probabilistic sample such as ours also could provide a suitable design for long-term monitoring of changes in pika distribution and enable the evaluation of hypotheses related to the effects of climate change on pika distribution.
Craters of the Moon might serve as an important long-term refugium for pikas due to the extent of the pahoehoe lava habitat there. The spatial configuration of habitat patches with relatively high predicted occurrence probabilities (e.g., >38%; Fig. 4) resembles the mainland–island configuration described by Harrison (1991) and Krohne (1997), and that paradigm has been used to describe at least 1 talus-dwelling pika population (Peacock and Smith 1997). Populations in mainland–island habitat complexes are generally considered likely to persist because of the high probability of recolonization from mainland patches (Krohne 1997; Peacock and Smith 1997). The gradual rather than discrete nature of patch boundaries in Craters of the Moon, and the potentially well-connected nature of these patches via rocky lava corridors, reinforce this possibility for Craters of the Moon.
Nevertheless, as noted above, accelerated climate change could erode the size and proximity of suitable habitat patches and dramatically change the future prospects for the species at Craters of the Moon. The species might be faced with an unpaid “extinction debt” (Tilman et al. 1994) in which observed declines in distribution and site occupancy lag behind declining habitat quality, as has been described for populations across the hydrographic Great Basin (Beever et al. 2003). This is particularly plausible given the extent of habitat at Craters of the Moon and its potentially large pika population and the relative longevity and sedentary nature of pikas (Smith and Weston 1990), all of which can increase the lag effect (Tilman et al. 1994). Historical factors generally are an important source of commission error (overprediction) in distribution models (Guisan and Thuiller 2005). If pika populations in Craters of the Moon are not in equilibrium, as was assumed here, but are already contracting in response to environmental change, this might explain some of the overprediction by our model. These questions underscore the importance of initiating long-term monitoring of pika populations in Craters of the Moon and elsewhere to identify population trends. Additionally, some low-elevation peripheral habitats with persistent pika populations could become particularly important for their potential to serve as refugia and as sources of genetic and phenotypic diversity (Channell and Lomolino 2000; Kark et al. 2005; Sexton et al. 2009), especially if climate change predictions prove accurate. Our results suggest that substrate, vegetation, and other habitat attributes, in addition to area and isolation metrics, should be considered when evaluating the conservation potential of these areas. Our results also suggest more refined hypotheses that should be evaluated through additional data collection and monitoring, such as the potential for certain lava types to support microclimates favorable for pikas, and that the distribution is contracting northward into the upper elevations of the Monument over time.
Acknowledgments
The order of authors, beyond the first author, is alphabetical and does not imply level of contribution. The project was funded by the National Park Service Upper Columbia Basin Network Inventory and Monitoring Program and Craters of the Moon National Monument and Preserve. J. Apel assisted in all aspects of the project and was instrumental to its success. D. Stucki, J. Bylund, students from the Oregon Museum of Science and Industry, and Student Conservation Association interns all assisted with field work. M. Lonnecker provided GIS support. C. Epps provided helpful comments during manuscript preparation. We thank associate editor J. Rachlow, A. Smith, and one anonymous reviewer for providing constructive comments that greatly improved the manuscript.
