Responses of sympatric canids to human development revealed through citizen science

Abstract Measuring wildlife responses to anthropogenic activities often requires long‐term, large‐scale datasets that are difficult to collect. This is particularly true for rare or cryptic species, which includes many mammalian carnivores. Citizen science, in which members of the public participate in scientific work, can facilitate collection of large datasets while increasing public awareness of wildlife research and conservation. Hunters provide unique benefits for citizen science given their knowledge and interest in outdoor activities. We examined how anthropogenic changes to land cover impacted relative abundance of two sympatric canids, coyote (Canis latrans), and red fox (Vulpes vulpes) at a large spatial scale. In order to assess how land cover affected canids at this scale, we used citizen science data from bow hunter sighting logs collected throughout New York State, USA, during 2004–2017. We found that the two species had contrasting responses to development, with red foxes positively correlated and coyotes negatively correlated with the percentage of low‐density development. Red foxes also responded positively to agriculture, but less so when agricultural habitat was fragmented. Agriculture provides food and denning resources for red foxes, whereas coyotes may select forested areas for denning. Though coyotes and red foxes compete in areas of sympatry, we did not find a relationship between species abundance, likely a consequence of the coarse spatial resolution used. Red foxes may be able to coexist with coyotes by altering their diets and habitat use, or by maintaining territories in small areas between coyote territories. Our study shows the value of citizen science, and particularly hunters, in collection of long‐term data across large areas (i.e., the entire state of New York) that otherwise would unlikely be obtained.


