Introduction

The ecological interactions between sympatric species through competition (interference or exploitation) are key phenomena that contribute to shaping the structure of ecological communities, as they can influence the abundance, distribution, habitat selection, and behavior of species within communities (^{Case and Gilpin 1974}; ^{Holt and Polis 1997}; ^{Caro and Stoner 2003}; ^{Hunter and Caro 2008}).

Interference competition is widely documented for mammals of the order Carnivora, being considered among the main factors that shape intraguild relationships between predators (^{Polis et al. 1989}; ^{Palomares and Caro 1999}; ^{Linell and Strand 2000}; ^{Donadio and Buskirk 2006}; ^{Palomares et al. 2016}). In fact, this type of competition between carnivores is generally higher when the species involved are morphologically similar and share similar diets (^{Morin 1999}). The strategy of species to achieve coexistence consists of minimizing competition through niche segregation in one or several dimensions, mainly spatial, trophic, or temporal (^{MacArthur and Levins 1967}; ^{Pianka 1969}; ^{Pianka 1973}; ^{Schoener 1974}). Within this guild, the potential for competition between sympatric species that use similar resources is largely determined by the spatial overlap between them (^{Kitchen et al. 1999}; ^{Palomares and Caro 1999}; ^{Grassel et al. 2015}; ^{Palomares et al. 2016}). To minimize interference competition, subordinate species display a range of ecological strategies: avoidance of encounters with individuals of dominant species, separation of their home ranges, and differences in habitat use (^{Case and Gilpin 1974}; ^{Palomares and Caro 1999}; ^{Linell and Strand 2000}; ^{Hampton 2004}; ^{Rosenheim 2004}; Donadio and Buskirk 2006; ^{Berger and Gese 2007}; ^{Hunter and Caro 2008}; ^{Chiang et al. 2012}; ^{Viota et al. 2012}; ^{Soto and Palomares 2015}; ^{Xia et al. 2015}; ^{Gompper et al. 2016}; ^{Palomares et al. 2016}).

The quantification of the size and overlap of the home ranges of carnivores, as well as the description of habitat use and selection, are essential for understanding the dynamics of ecological communities, as well as for species conservation and management (^{Bu et al. 2016}). However, these complex interactions between sympatric species are generally poorly known in the vast majority of the systems where they thrive (^{Melville et al. 2015}; ^{Gompper et al. 2016}).

In North America, the coyote (Canis latrans) and the gray fox (Urocyon cinereoargenteus) are mesocarnivorous species that are sympatric over large portions of their distribution ranges (^{Bekoff 1977}; ^{Fritzell and Haroldson 1982}; ^{Fuller and Cypher 2004}; ^{Servin et al. 2014}a; ^{Servin and Chacón 2014}). The spatial interactions and the coexistence process between coyotes and various species of foxes in the Americas have been extensively studied in northern areas of their geographic range (United States of America and Canada). Research on spatial dynamics between coyotes and red foxes (Vulpes vulpes) has shown marked spatial segregation and differentiated use of the local habitat between these species (^{Voigt and Earle 1983}; ^{Sargeant et al. 1987}; ^{Theberge and Wedeles 1989}; ^{Harrison et al. 1989}; ^{Sargeant and Allen 1989}; ^{Gese et al. 1996}; ^{Gosselink et al. 2003}; ^{Mueller et al. 2018}). In turn, research on the spatial dimension between coyotes and kit foxes (Vulpes macrotis) has shown the absence of spatial segregation; instead, a differential habitat use has been observed (^{White et al. 1994}; ^{White et al. 1995}; ^{Nelson et al. 2007}; ^{Moehrenschlager et al. 2007}; ^{Kozlowski et al. 2008}; ^{Kozlowski et al. 2012}; ^{Andrade-Ponce et al. 2020}). Most information on spatial interactions between coyotes and gray foxes has been recorded in the United States of America, mainly in coastal shrubland and xeric shrubland areas at low altitudes (<1000 m asl). Some studies reported no spatial segregation between the two species (^{Neale and Sacks 2001}; ^{Chamberlain and Leopold 2005}), while others evidenced that gray foxes avoid spatial coexistence with coyotes to reduce the risk of predation (^{Fedriani et al. 2000}; ^{Farias et al. 2012}). This topic has been scarcely studied in areas within their distribution range in México, and the details about the spatial dynamics between these canid species in their natural distribution range in the country remain unknown.

