Landscape of fear or landscape of food? Moose hunting triggers an antipredator response in brown bears

Abstract Hunters can affect the behavior of wildlife by inducing a landscape of fear, selecting individuals with specific traits, or altering resource availability across the landscape. Most research investigating the influence of hunting on wildlife resource selection has focused on target species and less attention has been devoted to nontarget species, such as scavengers that can be both attracted or repelled by hunting activities. We used resource selection functions to identify areas where hunters were most likely to kill moose (Alces alces) in south‐central Sweden during the fall. Then, we used step‐selection functions to determine whether female brown bears (Ursus arctos) selected or avoided these areas and specific resources during the moose hunting season. We found that, during both day and nighttime, female brown bears avoided areas where hunters were more likely to kill moose. We found evidence that resource selection by brown bears varied substantially during the fall and that some behavioral changes were consistent with disturbance associated with moose hunters. Brown bears were more likely to select concealed locations in young (i.e., regenerating) and coniferous forests and areas further away from roads during the moose hunting season. Our results suggest that brown bears react to both spatial and temporal variations in apparent risk during the fall: moose hunters create a landscape of fear and trigger an antipredator response in a large carnivore even if bears are not specifically targeted during the moose hunting season. Such antipredator responses might lead to indirect habitat loss and lower foraging efficiency and the resulting consequences should be considered when planning hunting seasons.


Section S2: Scaling experiment for resource selection in brown bears
As resource selection is scale-dependent in large carnivores (Pitman et al. 2017), we needed to identify the scale at which landcover types should be extracted.To choose the appropriate scale, we conducted a scaling experiment, whereby landcover types were extracted within buffers of various radii [0m (dummy coded), 50m, 100m, 150m, 250m, 350m] centered on the end location of each step.We used the fit_issf function from the amt package to fit conditional logistic regressions with step ID as strata to each bear-year individually (Signer et al. 2019).We used AICc, as described in the main document, to compare the performance of the following model structure: Used (0,1) ~ % Deciduous (5-15m) + % Conifer (5-15m) + % Young forest (˂ 5m, ˃7 years old) + % Bog + % Clearcut (˂ 7 years old) + Distance to Road + Terrain ruggedness + log(step length) + cos(turning angle), with datasets that differed only in the scale (i.e., buffer size) at which landcover types were extracted.We counted the number of bear-years for which each scale explained the most variance and our results suggest that most bears selected resources at the 50m scale (Fig. S2).We consequently used this dataset in further analyses.
The dummy model returned a warning regarding a potential infinite Beta for the Clearcut class in a single individual; however, that model was not the top-ranked within the set and it did not influence the results of the scaling experiment.S2).Model selection was conducted using second order Akaike Information Criterion (AICc) and its derived measures (ΔAICc and Akaike weight) from the AICcmodavg package (Mazerolle 2020).Models within ΔAICc < 2 were considered equivalent (Burnham and Anderson 2002).The model names reflect their underlying hypothesis.All models were fitted using generalized linear models with a binomial family and logit link function.The original Nationella Marktäckedata classes were 51 and 52 for the variable infrastructure.See table S1 for the other variables.
We also modelled resource selection at fine temporal scale in brown bears by using iSSF.To determine the best fixed effect structure, we first fitted conditional logistic regressions with step ID as strata by using the fit_issf function for each demographic group and for both day and night, separately (Signer et al. 2019).Each set contained three models built with a different set of variables and represented competing hypotheses (Table S3).In the Open model, brown bear habitat selection depends only on open habitats and included variables were clearcut, distance to road and open bog.In the second model (Natural model), brown bear habitat selection depends only on forest composition (i.e., young, deciduous, mixed and coniferous forests) and terrain ruggedness.The third model included the variables from both the Natural and Open models.We could not include the relative probability of moose kill in combination with landscape covariates due to lack of independence between this variable and the others.We added the log of step length and the cosine of the turning angle and an interaction between each covariate and the hunting period in all candidate models.Table S3.Structure of candidate models used to estimate habitat selection coefficients for brown bears (n = 53) in south-central Sweden, 2016-2019.
Notes: All models were fitted using conditional logistic regression with step ID as strata.
'/hunting' denotes an interaction between the variables and the hunting periods defined in Fig. 1 from the main manuscript.The best performing RSF model for moose hunter was the Full model, which had a ΔAICc of ≥11.48 to the next best model (Table S4) and was also attributed 100% of the weight within the model set (Table S4).The four best models all included the variables elevation, terrain ruggedness and distance to roads (Table S4).Notes: K = number of parameters, ΔAICc is the difference with the lowest AICc value, w is the model weight and LL is the log-likelihood of the model.

