The roles of habitat and intraguild predation by coyotes on the spatial dynamics of kit foxes

Intraguild predation (IGP) by a dominant predator can drive the spatial dynamics of a subordinate predator and may explain space-use patterns that deviate from theoretical predictions that species will use areas that maximize the availability of limited resources (resource availability hypothesis). Intraguild predation may suppress the distribution and abundance of mesopredators, but spatial resource partitioning may facilitate coexistence, with the subordinate carnivore utilizing suboptimal habitats. In arid systems, free-standing water was historically scarce, limiting the distribution of larger-bodied predators and offering large areas of refugia for smaller, arid-adapted species, such as the kit fox (Vulpes macrotis). In these systems, the development of artificial water sources may facilitate an increase in the distribution and abundance of larger carnivores (e.g., coyotes [Canis latrans]), perhaps to the detriment of kit foxes. We coupled noninvasive genetic sampling and dynamic occupancy models to evaluate the spatial dynamics of kit foxes and their intraguild predators, coyotes, in western Utah, United States. We evaluated the influence of habitat characteristics on coyote occupancy patterns, and then investigated the role of habitat and coyotes on kit fox space use at multiple scales. Coyote occupancy was unrelated to water availability, but was positively related to the proportion of shrubland and woodland cover, a pattern consistent with predictions of the resource availability hypothesis. Supporting predictions of IGP theory, kit fox occupancy was negatively related to shrubland and woodland cover, minimizing overlap with land-cover types favoring coyote occupancy. Furthermore, kit fox probability of local extinction was positively related to coyote activity. Interestingly, kit fox detection was positively related to coyote activity (i.e., kit fox detection was higher on spatial surveys with greater coyote sign), suggesting that at finer scales, kit foxes utilized riskier habitats to secure sufficient resources. Our results identified two alternative states predicted by IGP theory (i.e., intraguild predator dominated and coexistence of intraguild predator and intraguild prey) in a single system and elucidated the importance of considering dynamic processes and scale when investigating IGP.


