Management-mediated predation rate in the caribou – moose – wolf system: spatial con ﬁ guration of logging activities matters

. Landscape complexity can determine the population dynamics of interacting predators and prey. Yet, management plans are commonly developed from aspatial predictive models. This oversight may result in unexpected outcomes or the loss of opportunities to make spatial interventions that would increase a plan ’ s effectiveness. The management of the threatened woodland caribou ( Rangifer tarandus caribou ), boreal population, provides an example of such shortcomings when using an aspatial approach. Currently, the most in ﬂ uential management recommendation is to maintain at least 65% of undisturbed forests in areas occupied by caribou populations, regardless of the spatial con ﬁ guration of the forest cover. Using a spatially explicit individual-based model (IBM), we evaluated the effects of the spatial con ﬁ guration of cuts and roads on the mortality of boreal caribou living in sympatry with wolves ( Canis lupus ) and moose ( Alces alces ), an apparent competitor. Starting with a real forest landscape, we created forest management scenarios of the speci ﬁ c spatial distribution of cuts (mosaic, small, or large agglomeration) with increasing disturbance levels. We then ran the IBM with simulated agents, representing individuals of the three species, moving according to movement rules determined from radio-collared individuals. We found that movement responses to land cover types and roads differed among species. For example, caribou and moose generally avoided areas close to roads, contrary to wolves. Those differences in ﬂ uenced the mortality of caribou agents, which not only depended on the levels of disturbance but also depended on the spatial distribution of cuts and roads. After controlling for disturbance level, wolves were more successful when forest management required an extensive road network resulting in relatively high habitat fragmen-tation. Caribou agents experienced lower mortality in landscapes with low densities of road and disturbance-related edges. The effect remained much stronger, however, for the level than the spatial con- ﬁ guration of human disturbances. Still, our IBM demonstrated how landscape management could be used to manipulate species interactions, with the intent of either increasing or decreasing predation rates on speci ﬁ c populations, depending on management goals.

Abstract. Landscape complexity can determine the population dynamics of interacting predators and prey. Yet, management plans are commonly developed from aspatial predictive models. This oversight may result in unexpected outcomes or the loss of opportunities to make spatial interventions that would increase a plan's effectiveness. The management of the threatened woodland caribou (Rangifer tarandus caribou), boreal population, provides an example of such shortcomings when using an aspatial approach. Currently, the most influential management recommendation is to maintain at least 65% of undisturbed forests in areas occupied by caribou populations, regardless of the spatial configuration of the forest cover. Using a spatially explicit individual-based model (IBM), we evaluated the effects of the spatial configuration of cuts and roads on the mortality of boreal caribou living in sympatry with wolves (Canis lupus) and moose (Alces alces), an apparent competitor. Starting with a real forest landscape, we created forest management scenarios of the specific spatial distribution of cuts (mosaic, small, or large agglomeration) with increasing disturbance levels. We then ran the IBM with simulated agents, representing individuals of the three species, moving according to movement rules determined from radio-collared individuals. We found that movement responses to land cover types and roads differed among species. For example, caribou and moose generally avoided areas close to roads, contrary to wolves. Those differences influenced the mortality of caribou agents, which not only depended on the levels of disturbance but also depended on the spatial distribution of cuts and roads. After controlling for disturbance level, wolves were more successful when forest management required an extensive road network resulting in relatively high habitat fragmentation. Caribou agents experienced lower mortality in landscapes with low densities of road and disturbance-related edges. The effect remained much stronger, however, for the level than the spatial configuration of human disturbances. Still, our IBM demonstrated how landscape management could be used to manipulate species interactions, with the intent of either increasing or decreasing predation rates on specific populations, depending on management goals.