| INTRODUC TI ON
With a growing human population, it is increasingly important to determine how anthropogenic activities can influence wildlife populations and communities (McKinley et al., 2017;Nichols & Williams, 2006). However, detecting anthropogenic impacts is difficult in many systems due to a lack of appropriate data (Magurran et al., 2010). Long-term, large-scale data are particularly useful because it allows for the documentation of baselines and detection of broad-scale changes in ecological patterns that may be missed using sparse or small-scale data sets (Lindenmayer & Likens, 2010). Despite the obvious benefits, long-term studies remain relatively uncommon due to the costs and logistics required, especially for wide-ranging and cryptic species (Parsons, Goforth, Costello, & Kays, 2018).
One solution to overcome this problem is through citizen science, in which members of the public participate in scientific work (Bonney et al., 2009). Citizen science provides a means to achieve substantial data collection and data processing needs while also contributing to the awareness, knowledge, and understanding of the non-science public about the practice of science and the relevance of scientific outcomes (Parrish et al., 2019). Particularly in wildlife biology, it can increase the spatial or temporal scale of research projects while increasing the public's interest and knowledge toward wildlife and the environment (Frigerio et al., 2018).
Despite issues such as incomplete reporting and skewed sampling (Tye, McCleery, Fletcher, Greene, & Butryn, 2017), citizen science data can yield robust inferences when analyzed rigorously (e.g., Davies, Stevens, Meekan, Struve, & Rowcliffe, 2012;Kéry, Gardner, & Monnerat, 2010). Standardized monitoring schemes result in the highest data quality; however, opportunistic data collection is the next best option for long-term studies, overcoming uneven sampling with large amounts of data (Isaac, van Strien, August, de Zeeuw, & Roy, 2014). Studies based on field work are particularly challenging, due to potentially risky field conditions for researchers or the need for frequent sampling (McKinley et al., 2017). However, some subsets of the public (e.g., birders, hunters, hikers) are well-suited for this type of data collection because they regularly visit outdoor areas and have both the physical skills and the interest to participate in outdoor activities. Hunter indices, for example, have been used for decades to monitor wildlife population trends in the United States (e.g., Haskell, 2011) and Europe (e.g., Ericsson & Wallin, 1999).
North America has the highest abundance of hunters of any region in the world (Sharp & Wollscheid, 2009), providing a unique opportunity to use hunter-collected data to investigate broad ecological questions. The widespread occurrence of hunting facilitates use of such data to study species such as carnivores, which are generally difficult to detect given their low densities and secretive nature (Rich et al., 2013). While the geographic ranges of many North American carnivores have contracted over the past century, ranges of red foxes (Vulpes vulpes) and coyotes (Canis latrans) have increased from their historic geographic ranges by 16% and 40%, respectively (Laliberte & Ripple, 2004). Both species have a high degree of habitat and diet flexibility, which enables them to occupy habitats across a spectrum of human influences, ranging from undeveloped to highly urbanized (Bateman & Fleming, 2012).
Conversion of land to agriculture may further benefit coyotes due to increased prevalence of edge habitat, which can provide prey and facilitate travel (Theberge & Wedeles, 1989). Rural areas may thus support higher numbers of coyotes compared to forested landscapes (Richer et al., 2002). Coyotes also inhabit more developed areas where they can find a variety of anthropogenic food (Murray et al., 2015), but may be deterred from some such habitats due to conflict potential with humans (Gosselink, Van Deelen, Warner, & Joselyn, 2003).
Dietary overlap between coyotes and foxes can result in similarities regarding habitat selection (Theberge & Wedeles, 1989).
Like coyotes, red foxes may use urban and agricultural habitats (e.g., Gosselink et al., 2003;Harris & Smith, 1987). However, red fox habitat selection is also sometimes influenced by coyote predation risk in areas where they are sympatric (Holt & Polis, 1997;Sargeant & Allen, 1989). Increases in coyote abundance can depress red fox populations, either through competition or direct killing (Levi & Wilmers, 2012). Coyotes may also displace red foxes from preferred habitats, leading them to select habitats avoided by coyotes, such as areas of higher human density (Gosselink et al., 2003;Randa & Yunger, 2006).
Habitat selection occurs at multiple spatial scales, ranging from a microhabitat within a patch of an individual's home range (4th order selection) to the geographic distribution of the species (0 order selection; Meyer & Thuiller, 2006). Many studies of coyote and red fox habitat selection have examined individual home ranges, focusing on habitat selection at the 2nd or 3rd order scale (e.g., Gosselink et al., 2003;Theberge & Wedeles, 1989). Comparatively less work has examined habitat selection of multiple coyote and red fox populations on a broader spatial scale (1st order selection; e.g. Kays, Gompper, & Ray, 2008;Levi & Wilmers, 2012). Understanding these broad-scale patterns in canid habitat use is important in landscapes increasingly altered by humans, driving increased human-wildlife interactions and potentially conflict (Gompper, 2002;Weckel, Mack, Nagy, Christie, & Wincorn, 2010).
We examined how broad-scale landscape characteristics influenced coyote and red fox spatial distribution and relative abundance across rural areas of New York State, USA, using a citizen science dataset collected by volunteer bow hunters. Bow hunters are particularly suited to collect field data for scientific studies; they are camouflaged, observant, and quiet, typically remaining still in one place for several hours at a time, allowing relatively undisturbed observations of wildlife. We predicted that coyote and red fox relative abundance would show positive responses to the amount of human-dominated land cover (agriculture and low-density development), with a larger response by red foxes. We also predicted a positive correlation between coyote relative abundance and degree of landscape fragmentation. Finally, we predicted that red fox abundance would be negatively correlated with coyote abundance.

| Study area
Our study included legally hunted areas across New York State (total area 141,300 km 2 ), including state and private land. The vast major- The average WMU size was 1,364 km 2 (range 269-7,999 km 2 ).

| Canid relative abundance data
We obtained hunter counts of coyote and red fox from the New year. We excluded from the study three WMUs representing the urban areas of New York City and Buffalo, because no hunting was permitted (and thus no bowlog data collected) in these areas. Some records from 2012 were missing, so we excluded the entire year of data. Since the dates of the bow hunting season (and thus the dates of data collection) varied among years and between WMUs, we limited the dataset to logs recorded during 1-25 October each year to reduce potential temporal bias. Hunters did not record the time of day of their observations; we assumed hunter effort was concentrated at the beginning and end of daylight hours to coincide with greater deer activity (e.g., Higdon, Diggins, Cherry, & Ford, 2019), and that this pattern was similar over time and among WMUs. We found that there were some extreme outlier counts in the dataset, likely a result of data input error. To avoid bias introduced by these outliers, we removed from the dataset all counts greater than 10, which represented less than 0.01% of observations in the dataset. We also removed records with missing or unrealistic estimates (>12 hr in a single day) of effort (about 3% of observations).

