Derivation of Predator Functional Responses Using a Mechanistic Approach in a Natural System

1. Introduction

A long-standing problem in ecology is to measure how the acquisition rate of a predator varies with prey availability, namely the functional response. Functional response shapes are typically categorized as linear (type I), hyperbolic (type II), or sigmoidal (type III; Holling 1959b,a). This classification is commonly used by ecologists when incorporating predation into population and community models (Turchin and Hanski, 1997; Fryxell et al., 2007; Serrouya et al., 2015), and type II is the most widely applied model (Rall et al., 2012). The shape of the functional response can have major consequences on the outcomes of population and community models. For instance, a type III promotes stability or coexistence whereas a type II destabilizes predator-prey dynamics (Murdoch, 1973; Sinclair et al., 1998). Describing the functional response of pairwise trophic interactions is also important to understand higher-order interactions. For instance, the shape of the functional response alone can profoundly change predictions about the outcome of predator-mediated trophic interactions (Abrams et al., 1998; Holt and Bonsall, 2017).

Although functional responses are at the core of predator-prey theory (Solomon, 1949), most empirical research on functional responses has been conducted under controlled laboratory or field enclosure conditions (96%, n = 116 studies, reviewed by Pawar et al. 2012) where prey density is manipulated, predator consumption is recorded, and the functional response models are compared through statistical analysis (referred to as phenomenological models). These approaches are not well-suited for predators with relatively large home ranges and may fail to reproduce foraging conditions encountered in the wild. Determining the shape of functional responses in natural systems is however often limited by a combination of factors, including small sample sizes, a relatively narrow gradient of observed prey densities, the difficulty to observe predator-prey interactions directly, or the difficulty to estimate predator and prey numbers (Gilg et al., 2006; Therrien et al., 2014; Suryawanshi et al., 2017; Ellis et al., 2019). The large variability around predator acquisition rates observed in the wild can also prevent us from discriminating among functional response shapes, and hence limit our ability to accurately model predator-prey interactions in complex and natural ecosystems (O'Donoghue et al., 1998; Vucetich et al., 2002; Chan et al., 2017). Moreover, phenomenological models fail to identify the proximate mechanisms regulating predator acquisition rates. In this context, a mechanistically grounded approach to derive functional response is appealing.

Derivation of functional responses based on measurable features of species behavior (e.g., speed, attack and success probability) provides several advantages. Compared with phenomenological models, mechanistic models (1) allow assessing the shape of the functional response based on behavioral attributes of the predator, (2) are based on parameters with a direct biological interpretation, and hence have the potential to reinforce links between theory and data (Connolly et al., 2017). The number of mechanistic models of predator-prey interactions is growing, and most of them aim to predict trophic links based on species traits, especially body size (Gravel et al., 2013; Ho et al., 2019; Portalier et al., 2019). Mechanistic models of functional response further allow the integration of predator-prey pairs to describe trophic links, which can improve our ability to model complex ecological interactions. Despite their utility, mechanistic models based on the main components of predation in a natural system are, to our knowledge, virtually non-existent in vertebrate predators.

Using a combination of behavioral, demographic, and experimental data collected over 20 years in a natural system, we develop a mechanistic model to assess the functional response of a generalist mammalian predator to various prey species (4 predator-prey pairs). The originality of our approach is to assess functional response (i) by breaking down the components of predation (searching, chasing, capturing, and handling prey) and (ii) by using field experiments and detailed behavioral observations to parameterize each component included in the mechanistic model. We focused on the derivation of functional response of predator-prey pairs. We evaluated the coherence of our models using data from a long-term field study that estimated prey densities and predator acquisition rates. We also performed sensitivity analyses to identify the main proximate drivers of change in predator acquisition rates. Finally, we modeled the potential effects of density dependence in components of predation on the shape of the functional responses within the range of prey densities observed in the field.

The mechanistic model was developed for the arctic fox (Vulpes lagopus), a generalist predator of the tundra ecosystem, using highly detailed empirical observations from a long-term ecological monitoring program in the Arctic (Gauthier et al., 2013). This system offers several benefits to study predator-prey interactions among vertebrates, including a relatively simple food web, an open landscape and the continuous summer daylight allowing direct behavioral observations. The arctic fox is an active hunting predator that travels extensive daily distances within its territory in summer (Poulin et al., 2021). Lemmings and birds (mostly eggs and juveniles) are the main components of the summer diet of arctic foxes in most tundra ecosystems (Angerbjörn et al., 1999; Giroux et al., 2012). Lemmings exhibit population cycles with peak density every 3–5 years (Fauteux et al., 2015), and the arctic fox predation pressure on tundra ground-nesting birds is typically released at high lemming density (Summers et al., 1998; Bêty et al., 2002; McKinnon et al., 2014). Surprisingly, the exact mechanisms driving this well-known short-term apparent mutualism between lemmings and birds are still unclear, but they likely involve fox functional responses (Summers et al., 1998; Bêty et al., 2002).

A few studies attempted to quantify the functional responses of arctic fox using phenomenological models (Angerbjörn et al., 1999; Eide et al., 2005; Gilg et al., 2006). Relatively low sample sizes reduced the ability of previous studies to fully distinguish between different shapes of functional responses. Moreover, the hoarding behavior of arctic foxes was not considered in previous estimations of functional responses (Angerbjörn et al., 1999; Eide et al., 2005; Gilg et al., 2006). Like many other animals (Vander Wall, 1990), arctic foxes can predate more prey than they consume on the short-term, and such behavior can strongly increase prey acquisition rates, e.g., foxes foraging in goose colonies can hoard between 40 and 97% of eggs acquired during the bird nesting period (Samelius and Alisauskas, 2000; Careau et al., 2008). Although type III functional responses were previously used to model fox-prey population dynamics (Gilg et al., 2003, 2009), food hoarding may substantially reduce handling time and could therefore make the shape of the functional response linear or slightly convex (Oksanen et al., 1985).

2. Methods

2.1. Study System

During the summer, the southwest plain of Bylot Island, Nunavut, Canada (73° N; 80° W) harbors a large greater snow goose colony (Anser caerulescens atlanticus; ~20,000 pairs). Insectivorous migratory birds are also nesting in the study area and include the lapland longspur (Calcarius lapponicus), a passerine, and several species of shorebirds (primarily Calidris spp. and Pluvialis spp.). Two species of small mammals are present, the brown (Lemmus trimucronatus) and collared (Dicrostonyx groenlandicus) lemmings. The brown lemming has high-amplitude cycles of abundance with a 3–5-year periodicity, whereas the collared has low-amplitude cycles (Gruyer et al., 2008). The mammalian predator guild is dominated by the arctic fox and the ermine (Mustela erminea). The arctic fox is the main nest predator of geese (Bêty et al., 2002; Lecomte et al., 2008), sandpipers (McKinnon and Bêty, 2009; Royer-Boutin, 2015), and passerines (Royer-Boutin, 2015). Additional details on plant communities and general landscape can be found in Gauthier et al. (2013).

The model was parameterized and evaluated using data from Bylot Island, where foxes and their prey have been monitored since 1993. We observed foraging foxes using binoculars and spotting scopes (20 × 60x) from one or two blinds located in the middle of the goose colony during 10 summers between 1996 and 2019.

2.2. Mechanistic Model of Functional Responses

We used the Holling disk equation as a starting point to build the mechanistic model of functional response (Holling, 1959a) inspired by the general formalism of Pawar et al. (2012). Predation was broken down into four different components, which are searching, chasing, capturing, and handling of a prey item by a predator. Acquisition rate of a prey item (species i) by a predator [f(i)], namely the functional response, takes the following form:

where α_i is the capture efficiency (km²/h), N_i the prey density (number of prey/km²), and h_i the handling time of prey (h/prey). Capture efficiency is obtained by the product of predator speed (s; km/h), reaction distance (d_i; km), detection (z_i), and attack probability (k_i) of the prey by the predator, and the success probability (p_i) of an attack (Table 1):

Table 1. Definition of the parameters used in the functional response model.

The combination of the time spent chasing the prey once encountered ($\frac{T_{ci}}{p_i}$) and the time spent manipulating the prey once subdued (T_mi) define an overall prey handling time (h_i):

The time spent manipulating includes the time spent eating or hoarding the prey item.

α_i depends only on prey density, and we assumed that prey are randomly distributed. Satiety was not considered as a potential mechanism limiting acquisition rate. Indeed, foxes can predate more prey than they consume on the short-term; e.g., about 4% (n = 128) and 48% (n = 98) of predated eggs and lemmings are immediately eaten, respectively (Careau et al., 2007). Predator interference was not incorporated in the model as foxes rarely encounter and interact with other individuals while foraging within their summer territory (49 interactions, which represents 0.9% of the time over 118 h of direct observations of foxes foraging in the study area). The full model derivation is provided in Supplementary Material (section 1.1).

The general model of functional response (Equation 1) allows for a continuum between a linear (type 1) and a concave (type 2) functional response shape as a linear response can be elicited when handling time is negligible. In order to allow the mechanistic model to extend to a sigmoidal shape (type 3), we added density dependence in capture efficiency components that were expected to vary with prey density (i.e., reaction distance and detection and attack probabilities; see below).

2.3. Prey Specific Functional Responses

We adapted the general model (Equation 1) to each prey species based on their traits and anti-predator behavior (Figure 1). The specific models for each prey species are provided in Supplementary Table S2.

Figure 1. Conceptual mechanistic model of functional response of arctic fox to each prey species: attended (A) or unattended (B) goose nests, lemmings (C), and passerine and sandpiper nests (D). Predation was divided into four main components: search, detect, attack, and handle (which includes eating and hoarding). Arrows illustrate the probability that the predator reaches the next component. When there is no parameter indicated beside the arrow, the probability is 1. Parameters are as follows: d is the reaction distance, s the predator speed, z the detection probability, k the attack probability, p the success probability, P_c the probability of complete nest predation, T_C the time spent chasing, and T_m the time spent manipulating the prey.

For goose nests, the first modification was to add a component for complete and partial nest predation. This modification was necessary since a successful attack by the predator does not always result in complete clutch predation (Bêty et al., 2002), which can affect manipulation time and, ultimately, acquisition rates. The second modification was to split the general model into two components. A first component models acquisition rate of goose nests when the female is incubating or when one protecting adult is at <10 m from the nest (attended nest; Figure 1A). A second component models acquisition rate of goose nests during incubation recesses when both adults are at >10 m from the nest (unattended nest; Figure 1B). As geese can actively protect their nests against arctic foxes, their presence at the nest strongly influences fox foraging behavior (Samelius and Alisauskas, 2001; Bêty et al., 2002). This anti-predator behavior translates into changes in capture efficiency components. Thus, parameter values of capture rates were estimated separately for goose nests that were attended or unattended (Supplementary Table S1). When a nest is attended by a highly conspicuous snow goose, we assumed that nest detection probability is 1 within d (Figure 1A). For unattended nests, we used a detection probability function obtained from an artificial nests experiment (Supplementary Figure S1). Sometimes, unattended nests can be protected if parents detect a fox during an incubation recess and return quickly to their nest. Like attended nests, we thus estimated success probability (p) and complete clutch predation probability (P_c) for unattended nests (Figure 1B). The third and last modification was to introduce the nest attendance probability (w). We estimated this parameter by combining information on the average time spent on the nest by females and on the average distance between females and their nest during the goose incubation period (Reed et al. 1995; Poussart et al. 2000; see Supplementary Material, section 1.2).

The general model (Equation 1) was simplified for lemmings as we assumed that an attack is systematically initiated by the fox once a lemming is detected within d (Figure 1C). Sandpipers and passerines exhibit a variety of antipredator behaviors (such as distraction displays) to avoid nest detection by the predator (Smith and Edwards, 2018). However, sandpipers and passerines cannot protect their nest once detected by a fox. Thus, we assumed that once a nest is detected, it is always predated (attack probability is 1 and no chasing time is included in the model; Figure 1D).

We incorporated density dependence into the goose and the lemming models within the range of densities observed in our study system. For each parameter in which density dependence was incorporated, the minimum and the maximum parameter values were associated, respectively, with the minimum and the maximum prey density to calculate the slope and the intercept of the density-dependence relationship. In the goose model, we modified attack and success probabilities for attended nests, and reaction distance and detection probability for unattended nests. In the lemming model, we added density dependence in reaction distance, detection, and success probabilities. The rationale behind these additions is that predators may form search images for abundant prey, which can increase their ability to detect them (Nams, 1997; Ishii and Shimada, 2010). As predators could also increase their attack rate and success as prey density increases, we added density dependence in attack and success probabilities. We did not incorporate density dependence into the passerine and sandpiper nest models as the range of nest densities observed in our study system is likely too low to influence fox behavior (maximum of 12 nests/km² compared to a maximum of 926 goose nests and 414 lemmings per km²). See Supplementary Material (section 1.3) for more details on the incorporation of density dependence.

The model was implemented in R v. 3.6.0 (R Core Team, 2019).

2.4. Parameter Values

The model was parameterized mostly using data from Bylot Island but also from the literature when data were missing. Parameters were derived from field experiments using artificial nests or estimated using arctic fox GPS tracking data and direct observations of foraging foxes (Supplementary Table S1). See Supplementary Material (section 1.2) for a detailed description of the method used to extract each parameter.

2.5. Evaluating the Coherence Between the Mechanistic Model and Empirical Predator Acquisition Rates

Predator acquisition rates at different prey densities were assessed in the field annually using two independent methods. These data did not allow validation of the shape of the functional responses, but they provided a way to evaluate the performance of the mechanistic model in estimating prey acquisition rates at the various prey densities observed in our study system.

First, we obtained goose eggs and lemming acquisition rates by conducting direct observations of foraging foxes for 10 summers between 1996 and 2019 during the goose incubation period (details on behavioral observations can be found in Bêty et al., 2002; Careau et al., 2008). For each year, the acquisition rate was calculated as the total number of prey acquired (goose eggs or lemmings) divided by the total length of the observation bouts of individual foxes. The acquisition rate of a clutch of eggs was estimated by dividing the acquisition rate of goose eggs by the annual average clutch size. For the years where information was available, we also calculated the acquisition rate for attended and unattended nests. We estimated annual goose nest density either by visual counts of the nests located in the observation zone (range: 0.5–3 km²) during the incubation period (1996–1999, 2019) or over a fixed 0.2 km² plot within the intensively monitored core area of the goose colony (2004–2005, 2015–2016). We estimated lemming density annually with snap traps from 1994 to 2009 and with live traps from 2004 to 2019 (see Fauteux et al., 2018 for methods). We summed the density estimate of brown and collared lemming.

Second, we obtained passerine and sandpiper nest acquisition rates by monitoring annually (2005–2013) the fate of passerine and sandpiper nests (Gauthier et al., 2013; McKinnon et al., 2014). Nest density was estimated as the number of passerine and sandpiper nests found in a 8 km² plot systematically searched throughout the breeding season. We estimated acquisition rate of nest content (eggs or chicks) by using the daily survival rate of nests (dsr), the total number of nests found in the study plot (N_tot), the number of foxes foraging in the plot (N_fox), and the proportion of nests predated by foxes (P_fox). Since foxes establish territorial pairs on Bylot (Rioux et al., 2017), we assumed that 2 foxes were foraging in the study plot. We also considered that foxes were responsible for 100% (n = 19) and 81% (n = 25) of the failed sandpiper and passerine nests, respectively, as indicated by camera monitoring (McKinnon and Bêty, 2009; Royer-Boutin, 2015). An estimation of the acquisition rate is obtained by:

The daily nest survival rate was modeled using the logistic exposure method (Shaffer, 2004). Additional details on daily nest survival rate calculations and nest monitoring methods can be found in Royer-Boutin (2015). Density estimates for all prey species were standardized as the number of nests per km².

2.6. Uncertainty and Sensitivity Analysis

We quantified how uncertainty in parameter values affected estimation of predator acquisition rates by using the Latin hypercube sampling technique (an efficient implementation of the Monte Carlo methods; Marino et al., 2008). This analysis allowed us to investigate the uncertainty in the model output generated by the uncertainty and variability in parameter inputs. Each parameter was represented by a probability distribution (uniform or normal truncated) based on the distribution of empirical data (Supplementary Table S1). For some parameters, the biological information was limited, so we assigned a uniform distribution allowing for a large range bounded by minimum and maximum values. Latin hypercube sampling was then applied to each distribution (N = 1,000 iterations). This method involved dividing a probability distribution into N equal probability intervals that were then sampled without replacement, resulting in N iterations of the model using each combination of parameters values. This method allowed us to explore the entire range of each parameter and most of them encompass various environmental conditions (e.g., weather conditions, prey availability). We computed the median, the 90th, 95th, and 99th percentiles of the model output by using the empirical cumulative distribution.

We also conducted a local sensitivity analysis to identify key parameters of the mechanistic models within the range of prey densities observed in our study system. We modified each parameter value by ±100% while holding others constant, and we assessed how this variation affected the predator acquisition rate (expressed as % of change).

3. Results

From 1996 to 2019, we observed foraging foxes in the goose colony for 124 h. Average goose nest density was 409 nests/km² (range: 100–926 nests/km²) and lemming density was 193 ind./km² (range: 11–414 ind./km²; Supplementary Table S3). Based on direct observations of foraging foxes, average acquisition rates were 0.61 nest/fox/h (range: 0.19–1.82 nest/fox/h) for goose nests and 0.94 ind./fox/h (range: 0–2.85 ind./fox/h) for lemmings (Supplementary Table S3). The majority of eggs acquired (67%) were from unattended goose nests, while 33% were from attended nests (n = 218). Average passerine nest density was 7.7 nests/km² (range: 6.1–12.3 nests/km²), and sandpiper density was 2.5 nests/km² (range: 1.0–5.9 nests/km²; Supplementary Table S4). Based on nest monitoring, average acquisition rates were 0.10 nest/fox/h (range: 0.03–0.28 nest/fox/h) and 0.04 nest/fox/h (range: 0.002–0.169 nest/fox/h) for passerine and sandpiper nests, respectively (Supplementary Table S4).

The uncertainty analysis revealed that varying simultaneously all parameters used in the mechanistic model generated considerable variation in fox acquisition rates (Figure 2). Nonetheless, no parameter combinations resulted in a saturating functional response for all prey species within the range of prey densities observed in our study system: the acquisition rate at maximal prey density was below the saturation point in all simulations (see histograms in Figure 2, Supplementary Figure S5). Based on the value of the parameters estimated within the observed prey densities, acquisition rate at saturation was 8 nests/fox/h for goose nests, 17 ind./fox/h for lemmings, 166 nests/fox/h for passerine, and 26 nests/fox/h sandpiper nests. Depending on the prey species, most or all fox acquisition rates observed in the field fell within the 99th percentile of the values derived from the mechanistic models (Figure 2). The highest acquisition rates observed in the field were also much below the estimated saturation point in all prey species (Figure 2). As the goose nest model was split for attended and unattended goose nests, we also computed acquisition rates separately for each of these situations. Goose nest acquisition rate derived from the mechanistic model was higher for unattended nests than attended nests, which is consistent with empirical estimations (Supplementary Figure S6). Although most (66%) field estimates of acquisition rates fell within the 95th percentiles of the model output for unattented goose nests, all values were under the model median at nest densities above 200 nests/km² (Supplementary Figure S6). This may indicate a slight overestimation of the proportion of unattended nests at relatively high densities.

Figure 2. Functional response of arctic fox to density of goose nests (A), lemmings (B), passerine nests (C), and sandpiper nests (D). Black lines represent the median of the mechanistic model and the color bands represent the 90th, 95th, and 99th percentiles based on 1,000 simulations. Empirical data are represented by red and yellow dots respectively. Histograms in the inset show the distributions of acquisition rate at saturation for each simulation. Horizontal error bars in (A) indicate the range of nest density during the incubation period. Vertical errors bars in (A,B) represent standard errors calculated using bootstrapping. Errors bars in (C,D) represent 95% confidence intervals from daily survival rate estimates.

Sensitivity analyses indicated that predator speed was an influential parameter of the functional response of all prey species (Figure 3). Goose nest acquisition rate was generally more affected by parameters associated with unattended nests than attended nests (Figure 3A). The magnitude of change in goose nest acquisition rate related to the changes in manipulation time increased slightly with nest density. Lemming acquisition rate was not affected by chasing and manipulation time, whereas detection distance, and detection and success probability had an influence equivalent to predator speed (Figure 3B). Similarly, functional response models of passerine and sandpiper nests were not sensitive to change in manipulation time, whereas detection distance and detection probability had an influence equivalent to predator speed (Figure 3C).

Figure 3. Sensitivity of predator acquisition rates to changes in parameter values of the mechanistic models used to assess the functional response of arctic fox to goose nests [at 100 (A1), and 1,000 nests/km² (A2)], to lemmings at 250 ind./km² (B) and to passerine and sandpiper nests at 10 nests/km² (C). Sensitivity analyses are presented at intermediate densities for lemmings and passerine and sandpiper nests, since the results were very similar within the range of densities observed in our study system.

Adding density dependence into the goose and the lemming models had relatively minor effects on acquisition rates derived for low to moderate densities observed in our study system. The shape of the functional response changed slightly between models without or with density dependence in capture efficiency components (allowing for a gradient between type I and type III). At high densities, acquisition rates remained much below saturation points, and a maximum difference of 1.4 nests/fox/h at 1,000 goose nests/km² and 2.1 lemmings/fox/h at 450 lemmings/km² were found between models (Figure 4).

Figure 4. Functional response models of arctic fox to goose nests (A) and lemmings (B) with and without density dependence on capture efficiency components, within the range of densities observed in the field. Black lines represent the median of the mechanistic model and the color bands represent the 95th percentiles based on 1,000 simulations.

4. Discussion

Benefiting from a combination of behavioral, demographic, and experimental data collected over the past 20 years, we developed a mechanistic model of arctic fox functional response to four prey species. Our model derives the shape of the functional response of each predator-prey pair along a gradient from linear to sigmoidal. Predator acquisition rates derived from the mechanistic model were consistent with field observations, and the main proximate mechanisms driving predator acquisition rates were also identified. In all prey species, predator speed was an influential parameter, while handling time had a limited influence on acquisition rates. Although type III functional responses were previously used to model fox-prey population dynamics (Gilg et al., 2003, 2009), our simulations indicate that predator acquisition rate was not systematically limited at the highest prey densities observed in our study system. Our model allows for a mechanistic interpretation of the functional response of predator-prey pair and could be extended to more complex modules involving multiple predators and prey species.

Our results add to a growing body of research indicating that predators may not become systematically satiated or saturated at the highest densities of prey observed in nature (Novak, 2010; Chan et al., 2017; Preston et al., 2018). Holling's functional response models (type II and III), which are commonly used in population dynamics models (Turchin and Hanski, 1997; Gervasi et al., 2012; Serrouya et al., 2015), predict that predator acquisition rates should eventually saturate at high prey densities. Based on mechanistic models, which allowed us to vary simultaneously all components of predation, we found no evidence of arctic fox saturation at the highest prey densities observed in a natural system. Several factors may explain this result. First, the hoarding behavior of arctic foxes may substantially reduce handling time by limiting the constraints associated with digestion and satiety, which can make the functional response shape linear or slightly convex even at high prey densities (Oksanen et al., 1985). Second, while predator acquisition rates must theoretically become constrained by handling and/or digestion at high prey densities, the prey densities required to reach a saturation point could be rarely observed in natural systems. Indeed, empirical support for saturating functional response in the wild is relatively rare and comes mostly from controlled laboratory experiments in which the range of prey densities may exceed the range observed in nature [99% of all type II functional response were derived from controlled laboratory experiments (n = 61 studies); reviewed by Rall et al., 2012]. Such an issue can be avoided when mechanistic approaches are used to derive functional responses. One particularity of our system is the presence of a large goose colony where prey density can be quite high (up to ~900 nests/km²). Interestingly, even in this context, we found no evidence of predator saturation.

Historically, a categorical approach was adopted by ecologists to define functional responses. A linear functional response was traditionally attributed to filter feeders (Jeschke et al., 2004), a hyperbolic shape (type II) to invertebrates and a sigmoidal shape (type III) to vertebrate predators (Holling 1965, but see Hassell et al. 1977). Although, this categorization has some heuristic value in introductory texts and can be useful in some aspects of research where categorization is necessary, types I, II, and III should be considered simply as particular cases along a continuum. Instead of using a priori shapes to describe functional responses, our study illustrates how mechanistic models can generate functions linking prey density and predator acquisition rates that are specific, and hence more relevant, to the range of densities observed in a food web. Considering the strong effect of functional responses on the outcome of predator-prey models (Abrams et al., 1998; Sinclair et al., 1998), such specific functions should improve our ability to adequately simulate and quantify the strength of species interactions in natural communities.

We did not incorporate predator dependence in the functional response model, despite a growing body of studies indicating that some mechanisms (e.g., facilitation, interference) are likely to occur in functional responses (Novak et al., 2017). However, arctic foxes maintained summer territories (averaging 9.6 km²) with low overlap (Grenier-Potvin and Berteaux, submitted manuscript), which prevents potential interference within territories. We are thus confident that variation in predator density should not affect our main conclusions. Nonetheless, the mechanistic model could be extended to more complex predator-prey systems, including predator interference.

Habitat characteristics could affect several parameters of the mechanistic model, hence the functional response shape and magnitude could be modulated by the structural complexity of the landscape (Toscano and Griffen, 2013; Barrios-O'Neill et al., 2015). For instance, the detection distance of a nest by arctic foxes could be lower in dense vegetation (Flemming et al., 2016), the attack probability could be lower for nests located in wetlands and islets only accessible by swimming (Lecomte et al., 2008; Gauthier et al., 2015), and the success probability of an attack could be modulated by the presence of complex networks of lemming tunnels offering refuges. Exploration of the effects of structural complexity on functional responses remains rare (but see Lipcius and Hines 1986; Toscano and Griffen 2013; Barrios-O'Neill et al. 2015), and more empirical research is needed to integrate these sources of variation in mechanistic models.

The outputs of the mechanistic model were generally consistent with field observations. However, adding more complexity could improve its performance and our ability to identify the main drivers of predator acquisition rates. For instance, group defense and mutual vigilance are additional factors that may reduce predator acquisition rates at high prey density (Clark and Robertson, 1979). Although there is no evidence of group defense in geese (Bêty et al., 2001), the snow goose could benefit from the vigilance and early warning provided by neighbors nesting nearby (Samelius and Alisauskas, 2001). Beyond a threshold of goose nest density, such anti-predator behavior could reduce the proportion of unattended nests with increasing nest densities. Nest attendance probability was an influential parameter of the goose model and mutual vigilance may partly explain why acquisition rates observed in the field at moderate-high nest densities were under the model median (Figure 2A).

One mechanism often advanced for explaining the apparent mutualism between two prey sharing a common predator is predator saturation or satiation (Holt, 1977; Abrams and Matsuda, 1996). Our results showed that the arctic fox doesn't reach saturation at the highest lemming densities observed in our study system. This suggests that the underlying mechanism for the short-term positive effect of high lemming density on arctic bird reproductive success (Bêty et al., 2002; Blomqvist et al., 2002) is likely not predator satiation nor saturation. Instead, the apparent mutualism between birds and lemmings could arise from changes in other components of the functional response. For instance, the attack probability of an attended goose nest could be inversely dependent of lemming density, or daily distance traveled by the predator (speed) could be dependent of lemming density. As indicated by our sensitivity analyses, attack probability was not a strong driver of prey acquisition rates while predator speed was an influential parameter in all prey species. Hence, lemming-induced changes in predator speed through changes in predator activity budget (e.g., due to predator reproductive status or hunger level) could be an alternative hypothesis explaining the apparent mutualism between lemmings and arctic birds. Density dependent changes in components of the functional response have been observed in other systems and can generate nonlinearities in the functional response (Hassell et al., 1977; Abrams, 1982). The integration of all prey species into a mechanistic multi-species functional response model is the next step to fully identify the main proximate drivers of indirect interactions in natural communities.

5. Conclusion

Previous studies of functional responses typically tried to discriminate between predetermined shapes of functional responses. Our study illustrates how mechanistic models based on empirical estimates of the main components of predation can generate functional responses specific to a range of prey densities relevant to a given food web. Such mechanistically derived functional responses are needed to untangle proximate drivers of predator-prey population dynamics and to improve our understanding of predator-mediated interactions in natural communities. Although it would be unrealistic to resolve every pairwise interaction within ecological networks, our mechanistic model provides a starting point for studying higher-order effects such as indirect interactions that can emerge among prey species.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Ethics Statement

The animal study was reviewed and approved by Capture techniques and immobilization procedures were approved by the UQAR Animal Care Committee (CPA-64-16-169 R2).

Author Contributions

AB, DG, and JB conceived the study. AB, DG, DB, GG, JC, VC, NL, C-CJ, PR-B, and JB collected the data. AB, PR-B, and JC analyzed the data. AB wrote the manuscript with contribution from all authors. All authors contributed to the article and approved the submitted version.

Funding

The research relied on the logistic assistance of the Polar Continental Shelf Program (Natural Resources Canada) and of Sirmilik National Park of Canada. The research was funded by (alphabetical order): Arctic Goose Joint Venture, the Canada Foundation for Innovation, the Canada Research Chairs Program, the Canadian Wildlife Service, the Fonds de recherche du Québec-Nature et technologies, the International Polar Year program of Indian and Northern Affairs Canada, the Natural Sciences and Engineering Research Council of Canada, the ArcticNet Network of Centers of Excellence, the Northern Ecosystem Initiative Program (Environment Canada), the Northern Scientific Training Program, the Nunavut Wildlife Management Board, Polar Knowledge Canada, Université du Québec à Rimouski, Université Laval and the W. Garfield Weston Foundation.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We are especially grateful to the many people who helped us with field work over many years, the Mittimatalik Hunters and Trappers Organization and Park Canada's staff for their assistance. We sincerely thank both reviewers for their constructive and relevant comments. The article was published on BioRxiv (Beardsell et al., 2020).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2021.630944/full#supplementary-material

References