INTRODUCTION
The movement response of individuals to landscape features often vary among species, a difference sufficient to determine interaction rates between opponents and, ultimately, population and community dynamics (Holt 1984, Johnson et al. 1992, Morales et al. 2010, DeMars and Boutin 2018. Spatial heterogeneity in landscape attributes can govern food web interactions (Murrell 2005, Kauffman et al. 2007, Gorini et al. 2012. The key role of spatial heterogeneity has become increasingly recognized following Huffaker's (1958) landmark experiment. The study investigated the impact of numbers and spatial arrangement of food patches and movement barriers in a predator-prey system involving mites (Eotetranychus sexmaculatus as prey, Typhlodromus occidentalis as predator). Huffaker created different environments using oranges and rubber balls and followed changes in the interacting populations. His results showed that landscape complexity could stabilize predator-prey dynamics by reducing the rate of interactions among opponents and promote species coexistence (Huffaker 1958) under conditions where extinction would be expected in homogeneous environments (Gause et al. 1936). Thus, landscape composition and structure should be considered when studying population dynamics (Wiegand et al. 1999), but spatial patterns in vital rates remain largely overlooked in population studies.
Spatial heterogeneity has relevance beyond fundamental studies of population and community ecology. Wildlife management and conservation rely on the accurate forecast of biological systems (Wiegand et al. 1999, Fryxell et al. 2005, Sibly et al. 2009, Kareiva and Marvier 2015, Warwick-Evans et al. 2018), a challenge that may require spatially explicit approaches (Evans 2012, Stillman et al. 2015. For example, the spatial model of Mahoney et al. (2018) for coyote (Canis latrans) lethal control in a mule deer (Odocoileus hemionus) system demonstrated that predictions of the efficacy of predator management could be more reliable and management more effective by considering specific attributes of landscape heterogeneity (i.e., land covers, distance between patches, terrain elevation, forage quality for deer), which influence spatial overlap between coyote removal risk and deer space use. Implementing management plans developed from nonspatial models can have unexpected-even undesirable-consequences on the system because the result of management actions may vary with landscape composition and physiognomy (Wiegand et al. 1999, Mahoney et al. 2018. A persisting challenge of conservation planning lies in the preservation of organisms despite ongoing habitat changes due to the exploitation of natural resources (Margules andPressey 2000, Kareiva andMarvier 2015). Decision-makers can rely on scenario planning, which explores the uncertainty surrounding the future consequences of a decision by contrasting scenarios of a plausible future (Peterson et al. 2003). Scenarios can be tested with individual-based models (IBMs) that account for detailed knowledge of animal behaviors (e.g., movement, foraging, memory) in spatially explicit and dynamic environments Railsback 2005, McLane et al. 2011). Individual-based models are particularly useful when in situ experiments are difficult or impossible to conduct because of the complexity, size, and slow dynamics of ecological systems (Grimm and Railsback 2012). This situation represents the biological system involving the boreal population of woodland caribou (Rangifer tarandus caribou) in managed boreal forests.
The boreal population of woodland caribou (hereafter boreal caribou) is threatened in Canada (COSEWIC 2014). Boreal caribou are under top-down control (Seip 1992), and predation by gray wolves (Canis lupus) is considered the main driver of population decline (Rettie andMessier 1998, McLoughlin et al. 2003,Équipe de rétablissement du caribou forestier du Québec 2013). Anthropogenic disturbances, such as logging activities, increase habitat loss and fragmentation, which tends to intensify predation rate on caribou by triggering apparent competition (James et al. 2004, Wittmer et al. 2007). Deciduous vegetation invades cutovers few years after logging (Potvin et al. 2005, Gagné et al. 2016, which is the primary food source for moose (Alces americanus, Crête 1989). This increase in vegetation can trigger an increase in moose abundance, followed by a numerical response of the wolf population that can exacerbate predation risk of boreal caribou (Seip 1992, Wittmer et al. 2007). Moreover, linear features such as roads facilitate the movement of wolves across caribou ranges (James and Stuart-Smith 2000, Dickie et al. 2017, DeMars and Boutin 2018 and increase predation risk of caribou in proximity (Latham et al. 2011, Whittington et al. 2011. Because of such negative impact of humaninduced activities, Environment Canada (2011) recommended that critical habitat for boreal caribou should contain <35% of total disturbance, as defined by anthropogenic disturbances (plus a 500-m buffer) and by fires. With this threshold, a local population has a 60% probability of being self-sustaining (λ ≥ 1), given national averages in vital rates. However, the demographic consequence of maintaining at least 65% of undisturbed habitat varies among populations (Environment Canada 2012, Rudolph et al. 2017 because, for example, the same level of total disturbance can be associated with different compositions (e.g., a landscape with more linear features and fewer cutblocks than another) and that spatial configuration of landscape composition is likely to influence the spatial game of interactions between caribou, wolves, and moose differently . The three species do not react similarly to cuts, roads, and other landscapes features , 2015, Gagné et al. 2016, and those species-specific responses generate spatial patterns in predator-prey co-occurrence (Courbin et al. 2009 and prey mortality (Wittmer et al. 2007, Losier et al. 2015. The caribou-moose-wolf system thus provides an opportunity to use IBMs to inform conservation planning by revealing how different spatial patterns of human disturbances could influence predation rates.
In this study, we used a spatially explicit IBM to evaluate the relative effects of landscape configuration on caribou mortality for a population living in sympatry with wolves and moose in managed boreal forests. Specifically, we used realistic management scenarios differing in terms of spatial distribution of cuts and roads and levels of total disturbance to compare the relative rate of caribou mortality within a given season. Caribou, moose, and wolf agents moved in the virtual landscapes following rules that were empirically determined from the movement of GPS-collared individuals of each species. Our intent was not to simulate a specific system with all related model parameters but instead to parameterize the IBM with behaviors and vital rate outputs that are realistic for a caribou-moose-wolf system.

Forest management scenarios
We applied 12 virtual management scenarios to the undisturbed Montagne-Blanche region (50°33 0 N to 51°18 0 N, 69°46 0 W to 70°59 0 W) included in the Saguenay-Lac-St-Jean and the Côte-Nord area of Québec (Canada), representing a typical natural boreal landscape. We used a 1:20,000 ecoforest map provided by the Ministère des Forêts, de la Faune et des Parcs (hereafter MFFP) as our base map, which, among other information, characterized tree species composition and stand age of forested areas. Our scenarios (5065 km 2 , with a 25-m resolution) thus represented realistic landscape composition and configuration that individuals might encounter in nature. Furthermore, by using the same base map for all management scenarios, we were able to focus on the impact of different spatial distributions of cuts and roads on space use of caribou, moose, and wolves and prey mortality. We developed a total of 13 scenarios consisting of 12 management scenarios varying in terms of the levels of total disturbance and spatial distribution of cuts and roads, together with a reference scenario of human-undisturbed landscape (the original map of the region; Fig. 1). Following Environment Canada's (2011) approach, the levels of total disturbance were calculated as the percentage of the landscape of the nonoverlapping surface of burns, roads, and cuts, the latter two buffered by 500 m. Burned areas (≤50 yr old) covered 7% in all 13 landscapes. Cuts and roads were simulated in the 12 disturbed landscapes in order to obtain ranges of total disturbance of approximately 30%, 45%, 60%, and 75% of the landscape, which represent disturbance levels currently experienced by local populations of boreal caribou throughout its range (Environment and Climate Change Canada 2019). Cutblocks were organized in mosaic, small, or large agglomerations. Mosaic cutting consisted of dispersed clear-cuts with at least a 200-m buffer of residual forest around them and a 100-m buffer when cutblocks were <0.25 km 2 (MFFP 2018).
v www.esajournals.org Mosaic cutting thus had a checkerboard pattern of residual and logged forest patches of 0.85-1.00 km 2 each. Small-agglomeration cutting involved the grouping of adjacent cutblocks in sectors of 10-50 km 2 , whereas largeagglomeration cutting grouped cutblocks in larger sectors of 40-150 km 2 (Fig. 1). For a fixed level of disturbance, the total amount of cuts and roads differed between the three management strategies due to differences in the level of overlap between the 500-m buffers and differences in the roads required to connect the cutblock network (Fig. 2).
Forest management scenarios were generated using Remsoft's Woodstock and Stanley planning models 2017.11 (Walters et al. 1999, Remsoft, Fredericton, New Brunswick, Canada). The base map was split into spatial organization compartments (SOCs) of 50-150 km 2 , generally used by the Quebec government (MFFP 2018) to create spatial patterns of logging in line with current forest management harvesting guidelines. Each SOC was gradually available to harvest from south to north to consider the road network's deployment. Cuts in an available SOC occurred over 10 yr, and only stands that reached the exploitability age were harvested (as determined by annual growth rate curve; Pothier and Savard 1988). The targeted level of harvest was reached after a 50-yr period in each simulated landscape. Simulations of the IBM were then run on these resulting landscapes. Mosaic and smallagglomeration cutting were simulated using the Stanley planning model, whereas largeagglomeration cuttings were created with the AAUnit module of the Woodstock planning model. Road networks were constructed to connect centroids of harvested stands. First, we defined the same starting point for all scenarios, which was randomly selected within the south part of the landscape to create the first link of the road network. Second, we delimited a polygon for each set of adjacent cuts. For each of those polygons, we identified the centroid. Third, we connected each centroid to the nearest connected centroid using least-cost path algorithm to avoid crossing large rivers (i.e., water had a high value of cost and the rest had low value). Finally, roads were transformed into 25 m wide polygons. Arc-GIS 10.2.2 (ESRI 2014) and R software (R Development Core Team 2016) were used to create road networks and to incorporate them in forest management scenarios.

Movement rules derived from radio-tracking caribou, moose, and wolves
To identify species-specific movement rules implemented in the IBM, we radio-tracked caribou, moose, and wolves in landscapes under different levels of forest management. Radiocollared caribou were followed in the study area (48°N-52°N, 66°W-80°W) covering approximately 380,000 km 2 of the boreal forest in the Jamésie, Saguenay, and Côte-Nord regions of Québec, Canada, whereas radio-collared moose and wolves were followed over 11,500 km 2 in the Côte-Nord region. Forest harvesting is the main industrial activity occurring in those regions. Caribou densities vary between 1.6 and 2.5 individuals/100 km 2 (mean density throughout the area is 1.5 caribous/100 km 2 ;Équipe de rétablissement du caribou forestier du Québec 2013), and moose occur at densities ranging between 4.3 and 10.0 individuals/100 km 2 (Fortin et al. 2008, Lefort andMassé 2015). Although poorly documented, studies from nearby areas suggested wolf density range between 0.3 and 1.5 individuals/100 km 2 with a mean territory size of 1000 km 2 (Larivière et al. 1998, Jolicoeur andHénault 2002).
A total of 60 female caribou were equipped with GPS collars (i.e., 2200 or 3300 model, Lotek Engineering, Newmarket, Ontario, Canada, or TGW 3600 model, Telonics, Mesa, Arizona, USA, the total area disturbed (without the 500-m buffer, mosaic = 62.4 km 2 , small agglomeration = 60.0 km 2 , large agglomeration = 57.5 km 2 ; with the 500-m buffer, mosaic = 256 km 2 , small agglomeration = 155 km 2 , large agglomeration = 147 km 2 ). Other colors represent land cover types: red = burned area, blue = water body, dark green = closed-canopy conifer forest, green = open conifer forest without lichen, light green = open conifer with lichen, orange = mixed and deciduous and purple = other. Fig. 1 (2013), we defined three periods by merging caribou and moose biological seasons: pre-calving/calving season (mid-May through late July, which covers the calving season of both boreal caribou and moose), late summer (early August to the end of September), and winter (October to mid-May). Wolf locations were also separated according to the same periods.
We characterized land cover in the study area using Landsat Thematic Mapper image taken in 2000 with a 25-m resolution grid (Natural Resources Canada, Canadian Forest Service). The initial 48 land cover categories were regrouped into seven classes: closed-canopy conifer forest, open conifer forest with lichen, open conifer forest without lichen, mixed and deciduous forest, lichen heath community, water bodies, and other (including heath without lichen, open area, wetland, and unclassified areas, following Courbin et al. 2009). Land cover maps were updated every year by adding recent cuts (0-20 yr old), old cuts (21-50 yr old), roads, and fires (0-50 yr), which resulted in a total of 11 land cover types.

Statistical analysis of empirical movement rules
The empirical movement rules were determined from species-specific step selection functions (SSFs, sensu Fortin et al. 2005). SSFs Fig. 2. Relationship between road (a), edge densities (i.e., perimeter of 500-m buffer around cuts and roads) (b), and proportion of cuts (c), as a function of the proportion of total disturbance (burned areas and cuts and roads with 500-m buffer) for each of the three forest management strategies (i.e., mosaic, small agglomeration, and large agglomeration cutting). Overall relationships are represented by the dashed line.

Fig. 2 (Continued)
v www.esajournals.org provide the relative probability of selection among a set of options based on the comparison of observed and random steps (i.e., the linear segment between successive locations at 1-h interval) using conditional logistic regression. Accordingly, we drew 20 random steps with the same starting location as the observed step but with different lengths and turning angles (i.e., change in direction between the previous and the current step). Length and turning angle of random steps were both drawn from a uniform distribution within a disc of radius equal to the 99th percentile of step length distribution from all individuals of a given species (caribou = 1800 m, moose = 950 m, and wolf = 6800 m). For both observed and random steps, we then extracted  step length (SL), step turning angle (TA), land cover type (LC i ), and distance to the nearest road (DR j ) at the end of the step. We built speciesspecific SSFs taking the following form: where w(x) represents the SSF score of the step described by the vector x of variables SL, TA, LC i , and DR j with associated coefficient β LogSL , β SL , β TA , β LC i , and β DR j . The term LC i corresponds to dichotomous covariable i among a set of covariables (i.e., closed-canopy conifer forest, open conifer forest with lichen, mixed and deciduous forest, lichen heath community, water bodies, recent cuts, old cuts, roads, and fires, with open conifer forest without lichen being the reference category) representing the land cover type at the end of the step. The term DR j corresponds to dichotomous covariable j among a set of potential covariates representing classes of distance to the nearest road (i.e., 0-250, 250-500, 500-1000, and 1000-1500, with >1500 m being the reference category). The model includes both the natural logarithm of step length ln(SL) and step length SL (Table 1), as recommended by Nicosia et al. (2017). We parameterized the IBM's movement models with coefficients β LogSL , β SL , β TA , β LC , and β D providing some empirical evidence that their covariates triggered a movement response, that is, all β values with P values <0.1, a threshold selected to reduce the risk of type II error. Because each species responded differently to multiple habitat features, including roads and cuts, spatial changes in landscape configuration and composition should yield different outcomes in terms of mortality rates. We assessed the robustness of our models of all species with kfold cross-validation following Latombe et al. (2014). A detailed description of model fitting is available in Appendix S1.

Individual-based model
Overview. -Fig. 3 depicts the model flowchart during run time. At each time step, agents moved one by one, in random order, according to the empirically derived movement rules (see Statistical analysis of empirical movement rules). The movement process was executed first, and then, the predation process was computed. When a prey died, it was removed from the simulation, and no more information was recorded for that agent. For a given season, agents traveled in virtual landscapes with 1-h time step (as GPS-collared animals) for a total of 7320 consecutive steps. This large number of steps was used to obtain relatively stable mortality rates among replicates for a season, while considering computation time limitations. We ran 10 replicates for each scenario-season combination (13 scenarios × 3 seasons × 10 replicates), for a total of 390 simulations. Data collected at the end of simulation included individual identifier, location, date, time, and the number of kills made by each wolf. Computer simulations were conducted in virtual landscapes of 5065 km 2 with a 25-m resolution and were classified in the same 11 land cover types described above. Landscapes were superimposed with a raster of distance to road comprising five classes: (1) ≤250 m, (2) 251-500 m, (3) 501-1000 m, (4) 1001-1500 m, and (5) >1500 m.
Initialization of agents.-Caribou occur in groups for most of the year (Équipe de rétablissement du caribou forestier du Québec 2013, COSEWIC 2014) and spread themselves out during the calving season to reduce predation risk (Pinard et al. 2012). While we do not make the explicit link between group size and caribou behavior, our IBM implicitly considers the reaction of caribou within a group of typical size for each season because movement process of caribou was based on the behavior of 60 GPScollared individuals that were in groups. Each simulation was initialized with 95 caribou agents. When an agent died, it was removed from the virtual landscape. This approach was similar in all three landscapes, thereby allowing for the comparison of caribou mortality between management scenarios, for a given season.
For the three seasons, a moose agent represented a single individual, whereas a wolf agent represented a wolf pack (i.e., meta-individual). In our simulations, we accounted for the typical numerical response of moose that follows timber harvesting (Potvin et al. 2005, Anderson et al. 2018) and the concurrent increase in wolf density. We thus adjusted the number of moose in a simulation to the total disturbance and then adjusted the number of wolf packs to moose density based v www.esajournals.org on Messier (1994). Moose numbers varied between 200 and 250 and wolf packs between 3 and 8, depending on the levels of total disturbance.
Caribou and moose were randomly distributed across the landscape within a patch of their most strongly selected land cover type, that is, conifer with lichen forest for boreal caribou and mixed and deciduous forest for moose (Table 1). Wolves were also randomly distributed across mixed and deciduous forests because moose is a principal prey for wolves (Messier 1994, Tremblay et al. 2001. Details on movement and predation process. -At the beginning of a time step, agents had a 50% probability of leaving the land cover patch (i.e., set of adjacent landscape cells [25 × 25 m] of the same land cover type) they currently occupied. If they stayed, they randomly moved inside the patch (i.e., to another landscape cell), whereas if they left, they followed their movement rules parameterized with SSF coefficients (Table 1). When prey agents left the patch, 20 random steps were drawn within a buffer around their current location with a radius covering the 99th percentile of all observed step lengths for that species. The prey agent then moved to the location with the highest SSF score (w[x]). In the case of wolf agents, however, the buffer's radius (i.e., maximum step length) depended on whether the agent was hunting or not. Wolves display two fundamentals movement modes: the hunting mode, during which they actively search for prey, and the stationary mode, during which they consume prey and remain near the kill (Hebblewhite et al. 2003). Twenty locations of potential steps were drawn within a buffer radius of 6800 m in hunting mode and 200 m in stationary mode (Webb et al. 2008) to reproduce the two modes (Fig. 3). As prey agents, wolf agents moved to the destination of the step with the highest SSF score. Once every agent had moved, wolves were able to detect prey in their vicinity (≤1 km) and kill one when in hunting mode (wolves always started in hunting mode). This modeling approach is similar to Latombe et al. (2014) and is in line with the way our SSFs were estimated.
The predation process included two main components: (1) the probability of encounter and (2) the probability of death given an encounter (Lima and Dill 1990). First, we considered that an encounter occurred when a prey was within a 1km radius of a wolf, a distance within the detection range reported for wolves (Mech and Boitani 2003). This distance was also used in previous studies (DeMars et al. 2016, Spangenberg et al. 2019 and was similar to distances of 1.5 km chosen by Muhly et al. (2010) and of 1.3 km considered by Whittington et al. (2011). When both caribou and moose were within the detectable distance, wolves chose moose over caribou every time (Mech and Boitani 2003). Second, the probability of making a successful kill following an encounter reflected the combined probabilities of launching an attack given an encounter and making the kill following the attack (Lima and Dill 1990). For this latter probability, we considered that a prey faced a 20% risk of dying following an attack (as used by Latombe 2013). Previous studies have reported such relatively low killing success (~20%) of wolves hunting large herbivores , Wikenros et al. 2009). Given the 20% risk of dying following an attack, we adjusted the probability of an attack given an encounter so that the entire predation process resulted in 10% mortality over the 305 d that lasted a simulation (7320 time steps, see Overview) during which agents traveled with a season-specific, empirical movement rule. This mortality rate corresponds approximately to the annual mortality rates field reported in the study area (~10% for boreal caribou and~10% for moose; Crête and Courtois 1997,Équipe de rétablissement du caribou forestier du Québec 2013). Given the calibration process, the IBM should provide realistic prey mortality rates, despite some uncertainties in parameter values associated with the encounter rate (i.e., group behavior and detection distance) and kill rates (i.e., probability of kill given an attack). The calibration resulted in a 2.5% probability that wolves actually detect and then launch an attack when a prey is located ≤1 km (details of the calibration is provided in Appendix S2).
After a kill, wolves entered the stationary mode for 48 h, a period during which a pack could not kill another prey. Handling time varies largely among studies. Hayes et al. (2000) observed that wolf pack handled moose during a mean 2.9 d (range: 2.6-3.3 d) and caribou during 1.3 d. Large packs, however, are able to consume v www.esajournals.org  Messier and Crête (1985) reported that wolves may remain at a kill site for 8-23 d. Our assumption of a fixed handling time does not impact which prey will be killed next. Therefore, it should not influence our assessment of how the interplay between species-specific movement rules and complex changes in habitat features impacts the relative risk of mortality.

Analysis of IBM's outputs
Validation.-We validated our IBM with a bottom-up analysis strategy (Grimm and Railsback 2005) by verifying that individual-level behaviors of agents were consistent with the empirical data. To do so, we compared the proportion of use of land cover types and step length distributions between empirical data and the output of simulations conducted under the same landscape conditions. We ran the simulations in 5065km 2 landscapes (25-m resolution) that represented the study area of the Côte-Nord region from 2005 to 2009 (one landscape per year), that is, when and where GPS-collared individuals of all three species were tracked. We ran our simulations with 95 caribou, 210 moose, and five wolf packs, which corresponded to the mean densities of these three species observed over this study area. Simulations were run 10 times for each year and each season (10 replicates × 5 yr × 3 seasons = 150 simulations; see Appendix S3 for details). We then used relocation data from the radio-collared individuals to evaluate IBM's performance. Linear regression was used to compare patterns of land cover use by agents and real individuals. This process thus informs if the movement and predation rules imposed on agents were adequate and sufficient to reproduce the behaviors of the actual radiocollared animals.
Analysis of prey mortality.-First, we assessed whether the different spatial configuration of cuts and roads influenced prey mortalities. We used Poisson regression to relate the number of mortalities (caribou or moose) to the proportion of total area disturbed and its interaction with the three forest management types (i.e., either mosaic, small-, or large-agglomeration of cuts, with small agglomeration being the reference type). The model did not include the forest management type as the main effect and only included the different types as part of interaction terms. All management scenarios were implemented on the same basic landscape, and accordingly, mortalities in an unharvested landscape would necessarily be the same, that is, independent from the type of management that is later applied. In statistical terms, there should be a single intercept for the three management types for the relationship between mortalities and level of disturbance, but slopes could vary with management types. All assumptions were met for those models (e.g., variance inflation factors [VIFs] < 1.6).
Second, we quantified how wolf predation on caribou and moose varied with habitat features that mainly set apart the three forest management types, that is, road and edge densities (Fig. 2). This analysis thus provided an assessment that is less specific to the three forest management strategies and more generally linked to habitat structure. We used Poisson regression to relate the number of mortalities to the proportion of total disturbance and road and edge densities. The edge density for a scenario was estimated using the perimeter of all anthropogenic disturbance regions in the landscape. Following Environment Canada (2012), a disturbed region included the areas covered by cuts and roads together with their 500-m buffer (disturbed regions in Fig. 1). Road and edge densities increased with the level of total disturbance in the landscape (Pearson's r = 0.91 for roads and move one by one according to specific movement rules. For wolves, the mode in which they will influence the buffer radius within random points is drawn, but their movement rule remaining unchanged. Once all agents have moved, wolves in hunting mode can move to prey's location (i.e., they can move twice) and kill a prey depending on the success of the three stages (i.e., detection, attack, and kill). The cycle then starts over until the maximum number of steps is reached (i.e., 7320 steps).

Fig. 3 (Continued)
v www.esajournals.org 0.76 for edges). To evaluate the additive effect that road and edge densities had on prey mortality without facing multicollinearity issue, we used residual values of road and edge densities from the relationship they shared with the total amount of disturbance (Fig. 2). For example, a high residual value for road or edge density indicated a large deviation from the expected value (i.e., higher density than expected), given the proportion of total disturbance. Hence, it reflects differences in how a particular management practice requires a relatively extensive road network or generate relatively high edge density. The VIFs of all covariates were <2 in all final models.

Movements rules from radio-collared boreal caribou, moose, and wolves
The response of radio-collared individuals to forest cover types varied among species and seasons (Table 1). Boreal caribou and wolves selected mature conifer forests with lichen yearround, while moose avoided those forests. In contrast, mixed and deciduous forests were selected by moose year-round and by wolves in winter, whereas caribou always avoided those forests. The three species also avoided closedconifer forests, although those forests were used in proportion to availability by wolves during the calving season and by caribou in summer. Heaths with lichen were selected by caribou year-round and by wolves in summer, while moose avoided them in summer and winter. The species also selected burned areas during most of the year, except during calving for caribou and winter for moose.
The behavioral responses to anthropogenic features also varied among the three species and seasons (Table 1). Wolves avoided recent cuts year-round, caribou avoided those cuts in summer and calving, and moose selected them in summer. Although all species used old cuts in proportion to availability during most of the year, wolves selected them in summer, whereas caribou avoided them in summer and moose in winter. Wolves showed a strong selection for roads in all seasons. Unlike wolves, moose avoided roads during most of the year, except during summer. Caribou generally avoided areas within 1000 m of a road, except in summer, where road avoidance was within 500 m. Opposite to caribou, wolves selected areas closer than 250 m of roads all year, and up to 1000 and 1500 m during winter and calving, respectively. Moose selected areas between 250 and 500 m of the roads in summer, while they avoided being closer than 1500 and 250 m of roads in winter and in calving, respectively.

Model validation
Before using the IBM to assess the interaction between agents in virtual landscapes shaped by different management scenarios, we verified that our agents behaved in virtual landscapes as GPScollared individuals did in the actual landscapes with the same characteristics. We found that agent species were able to reproduce the behavior of GPS-collared individuals similarly under the same conditions (Appendix S3).
Step length distributions, however, were generally comprised of longer steps for agents than radiocollared individuals (Appendix S3). Differences between simulated and empirical distributions could be expected for two principal reasons: (1) When individuals stayed in a patch, they moved randomly without any constraint on distance (i.e., uniform step length distribution within the patch), and (2) when wolves attacked a prey, they moved to its location, which imposed a step length that can reach 1 km. Despite the longer steps, in the end, there was a strong link between the proportion of land cover used (i.e., number of observed steps in land cover/total of steps) by agents and by GPS-collared individuals (all R 2 > 0.70; Appendix S3: Fig. S1). Also, agents of each species used land cover types in relatively similar proportions as radio-collared individuals (Appendix S3: Figs. S2-S4).

Wolf predation on caribou and moose under different management scenarios
First, wolf predation on caribou and moose agents varied among the three forest management strategies (Fig. 4). As could be expected, prey mortalities consistently increased with the proportion of the total amount of disturbance in the landscape during any season (Table 2, Fig. 4). After controlling for disturbance level, however, the spatial configuration of cuts and roads resulting from the different forest management practices further impacted prey mortalities. During calving, both caribou and moose agents experienced a stronger increase in mortality in forests managed by mosaic cutting or large-agglomeration cutting than by smallagglomeration cutting (Fig. 4). In winter, caribou and moose experienced marginally higher mortality (P = 0.08 and P = 0.10, respectively; Fig. 4. Percentage of caribou and moose agents killed during simulations by wolves during calving, summer, and winter (left to right, respectively) as a function of the proportion of total disturbance (i.e., burnt areas, cuts, and roads with 500-m buffer) in virtual landscapes harvested either by mosaic (blue circle), small-agglomeration (yellow triangle), or large-agglomeration (gray square) cutting. Predictions of mortality number were computed accordingly to model 1 coefficients (Table 2) with increasing level of disturbance and then reported in percentage (represented by lines). Average mortalities (represented by points) and their standard errors of simulations (n = 10) are also represented for each species and spatial configuration of cuts. Percentage was computed based on the total number of agents at the beginning and the total number of agents killed at the end of the simulation. The caribou population included 95 individuals at the beginning of all simulations, whereas the initial number of moose increased from 200 to 250 individuals with increasing level of disturbance.
v www.esajournals.org Table 2) in forests harvested by mosaic cutting than by small-agglomeration cutting (Fig. 4). In summer, we did not detect differences in prey mortality among the three spatial arrangements of cuts (Fig. 4). Throughout the year, the most detrimental spatial configuration of cutover areas to the prey was mosaic cutting and, to a lesser extent, the large-agglomeration cutting strategy. These spatial arrangements of cutover areas were thus likely relatively unfavorable to caribou. Moreover, high levels of disturbance (which includes burned areas and a 500-m buffer around cuts and roads) were reached by harvesting a smaller portion of the landscape in largeagglomeration cutting, and even more severely, in mosaic cutting when compared to smallagglomeration cutting (Fig. 2).
Second, after accounting for disturbance levels, forest management requiring the deployment of a relatively dense road network (i.e., high residuals) resulted in higher mortality for moose agents during calving, and caribou agents during both calving and winter seasons (Table 2). For example, a 15% increase in road density (i.e., the difference between road density of mosaic cutting and small-agglomeration cutting for a proportion of total area disturbed of 0.60, Fig. 2a) would result in a 3% increase in mortality for caribou (i.e., a 24% increase in caribou mortality; [16.2-13.2%]/13.2%). Caribou agents also experienced higher mortality rates when living where forest management resulted in a relatively higher density of edges involving roads or cuts (high residuals) during calving and winter, while moose suffered this increase during summer and winter (Table 2). For example, a difference of 25% (e.g., the difference between edge density of mosaic cutting and large-agglomeration cutting for a proportion of total area disturbed of 0.72, Fig. 2b) would increase mortality by 3% for caribou during calving.

DISCUSSION
Our study reveals that spatially explicit approaches such as IBMs can enlighten our understanding of how the spatial organization of human disturbances can modulate predator-prey interactions. Specifically, our results show that, while disturbance levels have the upper hand, the spatial configuration of cutblocks and roads over forest management units can impact predation rates in multiprey systems. This information is relevant for habitat management of woodland caribou boreal populations, a species threatened with extinction (COSEWIC 2014). Table 2. Coefficients (β) and robust standard errors of generalized models relating the number of prey (i.e., caribou and moose) killed by wolves: (1) as a function of the proportion of disturbance in the virtual landscape and its interaction with categorical covariates representing forest management strategies (small-agglomeration cutting was the reference category), or (2)  Most studies guiding management of woodland caribou habitat focused on habitat loss by identifying habitat disturbance thresholds where selfsustaining caribou populations become at risk (e.g., Environment Canada 2008, 2012, Rudolph et al. 2017. Here, we demonstrate that the impact of human disturbances on boreal caribou populations is not only related to habitat loss linked with the disturbance levels but also related to how the spatial configuration of cuts and roads are deployed over the landscape, particularly under high disturbance levels. Admittedly, whereas the effect of spatial configuration is smaller than the effect of disturbance level, our study clearly shows that managing the spatial configuration of landscapes could be used to manipulate species interactions. Depending on management goals, the intent could either be to increase or decrease predation rates on specific populations. Our study is also relevant to habitat restoration planning given that all across their distribution range (Environment and Climate Change Canada 2017), boreal caribou commonly inhabit highly disturbed landscapes (above the 35% range disturbance threshold; Environment Canada 2011) where habitat restoration is necessary for population recovery. If managers have limited resources (e.g., human, financial) to carry out extensive habitat restoration and they need to choose among disturbed areas, caribou populations could benefit from restoration plans resulting in habitat configuration that favor caribou survival. Our study provides insights for such an endeavor.
In our IBM, the mortality of prey agents systematically increased with the levels of habitat disturbance, a typical trend reported in several empirical studies (Courtois et al. 2007, Environment Canada 2011, Rudolph et al. 2017. The relatively high mortality rate of caribou populations reported in harvested forest landscapes is often closely linked with the numerical response of wolves resulting from the increase in moose density, which in turn follows the growth of deciduous vegetation in clear-cuts during early post-logging succession (Seip 1992, Wittmer et al. 2007). Our IBM directly integrates the two numerical responses, with the density of wolves increasing from 0.6 to 1.6 pack/1000 km 2 and moose from 3.9 to 4.9 individuals/100 km 2 , as disturbance levels changed from 30% to 75% in our virtual landscapes. Accordingly, our study was not designed to assess the impact of forest disturbance levels on caribou mortality but, instead, to account for this negative relationship while assessing the impact of the spatial configuration of anthropogenic disturbances. Even though the disturbance levels (i.e., habitat loss) had the most impact on prey mortalities, our results indicate that survival can be manipulated by adjusting spatial configuration of logging in a way that influences the risk of predator-prey encounters. The effect of spatial configuration should be particularly manifested in landscape with high disturbance levels (e.g., >35%, Fig. 4), a situation representative of an increasing number of caribou populations (Environment and Climate Change Canada 2017, Fortin et al. 2017, Rudolph et al. 2017). In such a case, managers should first attempt to minimize the disturbance levels, and then select a spatial distribution of roads and cuts favorable to boreal caribou populations.
After controlling for disturbance levels, we found that mosaic cutting was the management strategy yielding the highest mortality rates for caribou during calving and (marginally) in winter. Our findings are consistent with previous studies that reported a negative impact of mosaic cutting on caribou populations Beange 1993, Hervieux et al. 1996). They also echo empirical temporal patterns of population declines that were recently documented on the three boreal caribou populations of the Jamésie region of Québec, where mosaic cutting has been extensively used over the two last decades (Rudolph et al. 2017). In these landscapes, mosaic cutting has benefitted moose (Jacqmain et al. 2008), and the construction of extensive road networks to connect cutblocks (Tittler et al. 2012, Rudolph et al. 2017) has likely facilitated wolf movements (James and Stuart-Smith 2000, Dickie et al. 2017, DeMars and Boutin 2018. Our IBM results show that mosaic cutting should be avoided in a caribou-moose-wolf system when the goal is to preserve caribou populations while minimizing the impact of conservation actions on resource extraction. We could have expected that the midway strategy between mosaic cutting and large agglomeration cutting-small agglomeration cutting-would have a moderate impact on caribou mortality. Our results, however, instead show that small clear-cut agglomerations had lower adult caribou mortality notably because this strategy tended to have fewer roads than large ones (Fig. 2a). Also, their spatial distribution created two contrasting patterns. The implementation of large aggregations of cutovers resulted in virtual landscapes characterized by a few large sectors with highly concentrated roads, which created few but highly risky areas. By contrast, small agglomerations of clear-cuts led to a more uniform distribution of roads. Predation risk was more broadly distributed across the virtual landscape in several areas of moderate risk of caribou. The higher mortality rate observed in a landscape managed with large agglomerations was an emerging property of the complex interplay between broadscale patterns of habitat disturbances and fine-scale responses to habitat features by the different species.
Our IBM was based on fine-scale movement rules derived from field observations of GPScollared individuals. The same rules (i.e., movement and predation processes) were used in all virtual landscapes for a given season; hence, differences in predator-prey encounter rates among scenarios were due to the deployment of cutblocks and roads that created increased predator-prey interactions. Given that we compared different spatial configuration scenarios while controlling for total disturbance level, we can determine for each spatial configuration how roads and edges impact predator-prey encounters. Indeed, wolf predation on caribou was better explained by these features than by categorical variables representing forest management types (i.e., ΔAIC between top model based on habitat features and top model based on management types = −2.7); our findings thus are broadly transferable to other potential management strategies. Roads facilitated wolf movements Stuart-Smith 2000, Dickie et al. 2017), their incursion into habitat refuges of caribou (Latham et al. 2011, DeMars andBoutin 2018), and their encounter with prey (Whittington et al. 2011, Kittle et al. 2017. As for edges, theoretical and empirical studies found that movement responses to cut and road edges create edge effects that affect predator-prey interactions across the landscape, especially in multi-prey systems . Our IBM shows that forest management that requires relatively extensive (i.e., high residuals, see Materials and Methods) road networks and create relatively high edge densities resulted in higher mortality rates for caribou. Consequently, a prey species may see one of its populations decline, while another increase in two landscapes impacted by the same disturbance level, simply because logging activities created different landscape configurations such that the declining population is more strongly exposed to predation.
A few caveats must be considered when interpreting the results. First, our study was specifically designed to assess relative mortality rates between landscapes shaped by different management strategies, for a given species in a given season. In other words, our study was not designed to contrast mortality rates between seasons nor between caribou and moose. This is because, for example, group dynamics of boreal caribou vary importantly among seasons, which can impact predator-prey encounter rates and, subsequently, mortality rates. Group dynamics, however, were not explicitly considered in our IBM. Grouping was still implicitly considered as our agents traveled following movement rules of a typical group-for a given season-because those rules were established under natural field conditions. Group size dynamics also vary between boreal caribou and moose, with moose being generally more solitary than caribou, especially outside of the calving season. Moreover, other factors may influence predation risk, such as the presence of calves (Pinard et al. 2012) and the number of alternative prey and predators (Latham et al. 2013, Wittmer et al. 2013). Our IBM's outputs are thus particularly informative when contrasting management scenarios, for a given season and a given prey species. Second, our study focuses on predation-related mortality of adult females because this factor is a major driver of population dynamics for long-lived species under top-down control (Gaillard et al. 1998(Gaillard et al. , 2000, such as boreal caribou . Given this focus, our study does not provide a full population viability analysis. To gain a more comprehensive understanding of how habitat configuration can impact caribou population dynamics, a next step would be to consider how forest management scenarios impact wolfv www.esajournals.org unrelated mortalities and recruitment to caribou populations. Overall, our study provides insights into habitat management planning for the boreal population of woodland caribou by revealing the influence of spatial configuration of cutover areas and road networks on species interactions. While multiple actions for caribou conservation have been proposed and tested , Environment and Climate Change Canada 2019, few actions involve mitigating the spatial configuration of human disturbances to minimize predator-prey encounter rates. For managing food webs involving large mammals such as woodland caribou, such an approach needs to be designed at broad spatial scales (>5000 km 2 ; Courtois et al. 2004). The use of individualbased modeling then becomes particularly useful to test various potential scenarios of spatial configuration of habitat conservation planning before implementing them over the landscape. Indeed, each application of a conservation plan can substantially impact endangered populations, with concurrent repercussions on regional biodiversity (Bichet et al. 2016). Our study exposes how the loss of woodland caribou critical habitat and its spatial configuration both matter in predator-prey relationships and how these trophic interactions can be dealt through land management to either increase or decrease biological populations.

ACKNOWLEDGMENTS
Simulations were run on the Colosse Supercomputer at the Université Laval, which is managed by Calcul Québec and Compute Canada. This project was funded by the Fonds de recherche du Québec-Nature et technologies (FRQNT) and Engineering Research Council of Canada, Silviculture and Wildlife Research Chair, Sentinelle Nord, the Ministère des Forêts, de la Faune et des Parcs du Québec and Centre d'étude de la forêt (CEF), a research center funded by the FRQNT strategic cluster program. We thank F. Bujold, M. Côté, and J. Rioux from the Ministère des forêts, de la Faune et des Parcs for their previous help with the conceptualization and the creation of the forest management scenarios. We also thank M. Boissonneault for his help with the programming part of the IBM and P. Racine for his suggestion on the construction of the road networks. Finally, we thank J.-P. Tremblay and Y.
Boulanger for their constructive comments on an early version of the paper.