For this reason, our objective was to evaluate the spatial ecological interactions between coyotes and gray foxes by analyzing the spatial segregation of the ecological niche under natural conditions in a temperate forest of the Sierra Madre Occidental, state of Durango, México. Our specific objectives were: 1) estimate the size and spatial overlap between the home ranges of both species and 2) evaluate habitat selection and use patterns to determine interspecific variations.

Our assumption was that the coexistence of these two species would be facilitated by the spatial segregation of their niches, exhibited by either a low or nil overlap of their home ranges or a pattern of differentiated habitat use. This is a case of an asymmetric interaction where coyotes display aggressive behavior against canids and other smaller species, which are displaced and even killed by coyotes, as reported for various fox species in North America (^{Sargeant and Allen 1989}; ^{Palomares and Caro 1999}; ^{Moehrenschlager and Sovada 2004}; ^{Moehrenschlager et al. 2007}). Thus, the gray fox (subordinate species) would be actively avoiding coyotes (dominant species) to reduce the risk of predation (^{Polis et al. 1989}; ^{Palomares and Caro 1999}; ^{Fedriani et al. 2000}; ^{Donadio and Buskirk 2006}; ^{Temple et al. 2010}; ^{Farias et al. 2012}).

Materials and Methods

Study area. This study was conducted in the buffer zone of “La Michilía” Biosphere Reserve (RBM), located in the municipality of Suchil, Durango, México, between coordinates 23° 21’to 23° 28’ N and -104° 09’ to -104° 21’ W (Figure 1). Physiographically, RBM is located in the transition zone between the Sierra Madre Occidental and the northern highlands of México (^{Halffter 1978}); besides, it covers part of the transition zone between the Nearctic and Neotropical biogeographic regions (^{Löwenberg-Neto 2014}; ^{Morrone 2014}; ^{Cuervo-Robayo et al. 2020}). Altitude in the study area ranges between 2,000 and 2,985 masl (^{Gadsden and Reyes-Castillo 1991}). To the north of the RBM, the climate is temperate and semi-dry (BS1k); in the rest of the zone, the dominant climate is temperate sub-humid (CW; ^{Garcia 2004}). The mean annual temperature is 12.6 °C, fluctuating between 2 °C (winter) and 22 °C (summer); the mean annual precipitation fluctuates between 600 and 900 mm (^{INEGI 2017}).

The main vegetation types are conifer forest (Pinus spp.) and oak forest (Quercus spp.); also present are natural grassland (Bouteloua spp.) and xeric shrubland (Arctostaphylos pungens, Acacia schaffneri). There are also transition zones between these types of vegetation, where the dominant species vary according to altitude, geomorphology, and microclimatic conditions, resulting in 22 different types of vegetation (^{González-Elizondo et al. 1993}; ^{Servín et al. 2014}b).

Figure 1 Geographic location of the study area in the buffer zone of La Michilía Biosphere Reserve (RBM), Durango, México. Home ranges of radio-collared coyotes (dotted red lines) and gray foxes (solid blue lines), derived from the minimum convex polygon (95 %), and the overlap between them, as well as the habitat types in the study area: Sv, disturbed vegetation (purple); QF, oak forest (light blue); MF, mixed forests (green); F-MS, forests with manzanita shrubland (pale pink); G, grassland areas (pale yellow). Areas in dark blue represent water bodies; the solid black line marks the border of the RBM and the dotted black line, the core zone of the RBM; gray lines are level curves (15 m). 

