The density of anthropogenic features explains seasonal and behaviour-based functional responses in selection of linear features by a social predator

Anthropogenic linear features facilitate access and travel efficiency for predators, and can influence predator distribution and encounter rates with prey. We used GPS collar data from eight wolf packs and characteristics of seismic lines to investigate whether ease-of-travel or access to areas presumed to be preferred by prey best explained seasonal selection patterns of wolves near seismic lines, and whether the density of anthropogenic features led to functional responses in habitat selection. At a broad scale, wolves showed evidence of habitat-driven functional responses by exhibiting greater selection for areas near low-vegetation height seismic lines in areas with low densities of anthropogenic features. We highlight the importance of considering landscape heterogeneity and habitat characteristics, and the functional response in habitat selection when investigating seasonal behaviour-based selection patterns. Our results support behaviour in line with search for primary prey during summer and fall, and ease-of-travel during spring, while patterns of selection during winter aligned best with ease-of-travel for the less-industrialized foothills landscape, and with search for primary prey in the more-industrialized boreal landscape. These results highlight that time-sensitive restoration actions on anthropogenic features can affect the probability of overlap between predators and threatened prey within different landscapes.

. Capture and handling protocols are described previously 37,38 , and were approved and carried out in accordance to university animal care protocols (University of Montana Animal Use Protocol 059-09MHWB-122209; University of Calgary Animal Use Protocol BI11R-17; University of Alberta Animal Care Committee Standards 99-69). Wolves were fitted with Lotek 2200/3300 (Lotek Engineering Systems, Newmarket, Ontario, Canada). GPS collars were originally programmed to acquire locations at a range of intervals between 15-min and 2-h, and we rarefied these locations to 2-h intervals to reduce autocorrelation and obtain more uniform sample sizes among individuals. We also partitioned wolf GPS www.nature.com/scientificreports/ data into "resting-feeding" and "travelling" locations because selection patterns likely differs between hunting, travelling, or searching behaviour and resting-feeding behaviour 11 . We defined resting-feeding locations as any location during which wolves spent at least six hours within a 300 m radius 39 . We also divided each year into three seasons (denning, rendezvous, and nomadic) based on wolf behaviour to account for variations in seasonal selection associated with life history requirements 38,40,41 . For analyses, we discarded individuals with fewer than 20 locations per landscape-season-behaviour classes per year, and accounted for within-pack correlation by removing one of every two GPS locations acquired from different individuals of the same pack when these locations occurred < 200 m apart during the same time interval 10 . The total number of GPS collar locations was 6,243 (less-industrialized foothills landscape: 3,731 locations from 8 individual-year; more-industrialized boreal landscape: 2,512 from 10 individual-year, Appendix S1).
Environmental variables. Using a 25-m digital elevation model, we derived topographic variables including slope (variable names in italics: Slope), elevation (Elev), topographic position index (TPI 42 ), and compound topographic index (CTI; terrain wetness 43 ). Predominant winds are from the south-west in this region and we therefore separated aspect into three binary variables (fFlat = 0°; fLee = from NW to E aspect; and fWind = from SE to W aspect, i.e., fLee represents a categorical variable describing the pixel as either NW-to E-facing (1), or not (0)). We also derived yearly landcover and percent canopy cover (%CC) variables from a combination of Moderate Resolution Imaging Spectroradiometer (MODIS) and Landsat imagery mapped at a 30-m resolution 44,45 . We grouped landcover classes into three categories (fMixed, fConifer, and fNon-forest), and recoded these categories into binary variables. We extracted timber harvest block locations and age from Alberta Vegetation Inventory (AVI) data provided by Forest Management Agreement (FMA) holders within the study area (Canadian Forest Products, West Fraser Mills Ltd., and Weyerhaeuser Co. Ltd.) to delineate harvest blocks < 25 years old. We then combined timber harvest data with road, seismic line, well site, and pipeline data provided by the Government of Alberta (GoA), FMA holders, and the Alberta Energy Regulator to calculate the density of anthropogenic features for each year of animal data (2003-2009) using a 70-m (A70-local) and 1-km (A1k-landscape) radii circular moving window average in ArcGIS 10.2 46 . We chose a 70-m radius to represent a local scale based on findings from DeCesare et al. 37 , and a 1-km radius as a conservative estimate of the influence of disturbance features at the landscape scale 20,47 . We also created a binary variable representing early successional forests (E. Seral) by merging vegetated, but non-forested landcover classes with AVI harvest blocks data < 25 year old. To represent the non-linear diminishing effect of small and large streams with increasing distances, we used an exponential decay function (1 − exp (−0.002 × distance (m)) ) as described by Nielsen et al. 48 (large streams: DistW1m and small streams: DistW20k). Finally, we also generated a raster surface representing the distance to the main highway intersecting the more-industrialized boreal landscape (DistHWY40), and a raster surface representing the distance to the nearest seismic line (Dist).

