Boreal predator co‐occurrences reveal shared use of seismic lines in a working landscape

Abstract Interspecific interactions are an integral aspect of ecosystem functioning that may be disrupted in an increasingly anthropocentric world. Industrial landscape change creates a novel playing field on which these interactions take place, and a key question for wildlife managers is whether and how species are able to coexist in such working landscapes. Using camera traps deployed in northern Alberta, we surveyed boreal predators to determine whether interspecific interactions affected occurrences of black bears (Ursus americanus), coyotes (Canis latrans), and lynx (Lynx canadensis) within a landscape disturbed by networks of seismic lines (corridors cut for seismic exploration of oil and gas reserves). We tested hypotheses of species interactions across one spatial‐only and two spatiotemporal (daily and weekly) scales. Specifically, we hypothesized that (1) predators avoid competition with the apex predator, gray wolf (Canis lupus), (2) they avoid competition with each other as intraguild competitors, and (3) they overlap with their prey. All three predators overlapped with wolves on at least one scale, although models at the daily and weekly scale had substantial unexplained variance. None of the predators showed avoidance of intraguild competitors or overlap with prey. These results show patterns in predator space use that are consistent with both facilitative interactions or shared responses to unmeasured ecological cues. Our study provides insight into how predator species use the working boreal landscape in relation to each other, and highlights that predator management may indirectly influence multiple species through their interactions.