Capture and Marking. We used Tomahawk^® live traps and jaw traps (Victor^® Soft Catch No. 3) to capture five gray foxes (females) and five adult coyotes (two females, three males), respectively. The ten individuals captured were sedated by intramuscular injection with a mixture of xylazine (xylazine hydrochloride) and ketamine (ketamine hydrochloride). The composite dose to induce anesthesia was 4 mg/kg ketamine plus 2 mg/kg xylazine for coyotes and 3 mg/kg ketamine plus 20 mg/kg xylazine for gray foxes (^{Servin and Huxley 1992}; ^{Kreeger and Arnemo 2018}).

While individuals were sedated, we recorded morphometrics, weight, and sex; age (pup, juvenile, and adult) was determined based on tooth wear. In individuals with weight and measurements of adult animals, we fitted a 150 MHz VHF radio transmitter collar (Telonics^®), weighing 120 g (model 200) for gray foxes and 170 g (model 300) for coyotes. The net weight of these radio collars accounted for 1.49 % of the mean weight of the coyotes captured (W = 11,400 ± 1418 g) and 3.88 % of the mean weight of gray foxes captured (W = 3,094 ± 205 g). After the radio collar was fitted, each individual was released and at the capture site on the same day.

The handling and physical and chemical containment of individuals were performed according to the guidelines recommended by the American Society of Mammalogy (^{Sikes et al. 2016}), under the scientific research collection license number SGPA/DGVS/12685/18 granted to Jorge Servin, issued by the Ministry of Environment and Natural Resources of México.

Radiotracking and Location Error. We gathered radiotelemetry data between September 2017 and August 2019 (Table 1). We located individual animals at any time of the day or night using portable receivers (Telonics^® Mod. TR-2) with “H”-type handheld antennas and fixed eight-element antennas known as zero-point systems (Wildlife Materials Inc.®). Animals fitted with radio collars were field-tracked using the “triangulation” method (^{Mech 1983}). This method consists of determining the location involving at least two directions (bearings or azimuths) using a compass from two different sites of known location separated from one another by at least one kilometer. A straight line was projected from each site to the bearings obtained so that the site where these lines crossed marked the location of the animal at that time. For the laboratory analysis of these measurements, we considered only those pairs of readings that were taken within 5 minutes and with a difference greater than 20° and less than 160°. To note, readings with differences less than 20° or greater than 160° produce triangles with very sharp vertices, which significantly increase location errors (^{White and Garrot 1990}). Prior to the start of the monitoring period for radio-fitted animals, we estimated the location error using reference transmitters placed at known sites, yielding an error of ± 3° (^{White and Garrot 1990}).

Home Range and Overlap. Using the location data recorded in the field, we constructed an Excel^® database, which was loaded into the LOAS^® program Location of a Signal, version 4.0.3.8 (^{ESS 2010a}); this returned a cloud of points in space and a database containing the georeferences of the locations of each radio-collared individual. With this database, we used the program Biotas^® version 2.0a 3.8 (^{ESS 2010b}) to calculate the size of the home range of each radio-collared individual, using the minimum convex polygon method set at 95 % (MPC; ^{Mech 1983}; ^{White and Garrot 1990}), while 50 % of sites were used to determine the core zone (i. e., the area with a high priority of use; ^{Powell 2000}). We used the MPC for its simplicity (^{White and Garrot 1990}) and to compare our results versus other studies addressing the species studied. To estimate the space shared between radio-collared animals, we measured the overlap of home ranges between pairs of individuals and then calculated the average of this overlap (^{Millspaugh and Marzluff 2001}).

We compared the size of the home ranges between the two species through a Student’s t-test for independent samples; in the case of coyotes, we compared the size of the home ranges between sexes through a Student’s t-test for a single sample (^{Sokal and Rohlf 1987}).

Habitat Use and Selection. We used a vegetation map of the RBM and its area of influence (1:50,000 scale) for the classification and assignment of habitat types according to the physiognomically dominant vegetation (sensu^{González-Elizondo et al. 1993}), which was digitized by the Laboratory of Wildlife Ecology and Conservation at Universidad Autónoma Metropolitana, campus Xochimilco. This map grouped habitat types into five categories (Figure 1): disturbed vegetation and agricultural areas (Sv); oak forest (QF), dominated by Quercus spp.; mixed forests (MF), with Pinus and Quercus as dominant or subordinate species; forests (pine, oak, or mixed) including patches of A. pungens shrubland (F-MS); and grasslands (G), areas where Bouteloua spp. occur as dominant or subordinate species.