| Covariate data
We obtained land cover rasters at 30-m resolution for the state of New York for years 2004,2006,2008,2011,2013, and 2016 from the National Land Cover Database (NLCD; Yang et al., 2018). For each of these years, we calculated the proportion of agricultural land use (classes 81-82), and the proportion of low-density developed land use (classes 21-22) for each WMU. We did not include medium-and high-density developed areas because they are typically not hunted and did not make up a large percentage of land area. We also excluded areas of open water from the total area in these calculations.

| Analysis
We fit generalized linear mixed models (GLMMs) to the canid count data in program R (R Core Team, 2019) using R package lme4 (Bates et al., 2019). We fit separate models for coyote and red fox. In each case, the response variable was the number of observations of the species in a given WMU in a given year. Given that the response variables were counts, we examined several possible generalized lin- we selected a Poisson model with an OLRE .

Accounting for overdispersion is important in Poisson regression
given that ignoring it can result in bias in parameter estimates, incorrect inference, or inflated estimates of R 2 (Harrison, 2014). Harrison (2014) reported that use of an OLRE in Poisson regression may not be suitable when overdispersion is due to zero-inflation. In our dataset, zero-inflation was likely not an issue; only 3.9% and 11.0% of total observations in a given WMU and year were zeros for coyotes and red foxes, respectively. Thus, we considered the OLRE modeling approach appropriate for our dataset. In both models, we also included log(observation hours) as an offset term in the linear predictor to account for varying effort across space and time, as well as random effects of WMU and year.
Fixed effects in each model included the proportion of WMU area in agricultural land-use types, the proportion of WMU area in the low-density developed land-use types, and the WMU agriculture land-use area-perimeter ratio. We also included second-order terms for each of these covariates based on a simple heuristic: If both the first-and second-order terms were statistically significant (p < .05), the second-order term was retained in the model; otherwise it was removed. In the red fox model, we also included a covariate representing coyote relative abundance in the WMU (number of coyotes observed per 100 observation hours). We tested for multicollinearity among these covariates using variance inflation factors (VIF), assuming multicollinearity was not an issue as long as all VIF < 3 (Zuur, Ieno, & Elphick, 2010). To assess model fit, we calculated marginal and conditional R 2 using R package MuMIn (Barton, 2019). We determined whether estimated fixed-effect coefficients were statistically different from 0 using Wald tests with α = 0.05.

| RE SULTS
Our final dataset included 342,710 hunter observations (i.e., counts of canids at a specific location on a single day). On average, there were 26,362 observations per year (range: 18,495-34,566).
Average observations per WMU (totaled across all years) was 3,851 (range: 377-12,530). Average hunter effort per observation was 3.61 ± 1.86 hr (mean ± SD). At least one coyote was counted in 12,047 (3.5%) hunter observations, and at least one red fox was counted in 10,356 (3.0%) hunter observations. Mean total effort per WMU-year was 1,070 ± 677 hunter-hours. Across all WMUs, mean count of coyotes per 100 hunter-hours was 1.62 ± 0.73, and mean count of red foxes per 100 hunter-hours was 0.82 ± 0.71. While observations varied among years, there were no consistent temporal trends in coyote and red fox counts over time (Figure 1). Coyote and red fox counts, as well as landscape variables, varied regionally across the state (Figure 2).
We found a significant negative effect of percent low-density developed land cover on coyote abundance, with a 1-standard deviation increase in low-density developed land cover (a change of about 11%) corresponding to a 21% decline in coyote counts per unit hunter effort (Table 1, Figure 3). Increasing fragmentation of agricultural habitat (i.e, an increase in the perimeter-area ratio of patches in the landscape) had a significant positive effect on coyote counts, with a 1-SD increase in perimeter-area ratio (about 0.008) increasing coyote counts by 14% (Table 1, Figure 3). Total amount of agricultural land cover did not significantly affect coyote counts, and no second-order terms were retained in the final model based on our heuristic. Marginal R 2 (variance explained by the fixed effects) for the coyote model was .09, and conditional R 2 (variance explained by both fixed and random factors) was .87.
In contrast to coyotes, both low-density developed and agricultural land cover had positive effects on counts of red fox per unit hunter effort (Table 2, Figure 3). For agricultural cover, both firstand second-order terms were significant and positive, corresponding to an increase in red fox counts with agricultural cover (Figure 3).
For low-density developed cover, only the linear term was included in the model; a 1-SD increase in low-density developed cover (about 11%) increased red fox counts by 33% (Table 2). Also in contrast to coyotes, perimeter-area ratio of agricultural land cover had a significant negative effect on red fox counts, with a 1-SD increase in the ratio decreasing expected red fox counts by 17% (Table 2, Figure 3).
There was no significant relationship between coyote and red fox counts per unit hunter effort (Table 2). Marginal R 2 for the red fox model was .28, and conditional R 2 was .87.

