Mesocarnivores in residential yards: influence of yard features on the occupancy, relative abundance, and overlap of coyotes, grey fox, and red fox

Study area

Our study was conducted in the north-western region of Arkansas, USA, which is a rapidly developing area with a population of approximately 350,000 people. Fayetteville is located within the Ozark Highlands Ecoregion, which is primarily forest intermixed with open areas and some pasture lands used for cattle and agriculture. Our study sites encompassed residential yards in and near the centres of Fayetteville, Springdale, Rogers, Eureka Springs, and Bella Vista and less developed areas up to 32 km from city centres (Fig. 1). The yards in which we conducted our study were owned by volunteers from the Arkansas Master Naturalist Program and University of Arkansas Department of Biological sciences. We restricted our study sites to yards occurring within an 80.5 km radius of downtown Fayetteville. We intentionally selected yards that represented a continuum of urban to rural setting as determined by surrounding housing-unit density (described below). We also ensured that we selected yards that contained (or did not contain) varying suites of a priori selected yard features associated with supplemental food, shelter, or water for wildlife (see below). Because volunteers came from the Master Naturalist Program and a Biology Department, not all yards contained features to attract wildlife and were representative of the yards available to wildlife in this region.

Fig. 1.

Map of study sites (residential yards) in north-western Arkansas, USA, during the summers of 2021 and 2022; study sites are indicated with a blue pin. (a) Large-scale view of the USA. (b) Finer-scale view of Arkansas, USA. (c) Finer-scale view of study sites within north-western region of Arkansas, USA.

Camera deployment

We deployed motion-triggered camera-traps (Browning StrikeForce or Spypoint ForceDark) in residential yards during the summers of 2021 and 2022. We deployed one camera per yard approximately 0.95 m above the ground on either a tripod or strapped to a tree. We placed cameras at least 5 m away from one another and no more than 100 m from a house. We positioned cameras to maximise wildlife detections, when possible, near features such as compost piles, water sources (natural or maintained), and fence lines or animal paths. We coordinated with homeowners to choose locations that would allow them to maintain their privacy and that also would not interfere with yard maintenance. If vegetation grew to impede the view of the camera, we would remove the vegetation. We set cameras to take a burst of three photos per trigger, with a 5 s reset time when triggered. We did not use any bait or lures at camera sites. Cameras were checked every 2 weeks to assess battery life and to download data. We moved cameras within yards up to three times within the season, to ensure that we captured the full range of wildlife present in each yard.

We surveyed each yard for features that we predicted would influence mesocarnivore occupancy (Table 1). We first recorded variables representing possible food sources, including the presence and number of bird feeders (both seed/suet and nectar), the volume of brush and firewood piles, area of compost piles, and the presence or absence of poultry. All area and volume calculations were undertaken by using a measuring tape to measure height, length, and width of the features. We reasoned that bird feeders could provide direct food to some of our focal species and could also serve as an attractant to prey species such as mice, rats, and birds (Saad et al. 2020). We considered the volume of firewood and brush piles as a food source because they could similarly be associated with high numbers of prey species (Goguen et al. 2015).

We also recorded the area of potential sheltering sites in a yard. Sheltering sites included the total available area of sheds and outbuildings under which animals could access as well as decking that was less than 0.3 m off the ground (Linduska 1947; Moll et al. 2018). We counted the number of supplemental water sources present in a yard, which included bird baths, garden ponds, and fountains that were actively maintained by the homeowner.

If a camera was deployed in an area that was surrounded by a fence, we scored the type of fence based on its ease of access to our focal species. We scored each yard with a binary score of zero or one. Yards without fences or with fences that presented minimal resistance to the movement of our focal species (barbed wire or those made of semi spaced slats 1 m tall that allowed for smaller-bodied species to pass through and most larger species to jump over) were given a zero. If a fence was fallen, had large holes in them, or had a gate that was always left open, the fence also received a score of zero. Privacy and chain-link fences that were 1.8 m or taller and were made of a solid material, that were likely to prevent the passage of focal species, apart from capable climbers, were given a one.

Landscape and environmental covariates