Maintaining these interactions is therefore vital to maintaining resilient ecological communities in an increasingly anthropocentric world (Valiente-Banuet et al., 2015).
The Canadian boreal forest is one biome experiencing dramatic landscape change due to natural resource extraction, which is creating a heterogeneous "working" landscape shared between industry and wildlife (Pickell, Andison, Coops, Gergel, & Marshall, 2015).
Logging alone has an area footprint of approximately 15 million hectares, while energy development has created over half a million kilometers of linear features across the region (Pasher, Seed, & Duffe, 2013). The cumulative effects of these extensive industrial footprints impact the distribution and abundance of a number of boreal mammals, though the strength and nature of influence vary by species Toews, Juanes, & Burton, 2018).
Behavioral and population changes in individual species can in turn affect interactions among species (Ritchie & Johnson, 2009), leading to broader indirect effects of industrial development in boreal ecosystems.
Like wolves, evidence suggests that black bears prefer linear features for ease of travel Tigner, Bayne, & Boutin, 2014). Mesocarnivores such as coyotes (Canis latrans) and Canada lynx (Lynx canadensis) may also benefit from linear features.
Human development, including linear corridors, is positively related to coyotes' northern expansion and persistence in the boreal forest (Hody & Kays, 2018) and lynx presence ; but see Bayne, Boutin, & Moses, 2008). With at least some members of the boreal predator community responding to extensive seismic line networks, we expect these networks to potentially influence the interactions among predators. More specifically, we predict that linear disturbances result in changes in the co-occurrence patterns of sympatric predator species. We hypothesize that these changes may be driven by three different types of intraguild interactions.
First, wolves may exert top-down influences on subordinate predators that influence how the latter use the landscape (Estes et al., 2011).
Second, subordinate predators may experience intraguild competition from one another. Although there is little evidence in the literature of intraguild competition between black bears and other subordinate predators, competition between coyotes and lynx has been insinuated in areas of human disturbance (Bayne et al., 2008) and may influence habitat selection (Murray, Boutin, & O'Donoghue, 1994). Third, all three subordinate predators should track availability of their prey. Black bears predate boreal ungulates (Linnell, Aanes, & Andersen, 1995;Pinard, Dussault, Ouellet, Fortin, & Courtois, 2012).
In this study, we investigated how apex predators -wolves-and subordinate predators-black bears, coyotes, and lynx-interact on a disturbed landscape by assessing spatiotemporal co-occurrences within a network of linear anthropogenic features. We developed predictions of co-occurrence patterns consistent with hypothesized interactions. Specifically, we predicted that (1) if black bears, lynx, and coyotes experience top-down pressures from wolves, they should spatiotemporally segregate from wolves; (2) subordinate predators should segregate from their intraguild competitors with whom they share resources (O'Donoghue et al., 2001;Guillaumet, Bowman, Thornton, & Murray, 2015); and (3) subordinate predators should overlap with their prey (Keim, DeWitt, & Lele, 2011;Theuerkauf, 2009). We further hypothesized that higher densities of disturbance would result in lower occurrences of black bears, but higher occurrences of lynx and coyotes , and that coyotes and lynx would occur less frequently in the winter (Pozzanghera, Sivy, Lindberg, & Prugh, 2016).
To test these hypotheses, we modeled black bear, coyote, and lynx spatiotemporal occurrences in relation to wolf and prey occurrences, and coyote and lynx spatiotemporal occurrences relative to one another (Table 1). We further compared the effects of interspecific interactions on predator occurrences with the effects of season and anthropogenic disturbance, allowing for interactions between the two (Table 1).

| Sampling design
Our study was located along the east side of the Athabasca River, approximately 70 km southwest of Fort McMurray, Alberta, Canada (56.2588 N, 112.6909 W; Figure 1). The 570-km 2 study area is bounded to the north and west by the Athabasca River and has linear feature density of 1.1 km/km 2 , including 523.6 km of seismic lines (Figure 1; Alberta Biodiversity Monitoring Institute, abmi.ca).
We deployed Reconyx PC900 HyperFire camera traps (Reconyx, Holman, WI) between November 2015 and 2016 as part of an ongoing monitoring project assessing wildlife responses to seismic line restoration. One of the objectives of this restoration was to deter predator movements, but the effect on predator line use was minimal (Tattersall, Burgar, Fisher, & Burton, 2019). We selected 60 camera sites (Figure 1) based on a stratified random design to sample along seismic lines that spanned a gradient of restoration strata. We installed one camera at each site at a height of approximately 1 m and facing across a seismic line. Cameras were set up at least 500 m apart to increase the probability of independent detections (Tigner et al., 2014). We set all cameras to take one image per trigger, with a one-second lag between triggers and no quiet periods.
We considered a survey period as 30 continuous months between November 2015 and April 2018 for coyotes and lynx, and the two eight-month periods between April and October 2016 and 2017 for black bears. We treated detection events of the same species as independent when occurring at least 30 min apart (Rovero & Zimmermann, 2016). All methods for wildlife monitoring were approved by the Canadian Council of Animal Care administered by the University of British Columbia (protocol A17-0035).

| Analytical framework
For large mammal species ranging across entire landscapes, interspecific interactions can be inferred using spatiotemporal relationships of species occurrences (Cusack et al.,  Note: Models were negative binomial GLMs at the spatial-only scale, and binomial GLMMs at the two spatiotemporal scales. Each model set corresponds to a hypothesized interspecific interaction. We tested models with co-occurring species as a predictor variable against three base models describing environmental effects. Candidate model sets for mesocarnivores (coyote and lynx) are identical, with mesocarnivore 1 describing the responding predator and mesocarnivore 2 describing the co-occurring intraguild competitor (e.g., when mesocarnivore 1 is coyote, mesocarnivore 2 is lynx and vice versa). At the spatial-only scale of analysis, we excluded season models for all species because the response variable aggregated detections across the entire survey period.
Correlation is not equivalent to causality, but examining these associations reveals whether patterns in predator co-occurrences are consistent with those predicted from interspecific interactions.
For interactions between predators, we assumed segregation of occurrences in space or time was suggestive of competition, while overlap in space and time was suggestive of facilitation (Cusack et al., 2016;Fahrig, 1992;Frey, Fisher, Burton, & Volpe, 2017;Karanth et al., 2017). To address multiple potential scales of interactions, we examined predator co-occurrences at the spatial-only scale (i.e., entire survey period) and at two finer spatiotemporal scales (weekly and daily). We predicted that spatial segregation was indicative of competition (Fuller & Keith, 1981), whereas spatial overlap alone indicated weak or indirect interactions. We assumed that overlap at finer spatiotemporal scaleswhere species co-occur at a given location within a given occasion length-suggested intentional proximity between two species, and thus was evidence of a facilitative interaction (Cusack et al., 2016;Swanson et al., 2016).

| Spatial relationships in species' occurrences
To examine spatial relationships between predators, we tested the degree to which spatial variation in species detections across the entire survey (i.e., relative abundance or frequency of use) was explained by the detections of potentially interacting species, using number of independent detections per camera trap site as the response variable. We included each camera's survey effort (number of days active) as a predictor variable, thus accounting for the effect of periods of camera inactivity on detections. This predictor variable was the only fixed effect in the null model and excluded from comparisons of effect size. We modeled detections in a generalized linear model (GLM) framework using a negative binomial distribution for overdispersed count variables (Bolker, 2009). We conducted all statistical analyses using the R package glmmTMB (Brooks, Kristensen, & Benthem, 2017).

| Fine and coarse spatiotemporal scales of occurrence
To test whether species co-occurring in space were also co-occurring in time, we examined co-occurrences at two temporal scales. At the finest temporal scale, we recorded the presence (1) or absence (0) of a species' at a given site within each day, producing a binary occurrence metric for each camera trap day.
Given our observed detection rates, a daily occasion length was the finest temporal scale we could reliably model. However, a single day occasion length results in low detection rates, leading to zero-inflated occurrence data (Rovero & Zimmermann, 2016). We therefore also discretized occurrence data into week-long periods to assess whether results were consistent across temporal scales.
For both spatiotemporal scales, we modeled occurrences using bi-

| Habitat modeling and scale analysis
We analyzed data using a two-step process: species-habitat modeling followed by co-occurrence modeling (Chow-Fraser, 2018; Cusack et al., 2016). In the first step, we modeled species response to habitat at multiple spatial scales (Fisher, Anholt, & Volpe, 2011;Levin, 1992). We ranked species-habitat models at multiple spatial scales using model selection with Akaike's information criterion (AIC, Burnham & Anderson, 2002)  the proportion of open forest as a measure of forest density (Murray et al., 1994).

| Co-occurrence modeling
In the second step of our two-step approach, we added heterospecific occurrences as variables to the best supported habitat model and weighed evidence for their ability to explain additional variation in predator occurrences (Table 1; Fisher et al., 2013). We created a model set for each of the interactions hypothesized to influence black bears, coyotes, and lynx, namely top-down influences from wolves, bottom-up influences from prey species, or competitive influences of lynx and coyotes on each other. We excluded models assessing interspecific interactions between black bears and coyotes or lynx because we found no supporting evidence for these interactions in the literature (Table 1). To assess influences of prey, we included a variable aggregating detections of all prey species for each of the target predators ( Table 2). As we further predicted that the strength of interspecific interactions could be influenced by season and level of anthropogenic disturbances, we also included additive and interaction models with variables for snow presence and linear density (  Figure 3b). We used data exploration techniques prior to modeling to assess all predictor variables for outliers, collinearities, and heterogeneity of variance (Zuur, Ieno, & Elphick, 2010). We also scaled all nonbinary variables by subtracting the mean and dividing by two standard deviations, thus improving model convergence and interpretation (Gelman, 2008).

| Model interpretation and model validation
Following co-occurrence modeling, we compared candidate models using AIC model selection to determine whether interspecific interactions influenced predator occurrences at each of the three spatiotemporal scales. We considered all models within 2ΔAIC of the top-ranked model as having similar explanatory power over the data (Burnham & Anderson, 2002). We consequently examined variables in these models for their influence on predator occurrence, with mean parameter estimates as measures of effect size and direction and 95% confidence intervals as measures of statistical significance (i.e., not overlapping zero). We used pseudo-R 2 to assess the proportion of variance explained by the top models for each species relative to the proportion of variance explained by a null model (McFadden, 1974). We considered R 2 values between 0.2 and 0.4 to be indicative of good model fit (McFadden, 1977).
In both top-ranked models for coyotes and lynx at the daily scale, standard errors around the estimates for interacting species were over three orders of magnitude larger than the estimates (−15.730 [−7943, 7912] for the top coyote model and −15.763 [−8613, 8582] for the top lynx model), and patterns in the model residuals indicated model misspecification. We therefore removed these models from all subsequent analyses.

| RE SULTS
From November 2015 to April 2018, total sampling effort was 14,054 camera-days for black bears (during their active season) and 32,436 camera-days for all other species (year-round). Of the four focal boreal predator species, black bears and wolves were detected most frequently, followed by coyotes and lynx (Table 3). At the daily spatiotemporal scale, both coyotes and lynx infrequently co-occurred with interacting species, and lynx infrequently co-occurred with their prey at the weekly scale (Table 3). For all species, models at the daily and weekly scales explained relatively little variance in occurrence data (pseudo-R 2 = 0.008-0.074), while models at the spatial-only scale performed much better (pseudo-R 2 = 0.166-0.283; Table 4).
Heterospecific occurrences were significant predictors of all three predators' occurrences, explaining variance in addition to that ex- Prey models were included in the top-ranked models for lynx at the both the spatial-only and weekly occurrence scales (Table 4).
Lynx were less likely to occur at high linear density sites where prey were present (β = −1.074 [−1.934, −0.214] at the spatial-only scale; Figure 2). Black bears and coyotes were not affected by prey occurrences at any scale (Table 4).
TA B L E 2 Full list of predictor variables used to model occurrence patterns of black bears, lynx, and coyotes

Predictor variables
Step of modeling For each species, we included all habitat variables in the first step of the modeling process, and retained habitat variables with confidence intervals that did not overlap zero to create a null model for the second step. We measured forest cover variables from the Alberta Vegetation Inventory (Alberta Vegetation Interpretation Standards, 2005) and linear feature data from the Alberta Biodiversity Monitoring Institute (ABMI, abmi.ca). We used camera trap data to extract all species occurrence and snow variables. Species' variables differed with modeling scale, as spatiotemporal scale affected occasion length and thus occurrence aggregation.
The occurrence of other mesocarnivores was not a strong predictor for either coyotes or lynx. For coyotes, models with lynx were included among top-ranked models at the daily and weekly scales, but neither main effects nor season and linear density interactions were significant predictors. The same was true for lynx, although coyote occurrences only helped explain lynx occurrences at the daily scale (Table 4).
All three predators responded to linear density on at least two

| Inferring interspecific interactions from predator co-occurrences
Black bears, coyotes, and lynx all co-occurred with wolves on at least one of the spatial and temporal scales observed in our study ( Figure 2). This suggests either that wolf spatial ecology is a determinant of space use by subordinate predators in working landscapes or that all predator species are cueing into the same resources in time and space in this complex environment (DeMars & Boutin, 2018 Relative to other species, wolves and black bears occur at high densities in the boreal forest (Burgar, Burton, & Fisher, 2018), and both select industrial linear features for easy travel, making them likely candidates for strong interspecific interactions in industrializing landscapes Latham, Latham, Boyce, & Boutin, 2011). Although anecdotal evidence describes aggressive interactions between individual wolves and bears (Palomares & Caro, 1999;Rogers & Mech, 1981), we propose that spatiotemporal overlap at the daily scale may be consistent with facilitative interaction between the two. Black bears were more likely to occur at a site even within a day of wolf occurrences, mirroring spatiotemporal patterns in lion-kill scavengers in Africa (Cusack et al., 2016;Swanson et al., 2016). Further, we found that black bears occurred more frequently at high linear density sites when wolves frequently occurred at those sites. Black bears are adept scavengers and may benefit considerably from carrion subsidies left by wolves (Allen et al., 2014;Wilmers et al., 2003). However, our results indicate that there are additional factors influencing black bear occurrences that were not included in our study (i.e., much unexplained variance in occurrences). Further research should investigate these sources of variance, as well as explore time elapsed between predator occurrences or analyze patterns in occurrences (i.e., which species occurs first) to observe these relationships at a finer temporal resolution (Schliep, Gelfand, Clark, & Kays, 2018;Swanson et al., 2016).
The positive association between coyotes and wolves at coarse spatiotemporal scales may also be a result of facilitation ( Figure 2).
However, as sympatric canid species with considerable niche overlap, coyotes and wolves are also likely to experience strong competition in which wolves frequently kill coyotes (Palomares & Caro, 1999;Paquet, 1991). We hypothesize this paradox could reflect density-dependent interactions: Where coyotes exist in low densities, they segregate themselves from wolves to reduce competition.
When coyotes exist in high densities, they may be able to reduce competition through behavioral mitigations such as increased group size or fine-scale temporal partitioning, thus increasing the benefits of scavenging (Atwood & Gese, 2008. Coyote density (2.64/100 km 2 ) eclipsed wolf density (0.77) by threefold south of our study area (Burgar et al., 2018 Note: Each value represents the total number of times both species were present at the same site and-for weekly and daily scales-within the same occasion. Rows represent response variables, and columns represent predictor variables. The total occurrences of wolves are given both within the summer-only sampling period for black bears (14,054 site-days) and the full sampling period for coyotes and lynx (32,436 site-days). Cells are marked with a dash where no interactions were hypothesized or tested.

TA B L E 3
Total co-occurrences of predator species across three spatiotemporal scales of analysis Note: For each species, top-ranked models are shown for each of the three spatiotemporal scales of analysis. Top-ranked models were those within 2ΔAIC of the highest weighted model. The column k is the number of parameters in the model, ΔAIC indicates the difference in AIC scores from the top model, and R 2 is a pseudo-R 2 measure describing the proportion of variance explained by each model relative to the variance explained in the null model. The top models for coyote and lynx at the daily scale were not included in subsequent analyses due to large confidence intervals.

TA B L E 4 Model selection tables of models of top-ranked models for black bears, coyotes, and lynx
detections, suggesting lower coyote relative abundance. Coyotes in this area may therefore be avoiding direct competition but benefitting indirectly via scavenging. This would be consistent with the observed co-occurrence at the weekly scale and lack thereof at the daily scale. However, as we did not observe segregation at the daily scale and these models had substantial unexplained variance, this hypothesis requires further inquiry.
Lynx and wolves spatially co-occurred, but did not temporally co-occur with wolves at the weekly or daily scale (Figure 2). This suggests that although lynx share a landscape with wolves, they may not interact at finer temporal scales. We also found an unexpected negative interaction between lynx prey and linear density at the weekly scale, where lynx were less likely to occur at high linear density sites where their prey were present. This result may be a spurious result arising from low co-occurrences between lynx and their prey (n = 2; Table 3). Given the poor fit of this model, we suggest that further research is needed on additional factors influencing lynx occurrences at fine temporal scales to elucidate this finding.
The negative effect of linear density on black bears, as well as the positive effect on coyotes and lynx, agrees with previous research conducted in northern Alberta Toews et al., 2018). We only found evidence of black bear and wolf co-occurrences changing as a function of anthropogenic disturbance at the spatial-only scale, indicating increasing overlap between these predators with increasing anthropogenic disturbance. This was also seen by other studies in tropical, semiurban, and mountain ecosystems (Chow-Fraser, 2018;Karanth et al., 2017;Wang et al., 2015). This system has less of a disturbance gradient than other landscapes in Alberta's boreal forest (Government of Alberta, 2017), which might explain the lack of response from other species or at finer temporal scales. Further, humans are largely absent from the study area, so direct human influence on interspecific interactions would be minimal. To better assess the influence of anthropogenic disturbance on interspecific interactions, a similar study could be conducted across a number of landscapes with varying levels and types of landscape change (Chow-Fraser, 2018).
An alternative interpretation of our results is that the overlapping occurrences we observed could be a result of predators responding to similar resource cues rather than responding to each other. Bears and wolves select linear features and linear feature density in similar ways (DeMars & Boutin, 2018;Finnegan et al., 2018), while coyotes and wolves may also sometimes select the same habitat (Latham et al., 2013). Like black bears and wolves, coyotes and lynx may also  . We suggest that further research into characteristics influencing linear feature use by predators would provide insight into species co-occurrences in a working landscape.

| Limitations
Fine-scale temporal analyses of interspecific interactions often require large sample sizes of detections to reveal patterns in activity and co-occurrence (Frey et al., 2017). In our study, low numbers of occurrences-and consequently, co-occurrences-may have limited our ability to reliably model effects of interacting species (Table 3).
Indeed, this is likely the cause of unexplained variance in our models at finer spatiotemporal scales. We urge caution when modeling co-occurrences from rare or elusive species at fine spatiotemporal scales. More robust methods to assess interactions for these species could use baited camera traps to increase detection probabilities (Stewart et al., 2016), include cameras in undisturbed habitat, or increase either spatial or temporal extent of the camera trap survey to increase sampling effort. Further, telemetry studies of interacting species could account for interspecific effects on movement patterns, thereby assessing how individual animals respond to space use by other species on a shared landscape (James, Boutin, & Hebert, 2004).

Caution is obviously necessary when inferring interactions and
their mechanisms from co-occurrence data. Although techniques exist to derive interaction strength and predict mechanisms, such approaches are either nontemporal (Dorresteijn et al., 2015) or require large amounts of data (Schliep et al., 2018;Swanson et al., 2016) or even direct observation (Atwood & Gese, 2008;Cusack et al., 2016). Whereas we assume that species co-occurrences indicate intentional proximity and thus suggest facilitative interactions, co-occurrences of predators with similar niches may be equally indicative of species competing for a shared resource (Chow-Fraser, 2018

| Management implications
In response to growing wolf populations, and out of concern of high predation rates on woodland caribou, the government of Alberta implements annual wolf reduction programs within some caribou herds at high risk of extirpation (Government of Alberta, 2017; Hervieux, Hebblewhite, Stepnisky, Bacon, & Boutin, 2014). Although effective in boosting caribou numbers in the short term, wolf removal is controversial and has direct consequences for the interactions structuring the boreal mammal community (Darimont, Paquet, Treves, Artelle, & Chapron, 2018;Sivy et al., 2017). As wolf reduction programs continue in caribou ranges in western Canada, we suggest that research should focus not only on caribou response, but also on responses of other species in the boreal mammal community.
Interspecific interactions arise from coexisting species partitioning space, time, and life-sustaining resources on a shared landscape where such resources are limited (Schoener, 1974). Understanding those interactions enables us to predict how they will respond when perturbed, empowering us to make informed and proactive management decisions. Here, we showed that nonapex predators exhibit spatiotemporal overlap with an apex predator on a working landscape. This overlap identifies patterns in how these four boreal predators use this landscape, which may indicate facilitative interactions or responses to the same ecological signals. These species additionally show individual responses to anthropogenic disturbances, though responses vary and further investigation is necessary to evaluate consequences for interactions. Results from this study highlight important considerations of the impact of predator management decisions, which may unintentionally alter the behavior of coexisting species (Burgar et al., 2018). The relationships observed in this study occur in the context of a landscape experiencing ongoing industrial development, offering insight into species coexistence patterns in the face of continuing anthropogenic landscape change.
To keep wildlife communities on such landscapes, we must commit to understanding the underlying relationships that allow them to thrive. Gareau (Canadian Natural Resources Ltd). At ITA, we thank L. Nolan,

F I G U R E 4
Effects of habitat features on predator occurrences in the habitat modeling step of analysis. Effect sizes are shown as parameter estimates (mean ± 95% confidence intervals) from negative binomial GLMs (spatial level) and binomial GLMMs (weekly and daily levels) of black bear, coyote, and lynx occurrences at three levels of analysis. Results are shown from habitat variables measured at the optimal spatial scale of influence: 250 m for black bears, 1750 m for coyotes, and 1500 m for lynx. Note that LowDecid is absent for black bears because lowland deciduous forest did not occur with 250 m of any camera stations. Significant habitat variables (with confidence intervals that did not overlap zero) were then included in the second step of the analysis to model effects of interspecific interactions on predators D. Pan, and A. Underwood for field data collection and GIS support.