2.1 Study area

The study area encompassed all six ranges of boreal caribou in north-east British Columbia (latitude: c. 56.620 to 60.000; longitude: c. −124.160 to −120.000). Landscape composition among these ranges was similar, consisting of a mosaic of deciduous and mixed-wood uplands, poorly drained low-lying peatlands and riparian areas. Terrain is predominantly flat to undulating and the climate is northern continental, characterized by long, cold winters and short summers. Forest fire is a frequent form of natural disturbance with a mean fire interval of c. 100 years. The study area is further notable because it encompasses the Horn River Basin, a geologic formation containing one of the largest deposits of natural gas shale in Canada. Consequently, oil and gas extraction activities are the dominant form of anthropogenic disturbance within caribou range. The density of LFs associated with these activities averaged 4.04 km/km² (SD: 3.23; Appendix S1) across all caribou ranges.

2.2 Caribou and predator spatial data

We used GPS data collected from caribou, wolves and black bears fitted with Iridium satellite GPS collars (Advanced Telemetry Systems, model #2110E; Lotek Wireless, model Iridium Track M 2D). During collar deployment, we targeted mature individuals occurring within or immediately adjacent to caribou range. All capture and handling procedures followed approved animal care protocols (BC RIC 1998; BC Wildlife Act Permits FJ12-76949 and FJ12-80090; University of Alberta Animal Use protocols # 748/02/13 and AUP00001309).

For caribou, we obtained data from 56 reproductive-aged females distributed among all six caribou ranges. Individual females were captured by aerial net-gunning during the winter months of 2011 (n = 25), 2012 (n = 2) and 2013 (n = 29). Collars were programmed to acquire one GPS location (or fix) every 2 or 4 hr during calving (15 April–15 July) and one to three times per day otherwise. Twenty females had functioning radiocollars through the 2014 calving season (i.e. beyond July 15, 2014).

For predators, we used data from 23 wolves and 19 black bears. All individuals were captured by aerial darting. Wolf collars were deployed during the winter months of 2012 (n = 3), 2013 (n = 10), 2014 (n = 4) and 2015 (n = 6) and were distributed among 13 packs with at least one pack being collared in each caribou range. Collars deployed in 2012 and 2013 had fix rates of every 15 min 1 May–30 June and once per day otherwise. Collars deployed in 2014 and 2015 had four 1-week intervals with fix rates of every 5 min and otherwise were either every 2, 3 or 4 hr. For all collars, one of the 1-week intervals occurred in June to capture fine-scale movements during the calving season of caribou. Bear collars were deployed in May 2012 (n = 4) and May 2013 (n = 15), and all had fix rates of every 30 min 1 May–30 June and once per day otherwise. Captured bears occurred in three of the six caribou ranges. All predator datasets spanned one calving season only.

Prior to analyses, we screened the raw data for potential errors (Appendix S2). Following these procedures, fix success rate per collar for caribou averaged 93% (range 68–100) during calving and 84% (range 28–99) outside of calving and per collar monitoring time averaged 541 days (range 112–745). For wolves, fix success rates averaged 90% (range 60–100) during the high-frequency monitoring intervals (i.e. every 5 or 15 min) and averaged 89% (range 56–98) otherwise. Monitoring times per wolf collar averaged 268 days (range 108–546) although one collar had a 4-month gap in monitoring post-calving. For bears, fix success rates averaged 81% (range 61–90) during calving and 52% (range 42–72) otherwise. Collar failures limited monitoring times for bears (mean time: 39 days; range 4–120) and three bears were excluded from further analyses as collar failures limited monitoring time to ≤1 day. Because only six bears had datasets extending beyond 30 June, we limited bear analyses to the calving season. To evaluate the potential influence of low fix success rates on inferences of resource selection (Frair et al., 2010), we conducted sensitivity analyses on population-level results for all three species by first excluding individuals with rates <80% then further excluding those with rates <90% (Appendix S3).

2.3 Predicting caribou parturition status and neonate survival

For caribou, we restricted analyses to GPS locations where a female was accompanied by a neonate calf (hereafter, calving locations). We identified calving locations using the two movement-based methods (MBMs) of DeMars, Auger-Méthé, Schlägel, and Boutin (2013) (Appendix S2). Briefly, MBMs predict the calving status of females (parturient vs. barren) and the survival status of neonate calves and further yield estimates of parturition date and calf loss date, where appropriate. For 2011–2013 data, MBM predictions of calf survival were further corroborated by data from aerial surveys conducted weekly during calving. These procedures yielded a final dataset of 44 female caribou predicted to have calved at least once during the 4-year study period. Nineteen females calved in 2 years, resulting in 63 caribou-calving seasons.

2.4 Resource selection analyses: General framework

To evaluate the spatial overlap hypothesis, we conducted resource selection analyses at two spatial scales for all three species. All analyses focused on evaluating selection patterns during the calving season of caribou. At the fine scale, we assessed selection of peatlands, LFs in peatlands and LFs in other land cover types for each species. At the larger scale, we assessed how LF density affected each species selection of areas both within and outside of peatlands.

For fine-scale analyses, we estimated step selection functions (SSFs), which compared resources associated with each observed movement “step” (i.e. an animal's movement linking successive GPS locations) to resources associated with a matched set of random steps sampled from a parameterized distribution of the observed movement process (Avgar, Potts, Lewis, & Boyce, 2016; Fortin et al., 2005; Appendix S4). For each observed step, we generated 20 random steps to adequately characterize resource availability (Northrup, Hooten, Anderson, & Wittemyer, 2013; Appendix S4). We compared resource values at the end of each step—rather than the average resource values along each step—as this approach is likely better for estimating selection of LFs (Thurfjell, Ciuti, & Boyce, 2014). To explicitly assess how LFs influence predator resource selection during movement, we constrained SSF analyses for bears and wolves to steps where the animal was actively moving and excluded those associated with resting or feeding (see Section 8 below). For caribou, we used all steps as our interest was in quantifying total exposure to LFs. All SSF analyses were confined to locations within the calving season when GPS collars were programmed for higher fix rates (caribou: every 2 hr, wolves: every 5 to 15 min; bears: every 30 min).

For large-scale analyses, we used two approaches, both focusing on a species’ response to LF density measured over two areal extents (200-m and 1,000-m radii). For caribou and wolves, we estimated third-order resource selection functions (Johnson, 1980), which compared resources associated with an animal's GPS locations to resources associated with available locations randomly sampled within the animal's annual home range. We characterized home ranges using the 90% isopleth of utilization distributions (UDs; Börger et al., 2006). For each UD, we specified the “reference bandwidth” as the smoothing parameter, which provided a good visual approximation of space use for each caribou and wolf (mean home range size for caribou: 802 km² [range 84–2,904]; for wolves: 1,222 km² [range 493–2,101]. For both species, we conducted sensitivity analyses to determine the time required for estimates of UD size to stabilize, excluding individuals falling below this minimum (Appendix S5). We also conducted sensitivity analyses to determine the minimum number of random points to adequately characterize availability within the UDs of each species (Appendix S5; Northrup et al., 2013). For bears, limited monitoring times prevented estimation of annual home ranges; therefore, we used SSFs with data rarefied to two locations per day, which allows sufficient time for bears to move among areas with potentially varying LF density. As with fine-scale analyses, we generated 20 random steps for each observed step.

2.5 Environmental covariates

We modelled resource selection with variables representing land cover and LFs. For land cover, we used Enhanced Wetlands Classification data from Ducks Unlimited Canada (Appendix S6), collapsing these data into seven biologically meaningful categories (Table 1). We specifically created a peatlands category by combining bogs and nutrient-poor fens to represent land cover highly selected by caribou during calving (DeMars, 2015). For LFs, we accessed provincial repositories to obtain year-specific data depicting roads, seismic lines and pipelines (Appendix S7) and we matched these data by year to the GPS radiocollar data. For each line type, we buffered the features by their average width plus the estimated GPS positional error of the radiocollars (±7.7 m, C. DeMars, unpublished data), resulting in buffered widths of 40 m for roads, 50 m for pipelines and 30 m for seismic lines. We considered any GPS radiocollar location or random location to be on a LF if it fell within this buffered footprint. Note that our objective was not to model animal response to specific LF types per se (e.g. Dickie et al., 2016), but rather to assess how LFs in general influenced spatial overlap among caribou, wolves and bears. For assessing effects of LF density, we merged data from all line types into one parsimonious layer and calculated LF densities in a moving window analysis using radii of 200 m and 1,000 m. We evaluated radii up to 10,000 m, but these larger measures began to approximate densities measured over the entire caribou range and obscured potential finer scale variation in LF density (Boyce, 2006).

2.6 Evaluating linear feature effects on neonate survival

We used a similar multiscale approach to assess LF effects on the survival of neonate calves. At a fine scale, we evaluated how a female's proportional use of LFs influenced calf survival. At the larger scale, we evaluated how female exposure to varying line density—measured in the same 200-m and 1,000-m radii—affected calf survival. We further evaluated bivariate models representing LF effects at both scales.

2.7 Statistical analyses

We estimated SSFs for each individual caribou, wolf and bear using conditional logistic regression (CLR). To constrain wolf and bear SSFs to actively moving steps, we used breakpoint regression to identify such steps as ≥150-m (Appendix S4). Across all three species, we estimated SSFs only for those individuals having ≥30 moving steps to allow for effective parameter estimation. For each individual, we estimated two forms of SSFs. The first incorporated covariates for land cover and LFs only (hereafter, SSF_Res indicating resource only models). The second form incorporated aspects of the animal's movement process into a mechanistic movement model that simultaneously estimated selection and movement parameters to potentially improve precision of parameter estimates (Avgar et al., 2016). With LFs likely affecting movement rates (Dickie et al., 2016), we incorporated the interaction of step length with LF variables into these models (hereafter, SSF_Res*SL; Table 2).

Variable specification within models differed among species (Table 2). Because caribou used predominantly peatlands during calving, land cover and LFs were specified as binary variables (i.e. peatlands/other; on/off lines). With wolves and bears being more habitat generalists, we incorporated all seven land cover categories into their models by coding six dummy variables and setting upland deciduous as the reference category. We specifically assessed the interaction of LFs with peatlands by further creating three dummy variables representing LF within a peatland, LF within other land cover types and peatland off a LF.

For each SSF, we used generalized estimating equations (GEEs) to calculate robust standard errors, a framework that accounts for potential autocorrelation among steps (Craiu, Duchesne, & Fortin, 2008; Appendix S4). To discriminate between SSF models within individual animals, we used the quasi-likelihood under independence criterion (QIC; Craiu et al., 2008). For population-level inferences on each species, we averaged parameter estimates across individuals, weighting each estimate by the inverse of its variance to give more weight to individuals with more precise estimates. We used bootstrapping to calculate 95% confidence intervals (CIs).

To evaluate model performance, we used k-fold cross-validation for CLR following Fortin et al. (2009). Briefly, this approach estimated SSFs using 75% of randomly selected strata with model outputs then used to generate predictions for observed and random steps within the withheld strata. Observed and random predictions were collectively ranked within each stratum with rankings tallied into bins across strata (nine bins for fine-scale wolf and bear SSFs; 21 bins for large-scale bear SSFs; four bins for caribou SSFs). We calculated a Spearman rank coefficient (r_s) to assess the correlation between bin rank and its associated frequency of observed predictions (r_{S_Obs}). We repeated this process to calculate r_s based on random expectation (r_{S_Ran}), using a randomly selected random prediction within each testing strata but excluding the observed prediction from stratum ranking. We repeated the entire procedure 30 times to calculate and for each SSF with increasing model performance equating to higher relative to . Because caribou SSFs yielded only four possible predictions, we also provide the frequency of observed steps falling within the top ranked category (TR_FreqObs) compared to the frequency of random steps (TR_FreqRan) to further assess model performance.

We estimated RSFs for caribou and wolves using generalized linear mixed-effects models, specifying individual caribou and individual wolf nested within pack as random grouping factors. We explicitly assessed the effect of LF density on peatland selection by assigning these variables—and for wolves, their interaction—as random slopes. This formulation yields variance estimates that appropriately reflect the animal as the sampling unit (Schielzeth & Forstmeier, 2009). Within each species-specific RSF, variables for land cover were specified as per SSFs and LF density variables were standardized to facilitate model convergence. We used Akaike's information criterion (AIC) to discriminate between models with density measured within 200-m and 1,000-m radii. We evaluated RSF performance using k-fold cross-validation (Boyce, Vernier, Nielsen, & Schmiegelow, 2002), randomly partitioning the data into fivefolds, estimating RSFs with fourfolds and then calculating r_s to assess how well model outputs predicted the withheld spatial locations. As with SSFs, we repeated this process 30 times. Note that for LF density, we considered both quadratic relationships and categorical transformations but these did not improve model performance and their results are not shown here.

For all resource selection analyses, we evaluated the robustness of inferences to functional responses, where selection may change as a function of resource availability (Mysterud & Ims, 1998). To do so, we assessed the strength of correlation between individual-level selection coefficients and individual-specific resource availability (Appendix S3). We considered a linear or quadratic relationship with an R² > 0.10 to be indicative of a functional response.

To assess LF effects on calf survival, we estimated mixed-effects Cox proportional hazard models, specifying individual caribou as a random grouping factor to account for females calving in multiple years. Cox models are time-to-event analyses, and in our formulation, the event is calf loss. Positive model coefficients are therefore interpreted as an increasing risk of calf mortality with an increase in the associated covariate. We considered three explanatory variables describing exposure to LFs: proportional LF use, calculated as the proportion of a female's GPS locations falling on a LF, and the mean density of LFs around each female's GPS locations measured within 200 m and 1,000 m radii. Because LF density data were highly skewed, we considered log-transformations of these variables. We fit univariate and bivariate models (proportional LF use and mean LF density at 200 or 1,000-m), and for the latter, we used variance inflation factors (VIFs) to assess for collinearity. We discriminated among models using AIC. For the top model, we evaluated relative fit using a LRT and tested the assumption of proportional hazards by assessing for linearity and a zero slope of the scaled Schoenfeld residuals.

All analyses were performed in r, version 3.3.0 (see Appendix S7 for the specific packages and functions used and associated references).

3.1 Fine-scale response to linear features

Female caribou with neonate calves generally avoided LFs when moving within their calving areas. Twenty-two of the 63 calving seasons had no GPS locations on LFs although these animals were also in areas with lower LF availability (p < .001 from a post hoc t test). Assessing selection across the remaining 41 calving seasons (mean steps/calving season = 195 [SD = 108]), of which five had SSF_Res*SL as the top model, avoidance of LFs was still evident (weighted  = −0.34 [95% CI −0.49, −0.20]). Eight calving seasons, however, had positive coefficients for LFs (Figure 1). Across individuals, SSFs had good predictive performance ( = 0.91 [range 0.30, 1.00],  = 0.91 [0.30, 1.00];  = 0.65 [0.21, 0.93],  = 0.60 [0.21, 0.89]; mean difference between and among caribou = 0.25 [0.5, 0.42]).

For wolves, LFs strongly affected selection of peatlands. We used 21 wolves (mean steps/individual = 792 [SD = 401]) to estimate population-level effects, removing two wolves that had similar movement paths—and therefore parameter estimates—to other wolves within their packs. Nineteen wolves had SSF_Res*SL as their top model. In general, wolves highly selected LFs in peatlands (weighted  = 0.85 [95% CI 0.68, 1.05]; Figure 2) yet avoided peatlands off LFs (weighted  = −0.21 [95% CI −0.38, −0.04]). For 12 wolves, LF in a peatland was the first or second ranked land cover type while peatland off LF was ranked lowest. Wolves also selected LFs in other land covers (weighted  = 0.93 [95% CI 0.79, 1.11]). Overall, model prediction was good across individual wolves ( = 0.95 [range 0.78, 0.99],  = 0.75 [0.55, 0.89], mean difference among wolves = 0.20 [0.04, 0.41).

LFs also influenced selection of peatlands by bears (mean steps/individual = 350 [SD = 201]). Whereas peatlands off LFs were selected at rates similar to the reference category (weighted  = 0.13 [95% CI −0.02, 0.28]), selection was higher for LFs in peatlands (weighted  = 0.51 [95% CI 0.39, 0.68]; Figure 2). LFs in other land covers were also highly selected (weighted  = 0.86 [95% CI 0.66, 1.06]). Two of 16 bears had SSF_Res*SL as their top model and predictive performance was good across all individual models ( = 0.91 [range 0.74, 0.97],  = 0.82 [0.64, 0.95, mean difference among bears = 0.08 [0.01, 0.19]).

3.2 Large-scale response to linear features

The three species had differing responses to LF density. For caribou, response was strongest at the 1,000-m scale (ΔAIC = 3,417 units lower;  = 0.82 [range 0.72, 0.88]) with females generally selecting calving locations with lower LF density (β = −2.00 [95% CI −2.83, −1.17]). Females, however, could not completely avoid exposure to LFs at the 1,000-m scale when calving (mean LF density per caribou year = 3.36 km/km²; range 0.33–12.69). Avoidance was also not consistent across all 63 females as 18 individuals had positive coefficients with 95% CIs that did not overlap zero.

Wolves were less affected by LF density. As with caribou, wolf response was strongest at the 1,000-m scale (ΔAIC = 800 units lower;  = 0.98 [range 0.96, 1.00) but wolves only showed weak avoidance of increasing LF density within peatlands (β = −0.16 [95% CI −0.33, 0.02]) and weak selection for increasing LF density within other land covers (β = 0.20 [95% CI −0.12, 0.51]). Specific to LF density within peatlands, only seven of 23 wolves had positive coefficients with 95% CIs that did not include zero.

LF density had variable effects on bear selection of land cover. Across the 11 individuals evaluated (mean steps/individual = 46 [range 32–65]), six responded most strongly to LF density measured within a 200-m radius and five bears had SSF_Res*SL as their top model. Estimating population-level parameters within both measured radii, bears generally avoided increasing LF density within peatlands (200-m: β = −0.32 [95% CI −0.53, −0.10]; 1,000 m: β = −0.20 [95% CI −0.43, 0.04] but were generally ambivalent to LF density in other land covers (200 m: β = 0.11 [95% CI −0.04, 0.22]; 1,000 m: β = 0.07 [95% CI −0.13, 0.25]. Only one bear selected for increasing LF density (1,000 m radius) within peatlands, having a positive coefficient with a 95% CI that did not overlap zero. Small sample sizes per bear prevented robust evaluation of model performance, and consequently, performance was low and variable across individuals (200 m:  = 0.39,  = −0.21; 1,000 m:  = 0.21,  = −0.09).

Across all resource selection analyses, inferences were generally robust to variation in fix rate and resource availability (Appendix S3). For the former, excluding animals with low fix rates did not change the direction of population-level parameters or substantially affect the overlap of 95% CIs with zero. For the latter, functional responses in selection were not evident for caribou (all R² < 0.10) while only weak negative relationships were noted for wolves at both scales (R² = 0.12–0.26) and for bears at a fine scale (R² = 0.11–0.16). These negative relationships indicate that the strength of selection decreases with increasing LF availability. Given the ubiquity of LFs within caribou range (Appendix S1), these relationships suggest that our results are conservative when inferring that LFs increase selection of peatlands by wolves and bears.

3.3 Linear feature effects on neonate survival

Of the 63 calves followed, 33 survived to 4 weeks old. Neonate survival was best predicted by a bivariate model describing a female's proportional use of LFs and the log₍₂₎ of mean LF density within a 1,000 m radius. This model was 2–3 AIC units lower than other competing models, and variable VIFs were 2.12, a value acceptable for valid inference (Zuur, Ieno, & Elphick, 2010). Model output suggested that the risk of neonate mortality increases as female use of LFs increases (β = 5.55 [95% CI 0.24, 10.85]) and decreases with every twofold increase in LF density (β = −0.43 [95% CI −0.80, −0.06]). The proportional-hazards assumption was generally supported (p ≥ .21 for nonzero linear trend in the scaled Schoenfeld residuals of each variable), and model fit was marginal (p = .08 from LRT), suggesting that, in general, LF effects are weak determinants of neonate survival.

4.1 Conservation implications

Mitigating effects of LFs has become a conservation priority because of their negative impacts on biodiversity (Laurance et al., 2009). Here, we demonstrated how LFs may reduce the efficacy of prey refugia, which suggests that mitigation efforts should be focused on limiting or restoring LFs that contribute to predator–prey spatial overlap. We caution, however, that spatial overlap is not the only mechanism by which LFs can affect predator–prey dynamics. For example, in our study system, LFs outside of peatlands may enhance wolf hunting efficiency of alternative prey, resulting in increased wolf numbers, which may then “spill over” into caribou habitat (Holt, 1984). Our results do not give insight into the relative importance of this mechanism vs. spatial overlap in terms of impact on caribou populations. Effective prioritization of mitigation efforts will therefore require further research to disentangle the importance of such effects, particularly given the substantial costs associated with LF restoration for wide ranging species inhabiting highly impacted forested systems (Schneider, Hauer, Adamowicz, & Boutin, 2010).

We further caution that our results do not give insights into the spatial requirements of prey refugia. Such requirements will likely be context-specific, depending on the life-history traits of both predator and prey. For caribou, previous research has suggested that disturbance effects are best measured at the scale of a herd's range and that these effects are largely negative (Environment Canada, 2012; Sorensen et al., 2008). Given this scale, the functional response of predators in selecting LFs and the sensory capabilities of predators to detect prey (e.g. for wolves, up to 2.5 km; Mech & Peterson, 2003), our results suggest that overall LF density will likely need to be low within caribou range (e.g. <1 km/km²; McCutchen, 2007) and that many peatland complexes will need to be essentially intact. Again, given the wide spatial extent of most caribou ranges, achieving such densities through LF restoration will not be a trivial task.

Finally, we acknowledge that different types of LFs may exert unequal impacts on predator–prey dynamics. For example, wolves have shown relative avoidance of LFs that were narrow (≤7 m wide) and sinuous, perhaps due to decreased line of sight or increased movement costs (Dickie et al., 2016; Wilson et al., 2013), and bears have shown decreased use of LFs ≤2 m wide (Tigner et al., 2014). Determining the relative impacts of specific LF types on predator–prey dynamics will therefore be necessary to better understand and quantify LF effects over large spatial scales.