A spatial theory for characterizing predator-multiprey interactions in heterogeneouslandscapes
by Daniel Fortin, Pietro-Luciano Buono, Oswald Schmitz, Nicolas Courbin, ChrystelLosier, Martin-HuguesSt-Laurent, Pierre Drapeau, Sandra Heppell,Claude Dussault, Vincent Brodeur, Julien Mainguy
SUPPORTING INFORMATION
Appendix S1.Detailed description of how the spatially explicit, mechanistic functional response was quantified
Parameters NC, NM, hM, hC, aandC
Values for the parameters NC, NM, hM, hC, aandCwere obtained from previous field studies that monitored the movement and interactions of wolves, moose and caribou in a northern boreal forest landscape. Overall density was set to NC = 0.020 individual/km2[1]for caribou, and NM = 0.045 individual/km2for moose [2]. While these densities are representative of the study area as a whole, local densities vary according to the response of each species to anthropogenic features. We also considered that hM= 0.35 (100 days/individual) and hC = 0.11 (100 days/individual), and ai = 17.84 (km2/100 days) [3]. We determined the relative use of caribou by wolves (C) in systems where caribou and moose co-occur but vary in density, based on the information provided in Table S2 of[3] . We calculated that the relative consumption of caribou by wolves covaried with caribou (NC) and moose (NM) density as: C= exp(Z)/(1+exp[Z]), where Z = 6.96(NC / [NM + NC]) – 4.11); R² = 0.37, F1,11 = 5.79, P = 0.04. This equation indicatesthat when both prey species are equally available to wolves,(i.e., NM = NC), for example, caribou should comprise 35% of the wolves’ diet. In this example, moose are more vulnerable than caribou.
ExpressionsvC(s)andwCi
The probability density function vC(s)was taken from Fig. 1 of[4], a trend that closely matches the spatial distribution of caribou over 162,920 km2 of boreal forest. Estimates of habitat selection coefficients for caribou were also obtained from [4], which is based on a comparison of habitat features at the locations of track networks and at locations randomly drawn within their 162,920-km2 study area in the boreal forest of Québec, Canada. The parameter wCi(s) was the exponentiated value of the selection coefficient for land cover iprovided in Table S2 in [4].
ExpressionsvM(s)andwMi
We estimated vM(s)byparameterizing the advection-diffusion model (Eq. 2) using field observations collected in the Côte-Nord region. The model was based on 0 and b(Eq. 2). In addition, solving the PDE requires quantifying, the distance at which the repulsive effect of the edge of the disturbance vanishes, and , the distance from the edge of the disturbance where the advection term starts decaying to 0 [4]. Data for this parameterization come from monitoring 15 moose with GPS collars or ARGOS/GPS collars scheduled to record a location every 4 hours during ca. 20 months between 2005 to 2009. According to [5], moose tended to move away from the nearest road (i.e., cosine of their orientation with respect to road [cos(ϴRoad)] < 0) up to 0.4km, at which point, they simply move randomly with respect to roads(see Fig. 6.3D of [5]). There was no distance-related taxis detected in response to cutovers. Because [5]did not report anystatistical model, we reanalysed the data with a piecewise regression,and detected a plateau at 0.37 km (95% confidence interval: 0.11 – 0.62 km) from the nearest road. The model follows: cos(ϴRoad) = 0.20 + 0.57sRoad, when sRoad0.37 km; otherwise cos(ϴRoad)= 0.0064 + 0.00013sRoad. Given that Φcorresponds to the distance at which the repulsive effect of the boundary vanishes,we set= 0.37 km, a reasonable approximation because of the strong link between distance to nearest road and to nearest cut. Given that radio-collared moose were located on average at 4 km from the nearest cut or road, we set = 4. Because the PDE should not be sensitive to [4]. The remaining parameters were set to 0 = 0.04 and b= 0.183[see 4]. The numerical computations to solve equation (2) of the main document for the and cases were performed using MEF++, a finite element method research software created, developed and maintained by the GroupeInterdisciplinaire de Recherche enÉléments Finis (GIREF) of the Université Laval in Québec, Canada (Public Communications: Accessed 25 June 2015). Details of the computational methods used can be found in the Supplementary Material of [4] with the time-step modified to Δt = 0.025. According to the PDE’s numerical solution,the moose density (density distribution function,) should gradually increase from the nearest road up to ca. 0.48 km, at which point density should level off (Fig. 3). We tested this prediction using independent data of moose spatial occurrences, i.e., with the data used for habitat selection analysis.
We estimatedhabitat selection coefficients for moose during winter in the Côte-Nord region (Québec, Canada) from snow track surveys. The Côte-Nord study area (50o00’ N to 51o55’ N, 67o40’W to 70o40’ W) is characterized by a large proportion of old-growth conifer forest dominated by black spruce (Piceamariana), balsam fir (Abiesbalsamea), and jack pine (Pinusbanksiana). Moose track networks were detected during aerial surveys carried out in late January 2006 from a fixed-wing aircraft by flying along equidistant north-south transects spaced 0.5 km apart in 72 60-km2 sample plots [see 6 for details]. Plots were spread over 21737 km2. Similar to [4]for caribou, we compared habitat features at the centroid locations of the 191 moose track networks with those at 720 random locations (i.e., 10 random points in each plot). For all locations, we determined the distance from the nearest cut and road, as well as the land cover type in which they occur: conifer forest (dominated by open or closed conifer stands with mosses), deciduous or mixed forest, burnt area, lake, lichen-heath community (dominated by lichens with 10–40% of conifers or by shrubs and lichens), or other type (dominated by bare ground, moss and shrubs, agricultural lands, or unclassified areas). The geographic information system (GIS) was based on a Landsat Thematic Mapper image (25-m resolution) taken in 2000, with clearcuts, roads and fires updated every year based on information provided by the local forest companies.We evaluated factors influencing the probability of moose occurrence using a binomial, semiparametric generalized additive model, contrasting centroid locations of track network (coded as 1) and random locations (coded as 0). Land cover types were independent parametric variables, whereas distance to the nearest road was an independent nonparametric variable fitted with a LOESS smoother. This analysis was restricted to observations < 15 km because most edge-effects of roads vanished farther than 5 km from the disturbance, and because we encountered model convergence issues when considering the entire data set.
We found that moose selected burned areas and mixed-wood forest stands, relative to cutblocks (Table S1). Moose did not respond to other land cover types differently than to cutblocks. After controlling for those habitat features, the relative probability of moose occurrence was lowest directly at the edge of the nearest road, and quickly reached a plateau a short distance from the edge (Fig. 3). The LOESS smootheridentified a local peak at 0.39 km, which is close to the plateau expected at 0.48 km from the steady-state solution of the PDE (Fig. 3), but the observed distance at which density reached the plateau was 2.6 km (Fig. 3). Consistently, a piecewise logistic regression identified a plateau at 2.98 km (95% confidence interval: 1.13 – 4.84 km) from the nearest road: logit(sRoad) = 1.72 + 0.31sRoad, when sRoad2.98 km; otherwise logit(sRoad) = -0.68 + 0.04sRoad. The 95% confidence intervals of all parameter values excluded 0, with the exception of the -0.04 coefficient. This indicates that the maximum probability corresponded to a plateau.
Table S1. Semiparametric generalized additive model of habitat selection by moose in winter. The model is based on a comparison between the characteristics of 720 random locations and 191 track networks of moose located from aerial surveys in the Côte-Nord region, Québec, Canada. Selection coefficients (s) are presented with standard errors (SE) and associated P-values.
Covariate* / / SE / t / PWater body / -0.765 / 0.588 / -1.30 / 0.19
Lichen-heathcommunity / -0.141 / 0.860 / -0.16 / 0.87
Conifer forest / 0.177 / 0.433 / 0.41 / 0.68
Peatland / 0.551 / 1.194 / 0.46 / 0.64
Burn / 0.928 / 0.482 / 1.92 / 0.05
Mixed forest / 0.894 / 0.424 / 2.11 / 0.03
Road / 1.261 / 0.768 / 1.64 / 0.10
Analysis of deviance Source / Df / Sum of squares / χ2 / P
Loess (SRoad) / 4 / 16.04 / 16.04 / 0.003
*Cuts (and other) were considered as the reference category
ParameterwPi
Nine wolves (3 males and 6 females) from four packs were monitored from 2005 and 2010 in winter (15 November – 31 March) in the Côte-Nord region. Wolves were captured with a net gun fired from a helicopter, and fitted with a Global Positioning System (GPS) collar (Lotek Engineering Inc., Newmarket, Ontario) scheduled to locate individuals every 4 h, for five days a week and every hour the rest of the week, or with an Argos/GPS collar (Telonics Inc., Mesa, Arizona) relocating individuals every 10 h (Argos/GPS). Each wolf was followed for 1.5±0.4 winters.We characterized all GPS locations as well as an equal number of locations randomly distributed within the 100% minimum convex polygon of the observed observations. This area exceeds the individual home ranges; hence, our evaluation of habitat selection is intermediate to Johnson’s [7] second and third order of selection. We used the same GIS representation of the landscape for wolves as for moose (see description above).
We built resource selection functions (RSFs) by contrasting habitat characteristics at observed and random locations. We estimated RSF parameters using a generalised linear mixed model (GLMM) with a binomial distribution, with the maximum likelihood being estimated using an adaptive Gaussian quadrature (4 quadrature points) procedure [8]. We applied three-level, mixed-effects RSF models to accommodate the hierarchical structure of wolves within packs with random intercepts at both the pack and individual levels [9]. Logging roads were a critical determinant of edge in the study area because they were closely associated with cutblocks, with a strong correlation between the distance to the nearest road and the nearest cut (r = 0.98, n = 77,860, P < 0.0001). To avoid multicollinearity issues, we thus considered distance either to the nearest cut or road (s).
We were specifically interested in understanding whether or not wolf selection of mixed and deciduous stands (the key land cover type for moose), lichen-heath communities (the key land cover type for caribou), and conifer stands varied as a function of the distance to human-disturbed sites. This was achieved by including interaction terms in the GLMM. We used an information-theoretic approach to identify the relevant independent variables, and to identify the most appropriate functional form of the relationship between the logit and the various interaction terms. As for moose and caribou, wPi(s)was the exponentiated value of selection coefficient of land cover i at distance s.
Radio-collared wolves selected mixed and deciduous forests most strongly (Fig. S1; Table S2). Near cuts and roads, they also focused their search on lichen-heath communities (Fig. S1), the vegetation type most strongly selected by caribou. The strength of wolf selection for lichen-heath communities declined, however, with distance to human disturbance. Hence, the probability of wolf occurrence was lowest at the farthest distances from the disturbances.
Table S2.Winter resource selection function for gray wolves in the Côte-Nord region of Québec, Canada. The model was estimated by comparing the characteristics of 15,572 locations gathered from 9 wolves that belonged to four packs versus 77,860 locations randomly placed within the 95% minimum convex polygon of individual packs. Covariates is the distance (km) to the nearest cut or road.Selection coefficients () are presented with robust standard errors (SE) and associated P-values.
Covariate / / SE / t / PWater body / -0.439 / 0.330 / -1.35 / 0.18
Lichen-heathcommunity / 0.527 / 0.180 / 2.92 / 0.004
Conifer forest / 0.376 / 0.067 / 5.64 / <0.0001
Peatland / 0.229 / 0.306 / 0.75 / 0.46
Burn / 0.186 / 0.090 / 2.06 / 0.04
Mixed forest / 1.120 / 0.111 / 10.03 / <0.0001
Road / 0.299 / 0.161 / 1.86 / 0.06
s2 / -2.0e-5 / 7.1e-5 / -0.25 / 0.80
Interaction term
s × Lichen-heath community / -0.017 / 0.006 / -3.02 / 0.003
s2×Conifer forest / -8.1e-4 / 1.2e-4 / -6.92 / <0.0001
s2 ×Mixed forest / -1.0e-4 / 5.0e-5 / -1.93 / 0.05
*Cuts (and other) were considered as the reference category
Fig. S1. Relative probability of occurrence of gray wolves among land cover types, as a function of the distance to the nearest clearcut or roadin the Côte-Nord region of Québec (Canada), as estimate fromTable S2.
Sensitivity of per capita mortality risk predictions for caribou
We tested the sensitivityof our predictions of per capita mortality risk for caribou by imposing a 15% change in the parameters NC, NM, hM, hC, a, C, wCi, and wMi. We found that the average per capita risk increased by more than 15% only when changing parameters NM, wC_lichen, wC_conifer, and wM_conifer. Even in these four cases, the shape of the spatial changes in mortality risk remained the same.In fact, the shape of the model’s predictionswas only sensitive to parameters (s)and ), which in turn largely depends on (see above).We thus determined how per capita predation risk should vary given the expected variations in for caribou and moose. We consideredthe average -value, together with the lower and upper limits of its 95% confidence interval, i.e., = 3.7, 0.35, and 6.00, respectively, for caribou(see Fig. B1 of [4]), and = 0.37, 0.11, and 0.62, respectively, for moose (see ExpressionsvM(s)andwMiabove). The PDE’s numerical solutionbased on these -values yieldedthree density distribution functions()for caribou (Fig. S2A) and for moose (Fig.S2B). Predictions of were much more sensitive to the range of -values observed for caribou than for moose. Nine combinations of species-specific -values (3 values for caribou × 3 values for moose) could then be considered to estimate per capita mortality risk. We found that mortality risk should increase for caribou at an intermediate distance from cuts and roads – wherecaribou occur at relatively high density –for six of the nine combinations of species-specific -values. Per capita risk was expected to simply decrease with increasing distance from cuts and roads when the predictions are based on the lower limitof the ’s 95% confidence intervalfor caribou ( = 0.35).This situation appears rather unlikely, however, because the prediction is based on a density distribution that does not reflect observed spatial patterns in caribou distribution[4]. Indeed, the prediction would emerge only if caribou were most abundant near cuts and roadsthan far away.By contrast, areas near these anthropogenic features have a relatively low probability of occurrenceof caribou [4, 10-13]. We thus conclude that grouping in the presence of moose is likely to increase the risk of predation for caribouunder most observed field situations.
Fig. S2. A) Expected density distribution as a function of distance from the nearest anthropogenic disturbances (), as predicted from an advection-diffusion movement model with parameters =7.4, 0=4, and b=0.183, together with = 0.35 (Lower), 3.7 (Average), or 6.00 (Upper) for forest-dwelling caribou, and B) withparameters =4, 0=0.04, and b=0.183, together with = 0.11 (Lower), 0.37 (Average), or 0.62 (Upper) for moose. C)Six expected distance-dependent trends in per capita probability of mortality forforest-dwelling caribou, as predicted by combining one of three -values for caribou and one of three -values for moose. Differences in expected mortality risk due to variation in -values for moose are visible at short distance from cuts and roads, but converge distance for each caribou’s -value at farther.
References
1.Équipe. de rétablissement du caribou forestier du Québec. 2013. Plan de rétablissement du caribou forestier (Rangifer tarandus caribou) au Québec — 2013-2023. produit pour le compte du ministère du Développement durable, de l’Environnement, de la Faune et des Parcs du Québec, Faune Québec, 110 p. (
2.Lefort S., Huot M. 2008 Plan de gestion de l’orignal 2004-2010 : bilan de la miplan. (Ministère des Ressources naturelles et de la Faune, Direction de l’expertise sur la faune et ses habitats, Service de la faune terrestre et avifaune, Québec. 38 p.
3.Joly D.O., Patterson B.R. 2003 Use of selection indices to model the functional response of predators. Ecology84, 1635-1639.
4.Fortin D., Buono P.-L., Fortin A., Courbin N., Gingras C.T., Moorcroft P.R., Courtois R., Dussault C. 2013 Movement responses of caribou to human-induced habitat edges lead to their aggregation near anthropogenic features. Am Nat181, 827–836.
5.Latombe G. 2013 Développement d’un modèle centré sur l’individu des déplacements du caribou, du loup et de l’orignal, et de leurs interactions, en forêt boréale aménagée. PhD thesis. Montréal, Québec, Canada, Université de Montréal.
6.Crête M., Rivest L.-P., Jolicoeur H., Brassard J.-M., Mesier F. 1986 Predicting and correcting helicopter counts of moose with observations made from fixed-wing aircraft in Southern Quebec. J Appl Ecol23, 751-761.
7.Johnson D.H. 1980 The comparison of usage and availability measurements for evaluating resource preference. Ecology61, 65-71. (doi:10.2307/1937156).
8.Bolker B.M., Brooks M.E., Clark C.J., Geange S.W., Poulsen J.R., Stevens M.H.H., White J.S.S. 2009 Generalized linear mixed models: a practical guide for ecology and evolution. Trends Ecol Evol24, 127-135. (doi:10.1016/j.tree.2008.10.008).
9.Hebblewhite M., Merrill E. 2008 Modelling wildlife–human relationships for social species with mixed-effects resource selection models. J Appl Ecol45, 834-844.
10.Courbin N., Fortin D., Dussault C., Courtois R. 2014 Logging-induced changes in habitat network connectivity shape behavioral interactions in the wolf-caribou-moose system. Ecol Monogr84, 265-285. (doi:10.1890/12-2118.1).
11.Courtois R., Gingras A., Fortin D., Sebbane A., Rochette B., Breton L. 2008 Demographic and behavioural response of woodland caribou to forest harvesting. Canadian Journal of Forest Research38, 2837-2849. (doi:10.1139/x08-119).
12.Leblond M., Frair J., Fortin D., Dussault C., Ouellet J.-P., Courtois R. 2011 Assessing the influence of resource covariates at multiple spatial scales: an application to forest-dwelling caribou faced with intensive human activity. Landsc Ecol26, 1433-1446. (doi:10.1007/s10980-011-9647-6).
13.Schaefer J.A. 2003 Long-term range recession and the persistence of caribou in the taiga. Conserv Biol17, 1435-1439.
1