The locations of the radio-collared individuals of both species were superimposed to the resulting map to quantify the frequency with which each individual was located in each habitat type within the RBM. We calculated the habitat selection coefficient by species, individual, and habitat type, as the ratio between observed habitat use and habitat availability (^{Manly et al. 2004}). The observed habitat use was determined from the radio-location points recorded for each individual by habitat type. Habitat availability was derived by multiplying the number of radio-location points of each individual in a particular habitat type by the observed proportion of that habitat type within its home range obtained through the MPC (^{Aebischer et al. 1993}; ^{Sankar et al. 2013}). This comparison is analog to Johnson’s third-order selection (^{Johnson 1980}). For each species, we calculated the habitat selection coefficient for the j-th individual and the i-th habitat type using the equation: ŵ _ij = u _ij / (π _i u _+j ), where ŵ _ij is the selection coefficient of individual j in habitat i; u _ij is the number of radiolocation points of individual j in habitat i; π _i is the relative availability of habitat i; and u _+j is the total number of individual radio location points of individual j (^{Manly et al. 2004}). We calculated a measure of the selection made by individuals of a given species as a group (taking into account the variation in the selection of habitats of each individual) with the following equation: ŵ _i = u _i+ / (π _i u ₊₊ ), where ŵ _i is the selection coefficient for habitat i; u _i+ is the total number of radiolocation points in habitat i; and u ₊₊ is the total number of radiolocation points for all individuals (^{Manly et al. 2004}). Under the assumption that individual j uses habitat type i randomly, the average value of the habitat selection coefficient is ŵ = 1 (use according to availability); thus, coefficients with values ŵ > 1 indicate a higher-than-expected use (i.e., preference), while ŵ values < 1 indicate lower-than-expected use (i. e., avoidance; ^{Manly et al. 2004}). To determine whether a value of habitat selection coefficient (ŵ _i ) was significantly different from 1, we generated and used the 95 % Bonferroni confidence intervals (sensu^{Manly et al. 2004}). We used a G-test or two-step log-likelihood ratio to test the null hypothesis that habitat was used according to habitat availability (^{Sokal and Rohlf 1987}). First, we performed the G-test for each individual; afterward, we added the values of the test statistics for all individuals of a given species to test the overall habitat selection of individuals (^{White and Garrot 1990}; ^{Manly et al. 2004}).

We calculated the habitat selection coefficient with the adehabitatHS package (^{Calenge 2006}) for R version 4.0.1 (^{R Core Team 2019}). All statistical analyses were performed with this software, considering a significance level α = 0.05. For those parameters that require so, we report the mean ± standard deviation.

Results

Although we captured and tracked ten individuals - five gray foxes and five coyotes -, a local inhabitant delivered one radio collar that we had fitted to a female coyote captured four weeks earlier, reporting that the collared coyote was found dead by gunshot. Therefore, below we report the data corresponding to nine individuals.

Home Range and Overlap. Between September 2017 and August 2019, we recorded a total of 945 radio location points, 429 corresponding to four coyotes (= 107.25 ± 17) and 516 to five gray foxes (= 103.20 ± 47; Table 1).

The average home range for coyotes was 12.20 ± 1.74 km² (n = 4; range 9.74-13.81 km²), with a mean core zone of 2.97 ± 1.03 km² (n = 4; range 1.98-4.05 km²); the home range of male coyotes (13.02 ± 0.70 km²; n = 3) was significantly larger (t = 8.07; d. f. = 2; P = 0.01) than the home range of the only female monitored (9.74 km²). The mean home range size for gray foxes was 5.37 ± 0.67 km² (n = 5; range 4.40-6.09 km²), with a mean core area of 1.01 ± 0.52 km² (n= 5; range 0.41-1.63 km²). A t-test showed that the mean home range size of coyote was significantly larger versus gray fox (t = 8.18; d. f. = 7; P = 0.001).

