Summer range occupancy modeling of non-native mountain goats in the greater Yellowstone area

Non-native species can have adverse impacts on native species. Predicting the potential extent of distributional expansion and abundance of an invading non-native species can inform appropriate conservation and management actions. Non-native mountain goats (Oreamnos americanus) in the greater Yellowstone area (GYA) have substantial potential to occupy similar habitats to native Rocky Mountain bighorn sheep (Ovis canadensis canadensis). To understand the potential for expansion of mountain goats in the GYA, this study evaluated detection-nondetection data derived from groundbased occupancy surveys of viewsheds partitioned into a grid of 100 3 100 m sampling units. Surveys were conducted over three summer seasons (2011–2013) in two study areas with well-established mountain goat populations. Relationships between scale-specific habitat covariates and mountain goat selection were evaluated to model occupancy and detection probabilities based on mountain goat detections in 505 of the 53,098 sampling units surveyed. Habitat selection was most strongly associated with terrain covariates, including mean slope and slope variance, at a spatial scale of 500 3 500 m, and canopy cover, heat load, and normalized difference vegetation index at a spatial scale of 100 3 100 m. These results provide new insight into multi-scale patterns of mountain goat habitat selection, as well as evidence that mean slope and slope variance are more informative terrain covariates than distance to escape terrain, which has been commonly used in published mountain goat habitat models. The model predicted 9,035 km of suitable habitat within the GYA, of which 57% is currently un-colonized. Seventyfive percent of all bighorn observations recorded in the GYA fall within predicted suitable mountain goat habitat. We also estimated that the GYA might have the potential to support 5,331–8,854 mountain goats when all predicted habitat is occupied, or approximately 2.5–4.2 times the most recent abundance


INTRODUCTION
Introductions and expanding distributions of non-native species are globally common and considered to be a significant component of human-caused global environmental change (Vitousek et al. 1997).The negative impacts of nonnative species to native biological diversity and ecosystem processes through competition for limited resources, predation, habitat alteration, hybridization, and disease transmission, have been documented across many taxa (Gordon 1998, Huxel 1999, Simberloff 2000, 2005, Pimentel et al. 2005, Peeler et al. 2011).However, effects of non-native species are not necessarily detrimental and may even be considered beneficial in certain conservation settings, for example, by providing habitat or food for native species, ecosystem services, or opportunities to preserve threatened species (Sagoff 2005, Schlaepfer et al. 2011).Understanding invasive species' patterns of habitat selection and predicting their range expansion is important for developing appropriate conservation and management actions (Lodge et al. 2006, Gormley et al. 2011).However, this can be challenging when species are rare on the landscape and inhabit heterogeneous, difficult-to-survey terrain.Traditional aerial survey methods are typically applicable only at a broad spatial scale and can lead to biased interpretations of habitat selection patterns due to imperfect detection (Gu andSwihart 2004, MacKenzie 2005).The more recent development of occupancy methods has provided the ability to incorporate estimates of detection into evaluations of habitat selection for rare species (MacKenzie et al. 2005(MacKenzie et al. , 2006)), and may also provide a means to understand selection at a fine spatial scale while effectively surveying large, heterogeneous landscapes.
The introduction and ongoing range expansion of non-native mountain goats (Oreamnos americanus) in the greater Yellowstone area (GYA) , a 56,000 km 2 mountainous region representing the largest intact temperate ecosystem in the world (Keiter and Boyce 1991), is among the leading conservation and management concerns of state and federal agencies (Schullery andWhittlesey 2001, Lemke 2004).Although agencies' specific management goals vary, the primary concern of regional biologists is the potential for non-native mountain goats to negatively impact native Rocky Mountain bighorn sheep (Ovis canadensis canadensis) through resource competition, behavioral displacement, and/or disease transmission (Pallister 1974, Stewart 1975, Dailey et al. 1984, Varley 1994).In western North America, bighorn sheep populations suffered catastrophic declines due to market hunting, disease die-offs, and habitat loss and degradation during European colonization (Buechner 1960, Montana Fish, Wildlife andParks 2010).Restoration and conservation efforts over the past 80 years have been modestly successful with the GYA representing one of the largest bighorn ranges in North America, supporting 5,000 to 9,000 bighorn sheep (unpublished data).Current conservation concerns are focused on increasing or maintaining viable bighorn sheep populations in the face of global climate change, continued habitat loss, and disease.
Mountain goats are expected to continue expanding their distribution throughout the steep and heterogeneous mountain terrain of the GYA that are currently inhabited by populations of bighorn sheep.Understanding patterns of habitat selection and the extent of range expansion of mountain goats is needed to guide management and conservation decision-making for both species.However, this information is lacking due to the remote and rugged sites often occupied by these species, as well as their general rarity and imperfect detectability on the landscape.Therefore, we used a novel field-sampling technique based on occupancy methods to assess habitat selection of mountain goats during the summer season.We fit suites of competing models to uncover the relative strength of association between habitat covariates and mountain goat occupancy, and extrapolated the habitat selection results based on the most supported model at the landscape level of the GYA to predict suitable summer habitat for mountain goats.Based on this prediction and recent abundance estimates of mountain goats in the colonized areas of the GYA, we estimated the abundance of mountain goats that may eventually be supported if all predicted habitat in the GYA becomes occupied.
Our general hypothesis is that mountain goats have considerably more suitable habitat in the GYA than currently occupied and thus distribu-tion and abundance can be expected to continue increasing for some time into the future.Two specific and novel hypotheses regarding habitat covariates were also evaluated.First, we expected that mountain goats select different habitat components at different spatial scales, as has been recognized for many species in investigations and discussions of scale-sensitive habitat selection (Johnson 1980, Wiens 1989, Hobbs 2003, Mayor et al. 2009, Laforge et al. 2015).To address this, we used model selection techniques to evaluate a range of possible spatial scales for each habitat covariate (grain), within Johnson's (1980) third-order selection (extent), to obtain a scale-specific habitat selection model (Hobbs 2003).Secondly, we sought to compare the performance of a predominant and well-supported habitat covariate, distance to escape terrain (DET), with new covariates capturing mean slope (SLP) and variance in slope (SLPv).DET, measuring proximity to the nearest area of a specified slope threshold (e.g., slopes !408), has been used to capture the close association of mountain goats with steep and rough terrain, often termed ''escape terrain'' in reference to sites used for predator avoidance (e.g., Gross et al. 2002, Poole and Heard 2003, Shafer et al. 2012).However, the choice of a threshold value has not only been variable (ranging from 25 to 508) and generally arbitrary, but also assumes mountain goat selection for slopes greater than that value is constant (e.g., for !258, everything from 25 to 908 is of equal quality), even though studies have indicated greater use of intermediate slopes (Gross et al. 2002, Poole andHeard 2003).DET also assumes that selection for steep slopes is the same regardless of the variability in slope.Therefore, we hypothesized that SLP and SLPv would provide a more flexible and biologically meaningful alternative than DET.