Open
For the bear iSSF, the combination of the natural and open models performed best for all demographic groups during both day and night, and this model structure was attributed 100% of the weight in all model sets (Table S5).Table S6.Coefficients of movement parameters with 95% confidence intervals estimated from integrated step-selection functions for female brown bears with dependent offspring (n = 18 bear-years), solitary females (n = 17 bear-years) and subadult females (n = 18 bear-years) for day and night in south-central Sweden, 2016-2019.The coefficients were estimated for each hunting period: before hunting (red), bear hunt (yellow), low intensity moose hunt (green), pause (blue) and high intensity moose hunt (purple).The coefficients presented in this table were estimated from the bear habitat iSSF models.
Step Following our habitat selection analyses, we decided to investigate the brown bears response to moose hunting.iSSF can also be used to conjointly investigate movement and resource selection (Avgar et al. 2016); however, we could not investigate habitat specific movement in our original iSSF and included the movement parameters as controls (Signer et al. 2019, Fieberg et al. 2021), which is essential because bears modify their activity patterns during the fall.Investigating movement in an iSSF framework requires extracting resources at the start of each step and adding interactions with movement parameters and habitat variables, whereas resource selection requires the extraction of resources at the end of each step (Signer et al. 2019, Fieberg et al. 2021).
In this article, it was not possible to build a model with both habitat-specific movement and resource selection because we split our analyses according to time of day and it becomes an issue for the steps the overlap the onset of legal hunting hours.Thus, the start and end of these steps were recorded in different time periods and including habitat specific movement and resource selection in the same model could have led to bias estimates.We could also not create separate models for each demographic groups and pooled all individuals due to convergence issues; however, this is not an issue because the movement response to a disturbance is similar in all demographic groups (Ordiz et al. 2013).Therefore, we estimated the movement response to the probability of moose kill at the population levels with separate models for day and night during each of the hunting period with the structure described in Table 1 in the main document.
We fitted tentative Gamma and Von mises distributions for step length and turning angle for moving steps only (moving step = hourly displacement of ˃ 40m) with all individuals combined with the fit_distr function [amt package; (Signer et al. 2019)].We used the coefficients for the log of step length as modifiers of the original Gamma shape, whereas the coefficients for the cos of turning angle were used as modifiers of the original Von mises concentration parameters (Fieberg et al. 2021).Where K0 is the original concentration parameter (i.e., kappa), β(cos turning angle), β(cos turning angle:RSFhunt) and β(cos turning angle:RSFhunt 2 ) are the model coefficient for the cos turning angle and the interactions between the cos turning angle and the probability of moose kill (RSFhunt) and its quadratic term (RSFhunt 2 ).
We calculated the expected speed (m/h) of female brown bears when traveling through low and high probability of moose kill (10% and 90% quantiles) by multiplying the scale and shape of updated gamma distributions for day and night during each hunting period (Appendix S1: Table S7).The updated Von mises distributions were generated using the dvonmises function from the circular package (Agostinelli and Lund 2017).We did not need to generate updated Gamma distributions of step length, since we could obtain the mean expected hourly displacement by multiplying the updated shapes with the original scales.
We similarly generated probability densities of turning angles when traveling through low and high probability of moose kill (10% and 90% quantiles) by using the updated Von mises distributions for day and night during each hunting periods (Appendix S1: Figure S5).when traveling through areas with low (10% quantile) and high (90% quantile) relative

Figure S1 .
Figure S1.Map of the moose management units in Sweden.The rectangle represents the study

Figure S2 .
Figure S2.The number of individuals (n = 53 bear-years, in south-central Sweden, during 2016-

Figure S3 .
Figure S3.Coefficient estimates for a) the selection of mixed forest, b) the selection of

Figure S4 .
Figure S4.Coefficient estimates for a) the log step length and b) the cosine of turning angle with

Table S1 .
Variables included in integrated step-selection functions of brown bears (n = 53) during the fall in south-central Sweden, 2016-2019.
Notes:The mean ± standard deviation (SD) represent raw landscape features extracted from buffers (radius = 50m) centered at the end of used and available steps.Step lengths and turning angles are log-and cosine-transformed, respectively.Original Nationella marktäckedata classes were 118 and 128 for young forest, 111 to 113 for coniferous forest, 114 for mixed forest, 115 to 117 for deciduous forest and 2 for open bog.

Table S2 .
Structure of candidate models used to estimate habitat selection coefficients for moose hunters in south-centralSweden, 2016Sweden,  -2019.   .

Table S4 .
Model selection by Akaike Information Criterion (AICc) for candidate generalized linear models used to estimate habitat selection coefficients for moose hunters in south-central Sweden, 2016-2019.

Table S5 .
Model selection by Akaike Information Criterion (AICc) for integrated step-selection models in brown bears during day and night in south-central Sweden, 2016-2019.ΔAICc is the difference with the lowest AICc value, w is the model weight and LL is the log-likelihood.