Preferential cannibalism as a key stabilizing mechanism of intraguild predation systems with trophic polymorphic predators

Theory predicts intraguild predation (IGP) to be unstable despite its ubiquity in nature, prompting exploration of stabilizing mechanisms of IGP. One of the many ways IGP manifests is through inducible trophic polymorphisms in the intraguild (IG) predator, where a resource-eating predator morph competes with the intraguild (IG) prey for the shared resource while a top predator morph consumes the IG prey. Cannibalism is common in this type of system due to the top predator morph’s specialization on the trophic level below it, which includes the resource-eating predator morph. Here, we explore the consequences of inducible trophic polymorphisms in cannibal predators for IGP stability using an IGP model with and without cannibalism. We employ linear stability analysis and identify regions of coexistence based on the top predator morph's preference for conspecifics vs. heterospecifics and the IG prey's competitive ability relative to the resource-eating morph. Our findings reveal preferential cannibalism (i.e. the preferential consumption of conspecifics) stabilizes the system when the IG prey and resource-eating morph have similar competitive abilities for the shared resource. Though original IGP theory finds the IG prey must be a superior resource competitor as a general criterion for coexistence, this is not typically the case when the predator has an inducible trophic polymorphism and the resource-eating morph is specialized in resource acquisition. Preferential cannibalism may therefore be a key stabilizing mechanism in IGP systems with a cannibalistic, trophic polymorphic IG predators, providing further insight into what general mechanisms stabilize the pervasive IGP interaction.

Methods
Model(s) overview: To examine the effect of preferential cannibalism in a trophic polymorphic predator on IGP system stability, we compared two models of varying complexity (Fig. 1) under two scenarios pertaining to IG prey competitive ability. The first model (referred to herein as the “base” model) is an extension of the original Lotka-Volterra IGP model first proposed by Holt and Polis (1997) with the separation of the IG predator into two states: a resource-eating morph that competes with the IG prey for the shared resource, and a top predator morph that consumes the IG prey. Biomass moves from one state to the other as a function of resource density, which is intentionally general to encompass changes in frequency of morphs across generations or individuals switching between morphs in a lifetime. In the second model (referred to herein as the “full” model), we build upon the first model to include cannibalism in the IG predator, where the top predator morph consumes both the IG prey and the resource-eating predator morph. We further include a preference parameter, s, that controls the top predator morph’s preference for conspecifics (resource-eating predator morph) or heterospecifics (IG prey). We explore preference over a range of s values, specifically three values of s which represent preference for conspecifics (s=0.7), preference for heterospecifics (s=0.3), and no preference (s=0.5).

Analysis: We explored stability in our models using linear stability analysis (Gurney and Nisbet 1998; Murdoch et al. 2003; See SI 1). In short, we linearized eqs. 1 and 2 around their interior equilibrium, and then examined the stability of these systems to small perturbations. We first solved for the equilibrium in which all species have positive, non-zero abundances using the “Solve” function in Mathematica (Wolfram Research, Inc., Mathematica, Version 12.0, Champaign, IL, 2019). We then evaluated the Jacobian matrix at this solution and numerically computed the eigenvalues using the “Eigenvalues” function, selecting the eigenvalue(s) in which the real part is negative and therefore stable. We further performed numerical simulations to show the effects of preference and competitive superiority on dynamics and equilibrium densities in SI 2. Dynamics were simulated in the R programming language (R Core Team 2023) using the package “deSolve.”