Study areas
To make inferences regarding potential mountain goat habitat selection at a landscape level, two study areas (Fig. 1) were selected that captured a range of variation in mountainous regions of the GYA and had well-established mountain goat populations.The Gallatin Crest (GC) study area (425 km 2 ) was located in the southern Gallatin Range of Montana and Wyoming between the Yellowstone and Gallatin Rivers, lying partially within the Gallatin National Forest and extending south into Yellowstone National Park (YNP).The North Absaroka (NA) study area (925 km 2 ) was located in the Absaroka Range of Montana and Wyoming, lying in the northeast corner of YNP and partially in the Gallatin and Shoshone National Forests.Elevations ranged from 1,800 to 3,300 m in the GC and from 2,000 to 3,475 m in NA.Average annual precipitation was 119 cm in the GC (Shower Falls SNOTEL at 2,469 m, 1981-2010 data) and 131 cm in the NA (Fisher Creek SNOTEL at 2,774 m, 1981-2010 data;public communication, http://www.wcc.nrcs.usda.gov/ snow).Vegetation in both study areas was dominated by dry-mesic spruce-fir forests (including Engelmann spruce, Picea engelmannii, and subalpine fir, Abies lasiocarpa) with Douglas-fir (Pseudostuga menziesii ) and lodgepole pine (Pinus contorta) forests at lower elevations.Subalpine fir and whitebark pine (Picea albicaulis) occurred in the subalpine (;2,400-2,900 m) giving way to turf and meadow graminoid and forb species in the alpine zone (;2,800-3,400 m).A full complement of potential mountain goat predators were present in both study areas including wolves (Canis lupus), coyotes (Canis latrans), mountain lions (Puma concolor), bobcat (Lynx rufus), wolverines (Gulo gulo) and eagles (Aquila chrysaetos, Haliaeetus leucocephalus).
At the time of this study, mountain goats had been present in both the GC and NA for approximately 50 years, sharing sympatric range with native bighorn sheep (Varley 1996, Lemke 2004).No mountain goat transplant history exists for these areas, suggesting the areas were colonized by dispersing individuals from releases in the 1940-50s in the Madison and Absaroka Ranges adjacent to the GC and the Beartooth Range to the northeast of the NA (Laundre ´1990, McCarthy 1996, Lemke 2004).

Sampling design
General mountain goat distribution data collected by management agencies during several decades of population surveys were used to delineate approximate areas within each study area for conducting occupancy surveys.Valleys and low-to mid-elevation, densely-forested slopes were excluded while all alpine and subalpine habitats and portions of forested slopes adjacent to subalpine zones were considered for sampling.All regions of each study area identified for sampling were overlaid with a 100 3 100 m grid in a GIS to delineate contiguous occupancy sampling units.Logistics and safety precluded randomized selection of sampling units.Instead, we used a systematic approach based on travel routes and groups of sampling units (viewsheds) that were visible from fixed vantage points (viewshed survey point).The travel routes (primarily along trails, drainage bottoms, and ridgelines) providing access to viewsheds were designed to accommodate 3-5 day backpacking trips and cover all accessible regions of the study areas during the summer season (June to October).The first viewshed survey point along each travel route was randomly selected within 3 km from the start, with subsequent viewshed survey points placed systematically every 3 km along the remaining travel route.A 500-m buffer around each viewshed survey point allowed observers to select a location that maximized visibility of viewsheds that were obstructed by immediate structures such as trees or rock outcrops.
Because vast areas can be visible from viewshed survey points, observers constrained the viewshed to areas that could be surveyed in a reasonable amount of time (typically less than 2 hours) and generally within 3 km of the observers, where optical performance was reliable (e.g., to avoid distortion from heat shimmer) and detection appeared practical.Observers surveyed sampling units within the viewshed with the aid of 103 binoculars and 20-603 spotting scopes to determine if each sampling unit was occupied by a mountain goat group (at least one mountain goat).At most viewshed survey points (;86%), two observers surveyed the same viewshed independently and simultaneously.Observers situated themselves such that neither was influenced by the behavior and activities of the other, typically by placing a barrier (i.e., rock outcrop, ridge, tree, or backpack) between each other, ceasing communication, and maintaining a distance of at least 5 m to limit visual cues of detections.These concurrent surveys were used to construct the detection history (presence ¼ 1, absence ¼ 0) for each sampling unit to account for imperfect detection of mountain goat groups.
Field computer tablets displaying the grid of sampling units layered over an aerial photograph of the area to be sampled and a custom program written in ArcPad 10.2 (ESRI 2013a) allowed observers to record occupied and unoccupied sampling units using the tablet's touch screen.Any sampling unit surveyed where a group was observed within its boundaries during the period of observation was recorded as occupied and all other cells surveyed were recorded as unoccupied.Group size, initial behavior of the majority of individuals (e.g., feeding, resting, traveling, or other), and sex and age class of individuals were recorded (Chadwick 1983, Côte ´and Festa-Bianchet 2003).Observers additionally recorded date, survey point ID and location, and survey start and end times.Daytime surveys were conducted between the third week of June and the second week of October 2011-2013.Sampling units were often revisited during subsequent backpacking trips, however, revisits were treated as independent sampling units for analysis.For a detailed description of the methods for data management and processing after field work, see O'Reilly (2012).

Model development
We used single-species, single-season occupancy modeling (MacKenzie et al. 2002) in package ''unmarked'' of program R version 3.1.0(Fiske and Chandler 2011, R Core Team 2013) to estimate occupancy (w) and detection ( p) probabilities of sampling units by mountain goat groups.w was defined as the probability of a group being in a sampling unit at the time of the survey, and p was defined as the probability of detecting a group given that it was present at the time of the survey.Sixteen habitat attributes were selected based on review of published mountain goat studies and hypotheses of habitat relationships, and included in models as site covariates on w (see DeVoe 2015 for covariate development).To address specific hypotheses as well as to simplify otherwise large and complex model lists, we arranged covariates into four model selection suites: sloped terrain, ruggedness, forage, and heat load (Table 1; see respective descriptions in Model suites).A range of possible scale-specific habitat selection was captured by considering each covariate at five spatial scales: 100 3 100 m (that of the sampling units), 300 3 300 m, 500 3 500 m, 700 3 700 m, and 900 3 900 m.These spatial scales are referenced as 100-to 900-m (e.g., ''500-m'' for 500 3 500 m) and v www.esajournals.orgindicated by subscript (e.g., SLP 500 ).Additionally, linear, pseudothreshold (natural log transform; after Franklin et al. 2000), and quadratic forms were considered for each covariate to allow flexibility in the hypothesized relationships that may better describe the data.Each functional form was evaluated at each of the five spatial scales.No interactions were modeled due to the limited amount of detections in the dataset, and a model that held w constant was also included in each suite for model comparison purposes.Pairs of covariates with Pearson correlation coefficients of jrj .0.60 were not considered together in the same model to avoid collinearity problems and reduce the number of possible models.To aid in model convergence and covariate comparisons, all covariates considered for linear and quadratic functional forms were mean-centered and divided by their standard deviation.If necessary, covariates considered for the pseudothreshold form were adjusted slightly by addition of 0.005 to avoid values less than or equal to zero, allowing log transformation.Although primary interest was in covariate relationships with w, we also modelled variation in p by including a pseudothreshold form of canopy cover measured at the scale of the sampling unit (100-m) for all models in each suite (Poole 2007, Rice et al. 2009).This form hypothesizes that increasing canopy cover decreases detection probability but at a declining rate (near-asymptotic).Post hoc, after arriving at a final inferential model for w, we evaluated other covariates on p including linear forms of canopy cover, as well as slope variance, ruggedness index, and vector ruggedness measure (see Model suites for covariate descriptions) to capture complex landscapes that may make detecting mountain goats difficult (Poole 2007, Rice et al. 2009).These covariates were measured at the 100-m scale and considered singly and in additive combinations of canopy cover and terrain covariates in the models.A model that held p constant was also included in the evaluation.
To avoid model convergence issues caused by modest sample size, we retained detections that were made subsequent to the initial detection as groups traveled across sampling units.For models that did not initially converge or had indications of poor convergence (condition num-bers .10 5 ), several methods were attempted: (1) the number of iterations were increased (from 1,000 to 5,000 to 10,000), (2) different optimization methods (including Nelder Mead, conjugate gradient, BFGS, and L-BFGS-B) were attempted, and (3) starting values from the optimization that had the lowest negative log likelihood estimate and no apparent problems with coefficient and standard error estimates were fed into subsequent optimizations (Fiske and Chandler 2011).Convergence was indicated by no apparent problems with coefficient and standard error estimates, no improvement to the negative log likelihood with further modifications, and condition numbers ,10 5 .
Because we suspected lack of independence of observations due to multiple detections of the same group, MacKenzie and Bailey's ( 2004) parametric bootstrapping technique for occupancy data was used to estimate an overdispersion parameter (ĉ).Due to the difficulty of determining a complex model from the numerous varieties and forms of covariates for estimating ĉ, a reasonable model was obtained by: (1) ranking candidate models within each suite using Akaike's information criterion (AIC) and model weights; (2) including the resulting models ,5 AIC units of the most-supported models in a combined model suite in all possible additive combinations (see Combined model suite); and (3) selecting the resulting top-ranked model.Based on the ĉ from this model, we used quasilikelihood (QAIC) theory to allow overdispersion in the model selection procedure and inflate the estimates of sampling variance of the parameters (White et al. 2002, MacKenzie andBailey 2004).Within each suite and among suites, candidate models were again ranked, and models ,5 QAIC units of the top models were compared in a combined model suite to identify well-supported models.From this final suite, noncompetitive models that were similar to well-supported models but containing uninformative parameters were discarded (Burnham andAnderson 2002, Arnold 2010).

Model suites
Sloped terrain.-Thissuite was developed to directly test the prediction that SLP and SLPv would outperform DET covariates.We created five versions of DET based on the literature v www.esajournals.org(Table 1) and compared them against SLP and SLPv.We included models with SLP alone, but never with SLPv alone because it did not make biological sense unless in combination with SLP.Models with both SLP and SLPv were considered at the same spatial scale while functional forms were allowed to vary.This suite included 66 models.
Ruggedness.-This suite was developed to evaluate other terrain features based on previous mountain ungulate studies.We selected two ruggedness indices, including terrain ruggedness (RUG; Poole et al. 2009), measuring the variability in the rate of change in slope, and vector ruggedness measure (VRM; Taylor et al. 2004, Sappington et al. 2007, Beus 2010), measuring the variation in the three-dimensional orientation of a neighborhood of elevation raster pixels.Mountain goats have also been found to be in close proximity (within ;400 m) to rock outcrops and substrates typically associated with steeper, broken terrain (Foster 1982, Haynes 1992, Taylor and Brunt 2007).Therefore, a rock (RCK) covariate, measuring the proportion of land cover data reclassified to rock, was evaluated in addition to RUG and VRM with predictions that mountain goat occupancy would be positively associated with all three covariates.Because RUG and VRM were considered different constructs of similar landscape features, they were not evaluated additively in the same model.Additive combinations of RUG with RCK and VRM with RCK were constrained to the same functional form and spatial scale within each model to maintain a model suite of manageable size.This suite included 76 models.
Forage.-Mountain goats depend on summer forage to replenish fat reserves lost during the winter period, support milk production and gestation in females, accumulate muscle mass for the rut in males, and enhance survival of offspring (Pettorelli et al. 2006, 2007, Festa-Bianchet and Côte ´2008).However, the majority of published mountain goat habitat models lack evaluation of covariates describing forage.Therefore, this suite was developed to evaluate four potential forage covariates.Because meadows may provide forage opportunities, we created a proportion meadow (MDW) covariate from land cover reclassification of meadow, grassland, shrubland, and alpine turf.From Landsat The-matic Mapper satellite imagery taken from July 2000 scenes, we derived a normalized difference vegetation index (NDVI) covariate, a measurement that correlates strongly with aboveground net primary productivity (Kerr andOstrovsky 2003, Pettorelli et al. 2005) and has been used to assess impacts to juvenile survival and condition (Pettorelli et al. 2006(Pettorelli et al. , 2007)), reproductive success (Pettorelli et al. 2006), and seasonal migration patterns (Sawyer and Kauffman 2011, Bischof et al. 2012, Lendrum et al. 2014) in ungulates.From the Landsat imagery, we also derived the tasseled cap wetness (TCAP) covariate, which correlates with vegetation biomass and soil moisture (Todd et al. 1998).Tasseled cap wetness has been evaluated in a limited number of habitat studies of other species (Carroll et al. 2001, Alexander et al. 2005), but has not been applied to mountain goats.For both NDVI and TCAP covariates, the effect of forested areas with greater than 15% canopy cover was masked out to capture potential foraging areas.We predict mountain goat occupancy will be positively associated with MDW, NDVI and TCAP.
Lastly, canopy cover has only been evaluated in studies during winter, with hypothesized relevance for access to understory forage and as thermal refuge from the cold (Poole and Heard 2003, Taylor et al. 2004, Turney and Roberts 2004, Taylor and Brunt 2007, Poole et al. 2009).Similarly, for interior mountainous regions during the summer, sites may be selected based on canopy cover to access forage species in the understory or in open canopy (Festa-Bianchet and Côte ´2008).Therefore, the percent canopy cover (COV) covariate was developed and included as a forage covariate with an expectation that mountain goat occupancy will decline at higher values of canopy cover.Canopy cover may also limit the ability to detect predators (Festa-Bianchet and Côte ´2008), contributing to predicted negative association.Alternatively, canopy cover may provide thermal refuge from the summer heat (see Heat load; Mysterud and Østbye 1999).The forage suite included all functional forms and spatial scales of COV, MDW, TCAP, and NDVI.Additive combinations of forage variables were also considered where possible (no collinearity) and biologically reasonable: COV with TCAP, COV with NDVI, MDW with TCAP, and MDW with NDVI.These were constrained to the same functional form and spatial scale within each model to maintain a suite of manageable size.TCAP and NDVI were considered as different constructs of similar forage ideas, so were not considered in the same models.This suite included 121 models.
Heat load.-Alteredbehavior for regulating body temperatures has been demonstrated in daily movement patterns of Alpine ibex (Aublet et al. 2009).Similarly, mountain goats may be sensitive to temperature during the summer, preferring sites with lower heat loads (e.g., sites with persistent snowbanks or deep shade).Several studies have shown preference for certain aspects, albeit mostly during winter (Foster 1982, Varley 1996, Gross et al. 2002, Poole and Heard 2003, Taylor et al. 2004, Poole et al. 2009).However, aspect measurements are difficult to interpret in relation to thermal tolerances given the effects of topographic shading and latitude.Heat load indices calculated from elevation raster data provide a more interpretable alternative that has not been evaluated in summer habitat models.This suite was developed to evaluate two heat load covariates that take into account topographic shading effects, the solar radiation index (SRI; after Ciarniello et al. 2005) and global radiation (GRAD; Fu and Rich 1999).The calculation of GRAD also takes into account effects of latitude.These two covariates were considered different constructs of similar heat load ideas, so were not considered in the same models.This suite included 31 models.
Combined.-Alladditive combinations of the resulting models ,5 (Q)AIC units of the mostsupported model and the uncombined top model from each suite were evaluated in this suite.For models combining the SLP with the RUG/RCK suite top models, only SLP or RUG was included due to collinearity between these variables.This suite included 90 and 448 models in AIC and QAIC selection procedures, respectively.While we had a priori predictions about covariates we were uncertain about scales and combinations.Thus, we ran a relatively large number of models and our analyses can be considered a mixture of a priori and exploratory modeling.Recent work on running all-possible combinations suggest that such approaches tend to find similar results to more restrictive a priori results (Murtaugh 2009, Doherty et al. 2012).Regardless, we consider our results to be important for helping to (1) refine thinking and (2) suggest important ideas for future confirmatory research.

Model validation
From the top model in the combined suite, we used the estimated coefficients for the covariates and their respective scale-specific rasters within a GIS (see DeVoe 2015 for development methods) to map a continuous surface of predicted occupancy probabilities for the northern GYA where established populations of mountain goats and independent location data (hereafter referred to as management data) from past aerial surveys performed by management agencies exist (unpublished data).This assigned each raster cell an occupancy probability between 0.0 and 0.419.To validate the model, we then classified the probabilities into very low (,0.0058),low (0.0058-0.023), medium (0.023-0.0458), and high (!0.0458) to represent suitability of habitat.Classifications were selected by first extracting predicted values for occupied sampling units and then designating a ''very low'' cutoff value that captured 5% of the occupied sampling units with the lowest predicted values.The remaining 95% was split evenly into the remaining categories (31.6% of occupied sampling units in each).This method was used because the distributions of probabilities for both the occupied and available sampling units were highly skewed right (Appendix A: Fig. A1a, b), which indicated that only a small proportion of the surveyed landscape (Appendix A: Fig. A1c) consisted of suitable habitat.We constrained the management data to the third week of June to the second week of October from 1968-2012 and to those recorded on 1:50,000 or finer scale maps, and overlaid them on the occupancy probability raster to extract values for each location.We then calculated the percentage of management data that fell within each of the habitat suitability classifications.Because errors in the spatial precision of the management data are inherent due to aerial survey methodology, we also calculated the distance from each management data point to the nearest cell containing at least low habitat suitability (i.e., !0.0058).The percentage of management data that fell within 100 and 200 m of these cells was calculated to indicate the prediction success of the model given the likelihood of a relatively small amount of spatial error in the management data.

GYA habitat suitability and abundance predictions
Using the same method for developing a prediction raster for the northern GYA, habitat suitability classifications for the entire GYA were derived, and the total area (km 2 ) of cells containing habitat greater than very low suitability was calculated.We then obtained recent estimates of abundance for two delineated areas periodically surveyed by management agencies that have well-established populations of mountain goats occurring at different densities (Appendix B: Fig. B1; unpublished manuscript).These survey areas included a northern (including portions of northern and north-eastern YNP) and south-western region of the GYA, where mountain goats have been present for about 50 and 40 years, respectively.Abundance estimates from each survey area were used to calculate densities of mountain goats per km 2 for any habitat with predicted suitability above very low.The lower density estimate from the northern and the higher density estimate from the southwestern survey area were used to project a bestguess low and high (respectively) abundance of mountain goats that could potentially be supported in the GYA given colonization of all habitat classified as at least low suitability.

Survey results
Between the third week of June and the second week of October 2011 to 2013, 551 total viewshed surveys documented detection and non-detection observations for 53,098 sampling units, of which 33,152 were surveyed by double-observers.The average number of sampling units surveyed from each viewshed survey point was 77 (ranging from 3 to 312), and the average viewshed survey duration was 53 minutes (ranging from 5 to 211 minutes).A mountain goat group was detected by at least one observer in 505 sampling units (0.92 observations per viewshed survey), of which 369 were recorded by one observer only (289 in double-observer and 80 in single-observer surveys) and 136 by both observers.Observations of the same group after initial detection were recorded in 195 sampling units.Group size averaged 3 individuals (ranging from 1 to 20 individuals), with 38.8% of groups observed resting, 38.6% feeding, 21.3% traveling, and 1.3% of other behaviors.

Model suites
The model from AIC model selection used for estimating ĉ included SLP 500 , SLP 500 2 , SLPv 500 ps , RCK 300 , RCK 300 2 , COV 100 , NDVI 100 , and GRAD 100 (Appendix C: Table C1), resulting in ĉ ¼ 3.46 from 1,200 simulations.This ĉ supported speculation of dependence in the data from multiple observations of the same group and was therefore used to adjust AIC scores and variances of model results for all within, among, and combined suite comparisons based on QAIC scores.
Sloped terrain.-Thetop four models supported similar relationships, including both SLP and SLPv, and were strongly supported (cumulative model weight ¼ 0.96; Table 2).The best-supported model included SLP 500 2 and SLPv 500 ps .The second best model was identical to the top model but included SLP 2 and SLPv ps at the 300-m instead of 500-m scale.The third and fourth best models were identical to the top model in scale but included quadratic and linear SLPv, respectively.Other models in the suite were not well supported by the data (DQAIC !5.8), and models with any DET covariate received little support.DET40 was ranked highest, at 10.3 DQAIC.Models containing covariates measured at scales other than 300-or 500-m received little support, with the nearest ranked model (covariates measured at the 700-m scale) at 14.6 DQAIC.
Ruggedness.-The top three models included both RCK and RUG, with the top model including both as pseudothresholds measured at the 300-m scale, and the second and third best models as quadratics measured at the 300-and 500-m scale, respectively.Other models in the suite were not well supported by the data (DQAIC !6.4).Models containing covariates measured at scales other than 300-or 500-m received little support, with the nearest ranked model (with covariates measured at the 100-m scale) at 8.2 DQAIC.The best ranked model with VRM had a DQAIC of 20.0.
Forage.-Only the model including linear COV 100 and NDVI 100 was supported by the data.
All other models in the suite were not well supported (DQAIC !5.6).The nearest ranked model containing covariates measured at scales other than 100-m was 9.3 QAIC worse than the top model.The nearest ranked models with TCAP and MDW were 24.9 and 25.4 DQAIC, respectively.
Heat load.-Theseven top models included GRAD at various scales and functional forms.The top three models contained GRAD 100 in the linear, pseudothreshold, and quadratic forms, respectively.The subsequent three supported models were identical but at the 300-m scale.The seventh best model included linear GRAD 500 .Other models in the suite were not well supported (DQAIC !5.6), with all GRAD models outranking SRI models.The best ranked model with SRI was 17.9 DQAIC.
Among suite comparison.-Thedata provided clear support for models from the sloped terrain suite, with the four top models (cumulative model weight ¼ 0.96) all being within 2.94 DQAIC.Top models from the ruggedness, heat load, and forage suites were greater than 52 DQAIC.
Combined.-Aftermodels with uninformative parameters were removed, four models within 5 QAIC units of the top model remained (Table 3).
Model selection uncertainty between these remaining models was a result of covariates considered at different functional forms; hence, there was clear consistency across the models for the covariates and scales supported.There was consistent support for quadratic SLP and linear COV and NDVI, as well as for the 500-m scale for SLP and SLPv and the 100-m scale for COV, NDVI, and GRAD.The top model included pseudothreshold SLPv and linear GRAD.The second and third next best models differed only by a pseudothreshold GRAD and a linear SLPv, respectively.The remaining model differed by inclusion of both these forms.The top model was considered adequate for inference due to the general similarities, as well as the consistency in coefficient estimates (Appendix C: Table C2), among competitive models.
In the top model, coefficient estimates parameterizing w followed expected predictions, indicating that w increased: (1) to a maximum at moderate levels of SLP (;36.88);(2) as SLPv and NDVI increased; and (3) as COV and GRAD decreased (Fig. 2; Appendix C: Table C2).The v www.esajournals.orgcoefficient estimate parameterizing p did not follow the expected prediction, indicating that p had no relationship with canopy cover.Estimates of p ranged from 0.43 to 0.52.In the post hoc evaluation of additional covariates on p, the intercept-only model was the most supported with all models within 3.5 DQAIC (Appendix C: Table C3).The top model estimated p at 0.47.Therefore, the exploratory modeling indicated that the top model used for inference was not ignoring important sources of heterogeneity in detection probability.

Model validation
Of 1,575 records in the management data, 27.9% (n ¼ 440), 28.2% (n ¼ 444), 21.9% (n ¼ 345), and 22% (n ¼ 347) were within areas identified as having very low, low, medium, and high habitat suitability, respectively (e.g., Appendix A: Fig. A2).Areas identified as having low suitability or better therefore contained 72.1% of the data, with 85.8% (n ¼ 1352) of the data falling within 100-m and 92.5% (n ¼ 1458) falling within 200 m of these areas.These values indicated that the majority of location data were in close proximity to predicted suitable habitat.

GYA habitat suitability and abundance predictions
The extent of the GYA in consideration (58,994 km 2 ) was estimated to contain 9,035 km 2 of habitat predicted as having at least low suitability (Fig. 3).The northern GYA survey region had an estimated mountain goat density of 0.59/km 2 based on approximately 1,942 km 2 of habitat Fig. 2. Predicted occupancy (w) and detection ( p) probabilities of sampling units with 95% confidence intervals for each covariate in the inferential model for habitat selection of mountain goats in the greater Yellowstone area, USA (each respective plot holding all other covariates constant at their mean value).
Table 3. Quasi-likelihood (QAIC) model selection results for combined models examining relationships of habitat covariates and occupancy (w) probability of sampling units for mountain goats in the greater Yellowstone area, USA.predicted as having at least low suitability and 1,154 estimated individuals.The southwestern GYA survey region had an estimated mountain goat density of 0.98/km 2 based on approximately 321 km 2 of habitat predicted as having at least low suitability and 316 estimated individuals.From these densities, we estimated the GYA might support 5,331-8,854 mountain goats given colonization of all habitat predicted at least low suitability.

DISCUSSION
This study used an occupancy design employing independent, double-observer surveys of a fine-scale grid for understanding spatial scalespecific habitat selection for mountain goats, a relatively rare species inhabiting a large, mountainous landscape.Through analysis of detection and non-detection observations collected from two mountain ranges across three summer seasons, we found evidence that mountain goat habitat selection occurred at two distinct spatial scales.Comparisons among model suites and effect sizes indicated that mountain goats more strongly selected for mean slope and slope variance terrain attributes at a broader scale of 500-m than for characteristics of canopy cover, heat load, and forage availability, which are selected at a finer scale of 100-m or less.Some uncertainty exists for selection of canopy cover, heat load, and forage at scales finer than the 100m spatial scale due to the size constraint of the sampling units.However, a clear distinction is evident between the spatial scales of the terrain covariates and the other covariates, supporting the hypothesis that mountain goats perceive different landscape attributes at different spatial scales.Similar results were obtained by Laforge et al. (2015) using GPS data from white-tailed deer in a resource selection analytical framework, suggesting evaluation of landscape covariates at multiple spatial scales may have general utility in habitat selection studies.The evidence for the importance of these two spatial scales likely indicates a hierarchal structuring of mountain goat habitat selection (i.e., mountain goats first select for slope attributes at a broader scale and then, within that slope habitat, secondarily select for finer scale habitat components).In this study we had no clear a priori ideas for the spatial scale mountain goats may perceive each landscape attribute, requiring exploration of a range of spatial scales for each covariate.Confirming the specific spatial scales of selection of each landscape covariate by mountain goats will therefore, require additional studies.
Our results support mountain goat selection for intermediate (optimal 378) and higher variability slopes.This highlights the advantage of using slope and slope variance over distance to escape terrain covariates that constrain selection to an arbitrary threshold value and assume equal selection of slopes greater than that threshold.The results indicate that not all steep slopes are alike; intermediate slopes and highly variable terrain are important components of terrain selected by mountain goats.The optimal slope value is similar to average slopes previously reported to be used by mountain goats, ranging from 33-418 (Foster 1982, Varley 1996, Poole and Heard 2003).Slope and slope variance covariates are strongly recommended in place of distance to escape terrain covariates for future modeling efforts and understanding of mountain goat habitat, and are likely applicable in habitat studies of other species associated with steep and variable terrain, such as bighorn sheep and other mountain ungulates.Distance-based covariates should be carefully considered to determine if less constraining and more biologically meaningful covariates exist (Martin and Fahrig 2012).
At the finer scale, the results support mountain goat selection for lower canopy cover, lower heat load, and higher values of NDVI.Based on effect sizes, selection occurred strongest for canopy cover and weakest for NDVI.Reduced canopy cover may be related to increased availability and productivity of shade-intolerant graminoids and forbs that are desirable as forage for mountain goats (Pallister 1974, Varley 1996).Alternatively, mountain goats may be minimizing the amount of immediate forest structure to increase visibility for predator detection, as demonstrated for nursery groups in interior Alberta, Canada (Festa-Bianchet and Côte ´2008).Mountain goats may also be attempting to thermoregulate by seeking cooler microclimates created by terrain features.Animals were commonly observed in deeply shaded aspects or couloirs or near persistent snow patches, as had been observed in other studies in the GYA (Varley 1996).These responses to warmer summer temperatures during the day appear to be similar to those of Alpine ibex (Capra ibex; Aublet et al. 2009) and desert ungulates in temperate systems (Cain et al. 2006), but contrast with Shafer et al. (2012) who found that mountain goats in coastal southeastern Alaska selected for sites with higher heat loads during the summer.This inconsistency may be due to regional climate differences.In the more northern, cooler climates, sites with higher heat loads may support the most vegetatively productive habitats and have less persistent snow, while in the warmer and drier interior temperate regions, lower heat loads provide cooler microsites for avoiding exposure to high temperatures and reducing metabolic water loss.These sites are also likely to retain snowpack longer, providing foraging opportunities as the snow slowly recedes over the summer and new plant growth becomes available (Varley 1996).Higher NDVI values are likely capturing these and other fine-scale sites containing more palatable vegetation irrigated by melting of persistent snowpack or by spring seeps.
Mountain goats are morphologically adapted to and tightly associated with steep and complex terrain.This adaptation has been primarily perceived as a strongly selected evolutionary strategy to minimize predation risk, presumably based on observations that mountain goats quickly move to steep and complex terrain in the presence of predators and hypotheses that mountain goats can be found further from steep and complex terrain in systems lacking predators (Adams et al. 1982, Fox andStreveler 1986).Steep and complex terrain undoubtedly provides escape opportunities from certain predators.However, mountain goats are likely adapted to these sites for a host of reasons, including, reduced effects of interference and exploitative competition with other ungulates like bighorn sheep, increased access to a diversity of microsites for foraging, or increased access to cool, shaded environments for thermal regulation.Terrain with high variance in slope at a broader spatial scale may be particularly important for providing sufficient foraging and shade opportunities at a finer scale.Therefore, it is reasonable that mountain goat association with steep and complex terrain is not simply the result of a predation risk trade-off strategy, especially con-sidering that predation risk by mountain lions and golden eagles may be elevated compared to gentler slopes and more homogenous terrain.Rather, this association is probably a function of an assemblage of environmental characteristics created by the terrain that provides the necessary resources for survival and reproduction to which mountain goats are specifically adapted (i.e., their ecological niche).Future interpretation and modeling of habitat selection for mountain goats and other animals associated with similar terrain may be improved by the de-emphasis of ''escape terrain'' terminology.
Adequately modeling detection probability was particularly challenging using the methods employed by this study.We found little evidence for an association of canopy cover or terrain complexity with detection probability, a relationship that has been previously found in studies from aerial surveys (Poole 2007, Rice et al. 2009).
Based on the inferential model, estimates of detection probability, ranging from 0.43 to 0.52, were roughly constant regardless of canopy cover.Additionally, the post hoc evaluation of covariates on detection resulted in considerable model selection uncertainty, with the mostsupported model predicting constant detection at approximately 0.47.This result may have been due to lack of relevant covariates or dependence among the data as a result of subsequent recording of detections as groups traveled across sampling units during the viewshed surveys.The subsequent detections may have eroded any signal of the effect of canopy cover or terrain complexity because observers were more likely to record a group in forested or more rugged areas if it had been previously detected.In addition, we suspect detection probabilities may have been biased due to two possible reasons: (1) non-simultaneous viewing of each sampling unit by double observers while groups traveled across the viewshed; and (2) unmodeled variation in detection probability from lack of relevant covariates (MacKenzie et al. 2006).As a consequence, it is likely that estimates of occupancy probabilities were also biased (MacKenzie et al. 2006, Rota et al. 2009).However, our model predicted mountain goat habitat corresponded well with an independent dataset describing mountain goat distributions in the northern GYA.If detections of groups were scaled up to v www.esajournals.orgthe viewshed, average detection probability was approximately 0.65, similar to estimates found in previous studies from aerial surveys in other regions (0.59-0.69;Houston et al. 1986, Poole et al. 2000, 2007, Gonzalez-Voyer et al. 2001, Pauley and Crenshaw 2006, Rice et al. 2009).
Future applications of these methods may consider increasing the size of the sampling units to reduce the chance of animal movement among sampling units during surveys, as well as enacting stricter survey protocols to increase the chance that observers look in the same sampling units at the same time.This may be possible with larger sampling units or smaller viewshed extents.Trade-offs of these strategies would include reduced sample sizes (with potential increased model convergence issues) as well as less discernment of habitat associations at finer scales.Although modeling detection was a limitation in this study, the inferential model performed well in predicting observations from an independent data source with a reasonably conservative spatial error of 200 m across a broad geographic area, lending support for a reasonable prediction of habitat selection.Further model validation is being conducted by other researchers in Grand Teton and Glacier National Parks.
Finally, we present a landscape-wide prediction for the spatial distribution and abundance of mountain goats for the GYA given colonization of all suitable summer habitats.The model depicted clear spatial delineation of areas projected to be suitable habitat.About 18% of the GYA extent in consideration contained predicted suitable habitat which was primarily associated with heterogeneous terrain and indicated varying extents and connectivity of suitable patches.The currently colonized areas of the GYA, including the northern GYA and portions of the Snake River Range (approximately the southwestern GYA survey region), contained about 43% of the total predicted suitable habitat.The remaining, un-colonized area was projected to be 6,131 km 2 , or about 1.3 times larger than the area currently colonized.The best-guess abundance of mountain goats potentially supportable by the suitable habitat in the GYA ranged from 5,331 to 8,854.Given the most recent abundance estimate for the GYA of 2,354 mountain goats (unpublished manu-script), these figures suggest that when mountain goats have expanded their range throughout the GYA, they may be 2.5 to 4.2 times more abundant than currently.
Suitable habitat overlaps considerably with the general summer distribution of bighorn sheep observations throughout the GYA with approximately 75% of all bighorn observations within the GYA occurring within areas defined as suitable mountain goat habitat (e.g., Appendix B: Fig. B2; unpublished data).Limited observations of interactions between the two species have suggested interference competition with mountain goats dominating (Reed 1986).However, the biological effects of any potential competition are largely unknown.Adams et al. (1982) and Varley (1996) suggest that bighorn sheep and mountain goats in native sympatric areas occupy separate ecological niches that overlap only marginally, alleviating potential competition for scarce resources.Adams et al. (1982) hypothesized, however, that non-native mountain goat populations occupied a broader niche during rapid population growth and range expansion in Colorado, where native bighorn sheep were absent.Conceivably, the bighorn sheep niche may also be broader in the absence of mountain goats due to lack of competition.This reasoning indicates that in situations where non-native mountain goats are expanding into occupied bighorn sheep ranges competitive overlap may be increased during the initial expansion phase of mountain goats, but may be alleviated as populations of both species adjust and stabilize, reflecting reduced niche overlap likely present in native sympatric populations.Given that the data for this study were collected from areas with well-established, sympatric populations, the prediction of suitable habitat likely reflects the latter condition in the GYA and inference may be stronger as mountain goats continue to expand into regions currently occupied exclusively by bighorn sheep.The results may be conservative in areas where mountain goats are predicted to be allopatric and a broader niche may be sustained due to lack of competition.
In addition, the realized occupancy and abundance of suitable summer habitat may be constrained due to competition for scarce resources in shared winter habitats, such as on high v www.esajournals.orgelevation, wind-swept ridges, limiting population sizes in certain regions of the GYA (Adams et al. 1982, Varley 1996).It is unknown to what extent such conditions existed in areas surveyed during this study, and there is likely large variation in the amount of shared winter habitats throughout the GYA.This study also did not take into account the potential effects of different regional climatic conditions on patterns of habitat selection and abundance of mountain goats.For example, the southern, more arid regions of the GYA may induce higher metabolic demand for thermoregulation and water (Cain et al. 2006), prompting greater preference for sites with lower heat loads, more canopy cover, or more succulent vegetation.Additionally, regions that experience harsher winter conditions (e.g., deep snow accumulations) and have less available winter habitat (steep, snow-shedding slopes or wind-swept ridges) may reduce the number of mountain goats occupying the predicted suitable summer habitat.Lastly, areas that are predicted as suitable habitat but have no historic observations of mountain goats from management agency surveys over 50 years after mountain goats were introduced into the northern GYA (e.g., the West Boulder Plateau in the northern GYA) may indicate that an important habitat covariate not considered in this study was not captured in the inferential model.Therefore, suitable habitat and abundance estimates of mountain goats should be interpreted carefully.
Future research focused on understanding the ecological niches and overlap of bighorn sheep and mountain goats in both summer and winter ranges in multiple regions of the GYA would provide valuable insight into the dynamics and effects of potential competition.Ongoing studies across the GYA using global positioning system (GPS) technologies are collecting such data and will provide a comparative method with new insights into patterns of habitat selection for both species, as well as further insights into potential competition between these two mountain ungulates.

Fig. 1 .
Fig. 1.The Gallatin Crest and North Absaroka study areas in the northern portion of the greater Yellowstone area, USA.YNP refers to Yellowstone National Park, while GTNP refers to Grand Teton National Park.

Fig. 3 .
Fig. 3. Predicted summer habitat suitabilities of mountain goats for the greater Yellowstone area, USA, overlaid by data collected from 1968-2012 by management agency personnel showing the general, current distribution of mountain goats.

Table 1 .
Covariates (developed at the 100-, 300-, 500-, 700-, and 900-m spatial scales excluding DET covariates) considered for modeling occupancy (w) probability arranged in model suites to evaluate mountain goat habitat selection in the greater Yellowstone area, USA.

Table 2 .
Quasi-likelihood (QAIC) model selection results for models arranged in model suites examining relationships of habitat covariates and occupancy (w) probability of sampling units for mountain goats in the greater Yellowstone area, USA.The top models within 5 DQAIC of the top-ranked model for each suite are presented.Covariate subscripts indicate spatial scale and superscripts indicate functional form (2 ¼ quadratic, ps ¼ pseudothreshold, [none] ¼ linear).Detection ( p) probability for each model was parameterized using COV 100 ps .Model variables defined in Model development.K ¼ no. of parameters, wt i ¼ QAIC weight.Quasi-likelihood for the top model from the across-suite comparison was À859.59.Results based on ĉ ¼ 3.46.
The top competitive models within 5 DQAIC of the top-ranked model are presented.Subscripts and superscripts are as in Table2.K ¼ number of parameters, wt i ¼ QAIC weight.Detection ( p) probability was parameterized using COV ps 100 .ID number for comparison of model coefficient estimates in Appendix C: TableC2.Quasi-likelihood for the top model was À848.634.Results based on ĉ ¼ 3.46.