INTRODUCTION
The resource availability hypothesis predicts that species will utilize areas that maximize limiting resources (Ernest et al. 2000, Blaum et al. 2007), but competition and predation can influence space use and may explain species distributions that deviate from expectations based on resources alone (Schoener 1974, Heithaus 2001, Thompson andGese 2007, Vanak et al. 2013).
Among carnivores, interactions are commonly characterized by intraguild predation (IGP), where dominant and subordinate predators (hereafter intraguild [IG] predator and IG prey, respectively) compete for shared resources, but IG predators also kill IG prey (Polis et al. 1989, Holt and Polis 1997, Verdy and Amarasekare 2010. Traditional IGP theory predicts that in systems involving an IG predator, an IG prey, and a shared resource, species persistence will be determined by the resource-ratio hypotheses (i.e., the R Ã rule; Polis 1997, Holt andHuxel 2007). Holt and Polis (1997) defined R Ã as the equilibrium density of the shared resource. Assuming the IG prey is a superior exploitative competitor, only the IG prey is predicted to persist at low resource densities (R Ã N ). Conversely, at high resource densities (R Ã P ), the IG predator is predicted to be sufficiently abundant to exclude the IG prey Polis 1997, Verdy andAmarasekare 2010). Alternative stable states of IG predator and IG prey coexistence are predicted at intermediate resource levels (R Ã NP ; where R Ã N < R Ã NP < R Ã P ), but generally require the IG prey be the superior competitor and the IG predator gain significantly from consuming the IG prey Polis 1997, Verdy andAmarasekare 2010).
Coexistence of IG predator and IG prey is common among mammalian carnivores (Palomares andCaro 1999, Kimbrell et al. 2007). Mammalian IGP is often unidirectional and characterized as an extreme form of interference competition, where the IG predator kills but does not regularly consume the IG prey, a characteristic in traditional IGP theory that destabilizes coexistence and leads to competitive exclusion (Holt and Polis 1997, Palomares and Caro 1999, Kimbrell et al. 2007, Lourenc ßo et al. 2014). Extensions to traditional IGP theory have explored how behavioral complexities (e.g., IG prey avoidance of IG predators, vigilance) and/or alternative prey resources influence coexistence (Heithaus 2001, Holt and Huxel 2007, Kimbrell et al. 2007). When alternative basal resources are available to the IG predator, theoretical predictions suggest that coexistence of the IG predator and IG prey can occur without the IG predator benefiting directly from the consumption of the IG prey (Heithaus 2001, Holt andHuxel 2007). Under asymmetrical IGP, IG predators are predicted to select habitats that align with resource availability, while IG prey are predicted to balance risk (i.e., behavioral avoidance of IG predators) and resource availability (Heithaus 2001). Similarly, vigilance of IG prey and basal resources may also facilitate coexistence in the absence of IG prey consumption. When the IG predator and IG prey share multiple basal resources (i.e., prey), relatively large attack rates on the IG prey are predicted to decrease IG prey foraging efficiency (via increased vigilance), reduce competition (via interspecific killing), and improve IG predator efficiency (via reduced basal prey vigilance from lowered predator densities; Kimbrell et al. 2007).
We investigated the spatial dynamics of the kit fox (Vulpes macrotis), a mesocarnivore native to North American deserts, at a site that has experienced a relatively recent increase in the distribution and abundance of a dominant IG predator, the coyote (Canis latrans; Arjo et al. 2007, Kozlowski et al. 2012. Kit foxes have declined and are of conservation concern across much of their range . Arid environments may naturally exclude or limit large mammalian predators with high water demands, offering refuge for arid-adapted mesocarnivores. However, water developments have become common within deserts and may facilitate increases in previously rare, larger-bodied carnivores, thereby altering the composition of carnivore communities and potentially increasing the frequency or magnitude of interactions between IG predators and IG prey (Kamler et al. 2003, Arjo et al. 2007, Atwood et al. 2011, Brawata and Neeman 2011. At our study site, the U.S. Army Dugway Proving Ground, Utah, United States (hereafter Dugway), artificial water sources have been developed since the 1960s (Arjo et al. 2007). Over this same period, coyotes have usurped kit foxes as the most abundant predator (Arjo et al. 2007, Lonsinger 2015. We coupled noninvasive genetic sampling and dynamic occupancy models to investigate the relative influences of habitat characteristics (e.g., water availability) and IGP (i.e., coyotes) on the spatial dynamics of kit foxes. Occupancy modeling allows researchers to estimate the probability a species occurs within a site, while accounting for imperfect detection, variation in environmental covariates (MacKenzie et al. 2002(MacKenzie et al. , 2003(MacKenzie et al. , 2006, and species interactions (Richmond et al. 2010). The application of dynamic models provides an opportunity to better understand how habitat characteristics and IGP influence the dynamic processes of colonization and local extinction that drive occupancy patterns (MacKenzie et al. 2003). We hypothesized that coyote occupancy would be influenced by water availability and shrubland and woodland cover (Kozlowski et al. 2012), aligning with the resource availability hypothesis and predictions of IGP theory that IG predators will select habitats to maximize resource availability (Ernest et al. 2000, Heithaus 2001). We also hypothesized that kit fox occupancy would be negatively associated with coyotes and habitat characteristics favored by coyotes (Nelson et al. 2007, Thompson and Gese 2007, Kozlowski et al. 2012. Consistent with predictions of IGP theory, we expected lower colonization and higher local extinction probabilities for kit foxes at sites also occupied by coyotes (Polis et al. 1989, Holt andPolis 1997).

Study area
The study area encompassed~3015 km 2 of Great Basin Desert in western North America, including portions of Dugway and surrounding federal lands (Fig. 1). Land cover was characterized by cold desert playa (primarily pickleweed Habitat classifications display the distribution of shrubland and woodland (SW) cover versus areas with lower (e.g., grasslands) or sparse (e.g., playa) vegetative cover. Water points indicate the locations of free-standing water sources (e.g., spring, tanks).
[Allenrolfea occidentalis]), chenopod shrubland (Atriplex confertifolia and Kochia americana dominated), and vegetated and unvegetated dunes at low elevations. Sagebrush (Artemisia spp.) shrubland and open juniper (Juniperus osteosperma) woodland were found at higher elevations. Greasewood (Sarcobatus vermiculatus) shrubland was distributed across elevations. Invasive grasslands, primarily cheatgrass (Bromus tectorum), were common at lower elevations. Winters were cold (January-February mean high = 2.6°C) and summers hot (July-August mean high = 34.9°C); mean annual precipitation was 17.4 cm. Elevation ranged from~1200 to 2154 m. Water was historically restricted to natural springs located primarily in the mountains (Arjo et al. 2007, Hall et al. 2013), but has since become widespread due to increases in artificial water sources (Fig. 1).

Field surveys and species identification
We randomly selected 103 sites across the study area using ArcGIS 10 (ESRI, Redlands, California, USA) for monitoring kit fox and coyote occupancy. Sites were 6.25 km 2 (2.5 9 2.5 km), an area similar to the average home ranges reported for both kit foxes (2.5-11.6 km 2 ; List and Cypher 2004) and coyotes (5.5-6.9 km 2 ; Gese et al. 1988, Nelson et al. 2007). Within each site, we established four 500-m transects, representing spatially replicated surveys, along dirt or gravel roads. While the use of spatial replicates, instead of temporal replicates, has been debated in the occupancy literature, utilizing this technique is unlikely to cause bias for highly mobile species (Kendall and White 2009, Guillera-Arroita 2011, Harris et al. 2014. We conducted surveys during four seasons, namely two winters (14 January-6 March 2013; 13 January-19 March 2014) and two summers (12 July-16 August 2013; 10 July-21 August 2014). We visited each site once per season and surveyed each transect for carnivore scats. When a scat was detected, we collected~0.7 mL of fecal material following procedures of Lonsinger et al. (2015).
We performed DNA extraction and polymerase chain reaction amplification in a laboratory dedicated to low-quality samples to minimize contamination risk. We determined species identification using mitochondrial DNA fragment analysis (mtDNA; De Barba et al. 2014) following methods detailed in Lonsinger et al. (2015). Samples that failed to amplify, mixed samples (i.e., amplified mtDNA for >1 species), or samples identified as a non-target species were excluded from subsequent analyses.

Occupancy modeling
We planned to use a dynamic, co-occurrence occupancy model (Richmond et al. 2010, Yackulic et al. 2014 to evaluate the relative roles of habitat and IGP on kit fox spatial dynamics. Assuming coyote spatial dynamics was independent of kit fox occurrence, we first employed a single-species dynamic model (MacKenzie et al. 2003) to estimate coyote occupancy (w) and identify factors influencing their spatial dynamics (probabilities of colonization [c] and local extinction [e]) and detection probabilities (p). Retaining the best-supported coyote model, we planned to explore the influence of coyote occurrence and environmental covariates on kit fox occupancy, detection, and dynamics via the co-occurrence model. Unfortunately, coyote occupancy estimates were very high across seasons (see Results and Appendix S1), effectively eliminating our ability to evaluate co-occurrence patterns (Richmond et al. 2010) and processes (Yackulic et al. 2014, Dugger et al. 2016. Coyotes are ubiquitous in our system; if kit foxes occurred at a site, they co-occurred with coyotes. Consequently, we used a dynamic single-species occupancy model to evaluate whether the spatial dynamics of kit foxes were influenced by environmental covariates or intensity of coyote activity, exploiting the variation in coyote activity at the site and transect levels. We interpreted variation in detection probability (i.e., among transects at occupied sites) of kit foxes as reflecting differences in fine-scale space use (i.e., a behavioral response). Because the occupancy dynamics analyses of coyotes did not inform our kit fox analyses, we limited the associated methods and results to Appendix S1.
We identified environmental covariates expected to influence detection, occupancy, colonization, and local extinction of kit foxes. Road characteristics can influence scat persistence (Lonsinger et al. 2016) and detection . We characterized the road type of each transect as (1) unmaintained two-track road, or maintained (2) single-lane or (3) twolane gravel road (sensu Lonsinger et al. 2016). Scat detection may also be influenced by the ❖ www.esajournals.org presence of snow, survey date, and/or survey time (Harris et al. 2014); we recorded these covariates during each survey. Snow may reduce detection if scats become covered. Date may further influence detection if canid activity changes throughout winter (e.g., during reproduction) or summer (e.g., increased juvenile activity; Ralls et al. 2010). Survey time may influence visibility and was standardized across seasons as time from solar noon.
We identified four environmental covariates that may directly or indirectly influence spatial dynamics of kit foxes at the site level. We expected soil to influence kit fox occupancy, as they utilize burrows year-round (Arjo et al. 2003, Kozlowski et al. 2008; soil layers were obtained from the Utah Automated Geographic Reference Center (http://gis.utah.gov/) and were reclassified into four categories (silt, fine sand, blocky loam, and gravel; sensu Dempsey et al. 2015). Water availability may influence canid space use (Arjo et al. 2007, Hall et al. 2013. Whereas kit foxes do not require free-standing water (Golightly and Ohmart 1984), coyotes may (Arjo et al. 2007), and therefore, kit fox occurrence may be indirectly influenced by water. Perennial water sources were identified using Dugway's Natural Resource Program GIS layers. Additionally, we utilized Google Earth imagery to identify the convergence of livestock and horse trails, then ground-truthed these points to locate additional water sources. For each site, we characterized water in three ways: (1) distance to nearest water, and the number of water sources within (2) 2.5 km and (3) 5 km from the site center. Data on prey densities and diversity were not available across sites, but are associated with habitats in our study area with shrubland and woodland habitats typically supporting higher abundance of small mammals and leporids (Arjo et al. 2007, Kozlowski et al. 2012. Although kit fox historically occupied shrubland habitats at our study area (e.g., greasewood flats, vegetated dunes; Egoscue 1962), recent research suggests that kit foxes may minimize overlap with coyotes and therefore avoid prey-rich shrubland and woodland habitats (Kozlowski et al. 2012). We used 2012 LANDFIRE vegetation data (http://la ndfire.cr.usgs.gov/) to calculate the proportion of shrubland and woodland cover (%SW) within in each site, with higher %SW presumably representing relatively higher prey availability (Kozlowski et al. 2012) and greater thermal cover (Blaum et al. 2007). Road density may also influence canid detection and/or occupancy. We obtained road layers from the Utah Automated Geographic Reference Center and calculated road density for each site. We processed all GIS layers with ArcGIS 10.
Finally, we expected kit foxes to minimize the potential for negative interactions with coyotes. To evaluate the influence of coyotes on kit fox spatial dynamics, we characterized coyote activity at the site and transect levels. Site-level coyote activity was characterized as (1) the total number of coyote scats detected and (2) the total number of transects on which coyotes were detected. Transect-level coyote activity was characterized as (1) the number of scats detected or simply (2) the detection or non-detection of coyotes.
We evaluated correlations among covariates with a Kendall's rank correlation test. Only the three characterizations of water availability were correlated with one another (r > |0.48|) and we never included >1 water characterization in a given model. We used a structured modeling approach, first identifying the best global model, then sequentially fitting simpler model structures for p, initial occupancy (w 1 ), and the dynamic parameters (e and c, together). We considered global structures for p containing road type, road density, presence of snow, date, sun (i.e., difference between survey time and solar noon), transect-level coyote activity, and temporal variation among seasons. We considered global structures for w 1 , e, and c that contained soil type, % SW, road density, water availability, and sitelevel coyote activity; for e and c, site-level coyote activity reflected that from the preceding season. Global structures for e and c also included temporal variation among seasons. To identify the best global structure, we fit models with all possible combinations of those covariates with >1 characterization (i.e., road type, water availability, and site-and transect-level indices of coyote activity). The most parsimonious characterization of each predictor was retained for subsequent analyses (Appendix S2: Table S1).
After identifying the best global model structure for kit foxes, we fit all possible combinations of predictors for p, while maintaining the global structures for w 1 , e, and c, to identify the best detection model. Next, using the most parsimonious structure for p and the global structures for ❖ www.esajournals.org e and c, we fit all possible combinations of predictors for w 1 to identify the best occupancy structure. Finally, we retained the most parsimonious structures for p and w 1 and simultaneously evaluated models for e and c, considering all possible combinations of predictors. We conducted all analyses with program MARK (White and Burnham 1999), employing an informationtheoretic approach (Akaike's Information Criterion with small sample size correction, AIC c ) to compare the relative fit of models and cumulative Akaike weights to evaluate predictor importance (Burnham and Anderson 2002).
To further explore the influence of coyotes (see Results) on kit foxes, we calculated the equilibrium occupancy (w Eq = c/[c + e]; MacKenzie et al. 2002) for kit foxes on sites with varying levels of coyote activity (i.e., 0-12, by 2). We based w Eq on estimates of colonization and local extinction probabilities for kit foxes from the final open period (winter 2014 to summer 2014), the most parsimonious model structure, mean values for other numeric covariates, and sites characterized by silty soils (the soil with the highest colonization probability; see Results).

Field sampling, species identification, and site and transect characteristics
Sampling effort was constant across seasons, with four 500 m transects being surveyed within each of 103 sites per season, resulting in 824 km of surveys. We collected 1702 scat samples, of which 64% were coyote and 18% were kit fox. We also detected domestic dog (<1%), red fox (V. vulpes; 1%), bobcat (Lynx rufus; 2%), and cougar (Puma concolor; <1%), and 15% of samples failed or were mixed (Table 1). Across seasons, na€ ıve estimates of occupancy were higher for coyotes than for kit foxes (Table 1, Fig. 1).

Patterns of occupancy and spatial dynamics
The best global model structure for kit foxes included the distance to nearest water and an ordinal relationship among road types (Table 2; Appendix S2: Table S1). The influence of coyote activity on occupancy dynamics and detection of kit foxes was best characterized as the total number of coyote scats detected at the site and transect levels ( Table 2; Appendix S2: Table S1). The best detection structure suggested that detection probability (p) of kit foxes was positively related to transect-level coyote activity (b ¼ 0:20, c SE ¼ 0:06, 95% CI = 0.08, 0.32) and varied by road type (b ¼ 0:23, c SE ¼ 0:15, 95% CI = À0.06, 0.51). Mean kit fox p  Occupancy was lower for kit foxes than for coyotes (Fig. 2). The best structure for kit fox initial occupancy suggested that %SW had a strong negative (b ¼ À13:46, c SE ¼ 3:98, 95% CI = À21.26, À5.66) influence (Fig. 3). In contrast, %SW was positively associated with coyote occupancy ( Fig. 3; Appendix S1). For kit fox w 1 , the cumulative weight for %SW was high and no other predictors carried substantial weight (Table 2; Appendix S2: Table S3). Dynamic parameters for kit foxes were influenced by coyote activity; the best model structure included a positive effect on both e (b ¼ 0:97, c SE ¼ 0:45, 95% CI = 0.09, 1.85) and c (b ¼ 0:23, c SE ¼ 0:15, 95% CI = À0.05, 0.52), though the effect on c was weaker (Table 2;  Appendix S2: Table S4). Additionally, e varied temporally among seasons and soil type influenced c (Appendix S2: Table S4). For c, sites characterized by silty soils had the highest probability of colonization, while sites with sandy or loamy soils had the lowest colonization probabilities. As predicted, kit fox occurrence was more stable on sites with lower coyote activity, with less change in occupancy status (i.e., turnover) between seasons (Fig. 4).
Among sites characterized by silty soils, estimates of e and c ranged from 0.12 to 0.99 and 0.23 to 0.83, respectively; both increased with increasing site-level coyote activity. For the range of coyote activity levels evaluated, kit fox w Eq was similar among sites with coyotes, ranging from 0.33 to 0.45. When no coyotes were detected, kit fox w Eq was substantially higher at 0.65.

DISCUSSION
Predators are elusive, wide ranging, and occur at low densities (Palomares andCaro 1999, Gompper et al. 2006). Mammalian IGP systems are challenging to investigate, as these attributes often apply to both the IG predator and IG prey. Indirect species sign (e.g., scat) is often conspicuous and noninvasive monitoring alleviates many of the challenges of detecting carnivores, facilitating multispecies monitoring (Gompper et al. 2006). Imperfect detection of sign can bias inferences regarding species occurrence (MacKenzie et al. 2002). Occupancy modeling explicitly addresses this concern and dynamic models improve our understanding of complex systems by providing insights into the processes of local extinction and colonization driving observed occupancy states (MacKenzie et al. 2002(MacKenzie et al. , 2003(MacKenzie et al. , 2006. Notes: Predictors: DistW, distance to nearest water; W2, number of water sources within 2.5 km of site center; W5, number of water sources within 5 km of site center; RTO, ordinal road type; RTC, categorical road type; RD, road density; snow, presence or absence; sun, difference between survey time and solar noon; date, days since surveys were initiated within season; SW, proportion of shrubland or woodland habitat within a site; soil, categorical classification of majority of soil types for a site (silt, fine sand, blocky loam, or gravel); CS, total number of coyote scats detected within a site; CT, total number of transects on which coyotes were detected within a site; CA, number of coyote scats detected at the transect level; CP, detection or non-detection of coyotes at the transect level; t, temporal variation among seasons. See Appendix S2 for the complete model sets used to evaluate each parameter. Bold indicates predictors in the best model. Global represents the evaluation of different characterizations for water, road type, and site-and transect-level coyote activity to identify a best global model. [S14]) occupancy probabilities with 95% confidence intervals for kit foxes (Vulpes macrotis) and coyotes (Canis latrans) in Utah, United States, 2013-2014. Both kit fox and coyote occupancy probabilities are plotted based on their best model structure, the median proportion of shrubland and woodland cover (13.0%), and mean values for other numeric covariates. Soil type (a categorical covariate) was present in the best kit fox model structure, and kit fox occupancy probability is plotted separately for each soil type. Among 103 sites, 44.7%, 34.9%, 11.7%, and 8.7% were characterized as predominantly silt, fine sand, blocky loam, and gravel, respectively. Fig. 3. Initial occupancy probabilities for coyotes (Canis latrans) and kit foxes (Vulpes macrotis) as a function of shrubland and woodland cover in Utah, United States, 2013-2014. Occupancy probability is plotted based on the best model structure for each species using the mean values for other numeric covariates. Soil type (a categorical covariate) was present in the best kit fox model structure and kit fox occupancy probability is plotted in the soil with the highest probability of colonization (i.e., silty soils).
Canid systems consisting of foxes (IG prey) and a larger IG predator (e.g., coyotes, jackals [C. mesomelas], dingos [C. lupus dingo]) have become model systems for investigating mammalian IGP (Nelson et al. 2007, Thompson and Gese 2007, Brawata and Neeman 2011, Kozlowski et al. 2012, Robinson et al. 2014. While co-occurrence modeling offers a framework to investigate IGP (e.g., Robinson et al. 2014), if either species is widely distributed, insufficient heterogeneity may exist in occupancy to effectively evaluate competitive exclusion. In our system, coyote occupancy was high and we demonstrated how variation in coyote sign can be exploited to investigate the influence of IGP on spatial dynamics of IG prey using singlespecies dynamic occupancy models. To our knowledge, our study was the first to couple noninvasive genetic sampling and dynamic occupancy models to explore a mammalian IGP system.

Patterns and processes influencing canid occurrence and space use
In arid environments, increased water can reduce physiological stress, increase survival, and facilitate persistence of large predators (Brawata and Neeman 2011), increasing the potential for conflict with IG prey (Atwood et al. 2011). Consequently, we expected coyote occupancy to be higher in sites with greater water availability and for water to influence kit fox occupancy through an indirect effect of influencing coyote space use. We did not detect a relationship between coyote or kit fox occupancy and water (Appendices S1 and S2). Instead, our results supported recent research suggesting that despite use of free-standing water, space use of coyotes at Dugway was not restricted by water (Hall et al. 2013, Kluever and.
Empirical observations suggest that canids commonly employ spatial partitioning to facilitate coexistence, with IG predators conforming to predictions of the resource availability hypothesis (Ernest et al. 2000, Blaum et al. 2007), while IG prey occupy habitats that minimize risk of IGP, aligning with expectations of IGP theory and mesopredator suppression (Soul e et al. 1988, Heithaus 2001. For example, swift foxes (V. velox) selected habitats that minimized risk of IGP by coyotes, which occupied resource-rich habitats (Thompson and Gese 2007). Similarly, coyotes displaced endangered San Joaquin kit foxes (V. m. mutica) from prey-rich shrublands Fig. 4. Mean turnover in occupancy state between seasons for kit foxes (Vulpes macrotis) as a function of sitelevel coyote (Canis latrans) activity in Utah, United States, 2013-2014. Mean turnover was plotted based on the best model structures and mean values for other numeric covariates. Soil type (a categorical covariate) was present in the best kit fox model structure and mean turnover is plotted for sites characterized predominantly as the soil with the highest probability of colonization (i.e., silty soils). (Nelson et al. 2007). Our results aligned with these findings and our predictions. At broad scales, coyote occupancy was positively associated with %SW cover. In contrast, kit fox occupancy was negatively associated with %SW, despite evidence that both historical (Egoscue 1962) and predicted (Kozlowski et al. 2012) distributions of kit foxes included shrubland habitats (e.g., greasewood shrublands) in this system. Shrubland and woodland habitats at Dugway supported relatively high prey resources (Arjo et al. 2007, Kozlowski et al. 2012, and the higher vegetative structure likely offered greater thermal cover for larger-bodied predators than other habitats (Blaum et al. 2007). In contrast, lower vegetation enhances visibility and predator detection for kit foxes (Arjo et al. 2003. Under traditional IGP theory, increases in resource enrichment (R*) tend to shift systems from IG prey dominated (at R Ã N ) to coexistence (at R Ã NP ) and IG predator dominated (at R Ã P ). Kit foxes and coyotes in New Mexico exhibited two alternative stable states in accordance with IGP theory: (1) only kit foxes in resource-poor nonshrubland habitats and (2) kit fox and coyote co-occurrence in more resource-abundant shrublands (Robinson et al. 2014). Robinson et al. (2014) suggested coyotes were likely precluded from selecting resource-poor habitats with insufficient resources in New Mexico. Our results suggested that no habitat had insufficient resources (e.g., water and prey) to support coyotes at Dugway. Instead, we observed kit fox and coyote co-occurrence in non-shrubland habitats with intermediate resources, and only coyotes in resource-rich shrubland habitats. The detection of alternative stable states in natural systems is rare, and coexistence is predicted to be unlikely in mammalian systems (Verdy and Amarasekare 2010). Yet, all three theoretical stable states have been documented between kit foxes and coyotes across these two systems.
Characteristics common to mammalian IGP (e.g., asymmetrical IGP, IG predators not consuming IG prey; Palomares and Caro 1999) are expected to destabilize coexistence and shift systems toward either IG prey or IG predator dominated under traditional IGP theory Polis 1997, Verdy andAmarasekare 2010). Despite rarely consuming kit foxes (Palomares and Caro 1999), coyotes are the leading cause of kit fox mortality (i.e., 56%; Kozlowski et al. 2008) in our system. Still, equilibrium occupancy estimates were comparable to observed occupancy, suggesting coyotes and kit foxes have likely reached stable coexistence.
Stable coexistence may be facilitated by alternative resources and behavioral proclivities (Heithaus 2001, Holt and Huxel 2007, Kimbrell et al. 2007. Although the IGP literature often portrays shared basal resources as prey (e.g., Kimbrell et al. 2007), IG predators and IG prey may compete for other limited resources as well. Water is commonly a limited resource in deserts and may be acquired through prey (as preformed or metabolic water) or free-standing (e.g., springs) water. We did not detect a relationship between coyote occupancy and water availability, a result that aligned with the findings of other recent studies in our system (Hall et al. 2013, Kluever and. Still, where artificial water sources were present at Dugway, coyotes were detected using these sources 231 times more often than kit foxes (Hall et al. 2013); these water sources may have served as an alternative resource, reducing prey requirements that would have otherwise been necessary to meet water demands (Golightly and Ohmart 1984). Additionally, although dietary overlap between coyotes and kit foxes was high, coyotes were capable of exploiting larger mammalian prey (ungulates) that kit foxes were incapable of killing (Kozlowski et al. 2008). Our results suggested coyotes were more likely to occupy shrublands and woodlands, which maximized the availability of preferred prey (i.e., water-rich mammals) and provided thermal relief from climatic conditions. This pattern of resource matching by the IG predator to habitats supporting alternative resources was in accordance with IGP theory extensions and is predicted to support stable co-occurrence (Heithaus 2001, Holt andHuxel 2007).
Behavioral avoidance and increased vigilance may also encourage stable coexistence in mammalian IGP systems (Heithaus 2001, Kimbrell et al. 2007. Coexistence among intraguild predators often requires behavioral adjustments by subordinate species in their activity patterns or fine-scale space use. Vanak et al. (2013) found predators were aware of competitors at various spatial scales and subordinate species adjusted ❖ www.esajournals.org movements in the presence of IG predators. Similarly, habitat selection of Cape foxes (V. chama) did not differ in the presence of jackals at broad scales, but they had atypically large home ranges in the presence of jackals, presumably to facilitate behavioral avoidance of jackals during foraging (Kamler et al. 2013).
We assumed the number of coyote scats detected along each transect was reflective of coyote activity. Our results demonstrated coyote occupancy was sufficiently high that wherever kit foxes occurred, they co-occurred with coyotes. At broad scales (among sites), kit foxes appeared to avoid habitats with increased IGP risk (i.e., safety matched), aligning with predictions of behavioral extensions of IGP theory (Heithaus 2001) and empirical observations (Arjo et al. 2007, Kozlowski et al. 2012. At finer scales (within sites), our results indicated kit foxes were more likely to use areas with greater coyote activity, presumably balancing predation risk and food availability to secure sufficient resources for coexistence. Finescale space-use patterns were consistent with predictions of IGP theory extensions incorporating avoidance and vigilance, which tend to stabilize coexistence (Heithaus 2001, Kimbrell et al. 2007). These patterns aligned with empirical observations of Hall et al. (2013), which indicated a lack of spatial separation between coyotes and kit foxes, and could result from either (1) both species congregating in prey-rich areas or (2) coyotes hunting kit foxes. While we cannot formally distinguish between these two hypotheses, kit fox remains have not been detected in coyote scats from Dugway, despite significant overlap in dietary content and nightly activity patterns (Kozlowski et al. 2008). Collectively, these patterns supported the premises that killing of kit foxes was competitive in nature and that fine-scale overlap in space use between kit foxes and coyotes was likely driven by resource matching. Kimbrell et al. (2007) suggested competitive killing can help stabilize coexistence by balancing the competitive inferiority of the IG predator through increased vigilance and decreased foraging efficiency of the IG prey, and decreased vigilance and increased susceptibility of shared prey resources to the IG predator.
Investigations of canid IGP systems have focused primarily on static occupancy patterns, but elucidating drivers of local extinction and colonization can improve our understanding of the relative roles of habitat features and IGP on space use by subordinate species (MacKenzie et al. 2003). As predicted, the probability of local extinction by kit foxes was elevated across sites experiencing higher coyote activity. Contrary to predictions, colonization by kit foxes was positively associated with site-level coyote activity; the positive relationship between coyote activity and dynamic parameters (e and c) for kit foxes suggested that sites with greater coyote activity experienced higher turnover in kit fox occupancy (Fig. 4) and may indicate sites experiencing higher IGP. Intraguild predation theory predicts local extinctions of IG prey may be regulated by an IG predator, and that IGP effects will be more acute when the two have high dietary overlap (Polis et al. 1989, Holt andPolis 1997), as observed between kit foxes and coyotes (Kozlowski et al. 2008). When sympatric, coyote predation can account for a significant proportion of kit fox mortalities (up to 78%; Ralls and White 1995, Nelson et al. 2007, Kozlowski et al. 2008. Consequently, local extinction may result from a decreased ability to avoid IGP at sites with higher coyote activity, which once unoccupied by kit foxes, become sites available for colonization. Among available sites, those characterized by silty soils, which facilitate burrow excavation and therefore escape cover (Egoscue 1962, promoted colonization and had the highest kit fox occupancy probabilities (Fig. 2).

CONCLUSIONS
Our results identified two alternative states predicted by IGP theory (i.e., IG predator dominated and coexistence of IG predator and IG prey) in a single system. The stability of mammalian IGP systems is increased in the presence of alternative resources and behavioral responses, but the shift in dominance (e.g., between IG prey, coexistence, and IG predator) with changing productivity predicted by traditional IGP theory is often maintained Polis 1997, Holt andHuxel 2007). Equilibrium occupancy estimates suggested the distribution of kit foxes would increase greatly if coyotes were more spatially limited, but also that kit foxes and coyotes may have reached a stable state of coexistence. Observed stable coexistence may be enabled by a combination of kit foxes employing broadscale safety matching (e.g., behavioral avoidance of coyotes) and fine-scale resource matching (i.e., balancing predation risk and prey acquisition through increased vigilance). Our study elucidates the importance of considering dynamic processes and scale when investigating IGP systems. While we were unable to detect a direct influence of coyotes on kit foxes when considering occupancy, dynamic parameters suggested coyotes influenced the stability of kit fox occupancy. In particular, despite evidence of broad-scale avoidance of coyotes by kit foxes, elevated rates of local extinction (and turnover in occupancy) by kit foxes at sites with greater coyote activity suggested coyotes may still exclude kit foxes from some areas; thus, broad-scale patterns of space use may have resulted from competitive exclusion of the IG prey by the IG predator coupled with avoidance of the IG predator by the IG prey.

ACKNOWLEDGMENTS
Funding was provided by the National Geographic Society's Conservation Trust (C248-13) and the U.S. Army Research Laboratory and the U.S. Army Research Office through the U.S. Department of Defense's Environmental Security Technology Certification (12 EB-RC5-006) and Legacy Resource Management (W9132T-12-2-0050) Programs. The Dugway Natural Resource Management Program provided logistical support. The Utah Division of Wildlife Resources provided housing. We thank B Kluever, M Smith, T Edwards, M Melham, J Decotis, C Perkins, M Richmond, E Burke, and K Cobb for assistance with field work. J Adams provided laboratory guidance and Waits lab group technicians assisted with laboratory analyses. P Lukacs and J Rachlow reviewed earlier versions of this manuscript and provided helpful improvements. We thank G Roemer, J Cain, and one anonymous reviewer for insightful comments and suggestions that greatly improved this manuscript.