Seismic line variables.
We used point cloud LiDAR data collected between 2003 and 2008 to extract the maximum vegetation heights, at a 1-m resolution, for 100 m sections of seismic lines spanning the study area (VegHT). Details of LiDAR-based extractions of vegetation heights on seismic lines are described previously 10 . We derived seismic line wetness (WAM) under each 100 m section of seismic line from the average depth-towater values extracted from wet areas mapping 49 , and used an exponential decay function (1 − exp −1.55×WAM(m) ) to rapidly decrease the effect of depth to water at depths > 2 m, and to set values > 3 m as constant because the mean root depth of boreal forest vegetation is 2 ± 0.3 m 50 . Finally, we used Geospatial Modeling Environment (GME) 51 to determine the landcover category that intersected each 100 m segment of seismic line. When seismic line segments fell within two or more landcover types, we used the majority landcover type along the seismic line section.

Statistical analyses.
Our goal was to understand how wolves, on average, respond to attributes of seismic lines, and whether the density of anthropogenic features surrounding seismic lines leads to scale-specific functional responses in habitat selection. For this, we developed within-home range (3rd order) resource selection functions (RSFs) for each individual and used the inverse of the variance associated with each coefficient to calculate weighted averages-i.e., we calculated population averages from individual-based models giving more weight to coefficients derived from animals with more precise estimates 52 . Before generating individual-based models, we first built population-level baseline generalized linear models (GLMs) for 'resting-feeding' and 'travelling' locations in each landscape. We generated these baseline GLMs as a tool to identify environmental variables that should be included in all individual-based models to avoid the potential for unreliable interpretations associated with averaging across models built with different sets of baseline variables. Using this approach, we included non-informative variables in individual models but performed model-averaging on a consistent set of variables for each individual within each landscape-season-behaviour class. To optimize model fit, we only retained variables that were influential at the population-level (coefficients from weighted averages that do not overlap zero) for each landscape-behaviour dataset. Collinearity and correlation between variables differed per season and resulted in different baseline models per season within each landscape (Appendix S2). Using these baseline models (M1) and null models (M0) as starting points, we generated 13 population-level GLMs based on a priori consideration of multiple working hypotheses to consider as candidate models for our final individualbased models 53 (Table 1). We used interaction terms between densities of anthropogenic features at the local (A70) and landscape (A1k) scales (i.e., measures of changing habitat availability) and the distance to the nearest seismic lines (Dist) to assess the potential for functional responses in habitat selection driven by changes in availability of anthropogenic habitat. We expected the influence of anthropogenic features on habitat selection patterns to decrease rapidly and non-linearly with increasing distance between individuals and targeted anthro- www.nature.com/scientificreports/ pogenic features 54 . We therefore evaluated two decay functions for the distance-to-the nearest seismic line variable (Dist) derived from methods described in Nielsen et al. 48 (1 − exp (−0.002 × distance (m)) , and 1 − exp ), compared univariate models of the Dist variable without any decay function to these two decay functions using AIC, and selected the Dist variable with the lowest AIC per dataset for subsequent analyses (Appendix S3).
There is currently no consensus on how to best approach model selection using individual-based models to infer population-level behaviours, however, model selection is straightforward when using a population-level approach 53 . Parameters of a variable of interest can only be estimated from averaging individual-based models if each individual is exposed to a range of environmental conditions associated with that variable, and each individual is exposed to the same variables themselves 55 . When using individual-based models: (1) there is currently no clear consensus on how to average across models with different sets of variables, which can yield ambiguous results, and (2) eliminating individuals that lack exposure to the full set of variables of interest inflates the influence of individuals exposed to the full set of variables when identifying the 'top models' 56 . Moreover, the choice of adopting individual-or population-level inference depends on specific ecological questions, which in our study were inferences at the population-level 55,57 . We therefore chose to perform model selection using population-level models, but used the inverse-weighted coefficients 52 of individual-based models for populationlevel inferences (sensu 58,59 ). Appendix S4 further explains our rationale, and Supplementary Table S5 compares the results of model selection performed on population-level models to results from model selection performed on individual-level models.
We derived population averages for animals across each landscape-season-behaviour classes (two landscapes: less-industrialized foothills landscape vs. more-industrialized boreal landscape, three seasons: denning, rendezvous, and nomadic, and two behaviours: resting-feeding, and travelling, for a total of 12 distinct datasets) following a 'design III' use-availability approach 53,60 . We used GME 51 and ArcGIS 10.2.2 46 to generate 20 random 'available' point locations for every used GPS collar location per animal-year-season minimum convex polygon (MCP). Fix rates from GPS collars were < 80% and the number of GPS locations per individual in each dataset varied largely due to missed GPS fixes that over-or under-represented certain individuals (Appendix S1). We therefore accounted for unequal probabilities of obtaining successful GPS collar locations (P fix ) with changing terrain and habitat by using the weight argument of the svyglm function within the survey package in R to specify observation-specific inverse probability weights for used locations of each individual and held the value of P fix for available locations constant within general linear models [61][62][63][64] . We restricted available locations to at least 30 m from one another; the size of the raster pixel used for analyses. We excluded one of two correlated variables (r ≥ 0.5) from any model and also excluded variables with Variance Inflation Factors (VIFs) > 3. We standardized all continuous variables to improve model convergence and used a population-level informationtheoretic (IT) approach with Akaike's Information Criterion (AIC) to assess these variables 53,55 . We carried out data exploration and statistical analyses in R 64,65 . We evaluated model performance using leave-one-out (LOO) cross sample validation; iteratively refitting models to subsets of the data following Matthiopoulos et al. 66 , where mean correlations indicate average model performances, with values closer to 1 (range from 0 to 1) indicating stronger fit. We present model results as beta coefficients (β) ± 95% Confidence Intervals (CI) unless otherwise noted, and to address model selection uncertainty, we report influential interactions for all models with weights of evidence (ω i ) ≥ 0.1 53 . Table 1. Candidate models and associated working hypotheses proposed to explain seasonal (denning, rendezvous, and nomadic) selection of areas near regenerating seismic lines for wolves travelling and restingfeeding in a more-industrialized boreal and less-industrialized foothills landscape of west-central Alberta, Canada between 2003 and 2009. VegHT vegetation, WAM wet areas, E.Seral early successional forests, Elev elevation, Dist distance, densities of anthropogenic features at the local (A70) and landscape (A1k) scales, and topographic position index (TPI). Variables are fully described in "Materials and methods" section. 'Base' refers to the suite of variables included in the respective landscape-season-behaviour baseline models for each dataset (See Appendix S2 for details of baseline models).

