Accounting for food web dynamics when assessing the impact of mesopredator control on declining prey populations

1. Increasing populations of mesopredators are suspected to cause declines in vulnerable wildlife to the extent that mesopredator decimation actions (culling) have become commonplace. Design constraints, especially a lack of spatial replication, often hamper the assessment of the impact of such actions. However, extensive temporal replication (i.e. time series) and accounting for potentially confounding variables may alleviate this problem. 2. In alpine-arctic tundra, the red fox Vulpes vulpes is increasing, while many bird species are declining, likely due to increased predation. Here, we assessed the impact of a long-term (12-year) and spatially extensive ( ~ 3,500 km 2 ) red fox culling action on the red-listed willow ptarmigan Lagopus lagopus in the Norwegian Arctic. Ptarmigan populations were monitored annually in the impact area and in an adjacent no-action area, including a 5-year period before the action commenced. While logistical constraints prohibited monitoring of red fox population densities, the number of culled foxes and three influential food web covariates were monitored after the onset of the culling action. 3. A Before-After-Control-Impact-Paired-Series (BACIPS) analysis without food web covariates indicated that red fox culling curbed the decline of the population in the impact area, and that ptarmigan population density became ~


| INTRODUC TI ON
Anthropogenic landscape-use and climate change in combination with extirpation of apex predators have changed predator communities in many ecosystems (Estes et al., 2011;Prugh et al., 2009;Ritchie & Johnson, 2009). In particular, generalist mesopredators have increased in abundance and expanded their distribution ranges worldwide in recent decades (Gregory & Marchant, 1996;Tapper, 1992).
Anthropogenic-induced release of mesopredators enhances predation pressure on their prey (Conner & Morris, 2015;Prugh et al., 2009;Roos et al., 2018), and ground nesting birds are particularly sensitive as predation is often the most frequent cause of avian reproductive failure (Martin, 1993(Martin, , 1995. Increased nest predation from mesopredators has become a genuine concern regarding the conservation of many bird populations (Côté & Sutherland, 1997;Crooks & Soule, 1999;Kubelka et al., 2018;Marolla et al., 2019;Roos et al., 2018). Increased nest predation appears presently to be particularly acute in the Arctic Kubelka et al., 2018), although no studies have so far been able to attribute this to any particular mesopredator species.
To mitigate predation-induced declines of wildlife populations, mesopredator culling actions have been implemented in many places, however, with quite variable success (Conner & Morris, 2015;Kämmerle & Storch, 2019;Salo et al., 2010). The high mobility and demographic resilience of mesopredator species require that culling actions are intense, long term and large scale. This implies logistic and financial constraints on the design of actions, in the sense that the spatial replication of action and controls that facilitate robust impact assessments is often lacking (Conner & Morris, 2015;Kämmerle & Storch, 2019). Such design constraints represent a common challenge to impact assessments of largescale management interventions (De Palma et al., 2018;Saterson et al., 2004;Stewart-Oaten & Bence, 2001;Sutherland et al., 2004), especially when the intervention is spatially or temporally confounded with other drivers of the management targets (Hewitt et al., 2001;Underwood, 1992). Taking into account such drivers of natural dynamics as covariates in analyses may improve the accuracy of impact assessments when management interventions lack randomization and replication (Stewart-Oaten et al., 1986. The red fox Vulpes vulpes is a generalist mesopredator with global relevance for conservation and wildlife management, and now subject to management actions (i.e. culling) in many places in the world (Fletcher et al., 2010(Fletcher et al., , 2013Kämmerle & Storch, 2019;Marolla et al., 2019). Over the last century, the red fox has expanded into Arctic and Alpine tundra (Gallant et al., 2020). In Scandinavia, this is considered the main threat to critically endangered populations of the arctic fox Vulpes lagopus (through interspecific competition; Angerbjörn et al., 2013 and the lesser white-fronted goose Anser erythropus (through predation; Marolla et al., 2019). For more than a decade, red fox culling has been implemented to conserve these populations. However, increased abundance of red fox in tundra has wider implications (e.g. Elmhagen et al., 2015), such as the recent community-wide decline in arctic-alpine birds in northern Europe (Lehikoinen et al., 2014).
With its circumpolar distribution and high value as a game bird, the willow ptarmigan Lagopus lagopus is a prominent member of the alpine-arctic bird community. Many ptarmigan populations are presently declining in the Arctic (Fuglei et al., 2019). In Norway, where the willow ptarmigan is probably the most comprehensively surveyed terrestrial bird, it was placed on the national red-list in 2015 as 'near threatened' (Henriksen & Hilmo, 2015). While predation is a key driver of alpine-arctic ptarmigan population dynamics (Angelstam et al., 1984;Breisjøberget et al., 2018;Henden et al., 2017;Kausrud et al., 2008;Marcström et al., 1988), the role of different predator species (including the red fox) has been difficult to disentangle. This is because a diverse community of predator species prey on ptarmigan  and the overall predation pressure is also driven by the dynamics of other prey species (e.g. small rodents) in the food web (Moss & Watson, 2001). Yet the red fox is often highlighted as the most influential ptarmigan predator (e.g. Breisjøberget et al., 2018), and red fox culling to mitigate the decline in ptarmigan and other alpine birds is increasingly discussed among stakeholders.
Here, we take advantage of a large-scale and long-term red fox culling action aimed to conserve the endangered arctic fox population in northern Scandinavia. Annual population surveys of willow ptarmigan were conducted both within the red fox culling area (impact) and in an adjacent no-action reference area (control), during 5 years before and 12 years after the action commenced. This allowed us, as a first step, to assess the impact of red fox culling on the ptarmigan populations by means of a Before-After-Control-Impact-Paired-Series (BACIPS) analysis (De Palma et al., 2018;Stewart-Oaten & Bence, 2001). Monitoring data of drivers of food web dynamics known to influence ptarmigan population dynamics were available for the 12 years after the onset of the red fox action. This allowed us in two following steps to assess whether the inclusion of such drivers as covariates in Control-Impact-Paired-Series (CIPS) analyses increased the accuracy of the impact assessment. the target species so that such drivers can be included as covariates in the analysis. This applies in particular to declining bird populations in boreal and arctic food webs ruled by strong multi-annual interaction cycles.

K E Y W O R D S
carrion, culling, generalist predators, management evaluation, red fox, rodents, tundra, willow ptarmigan 2 | MATERIAL S AND ME THODS

| Study system
Our study was conducted across an area of ~6,700 km 2 at 70-71°N in the northeastern part of Finnmark County, Norway. This includes two study regions: Varanger Peninsula (hereafter Varanger), where the red fox culling action took place (i.e. impact region), and the adjacent no-action (control) region Ifjord-Tana (Figure 1). To account for spatio-temporally varying food web dynamics, three subregions corresponding to two separate river basins within Varanger (Komagdalen (Ko) and Vestre Jakobselv (VJ)) and one in Ifjord-Tana (Ifjord) were delimited (Figure 1).
Alpine and low-arctic tundra cover most of the study area. In highlands (>350-400 m a.s.l.) the tundra is barren without continuous cover of vascular plants (Figure 1), while lower altitudes contain more productive shrub-tundra (Killengreen et al., 2007). Optimal willow ptarmigan habitats are located in the lushest parts of the shrub-tundra (in particular riparian tall-shrub; Ehrich et al., 2011) and in the forest-tundra ecotone (Figure 1; Pedersen et al., 2012).
The structure of the tundra food web is qualitatively similar across the study area Killengreen et al., 2007).
However, the study area exhibits spatial variation in the dynamics of important food web components (e.g. Soininen et al., 2018). We here emphasize four food web processes that potentially constitute important drivers of ptarmigan population dynamics in the study area (cf.   Figure 2); small rodents (multi-annual interaction cycles), reindeer (carrion subsidies to predators), red fox (predation and control actions) and humans (ptarmigan hunting).
Owing to their multi-annual population cycles that strongly influence food web dynamics, small rodents (voles and lemmings) constitute the functionally most important component of the herbivore F I G U R E 1 Study design and annual cycle of data sampling. (a) Study areas in northeastern Finnmark. The shades of grey denote tundra in different elevation zone (100 m altitude intervals), where the two light grey zones are areas >400 m a.s.l. with very sparse vegetation. Green denotes areas with mountain birch forest. Large red square denote the Varanger region with the red fox action, while the green square denote the no-action (control) Ifjord-Tana region. Broken lines delineate the two sub-regions within the Varanger action region where the number of red foxes culled was enumerated for the step 3 analysis. The ptarmigan line transects included in the sub-regional analyses are shown as black line segments, while blue transect lines at Ifjord-Tana denote lines included in the regional analysis. Squared ellipses denote areas with monitoring data on small rodent abundance. (b) Annual cycle of data acquisition and the management action (red fox culling)

Iłord-Tana
community in tundra food webs (Ims & Fuglei, 2005). Due to predator functional and numerical responses (cf. alternative prey mechanism; Lack, 1954), the rodent cycle is an important driver of the population dynamics of willow ptarmigan (Kausrud et al., 2008;Myrberget, 1984) and the community of tundra-dwelling bird species (Blomqvist et al., 2002;Ims et al., 2013;Lehikoinen et al., 2014;McKinnon et al., 2014). In our study area, three rodent species exhibit spatially and interspecifically synchronous 4-year cycles Kleiven et al., 2018). However, the amplitude of the rodent cycle varies among the sub-regions (Soininen et al., 2018), which implies spatial variability in predation pressure on ground nesting birds .
Semi-domestic reindeer Rangifer tarandus are abundant and functionally important herbivores in northern Scandinavia (Henden et al., 2014;Ims & Henden, 2012;Ims et al., 2007). Reindeer herds in the study area display density-dependent and climate-driven mortality rates that differ between regions . Resultant reindeer carrion constitute a significant, but spatio-temporally variable, subsidy to predators  that may have predation-mediated knock-on effects on bird populations in tundra Marolla et al., 2019).
The red fox is an omnipresent, abundant generalist mesopredator in the study regions Killengreen et al., 2012). It responds both functionally and numerically to the abundance of rodents and reindeer carrion (Henden et al., 2009(Henden et al., , 2014Killengreen et al., 2012). Over the last century, the red fox has increased in numbers and expanded its range in northern Scandinavia, and now poses a serious threat to several wildlife species (Elmhagen et al., 2015(Elmhagen et al., , 2017Henden et al., 2009). Consequently, red fox culling campaigns to reduce these negative impacts are ongoing in several places (Angerbjörn et al., 2013;Marolla et al., 2019). In northern Scandinavia, the most intensive and spatially extensive campaign is conducted on Varanger Peninsula as an effort to conserve the endangered arctic fox population . Over the years 2005-2016, field inspectors from the Norwegian Environment Agency and local hunters have culled ~2,550 red foxes. The local hunters have been encouraged by the payment of a financial incentive for each culled fox.
Ptarmigan hunting affects ptarmigan population dynamics proportionally to the harvesting rate Pedersen et al., 2004). Considerable changes in hunting regulations in the study area over the last two decades have likely caused temporally variable impacts of hunting on ptarmigan populations. In particular, a bag-limit implemented since 2010 has probably contributed to a lower harvesting rate. The hunting pressure is also likely to vary spatially, for instance, depending on the distance from roads and cabins to hunting grounds.

F I G U R E 2
Conceptual model denoting the main food web and management drivers of willow ptarmigan density. Solid lines denote direct effects, while stippled lines denote indirect effects of different drivers on ptarmigan population density. Black filled circles denote the expected signs (±) of driver effects. Boxes with grey perimeter lines denote predictor/covariates and response variables included in the CIPS model. Values with red perimeter lines denote estimated standardized coefficients of the covariates/predictor with 95% CI

| Data sources and variables
Here we provide a brief overview of the willow ptarmigan data.
Detailed descriptions of the variables used as food web covariates in the analyses (ptarmigan harvest, small rodent abundance, carrion abundance and red fox culling) are given in Appendix S1.
Time series of willow ptarmigan population counts during 2000-2016 were obtained from Hønsefuglportalen (http://honse fugl.nina.no/ Innsy n/) and are based on annual surveys of transect lines distributed across Finnmark County Nilsen et al., 2018).
Transect lines are surveyed annually from the 5th-20th of August by trained personnel with pointing dogs according to a distance sampling protocol (Buckland et al., 2001). A subset of 28 transect lines are located within our study area (Figure 1), 14 in the Ifjord-Tana region and 14 in the Varanger region. In the two sub-regions Ko and VJ in Varanger, we used eight and six transect lines, respectively, while in the sub-region Ifjord we used the eight transect lines closest to the focal watershed ( Figure 1). We obtained regional and sub-regional ptarmigan population density estimates from the state-space modelling approach described in Appendix S1.
We did not include food web covariates in this analysis, because the data on small rodents were not available before the onset of the red fox culling action. The data on ptarmigan harvest rate were also of poorer quality in the first part of the time series, due to low reporting rates of harvest from hunters . Hence, we conducted the BACIPS analysis at the spatial scale of the two regions Ifjord-Tana (control) and Varanger (impact). To correct for large-scale synchrony in ptarmigan population dynamics between the two regions and thus to partial out the within-year differences in population densities between the two regions, year and region were specified as nested random factors. To obtain region-specific estimates (i.e. contrasts) of change in population density between the two periods before and after the red fox action and between the impact and control region for both these period, we specified a fixed factor with four levels (control_before, control_after, impact_before and impact_after). Because potentially influential covariates were not included in this analysis, we corrected for serial correlation in the residuals by including a first-order autoregressive (AR1) term as a random effect in the model.
In steps 2 and 3, the continuous covariates were standardized (i.e. scaled with mean = 0 and SD = 1) to ease comparison of their relative effects. To correct for attributes of rodent cycles that were not accounted for by annual rodent abundance index such as different species composition and seasonal dynamics in the different cycles, the three consecutive rodent cycles were specified as random effects (Pinheiro & Bates, 2000). Moreover, we corrected for serial correlation in the residuals by including a first-order autoregressive (AR1) term as a random effect in the models.
In all analyses, we applied linear mixed-effects models (LMM, package nlme) in the software R (R Core Team, 2019) fitted to log-transformed ptarmigan population density as the response variable. We assessed models for constant variance of the residuals, presence of outliers and approximate normality of the random effects.

| Accounting for sub-regional scale food web dynamics (CIPS analyses)
The time series of the three food web covariates and the ptarmi-

| D ISCUSS I ON
Assessments of large-scale management actions aimed at conserving vulnerable wildlife will seldom attain the strength of inference typically obtained by properly designed experiments. When replication and randomization of actions and controls are impossible due to the logistic/financial constraints implied by large-scale management actions, the best option for making robust inferences is the Before-After-Control-Impact-Paired-Series (BACIPS) design (Christie et al., 2019;De Palma et al., 2018;Stewart-Oaten & Bence, 2001). The lack of spatial replication of most large-scale BACIPS designs is alleviated by extensive temporal replication (i.e. paired time series).
However, the merit of this approach hinges on whether other drivers of temporal dynamics are confounded with BACIPS design. In case of confounding drivers, assessments should be model-based, in the sense that influential external drivers are included as temporal covariates in the statistical model applied to estimates of the effect of the management action. This is rarely done, however, as time series in BACIPS designs are usually too short (cf. De Palma et al., 2018) to allow robust parameterization of covariates. Also, the design of ecological monitoring programs associated with management actions is often inadequate, which means that data on influential covariates are lacking (Lindenmayer & Likens, 2010). Both of these shortcomings apply to assessments of predator control actions aimed to conserve ptarmigan and other grouse species (Kämmerle & Storch, 2019).
Using data from an exceptionally large-scale and long-term mesopredator culling action in the Norwegian Arctic, we significantly improved the accuracy of the impact assessment by including influential drivers of food web dynamics as covariates in the analysis.
Although the BACIPS analysis also provided evidence for an effect of the culling action, the effect appeared to be profoundly underestimated ( Figure 5). Adjusting for food web covariates after the onset of the action yielded an estimated difference (contrast) in ptarmigan population density between the action and control region (estimate (log): 1.45, 95% CI: [0.91, 1.99]) that was 2.3 times larger than the unadjusted contrast estimate (BACIPS model estimate: 0.63, 95% CI: [0.29, 1.45)]. The reason for the bias in the unadjusted estimate appears evident from the sub-regional dynamics of two covariates (see Figure 4). Compared to the control sub-region (Ifjord), the peak of the last rodent cycle was substantially lower and the ptarmigan harvest rate was generally higher in the action sub-regions (Varanger).
These differences in two influential drivers of sub-regional ptarmigan population dynamics appeared to diminish the unadjusted impact estimate of the management action.
Compared to the unadjusted contrast estimate before the action (0.24, 95% CI: [−0.54, 1.01]), the covariate-adjusted contrast estimate after the action (estimate: 1.45, 95% CI: [0.91, 1.99]) was six times larger. While we could not provide a covariate-adjusted contrast estimate before the management action due to lack of pre-action monitoring data, the difference between the (unadjusted) before and (adjusted) after contrast estimates is so large that it is likely that most of the difference can be attributed to the impact of the red fox culling action. Also, notice that coefficients of the standardized covariates were much smaller than the estimated effect of the red fox culling action ( Figure 2). Finally, support for a true effect of culling action was also evident from the estimated effect of the spatio-temporal variation in number of culled foxes within the action region (i.e. Step 3 of the analysis).
Meta-analyses of predator removal studies (Kämmerle & Storch, 2019;Salo et al., 2010) show that the evidence for effects of removal on prey species varies among studies. Salo et al., (2010) found that studies that often failed to find significant effects; (a) were based on un-replicated/non-experimental management actions, (b) targeted prey species with complex population dynamics (e.g. population cycles) and/or (c) involved the removal of a single predator species in complex food web settings. Additionally (d), Côté and Sutherland (1997) noted that predator removals mostly failed to halt the decline of populations that were already declining. It is interesting to note that this study possessed all of these four characteristics and yet was still able to find evidence for a strong impact of the red fox culling action.
Besides strengthening the impact assessment by accounting for influential covariates in the analysis, our case study benefitted from the large spatial and temporal scale of the management action, which have been identified as important success factors (e.g. Conner & Morris, 2015). Regarding spatial scale, over 90% of studies reviewed in Salo et al. (2010) took place on areas <100 km 2 , while our management region was ~3,500 km 2 . Spatially extensive action areas are likely particularly important in the context of wide ranging mesopredators, such as boreal red foxes (Barton & Zalewski, 2007;Norén et al., 2017;Walton et al., 2018). Spatial scale becomes an issue also when density-dependent dispersal of the management targets cancel the density contrast between small control and action areas (Côté & Sutherland, 1997).
Regarding temporal scale, we had access to 17-year ptarmigan time series (5-years before and 12 years after the management action), while most BACIPS assessments last only 2-5 years (De Palma et al., 2018). Long time series, in particular after the onset of actions (Thiault et al., 2017), are generally important for robust impact assessments, but are particularly crucial in the context of populations and food webs with complex dynamics. The 12-year action period in this case study included three 4-year rodent cycles. The amplitude of these cycles, which is known to affect ptarmigan dynamics mediated by predator functional and numerical response to rodent abundance (Angelstam et al., 1984;Lack, 1954;Steen et al., 1988), varied in both time and space. Having sufficiently long time series at the scale of sub-regions allowed us to account for this source of complexity and thereby improve the assessment of the culling impact.
Although the red fox is only one of the predators that link the population dynamics of ptarmigan and rodents (Henden et al., 2009), there was still a significant signal of the rodent cycle in the ptarmigan dynamics despite a large number of culled red foxes in the action areas. Marolla et al. (2019) found that a similar red fox culling action, that aimed to conserve the lesser whitefronted goose in sub-arctic tundra ~200 km west of our study area, did not diminish the strong signal of the rodent cycle in goose breeding success. As pointed out by Salo et al. (2010) the effect of removing a single predator species can be expected to be at least partly compensated for by other predators, especially in predator-rich tundra food webs . While we conclude that the red fox culling action had a clear positive effect by curbing the population decline seen in most places in Finnmark , it did not lead to a population increase.
A likely reason is that several other generalist predators like corvids exert increased impacts on alpine-arctic birds in a warming climate . An important challenge for future studies is to obtain direct estimates of how populations of red fox and other members of arctic predator communities are impacted by climate change and management interventions.

| Management implications
Our study provides the first evidence that the expansion of red fox in the Arctic, at least regionally, is playing a role in the decline of bird populations belonging to the tundra biome. We also show that actions to decimate the red fox may be able to curb the population decline of an ecologically and socio-economically important species-i.e. the willow ptarmigan-given that such management actions are large scale and long term.
More generally, this study provides a strong case for the value of ecosystem-based monitoring with a food web approach (cf. Henden et al., 2017; in the context of management actions aimed to conserve wildlife. Species of conservation concern are often subjected to strong biotic interactions because of their central position within the food web. Such interactions vary in space and time in a manner that may cause confounding between management actions and other drivers of change, especially when management interactions are large scale and lack spatial replication. We suspect that the variable outcome of predator removal assessments may be due to the lack of adequate monitoring (sensu Lindenmayer & Likens, 2010). Monitoring drivers of food web interactions, so they can be accounted for when making model-based assessments about the impact of management actions, can profoundly increase the accuracy of management impact assessments.

ACK N OWLED G EM ENTS
This study was supported by the RCN funded project SUSTAIN and COAT-Climate-ecological Observatory for Arctic Tundra. The red fox culling and data collection was financed by the Norwegian Environmental Agency. The willow ptarmigan monitoring is organized by the Finnmark Estate (FeFo) and Hønsefuglportalen (http:// honse fugl.nina.no/Innsy n/) and conducted by volunteer hunters with pointing dogs. We thank Einar Asbjørnsen (FeFo) for access to harvest statistics and Hønsefuglportalen for access to the ptarmigan line transect data for Finnmark County, and all COAT field assistants who have participated in collecting small rodent data. We also thank Jarad Pope Mellard for constructive comments and grammatical proof-reading.

AUTH O R S ' CO NTR I B UTI O N S
All authors contributed to conceive the ideas and collect the food

DATA AVA I L A B I L I T Y S TAT E M E N T
Data available via DataverseNO, UiT Open Research Data https:// doi.org/10.18710/ 2ZG3EI .