To acquire landscape-level covariates to explain mesocarnivore occupancy and relative abundance, we calculated the area of several landcover types surrounding each focal yard. First, we imported all camera-location points to GIS (ArcGIS Pro 10.2; ESRI, Inc. Redlands Inc.) and created 1.5 km buffers around each camera on the basis of a reasonable approximation of mesocarnivore home ranges (Šálek et al. 2015). Within each buffer, we then calculated the area of forest cover, developed open land (e.g. cemeteries, parks, and grass lawns), and agricultural land (grasslands, cultivated crops, and pastureland) by using the 2019 National Land Cover Database (Dewitz and US Geological Survey 2021). We used two different development categories for analyses because we predicted that our focal species might respond differently to development of varying intensities. We used low-intensity development, which was represented by areas containing 20–49% of impervious surface. We then grouped the National Land Cover Database (NLCD) categories of medium- and high-intensity development (which we refer to as high intensity development) that contained areas with over 49% impervious surface. We also quantified the maximum housing-unit density (HUD) around each camera by using the SILVIS Housing Data Layer (Radeloff et al. 2018).

Mesocarnivore occupancy

To explore the effects of covariates on the detection and occupancy of our three focal mesocarnivores (coyotes, red and grey foxes), we used a single-season multi-species occupancy model (MacKenzie et al. 2002). Our occupancy models were used to estimate the probability that a mesocarnivore occurred in a yard and the influence of covariates on occupancy probability while accounting for imperfect detection (MacKenzie et al. 2004). We chose multispecies models to account for avoidance or attraction between our focal species. Multispecies occupancy allowed us to simultaneously model environmental covariates, while letting the occupancy of one species vary on the basis of the presence of another. We used individual yards as our survey locations and within each survey location, we used 1-week sampling periods. For each sampling period, we assigned a one or a zero for each species at each survey location. If the focal species was detected during a survey period at a site, a one was assigned, and if none was detected, a zero was assigned. We used a 1-week sampling period because this represented an appropriate amount of time as to not over- or under-compress statistical power and is consistent with numerous occupancy analyses conducted using camera-traps (Trolle and Kéry 2003; MacDougall and Sander 2022). During the summer of 2022, we resampled 41 of our study yards from 2021. For these sites, we created year–site combinations and treated them as independent sites, as in previous studies (Hines et al. 2014; Devenish-Nelson and Nelson 2021; Murray et al. 2021). We used year as an occupancy covariate (Linden and Roloff 2013) to explore whether patterns differed between the 2 years of the study. Not all cameras in yards were actively functioning for the duration of the entire season because of camera malfunctions or staggered deployment and pick-up dates. To correct for this, we censored weeks in which the camera was not operating continuously for all 7 days. We censored these camera weeks by assigning a ‘NA’ value as opposed to a one or a zero.

We also compiled a number of environmental covariates that we predicted would influence the behaviour or activity of our three focal species. Both rainfall and temperature can influence how active our species were and therefore could affect detection probability (Madsen et al. 2020). We defined detection probability as the chance that a species would be detected on a given day. To gather environmental variables, we used publicly available data from a National Oceanic and Atmospheric Association (NOAA) weather station located at the Fayetteville Experimental Station. This weather station was 75 km from the furthest yard studied. We used this data to calculate the average weekly air temperature and average weekly precipitation for each camera site.

Before model fitting, we assessed collinearity of covariates and considered any two covariates to be collinear if they had correlation coefficients ≥|0.6| (Gilhooly et al. 2019). From those, we would then decide which of the two variables were predicted to be more meaningful and only included that variable in subsequent analyses (Gilhooly et al. 2019). For detection covariates in occupancy analyses, we found that temperature was correlated with both precipitation (r = −0.61) and week of the survey (r = 0.72) and was therefore subsequently removed. When assessing occupancy covariates, high-intensity development was correlated with area of low intensity development (r = 0.75) as well as forested area (r = −0.74). Area of low development was also found to be correlated with forested areas (r = 0.73). Both development levels were subsequently removed from analyses. So as to facilitate comparison between variables measured on different scales, we centered all continuous covariates by subtracting the variable mean value and then dividing by one standard deviation (Schielzeth 2010).

To avoid model over-fitting, we used a multi-stage fitting approach (Fuller et al. 2016) to select for the best detection covariate. We modelled mean precipitation, week of the year, an additive model of week + precipitation, and a null model as covariates for detection against null occupancy parameters and selected the top covariate model using Akaike’s information criterion corrected for small sample size (AICc). The covariate(s) in the top model were then used as the detection covariates in all subsequent analyses of occupancy probability.