Results
Baseline models. During all seasons, wolves in both landscapes showed similar baseline selection pat-  Table S10 Leave-one-out model validation indicated moderate to high model fit with mean correlations for individuals ranging from 0.4 to 0.7; Appendix S5).  Table S7). During the rendezvous season, wolves selected areas near relatively wet and low vegetation height seismic lines more than dry, low vegetation height seismic lines, and selected areas near high vegetation seismic lines, regardless of seismic line wetness (WAM*Dist*VegHT = -0.3 ± 0.07, Fig. 2B, Appendix S5 Table S7). Finally, during the nomadic season, we also observed evidence of a functional response in habitat selection driven by the density of anthropogenic features. Wolves selected areas near relatively low vegetation height seismic lines more within areas of low density of anthropogenic features at the landscape scale, and decreased selection for areas near low vegetation height seismic lines more with increasing density of anthropogenic features (A1k*Dist*VegHT = 0.2 ± 0.07, Fig. 4B, Appendix S5 Table S7). Based on the second best model (M6: Ease-of-Travel & Elevation, ω i = 0.4), wolves also selected areas near relatively low vegetation height seismic lines compared to high vegetation height seismic lines, and especially selected areas near low and moderate vegetation height seismic lines at low elevations (Elev × VegHT × Dist = 0.2 ± 0.04, Appendix S7 Table S10). Leave-one-out model validation indicated low to moderate model fit with mean correlations for individuals ranging from 0.2 to 0.5; Appendix S5). During the denning season, wolves selected areas near seismic lines at low elevation (Elev*Dist = 0.3 ± 0.08, Appendix S1 Table S3). During the rendezvous season, wolves selected areas near relatively wet seismic lines more than dry seismic lines (WAM*Dist = 0.1 ± 0.08, Appendix S5 Table S8). Based on the next best models (M6: Ease-of-Travel & Elevation, ω i = 0.3), these wolves also selected areas near relatively low vegetation height seismic lines compared to high vegetation height seismic lines, and especially selected areas near low and moderate vegetation height seismic lines at low elevations (Elev*VegHT*Dist = 0.1 ± 0.08, Appendix S7 Table S10). During the nomadic season, wolves generally selected areas near relatively dry seismic lines compared to wet seismic lines (WAM*Dist = − 0.1 ± 0.02, Appendix S1 Table S3), and this result was also supported by the next best model (M10: Wet areas & early seral forest for the rendezvous season, ω i = 0.3; WAM*Dist = − 0.1 ± 0.03, Appendix S7 Table S10). We could not assess model fit for the denning season because of low sample size but leave-one-out  Table S9). However, based  Table S10). During the rendezvous season, wolves selected areas near wet seismic lines more than dry seismic lines (WAM*Dist = 0.2 ± 0.09, Appendix S5 Table S9). During the nomadic season, wolves selected areas near wet seismic lines in early successional forests more than areas near dry seismic lines, and also selected areas near wet seismic lines in early successional forests more than areas near seismic lines in mature forests regardless of seismic line wetness (WAM*Dist*E.Seral = 0.7 ± 0.5, Appendix S5 Table S9). Finally, based on the next best model (M4: Distance & local functional response, ω i = 0.2), these wolves also selected areas near seismic lines when in areas with high densities of anthropogenic features at the local scale (A70*Dist = − 0.3 ± 0.2, Appendix S7 Table S10). Leave-one-out model validation indicated moderate model fit with mean correlations for individuals ranging from 0.4 to 0.5; Appendix S5).