The mean overlap of home ranges between coyotes (intraspecific overlap) was 43.7 ± 21 % (n = 12; range 18-77 %), whereas for gray foxes, the mean overlap was 6.6 ± 5 % (n = 8; range 1-14 %). The overlap of home ranges was significantly greater between coyotes than between gray foxes (t = 4.87; d. f. = 18; P = 0.001). The mean overlap of home ranges between coyotes and gray foxes (interspecific overlap) was 42.1 ± 27 % (n = 21; range 13-98 %; Figure 1).

Habitat Use and Selection. We found that coyotes use the different habitat types according to their availability, both as a group (G = 18.36; d. f. = 13, P = 0.14) and as individuals (P > 0.05; Table 2), although the highest habitat selection coefficient was obtained for grassland areas (G; ŵ _i = 1.17) and the lowest for forests with manzanita shrubland (Arctostaphilos pungens; F-MS; ŵ _i = 0.77; Table 2).

On the other hand, although individual variations were observed, gray foxes showed selective habitat use as a group (G = 113.93; d. f. = 14; P < 0.001; Table 2). The gray fox preferred mixed forests (MF; ŵ _i = 1.43 ± 0.17) and avoided disturbed vegetation (Sv; ŵ _i = 0.17 ± 0.0), grassland areas (G; ŵ _i = 0.27 ± 0.07), and forests with manzanita schrubland (F-MS; ŵ _i = 0.61 ± 0.11). Separately, selection coefficient values and their confidence intervals indicated that oak forest (QF; ŵ _i = 0.96; CI_B 0.17-1.74) was used according to its availability (Table 2).

Discussion

In a previous work carried out in the study area, ^{Servin (2000}) radio-tracked fifteen coyotes (eight males and seven females) over two years, reporting a mean home range size of C = 11.8 ± 2.71 km² for coyotes in general, C_M = 13.1 ± 2.5 km² for males, and C_H = 9.9 ± 3.3 km² for females. These values are similar to the ones obtained in the present study. The home range size of coyotes is a highly dynamic variable influenced by climate, prey availability, and habitats suitable for reproduction, as well as by population density and mortality rate (^{Danner and Smith 1980}; ^{Laundré and Keller 1984}; ^{Gese et al. 1988}; ^{Servín and Huxley 1995}; ^{Servín et al. 2014b}). While the home range size of a species varies geographically (^{Holzman et al. 1992}; ^{Chamberlain et al. 2000}), our results indicate that home range size in the study area lies within the range of values reported for coyotes in different habitats across its range (^{Bekoff 1977}; ^{Andelt and Gipson 1979}; ^{Young et al. 2006}), consistent with most of the studies carried out in temperate zones of North America (11.6-35.8 km²; ^{Servín and Huxley 1995}; ^{Servín 2000}).

In the case of gray fox females, the mean home range size reported here was 5.37 ± 0.67 km², an area 2.4 times larger than the one reported for females by ^{Servin et al. (2014}b), which was 2.24 km², in the same study area between 1991-1993. In this regard, some studies have reported that gray fox females tend to display a larger home range than males (^{Trapp and Hallberg 1975}; ^{Servín et al. 2014b}) and that the home range size of this species may vary depending on habitat quality and resource availability (^{Fuller and Cypher 2004}). Our results fall within the range of variation reported elsewhere for this species (^{Fritzell and Haroldson 1982}; ^{Fuller and Cypher 2004}; ^{Macdonald and Sillero-Zubiri 2004}).

In the study area, the average spatial overlap between the home ranges of coyotes and gray foxes was moderate (42.1 %). However, one-third of the interspecific pairs (diads) analyzed (n = 21) to derive this data showed high overlap values (> 60 %), as reported in other studies (^{Neale and Sacks 2001}; ^{Chamberlain and Leopold 2005}). Our results suggested that, since there are no apparent patterns of spatial avoidance of the gray fox toward the coyote through spatial segregation of the ecological niche, the spatial dynamics between these species is not fully explained by interference competition, as reported for these canid species in other areas where they display a sympatric distribution (^{Fedriani et al. 2000}; ^{Farias et al. 2012}).