For occupancy covariates, we included the area of compost pile, supplemental water sources, poultry presence, number of bird feeders, fence type, as well as the volumes of denning sites, firewood piles, and brush piles (Table 1). We also included four landscape-level occupancy covariates, including the areas of open development, agriculture, forest, and HUD.

Our candidate model set included all possible two-way additive combinations of occupancy covariates with our top identified detection covariate(s) as well as a global and a null model. All model fitting was performed in R (R Core Team 2022, available at https://www.r-project.org/) with the unmarked package (Fiske and Chandler 2011, available at https://cran.r-project.org/web/packages/unmarked/index.html). To improve clarity in presenting model selection tables, we display only models that were competitive within two ΔAICc. Because there was model uncertainty, we used full model averaging (Lukacs et al. 2009) for all the model estimates to generate an average of unconditional occupancy estimates (Cade 2015).

Mesocarnivore relative abundance

In addition to mesocarnivore occupancy, we also examined how yard and landcover covariates influenced the relative abundance of each focal species. We defined relative abundance as the number of focal species detections at each site (yard) divided by the total number of trap-nights. Although this metric often correlates with true abundance (Campbell et al. 2015), it is more often used as an index of intensity of use within an area and thus could be a good complement to occupancy, which is a coarser metric (Martin et al. 2010). Because camera-traps take multiple photographs at once, we combined all photos taken within a 5 min period into one ‘detection’ to minimise double-counting of individuals. To evaluate which landcover and yard variables most influenced the relative abundance of our mesocarnivore species in each yard, we used a generalised linear model (GLM) analysis. We used separate GLM analyses for each species to explore how landscape and yard features affected how frequently a mesocarnivore was detected in a yard. For each analysis, we evaluated the same candidate models as those used in the occupancy analyses. The candidate model set for each analysis consisted of simple one-way variable models and all additive two-way combinations as well as a global and a null model.

For each analysis, we ranked candidate models by using an information-theoretic approach with AICc. When appropriate, we derived parameter estimates for candidate models by model averaging all models within two ΔAICc (Burnham and Anderson 2002) in R (R Core Team 2022) with the AICcmodavg package (R package version 2.3.2, available at https://cran.rproject.org/package=AICcmodavg).

To improve clarity in presenting model selection tables, we display only models that were competitive within two ΔAICc for each analysis. Initial exploratory analyses indicated that relationships between predictor variables and response variables were linear and thus models were not corrected. Full model averaging was used to avoid bias away from zero (Lukacs et al. 2009). Model goodness-of-fit was assessed using residual plots and examining the R² value (R² = 0.12).

Interspecific temporal overlap

To better understand interspecific temporal overlap, we calculated the activity overlap between each species pair by using the R package overlap (Meredith and Ridout 2021; available at https://cran.r-project.org/web/packages/overlap/index.html). We created kernels of activity for each species to quantify when during the diel period they were most active. We first converted the time that each species was detected to radians. We then calculated a species-specific non-parametric kernel around the detection times to quantify when each species was active. For each species pair, we calculated the overlapped area under the kernel curve for the species pair. This number was quantified along a 0–1 index, with zero representing no overlap and one representing complete overlap (Schmid and Schmidt 2006). For our purposes, we used the Δ4 value (a vector of densities estimated at the time of observation) because we had over 75 detections of each species (Schmid and Schmidt 2006; Ridout and Linkie 2009). We used detections from across all of our study sites collectively to create activity kernels.

To assess the time between interspecific detections, we calculated the time elapsed between detections of one species in a yard and the next chronological detection of the next focal species in the same yard. We first did this by days between, and then for anything with less than 24 h between detections, we calculated hours between the detections.

Mesocarnivore occupancy

Cumulatively, we conducted 1456 week-long surveys for a total of 10,192 trap-nights and accumulated 1,494 focal mesocarnivore detections (507 coyote, 157 grey fox, and 830 red fox). We detected at least one mesocarnivore (coyote, grey fox, or red fox) in 74% of yards. Two species were present at 25 yards, whereas only six yards had all three species present.

Coyotes were detected at 51% of yards during 191 survey periods, grey fox at 14% of yards in 52 surveys, and red fox at 40% of yards during 160 surveys. The occupancy probability (ψ) for at least one mesocarnivore in a yard was 0.74 (95% credibility interval (CI): 0.66–0.80). The ψ of coyotes was 0.54 (95% CI: 0.43–0.63), ψ of grey fox 0.14 (95% CI: 0.09–0.26), and ψ of red fox was 0.41 (95% CI: 0.31–0.53). The detection probability (p_i) was similar among the three species, with p_i of coyotes averaging 0.18 (95% CI: 0.16–0.22), p_i of grey fox averaging 0.19 (95% CI: 0.14–0.25), and p_i of red fox averaging 0.20(95% CI: 0.17–0.24). Preliminary modelling for detection covariates suggested that the best predictors for mesocarnivore detection at a yard were the additive effect of week and average precipitation, and this model accounted for 63% of the weight of evidence (Table 2). Coyotes and red fox were less likely to be detected later in the summer and coyotes and grey fox were more likely to be detected during rainy weeks. Therefore, we used week and average precipitation as the p_i covariates for all analyses of ψ.

The three candidate models predicting mesocarnivore occupancy that were within 2 ΔAICc units of one another all assumed pairwise dependence between species. These top models accounted for 36% of the weight of evidence (Table 3). Two of the top models included fence type, one of which also included compost piles in the model; however, credibility intervals overlapped for all species for effect of compost pile. The other model that received some weight of evidence was the area of shelter sites (Table 3 and Supplementary material Table S1). Fences reduced coyote ψ (β = −0.92, 95% CI: −0.93, −0.07) but coyote ψ increased when sheltering area was abundant (β = 0.01, 95% CI:0.00, 0.01) (Fig. 2). We found no evidence that any of the measured covariates influenced ψ of red or grey fox (Table 4).

Fig. 2.

95% credibility intervals for covariates from our top a priori models evaluating detection (p_i) and occupancy probability (ψ) of coyote (Canis latrans), grey fox (Urocyon cinereoargenteus), and red fox (Vulpes vulpes) occurring in residential yards during the summers of 2021 and 2022 in north-western Arkansas, USA.

We found that Ψ of all our study species increased if another species was also present in that yard (Table 5). This relationship was strongest for coyote, given that a red fox was present, with grey fox having the weakest influence from coyote presence.

Relative abundance

We found little evidence for the effects of covariates influencing the relative abundance of coyote or red fox in a yard and the null model was the top-ranked model for both species. We found two top models for grey fox relative abundance within two ΔAICc that accounted for 42% of the weight of evidence (Table 6). Both top models included poultry presence and the top model also included forested area (Table 6). These models suggest that relative abundance of grey fox was highest in yards with a chicken coop (model-averaged β = 0.06, 95% CI: 0.01–0.3) (Fig. 3).

Fig. 3.

Influence of poultry presence on relative abundance of grey fox (Urocyon cinereoargenteus) in residential yards on the basis of a 2-year occupancy study conducted in 2021 and 2022 in north-western Arkansas, USA. Values are means with 95% credibility intervals.

Species temporal overlap

When assessing the activity times of our focal species, we found that they were all predominantly nocturnal. All species had peak activity levels approximately at midnight. Coyotes and grey fox had an 87.8% temporal overlap (95% CI: 79.0–95.3) (Fig. 4) and coyotes and red fox had an 89.6% temporal overlap (95% CI: 84.7–94.4) (Fig. 4). Grey and red fox had a temporal overlap of 78.1% (95% CI: 68.5–87.0) (Fig. 4).

Fig. 4.

Kernel temporal overlap between species across all sites within a 2-year study conducted in 2021 and 2022 in residential yards in north-western Arkansas, USA. (a) Grey fox (Urocyon cinereoargenteus) and coyote (Canis latrans). (b) Red fox (Vulpes vulpes) and coyote (Canis latrans). (c) Grey fox (Urocyon cinereoargenteus) and red fox (Vulpes vulpes).

Most of our focal mesocarnivores were not detected during the same night in the same yards. For example, grey fox and coyote detections were, on average, 4 (±5.5 s.d.) days apart (min: 0; max: 22). Coyote and red fox had the longest separation between detections, with an average of 12 (±19) days (min: 0; max: 88). Grey fox and red fox were separated by an average of 7 days (±4.4) days (min: 0; max: 33). If the species were detected within the same 24-h period, we calculated the average time between detections and found that coyote and grey fox had an average of 11 (±9.7) h between detections (min: 0.4; max: 23.5) at nine sites where they were both detected within the same 24-h period. Coyotes and red fox had an average of 7.6 (±7.6) h between detections at 15 sites (min: 0.96; max: 22). Red and grey fox were detected within the same 24-h period at only two sites and were detected an average of 14 (±5.9) h apart (min: 10.14; max: 18.5).