Discussion
Using GPS collar data spanning 7 years and collected from eight resident wolf packs in west-central Alberta, Canada, we demonstrated evidence of consistent landscape-scale functional response in habitat selection driven by industrial activity, where wolves were more likely to rest or feed in areas near relatively low vegetation height seismic lines within areas of low densities of anthropogenic features. Overall, we confirmed seasonal and landscape-based behavioural differences that (a) support evidence that seismic lines facilitate ease-of-travel 9,10 , especially in landscapes with comparatively low densities of anthropogenic features, and during the denning season when adults travel long distances to hunt, patrol the territory, and return to dens, and (b) that seismic lines potentially facilitate search for prey species while pup rearing during the snow-free rendezvous season when wolves typically select wet meadows 10,11,31 , especially within landscapes with comparatively high densities of anthropogenic features. Overall, our results suggest that effects of industrial activity drive a functional response in habitat selection for wolves selecting areas near seismic lines. Although seasonal patterns of selection were generally consistent for wolves across west-central Alberta, Canada, landscape heterogeneity highlights the usefulness of seismic lines for travel within landscapes with low densities of anthropogenic features, and the potential benefits of seismic lines towards search for primary prey within landscapes with comparatively high densities of anthropogenic features.
Our study is unique because although the cumulative effects of industrial activity have been investigated previously for wolves (e.g. 23,24 ), our study is the first to specifically observe evidence of (1) behaviour-specific (i.e., resting-feeding vs. travelling) responses to variation in vegetation height on seismic lines, and (2) yearround landscape-scale functional response in habitat selection in response to variation in vegetation height on seismic lines. In two distinct landscapes with largely dissimilar densities of anthropogenic features, we found that wolves consistently selected for seismic lines with relatively low vegetation height when in areas of low densities of anthropogenic features. These results support evidence that seismic lines covered with lower vegetation are likely more useful to wolves for travel and search for prey compared to seismic lines with high vegetation 9,10,19 , www.nature.com/scientificreports/ and especially in areas where fewer of these linear features are available on the landscape. Although seismic lines in both landscapes are utilized by humans using motorized Off-Highway Vehicles (OHV 15 ), human activity on seismic lines is likely not a major deterrent for wolves using seismic lines in either landscapes because unlike frequent vehicular traffic observed on roads 11,24 , the frequency of OHV traffic on seismic lines is minimal. We therefore propose that the functional response in habitat selection observed here is likely a result of seismic lines enhancing ease-of-travel or potentially, search for prey, and that with increasing densities of anthropogenic features, the enhanced value of individual linear features to wolves decreases, and therefore dilutes the overall wolf response, especially within areas where seismic lines are pervasive. www.nature.com/scientificreports/ In accordance with previous studies, e.g. 24,67 , we observed evidence of functional response in habitat selection driven by densities of anthropogenic features primarily during the snow-free seasons, although wolves in the more-industrialized boreal landscape were also more likely to select areas near relatively low vegetation height seismic lines during the nomadic season. Although we were unable to evaluate the effect of snow depth and snow compaction on wolf behaviour in our study, Droghini and Boutin 68 recently observed that snow compaction from motorized vehicle use of linear features reduced wolf sinking depth, and therefore likely favors travel efficiency. Results from Droghini and Boutin 68 therefore offer a potential explanation for the selection of areas near relatively low vegetation heights seismic lines during winter in the more-industrialized boreal landscape, and for the generally low model fit that we observed with the nomadic season in this landscape.
Within the less-industrialized foothills landscape, wolf selection patterns relative to seismic lines were generally in line with previous evidence that seismic lines facilitate travel, potentially because the relatively low densities of anthropogenic features in this landscape amplifies the usefulness of individual linear features. Whittington et al. 11 observed that wolves increasingly selected for linear features at high elevation, and here, we observed that wolves travelling in the less-industrialized foothills landscape increasingly selected for relatively low vegetation height seismic lines when travelling at low elevation. It is likely that relatively low vegetation seismic lines at low elevation and in more rugged terrain below treeline facilitate wolf travel compared to seismic lines with high vegetation or compared to seismic lines above treeline where vegetation is generally sparse. Wolves resting or feeding in the less-industrialized foothills landscape generally selected dry seismic lines, and these results are consistent with selection of seismic lines to improve ease-of-travel when prey are more diffuse, and when adults periodically return to dens or central locations 28,40 . Looking at broad-scale travelling locations within west-central Alberta, Finnegan et al. 10 also observed that wolves moved towards seismic lines with low vegetation heights more during the denning and rendezvous seasons, and especially towards low vegetation seismic lines in non-mature forests compared to higher vegetation seismic lines in mixed and conifer forests. Finnegan et al. 10 also observed that wolves stepped towards high vegetation height seismic lines during the denning season, and our analysis investigating functional responses in selection of seismic lines while accounting for landscape heterogeneity, further explains the findings of Finnegan et al. 10 . Specifically, we observed higher selection for low vegetation height seismic lines within landscapes with low and moderate densities of anthropogenic features, and increased selection for areas near high vegetation height seismic lines with increasing densities of anthropogenic features. We therefore highlight the need to consider landscape heterogeneity and functional response when investigating animal behaviour. In addition, consistent with patterns observed in the more-industrialized boreal landscape and habitat favoring high-quality browse species preferred by primary prey 2,10,69 , we observed that wolves travelling within the less-industrialized landscape during the rendezvous season also selected areas near wet seismic lines in early successional forests, indicating that although seismic lines are likely used to improve ease-of-travel, wolves in the less-industrialized landscape may also be targeting seismic lines in areas potentially preferred by moose, deer, and elk while pups aren't able to travel long distances.
Within the more-industrialized boreal landscape, wolf selection patterns relative to seismic lines generally supported behaviour more in line with searching for primary prey, potentially because high densities of anthropogenic features available in this landscape dilutes the usefulness of linear features for ease-of-travel. Consistent with habitats likely preferred by primary prey 35,36 , and with previous observations that wolves may opt to increase hiding cover and prey availability for pups during the rendezvous season 31 , wolves resting, feeding, or travelling in the more-industrialized boreal landscape during the rendezvous season selected relatively wet and low vegetation height seismic lines, and also selected relatively high vegetation height seismic lines regardless of seismic line wetness. Interestingly, selection patterns during the denning season, when adults may travel long distances to find prey while having to return to dens to feed pups, were more consistent with observations from the less-industrialized landscape where wolf selection supported ease-of-travel. It is plausible that because the more-industrialized boreal landscape has a high overall density of anthropogenic features, the usefulness of seismic lines to improve ease-of-travel is only apparent when travel efficiency is most needed, such as during the denning season. During the nomadic season, selection patterns for wolves resting or feeding were also more consistent with observations from the less-industrialized landscapes where wolf selection supported ease-oftravel. However, wolves travelling in the more-industrialized boreal landscape during the nomadic season were again consistent with areas associated with high-quality browse species for primary prey as we observed selection for areas near relatively wet seismic lines in early successional forests, and no selection for areas near seismic lines within mature stands 29,36 .
Our results highlight the importance of considering landscape heterogeneity and habitat characteristics, effects of industrial activity, and the functional response in habitat selection when investigating seasonal behaviour-based selection patterns. This is especially true in landscapes where restoration activity could influence predator-prey dynamics that have implications for management and recovery action affecting threatened species. We propose that by considering the spatial composition and arrangement of landscape features (i.e., spatial heterogeneity), when developing time-sensitive restoration actions, land managers have the ability to alter predator-prey dynamics. For example, to potentially avoid heightening wolf selection for linear features that remain unrestored within an otherwise restored area, land managers would need to consider the density of remaining anthropogenic features, the vegetation height on remaining linear features, and the potential for additional anthropogenic features with future resource extraction. By ensuring that restored areas do not confine and redirect wolves onto remaining attractive linear features that could improve travel efficiency for wolves (i.e., low vegetation seismic lines within areas of low densities of anthropogenic features at the landscape scale), land managers could potentially reduce predation risk for threatened prey species.
The results of this research offer unique insights into the effects of specific management actions on the functional response in habitat selection, and ultimately, on predator-prey dynamics in relation to landscape-specific spatial heterogeneity. Our results can help managers prioritize actions aimed towards the recovery of threatened Scientific RepoRtS | (2020) 10:11437 | https://doi.org/10.1038/s41598-020-68151-7 www.nature.com/scientificreports/ species such as caribou by ensuring that restored areas do not heighten wolf selection of remaining and attractive linear features such as low vegetation seismic lines in areas with low densities of anthropogenic features.
Restoring anthropogenic features that are most selected by predators, and that therefore likely facilitate hunting efficiency 11,70 , can decrease the potential overlap between predator and prey, but may also increase the selection of remaining seismic lines by predators (i.e., increase the strength of habitat-driven functional responses for predators on remaining and attractive anthropogenic features), therefore potentially heightening predation risk in unrestored areas. With the goal of effective landscape restoration in mind, we propose the need for management actions that consider the potential ripple effects of specific human activity (restoration or continued industrial activity) on predator-prey dynamics.

Data availability
The datasets generated during and/or analysed during the current study are available in the Dryad repository (https ://doi.org/10.5061/dryad .djh9w 0vxd).