Our results also suggest that the spatial coexistence dynamics between coyotes and gray foxes in the study area is governed by space-use mechanisms at a fine scale (^{Lonsinger et al. 2017}) mediated by differential habitat use. The gray fox used oak forest (QF) according to its availability and showed preferences for mixed forests (MF), as already reported for this species in the study area (^{Servín et al. 2014}b), as well as in other areas over its geographic range (^{Haroldson and Fritzell 1984}; ^{Chamberlain and Leopold 2000}). These forests offer vast areas that provide protection and shelter for gray foxes (^{Servín et al. 2014b}), being an important element within the home range of this species (^{Fritzell and Haroldson 1982}). The complex architecture of mixed oak-pine and pine-oak forests in the study area provide natural structures that can be used as resting sites and shelters; at the same time, these forests serve as escape routes and, therefore, are useful to avoid the risk of predation, as gray foxes are able to climb trees and even jump between tree branches (^{Fritzell and Haroldson 1982}; ^{Fuller and Cypher 2004}). In addition, foxes can use the tree stratum as a foraging zone, as its branches are habitats for potential prey that are part of their diet, such as passerine birds, squirrels (Sciurus nayaritensis, Tamias bulleri, and T. durangae), small rodents (Peromyscus spp. and Reithrodontomys spp.), lacertids (Sceloporus spp.). and insects. On the other hand, gray foxes avoided disturbed vegetation (Sv) and grassland areas (G). It has been shown that the risk of predation by coyotes can influence resource use by gray foxes (^{Fedriani et al. 2000}; ^{Chamberlain and Leopold 2005}). Thus, foxes are likely to be avoiding these open areas as these are devoid of shelters, hence offering lower evasion opportunities against the potential chase by coyotes (which use these habitats according to their availability) to avoid intraguild predation (^{Temple et al. 2010}).

In the case of coyotes, although no apparent preference for or avoidance of any particular habitat type was observed, a certain trend towards the preferential use of pasture areas (G) was noted since it attained the highest habitat selection coefficient (ŵ _i = 1.17). This trend of preferential use is consistent with data reported for coyotes in the RBM, as this species forage and catch their main prey (rodents and lagomorphs) preferentially in areas with open vegetation, such as grasslands (P; ^{Servin and Huxley 1991}; ^{Servin et al. 2003}). These open vegetation areas in the RBM are also home to the checker bark juniper or táscate (Juniperus deppeana) with varying abundances in different areas (^{González-Elizondo et al. 1993}). Juniper fruits are an important element in the diet of coyotes, being the plant food most frequently consumed by coyotes in the study area (^{Delibes et al. 1989}; ^{Servin and Huxley 1991}); this food category is actively sought and consumed by coyotes in open and grassland areas.

An aspect worth highlighting is the role of forests with manzanita shrubland (F-MS) in habitat selection and use by both species. On the one hand, coyotes used this habitat as expected, while gray foxes avoided it. One potential explanation lies in the different frequency of consumption of manzanita fruit by both species. These fruits represent a food resource highly consumed by coyotes, especially in the dry season (February-May; ^{Servin and Huxley 1991}), while it is consumed to a lesser extent by gray foxes (^{Delibes et al. 1989}).

In the present study, we showed that the home ranges of coyotes and female gray foxes showed a moderate interspecific overlap, so no spatial segregation occurred. However, differential use of habitat was observed, which explains the coexistence of these canids in the same area because their antagonistic behavioral interactions decrease through a trend towards the differential use of resources (^{MacArthur and Levins 1967}; ^{Tilman 1982}; ^{Holt 2001}). Our results are consistent with the theoretical hypothesis on intraguild predation ^{(Holt and Polis 1997}; ^{Polis et al. 1989}), which suggests that the coexistence between species in the same guild sharing basic resources requires that the subordinated species (gray fox) be better at exploiting the resources shared with the dominant species (coyote).