| D ISCUSS I ON
Using citizen science data, we found differences in the responses of coyotes and red foxes to anthropogenic habitats across New York State. We found mixed support for our prediction that both species would show positive responses to human-dominated land cover, as development was negatively associated with coyote abundance, but positively associated with red fox abundance. Similar Coyotes also did not show the predicted response to agricultural habitat. Suitability of agricultural habitat for coyotes varies extensively, which may have contributed to its lack of influence on coyote counts. Coyotes select crops that provide cover and food, such as corn, while avoiding those that do not, such as sorghum and soybeans (Andelt & Andelt, 1981;Gosselink et al., 2003). Row crop fields are dynamic environments in which harvest results in a substantial loss of vegetation over a short time period. This influences the amount of food and cover provided by the crop and can lead to seasonal changes in agricultural habitat use (Andelt & Andelt, 1981).
Abrupt losses in vegetation cover can also put coyotes at risk of being killed by humans due to increased visibility and lack of refuge (Van Deelen & Gosselink, 2006). Thus, benefits provided by increased foraging opportunities in agricultural habitats may be offset by increased mortality risk. This fine-scale selection of agricultural habitat and temporal shifts in its suitability could contribute to lack of coyote response to agricultural habitat across the broad spatial and temporal scales we assessed. These differences in spatial and temporal scales may also be reflected in the overall poor fit (i.e., low marginal R 2 value) of the coyote model. Coyote count data collected at finer scales (e.g., via camera traps) are likely required to better understand relationships with land-use variables.
The patterns of coyote abundance we observed may also reflect their population expansion into the state. Coyotes initially colonized New York from the north and circled the Adirondack region before  spreading south and west (Fener et al., 2005). This largely corresponds with the densities we documented with highest densities in the north and decreasing density moving south and west (Figure 2).
Current coyote abundances across New York may in part be a relic of their historic expansion route, with highest densities in the regions into which they first expanded. Observed lower coyote abundance in low-density developed areas may also be a function of changes in coyote behavior. Coyotes appear to become more nocturnal in urban and suburban areas (Grinder & Krausman, 2001;Jantz, 2011), which would make them less detectable by bowhunters active between sunrise and sunset, relative to coyotes in rural areas.
In contrast to coyotes, red fox counts were positively associated with both agricultural and low-density developed land use, supporting our predictions. While it may not be ideal habitat, particularly following harvest (Gosselink et al., 2003), agricultural habitats in New York provide a variety of food for red foxes (Cook & Hamilton, 1944) and open fields are frequently selected for denning because the loose soil facilitates easier digging than in wooded areas (Sheldon, 1950).
Our findings match previous work showing red foxes are dependent on human-dominated habitat types in areas where they overlap with coyotes (Gosselink et al., 2003).
In accordance with our predictions, coyote abundance increased with agricultural perimeter ratio, while red fox abundance declined.
Though edge effects are often considered beneficial to carnivores, there is considerable debate regarding this generality. Theberge and Wedeles (1989) also found coyotes to have higher abundances along edges, while others have reported no effect or a negative response to edges (Hinton et al., 2015). Within the same population, resident and transient individuals may have divergent responses to edge habitat (Hinton et al., 2015). Similarly, red foxes have been shown to have higher abundances (Heske, 1995) or no response to edges (Villaseñor, Blanchard, Driscoll, Gibbons, & Lindenmayer, 2015), although we documented a negative response. These disparate responses to habitat fragmentation across their range reiterate the intraspecific habitat selection variation that exists for these species, which has contributed to their widespread geographic ranges. Our results suggest fragmentation by agriculture across New York could augment populations of coyotes but may have the opposite impact on red foxes.
Several mechanisms facilitate coexistence of coyotes and red foxes, which may account for the lack of the predicted inverse relationship between species abundances. Habitat flexibility allows red foxes to shift habitat use to areas avoided by coyotes (Gosselink et al., 2003;Mueller et al., 2018). Red foxes may also alter their diets to avoid competition and potential interactions with coyotes (Theberge & Wedeles, 1989). Additionally, red foxes have smaller home ranges than coyotes, enabling them to persist in relatively small areas between coyote territories (Harrison, Bissonette, & Sherburne, 1989). Coexistence and co-occurrence at the larger, regional scale (sympatry) may be facilitated by spatial segregation at smaller, local scales (Amarasekare, 2003).
Our study adds to the growing research that employs citizen science data to understand canid ecology. Other studies have used citizen science data to examine the distribution and density of coyotes and red foxes (Mueller, Drake, & Allen, 2019;Scott et al., 2018). Citizen science data have also provided insights into the effects of roads on canids and the factors contributing to humancanid conflicts (Vercayie & Herremans, 2015;Wine, Gagne, & Meentemeyer, 2015). Pack dynamics of a recovering wolf population were studied using scat collected by volunteers (Granroth-Wildling et al., 2017) and photographs captured by the public provided data on disease spread of urban foxes (Scott et al., 2020). Collectively, this work suggests the potential for substantial advances in understanding human influences on canids using citizen science data. Such potential warrants the investment of resources into refining data collection procedures for citizen science data to glean the maximum insights from this resource.
Due to the characteristics of the surveys used in this study, our data closely match Kamp, Oppel, Heldbjerg, Nyegaard, and Donald (2016) definition of structured citizen science data or monitoring program: All surveys occurred at a specific time of year (October), we accounted for observation effort (hours in the field), and all species observed within a specific list were documented, including common ones (no rarity bias). However, like unstructured data, our data were comprised of opportunistic field observations not randomly located throughout the state; the counts represent relative abundance in areas used by bowhunters and not a random spatial sample within a WMU. Having volunteers document hours spent per survey was key to modeling the sampling process, providing the greatest potential for delivering robust and timely trends based on opportunistic data (Isaac et al., 2014). Population indices derived from chance observations of wildlife per unit of effort are not only of low cost to state wildlife agencies but can provide more accurate assessments than those derived from other types of data, such as spotlighting (Winchcombe & Ostfeld, 2001) or furbearer indices (Obbard et al., 1987), which are influenced by variations in pelt prices (e.g., Conlee & Johnston, 2018). The New York State DEC seeks a more uniform distribution of survey respondents by actively recruiting additional bowhunter volunteers in areas with lower participation (Fies & Norman, 2003).
Voluntary work by hunters and the general public is the foundation of this and many other wildlife monitoring programs around the United States and the globe. Many eastern and midwestern US states monitor target species through hunter observation data that are used for quantifying trends and making management decisions (Fies & Norman, 2003 Despite its shortcomings, opportunistically gathered data have the potential to make meaningful contributions in biodiversity science and policy-making (Isaac et al., 2014;Tulloch, Possingham, Joseph, Szabo, & Martin, 2013). Particularly, large citizen science data sets collected at broad spatiotemporal extents can reduce the amount of unsampled variation and lower the risk of false inferences (Bain, Wayne, & Bencini, 2015), and have potential to be integrated with other data sources to model population level dynamics and improve statistical inferences (Sun, Fuller, & Hurst, 2018). Data collected by citizen scientists can be equivalent in reliability and precision to that produced by traditional monitoring approaches (Rafiq et al., 2019). Citizen scientists can provide data used to investigate an array of ecological processes such as the spread of invasive species, wildlife disease dynamics, and the impact of climate change (Dickinson et al., 2012). Participants in citizen science projects often emerge with a greater appreciation of wildlife and nature (Toomey & Domroese, 2013). The engagement of citizens in scientific processes has the potential to combine the collection of useful data with outreach and education, helping to close the knowledge gap between academia and the general public (Jenkins, 2011).

ACK N OWLED G M ENTS
We thank the many bow hunters in New York State who collected data used in this study. Two anonymous reviewers provided helpful comments on an earlier version of this manuscript. Support for this research was provided by the Federal Aid in Wildlife Restoration Act (Grant W-173-G) and the Camp Fire Conservation Fund.

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data and code used in the paper are available on Data Dryad with https://doi.org/10.5061/dryad.0cfxp nvzs.