Strengthening intraguild predation increases the temporal variability of biomass across all trophic levels in model food webs

Strengthening intraguild predation increases the temporal variability of biomass across all trophic levels in model food webs | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Strengthening intraguild predation increases the temporal variability of biomass across all trophic levels in model food webs Chase J. Rakowski , Mathew A. Leibold , Caroline E. Farrior doi: https://doi.org/10.1101/2025.06.25.661600 Chase J. Rakowski 1 Department of Integrative Biology University of Texas at Austin Austin , TX, USA 2 Institute of Arctic and Alpine Research University of Colorado at Boulder Boulder , CO, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: chase.rakowski{at}colorado.edu Mathew A. Leibold 3 Department of Biology University of Florida Gainesville , FL, USA 4 Institut Natura e Teoria en Pirineus 09400 , Surba, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site Caroline E. Farrior 1 Department of Integrative Biology University of Texas at Austin Austin , TX, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Abstract Full Text Info/History Metrics Supplementary material Data/Code Preview PDF ABSTRACT Multiple global-change forces, from habitat alterations to warming, are altering food webs and trophic interaction strengths. Such changes in trophic interactions have important implications, as it is a tenet of ecology that trophic interactions are linked to the functioning and stability of ecosystems. For example, changes in the presence or strength of intraguild predation (IGP), the consumption of a predator by another predator that competes for shared prey, can have cascading effects on the biomasses of species and trophic levels. For this reason, IGP can affect key ecosystem functions at the base of the food web and is of special interest to practitioners of biological pest control. However, the relationship between IGP and ecosystem stability is not yet well understood, especially whether and how IGP might affect the stability of non-adjacent lower trophic levels including primary producers. In this study we simulate the dynamics of a six-species, four-trophic-level food web plus a limiting nutrient to explore the relationship between IGP strength and the temporal variability of species- and trophic group-biomass. By varying the IGP rate given the abundance of the eaten predator, we find that the model food web abruptly shifts between equilibria in which all species maintain either constant biomass or stable limit cycles where all trophic levels exhibit sustained and significant oscillations. While complex feedback in the model creates a divergence between the IGP functional response and the resulting realized IGP strength, both stronger IGP functional responses and stronger realized IGP are associated with a higher likelihood of oscillations. Furthermore, analyses indicate that the strongest consumptive interaction induces the oscillating behavior in an indirect effect initiated by the change in IGP. Overall, these results suggest that as food web structure changes in ecosystems worldwide, strengthening IGP runs the risk of inducing destabilizing effects that extend to the base of food webs, while weakening IGP could confer stability to ecosystem functions such as primary production. Finally, we discuss relevance to management, including the implication that IGP among biological control agents should be minimized to maintain stable crop production. INTRODUCTION Global change forces are altering community composition and predation rates in ecosystems around the world, thereby changing food web structure in local communities worldwide. Predators face fast extinction rates yet are sometimes replaced by introduced predator species, while other aspects of global change are modifying relative species abundances, leading to changes in predator and prey composition ( Purvis et al. 2000 , Cucherousset et al. 2012 , Bond and Lavers 2014 , Hu et al. 2022 ). Furthermore, aspects of global change such as warming can directly and indirectly affect predation rates independent of changes in community composition ( Hammill et al. 2018 , Frances and McCauley 2018 ). For these reasons trophic interaction strengths are likely to change in many ecosystems ( Frances and McCauley 2018 ). One such trophic interaction prone to change is intraguild predation (IGP), a very common aspect of food webs defined as the consumption of a predator by another predator which competes for shared prey ( Polis et al. 1989 ). Food web structure is generally an important determinant of ecosystem functioning, and IGP appears to be no exception. Species diversity, the prevalence of weak trophic interactions and omnivory, and the number of trophic levels can all affect ecosystem functions and stability ( McCann and Hastings 1997 , McCann et al. 1998 , Thébault and Loreau 2005 , 2006 , Jiang and Pu 2009 , Tilman et al. 2014 ). IGP in particular has been shown to affect the biomass of various trophic levels. Specifically, empirical and theoretical work has shown that IGP can often weaken total top-down control by multiple predators, and this effect can reverberate through trophic levels to alter trophic cascades ( Finke and Denno 2004 , Straub et al. 2008 ). Additionally, IGP can have varying effects on trophic level biomasses depending on its strength ( Chang et al. 2020 , Chang and Cardinale 2020 ). The nature of the relationship between IGP and ecosystem stability is murkier. Relevant empirical investigations are rare, and theoretical studies offer a mixture of conclusions. Some papers suggest that omnivory – of which IGP is a specific category – is generally destabilizing, both in the sense of constraining coexistence and in the sense of creating oscillations and chaos ( Pimm and Lawton 1978 , Holt and Polis 1997 , Tanabe and Namba 2005 ). However, other studies emphasize that omnivory can either destabilize or stabilize depending on omnivory strength and other conditions ( McCann and Hastings 1997 , Vandermeer 2006 ), and Křivan (2000) found that IGP strength is positively correlated with stability. Importantly, these studies analyzed models with only three or four species, limiting the complex food web feedbacks and prey switching that could underlie any relationship between trophic interaction strength and stability in more diverse food webs ( Post et al. 2000 ). Furthermore, it is not yet known whether IGP strength might influence the variability of biomass in non-adjacent lower trophic levels. Such a relationship would have important implications both for natural ecosystems and for biological control of agroecosystems, because the stability of basal trophic levels provides the foundation of ecosystems, climate regulation, and reliable crop production. Here we use a six-species, four-trophic level food web model (two primary producers, two herbivores, and a smaller predator consumed by a top predator) to test the effects of IGP on the central tendency of biomass and its variability throughout the food web. To explore a range of IGP strengths, we manipulate the maximum IGP rate, a parameter that defines the rate of IGP given density of the smaller predator, rather than adding or removing IGP altogether. The maximum IGP rate captures an aspect of the top predator’s behavior akin to its handling time or preference for the smaller predator relative to other prey species. However, this parameter is not the same as the realized IGP strength, defined here as biomass of the smaller predator consumed per unit time per unit biomass of the top predator, which depends on the maximum IGP rate as well as the densities of all the predator and prey species. We predicted that weaker IGP (of either definition) would allow greater total top-down control of the herbivorous shared prey, leading to reduced herbivore biomass and that the effect would extend to a trophic cascade, with reduced herbivory allowing producers to maintain higher and more stable biomass. Interestingly, we found that the maximum IGP rate has a complex relationship with realized IGP strength due to the feedback among species. Changes in the maximum IGP rate cause state changes and complex rearrangements in biomasses across the food web, thereby altering interaction strengths including IGP strength in a nonlinear manner. However, in general both lower maximum IGP rates and weaker realized IGP are associated with stable biomass of all species and trophic levels. METHODS We simulate temporal dynamics of food webs with a single limiting nutrient 𝑁, two primary producer species A 1 and A 2 , two herbivorous species H 1 and H 2 which consume both A 1 and A 2 , a predator species 𝑃 1 that consumes both H 1 and H 2 , and a top predator 𝑃 2 that consumes both H 1 and H 2 as well as 𝑃 1 , representing IGP ( Fig. 1a ). Within a trophic group, the species with subscript 1 is smaller-bodied than the species with subscript 2, which affects their growth rates and other relevant properties via allometric scaling as explained below. Our model is based on the work by Ceulemans et al. (2019) who study relatively complex food web dynamics in a plankton chemostat framework but exclude IGP. Besides adding IGP and making some relevant parameter changes (see below), we also reduce the model to two rather than four herbivore species, and alter the model to represent natural food webs by allowing mortality rates to vary across species. Download figure Open in new tab Figure 1: a) Diagram representing the model food web structure. Circles represent species; relative circle size represents relative body size within trophic levels. Arrows represent energy flow, with thickness representing qualitative differences in relative per capita consumption rates (determined by half-saturation constants 𝑐 X i Y j and maximum consumption rates 𝑟 𝑌 ) within trophic levels. The red arrow represents intraguild predation (IGP). b-c: Example functional response curves showing IGP rate over density (biomass) of the small predator P 1 , which in this context is the prey. b) 𝑟 P 1 P 2 = 0.04 (low IGP rate at high P 1 densities); c) 𝑟 P 1 P 2 = 0.1 (high IGP rate at high P 1 densities). We simulate biomass changes in each species from one time point to the next as a function of its own abundance, the abundance of its predators, and the abundance of its resources through time. The concentration of available N , n , is replenished at a rate 𝛿 proportional to its current deviation from the loading N concentration n 0 . N is taken up by the producers ( A 1 , A 2 ) following a Monod function: where is the N :carbon weight ratio of producers (parameterized for nitrogen), 𝑟 𝐴 𝑖 is the intrinsic rate of increase of A i , a i is the biomass of A i , and u i is the half-saturation constant of N uptake by A i . We model the consumption of a prey species at any trophic level X by a predator species, Y , with the functional response: where 𝑟 𝑌 is the maximum consumption rate of the predator Y , equivalent to the inverse of handling time, y is the biomass of Y , x is the biomass of X , and where l controls the curvature of the functional response at low prey densities ( l =1 yields a Holling type II while l =2 yields a Holling type III), and 𝑐 𝑋,𝑌 is the half-saturation constant of the consumption of X by Y . Here, l is set to 1.1 so that the functional response is saturating but is less destabilizing than Holling type II ( Fig. 1b-c ; Williams and Martinez 2004 ). When 𝑐 𝑋,𝑌 is high, Y consumes X at a low rate at a given biomass density of X ( x ), and vice versa. Conversely, when 𝑟 𝑌 is high, Y consumes X at a high rate at a given biomass density of X ( x ), and vice versa. Note that consumption of one prey species, 𝑋 1 is influenced by the other prey species, 𝑋 2 indirectly here due to predator satiation captured by the sum over all prey species in eqn. 2 . The change in biomass of producer species A i ( a i ) is governed by: where h j represents the biomass of H j , and the change in biomass of a heterotroph species (𝑌 𝑗 ) that consumes X species and is consumed by Z species is: where 𝑏 X i Y j is the conversion efficiency of 𝑋 𝑖 to 𝑌 𝑗 and 𝑚 𝑌 𝑗 is the baseline mortality rate of 𝑌 𝑗 . The functional responses of predators for given prey items, f , is given by eqn. 2 . For herbivores ( H species), A species take the place of the X ’s and P’s take the place of the Z’s . For 𝑃 1 , X = { H 1 , H 2 } and Z = 𝑃 2 ; and for 𝑃 2 , X = { H 1 , H 2 , 𝑃 1 } and there is no Z . Parameters r and m are different for each 𝑌 𝑗 , and 𝑐 is different for each combination of 𝑋 𝑖 and 𝑌 𝑗 . The conversion efficiency b is kept the same except for the case of IGP, where it is increased (see below). All differential equations that govern the model are written out in Appendix S1 (eqns. S1-S5). We primarily use parameters from Ceulemans et al. (2019) but with several exceptions. While Ceulemans et al. (2019) varied within-trophic level species trait diversity, we only use the parameters representing 75% of the maximal species trait diversity used in their study for pairs of species in a trophic group. The parameters vary based on allometric scaling ( Table 1 ): where, as in Ceulemans et al. (2019) , we assume body mass ratios of 10 3 across adjacent trophic levels and an allometric scaling exponent of -0.15. View this table: View inline View popup Download powerpoint Table 1: Central parameter values used in the model. However, rather than setting all mortality rates equal to each other to simulate a chemostat experiment, we use the same allometric scaling exponent to determine mortality rates across species and better represent natural food webs. Species in lower trophic levels have lower body masses, and one species is smaller than the other within each trophic level. Smaller species have higher maximum growth rates (𝑟 𝐴𝑖 for producers and 𝑏 X i Y j 𝑟 𝑌 𝑗 for consumers), higher rates of nutrient uptake or prey consumption (lower 𝑢 𝐴 𝑖 or lower 𝑐 X i Y j : see eqn. 1 and 3, respectively), and higher baseline mortality rates (𝑚 𝑌 𝑗 ), but they are less defended – that is, they are consumed faster. Within trophic levels, we assume the smaller consumer species consumes the smaller prey species at a higher rate (lower 𝑐 X i Y j ) than it consumes the larger prey species, whereas the larger herbivore is a generalist, which is to say 𝑐 𝐴 1 H 2 = 𝑐 𝐴 2 H 2 . While Ceulemans et al. (2019) modeled four herbivorous species representing all possible combinations of the defended/undefended and selective/generalist traits, we model only an undefended selective ( H 1 ) and a defended generalist herbivore ( H 2 ) to simplify the analysis and focus on the case where the herbivores have a tradeoff with the best chance of enabling coexistence. To allow exploration of a range of IGP strength, we alter some of the predators’ parameters from their allometrically-determined values in a way that enables the predators to coexist in the presence of IGP. First, we assume the predators have a tradeoff such that IGP is especially beneficial for 𝑃 2 (higher conversion efficiency, 𝑏 P 1 P 2 , than for other interactions), but that the smaller predator 𝑃 1 is a superior competitor for their shared herbivore prey (lower 𝑐 H 𝑖𝑃𝑗 value; Holt and Polis 1997 ). Additionally, we lowered 𝑃 2 ’s maximum herbivore predation rate 𝑟 𝑃2 . Finally, the predators’ baseline mortality rates are set equal to each other rather than setting P 1 ’s mortality rate higher in line with its smaller size ( Table 1 ). Note that the smaller predator’s effective mortality rate, including both baseline mortality and mortality due to predation, is still lower than that of the top predator. Also note that, as with other parameters, the predators’ mortality rates were separately varied and therefore no longer equal in the sensitivity analysis described below. We ran simulations across the entire range of values of the maximum IGP rate 𝑟 P 1 P 2 for which all 6 species coexist. The parameter 𝑟 P 1 P 2 is equivalent to the inverse of handling time for the IGP interaction and along with other parameters, especially 𝑐 P 1 P 2 , it determines the degree to which 𝑃 2 exploits 𝑃 1 at a given density of 𝑃 1 . Species extinction was defined by the biomass of a species falling below 0.01 µgC/L after removing transient dynamics. To test the robustness of our results, we re-ran this parameter sweep of 𝑟 P 1 P 2 for different sets of parameter values for nearly all other parameters ( Table 1 ). Each set of parameters was developed by randomly and separately varying parameters by up to 10% above and 10% below the central value, then retaining parameter combinations that adhered to the allometric assumptions. These assumptions include that compared to larger species, smaller species grow faster, take up nutrients or consume prey faster, and experience faster baseline mortality rates (except the smaller vs. top predator), and that within trophic groups, the smaller species consume smaller prey faster than they consume larger prey (see Appendix S1: Table S1 for specific parameter constraints). In all cases, we ran simulations for each parameter combination until equilibrium, either a fixed point or a stable limit cycle, was reached. Then we removed transient dynamics and recorded available nutrient masses, species biomasses, and trophic interaction strengths (nutrient uptake rates per unit biomass of the producer, and consumption rates per unit biomass of the consumer) for at least 2,500 days (6.85 years). Importantly, the realized trophic interaction strengths are different from the parameters defining trophic interactions. For example, while we manipulated the maximum IGP rate 𝑟 P 1 P 2 , we measured the realized IGP strength as the biomass of P 1 consumed per day per unit biomass of P 2 as defined in the model (eqns. 2 and 5 , Appendix S1 eqn. S4): where p 1 , h 1 , and h 2 , are results of the simulation. Under the central parameters and first sets of randomized parameters we thoroughly tested the effects of the initial biomasses of all species on long-run biomass dynamics. Upon finding that all possible alternative stable states were reliably reached across initial predator biomasses of { P 1 = 5, P 2 = 2500}, { P 1 = 500, P 2 = 500}, and { P 1 = 2500, P 2 = 5}, with initial biomass of all other species set to 500, we used these starting conditions in analyzing all subsequent randomized parameter sets. We further conducted analyses to better understand which species induce biomass oscillations. To do so, we first ran simulations with 𝑟 P 1 P 2 values that cause an oscillating state, except that we forced the biomass of one individual species to remain constant at its median. Then we repeated this exercise with every other species, and in both oscillating states (1 and 3), allowing us to identify which species are responsible for the oscillations. These analyses were conducted using the central parameter values ( Table 1 ). All simulations and analyses were conducted in R v. 4.4.0 ( R Core Team 2024 ). RESULTS Our results show that there is a complex, three-phase effect of varying the maximum IGP rate 𝑟 P 1 P 2 – which determines the rate of IGP at given densities of the predators – in which the food web abruptly transitions between sustained oscillations and point equilibria ( Fig. 2 ). Note that under randomized parameter sets, simpler two-phase relationships with an oscillating equilibrium and a point equilibrium were also common (Appendix S1: Fig. S1). The point equilibria are generally associated with low 𝑟 P 1 P 2 and low realized IGP strength ( eqn. 7 ). However, we also find qualitatively distinct oscillations at very low versus high 𝑟 P 1 P 2 , both of which are associated with high realized IGP strength. Therefore, the maximum IGP rate 𝑟 P 1 P 2 and the realized IGP strength are often unaligned. Furthermore, the transitions from stability to both oscillation types often depend on initial predator biomasses. Adding to the complexity of results, we find that for some values of low and high 𝑟 P 1 P 2 , one or more species go extinct under some initial predator biomasses while a stable coexisting state is reached under different initial predator biomasses ( Fig. 3 , Appendix S1: Fig. S1). Download figure Open in new tab Figure 2: Example time series demonstrating long-run behavior of masses of each species and of the available nutrient for each of the three qualitative states. a) State 1, slow oscillations (𝑟 P 1 P 2 = 0.045); b) State 2, static (𝑟 P 1 P 2 = 0.06); c) State 3, fast oscillations (𝑟 P 1 P 2 = 0.08). All simulations were run with the central parameters ( Table 1 ). Download figure Open in new tab Figure 3: Stacked bar chart representing the frequency of qualitative equilibrium states across maximum IGP rates 𝑟 P 1 P 2 , summarized across all randomized parameter sets. Possible states include State 1: slow oscillations, State 2: point equilibrium, and State 3: fast oscillations, as well as alternative stable states where State 1 or 2 are possible or where State 2 or 3 are possible, depending on initial conditions, and finally the outcome where all six species fail to coexist (“no coexistence”). Qualitative states State 1 is characterized by a stable limit cycle of slow biomass oscillations and high producer biomass, specifically a very high R 1 biomass despite R 2 biomass being very low. The available nutrient is quickly taken up by the large R 1 population and does not often reach a high concentration ( Fig. 4a , Appendix S1: Fig. S2a). Herbivore and predator biomass is also somewhat low, although P 1 stays at a higher biomass than in the other states, uniquely resulting in similar biomass of the two predators. State 1 is the equilibrium for a small range of very low maximum IGP rates: 𝑟 P 1 P 2 between 0.043 and 0.045 with the central parameters ( Fig. 4 ), and between 0.036 and 0.051 among all randomized parameter combinations studied ( Fig. 3 , Appendix S1: Fig. S1). Under some parameter combinations State 1 is never possible ( Fig. 3 , Appendix S1: Fig. S1). Under the central parameters, the system only reaches State 1 if both predators begin with moderate biomasses, e.g. 500 each; if one predator begins with a much higher biomass than the other, it reaches State 2. Download figure Open in new tab Figure 4: Long-run mass of a) the available nutrient N and total biomass of b) the producers R i , c) the herbivores H i , and d) the predators P i , for simulations that differ in the maximum IGP rate 𝑟 P 1 P 2 and are run with the central parameters. Many initial conditions were tested, but three initial conditions of P 1 and P 2 biomass were sufficient to reproduce all possible equilibrium states across 𝑟 P 1 P 2 values. All outputs that result in coexistence of all six species are shown. In State 2, the species reach a point equilibrium with no oscillations. The top predator P 2 has a very high biomass, resulting in a high total predator biomass, but P 1 has a very low biomass. The herbivores and producers have moderate biomasses, with R 2 maintaining a somewhat high biomass but R 1 maintaining a very low biomass, keeping the available nutrient at a moderately low concentration ( Fig. 4 a-c, Appendix S1: Fig. S2). This state occurs in a range of moderate to high maximum IGP rates: 𝑟 P 1 P 2 between 0.043 and 0.082 with the central parameters ( Fig. 2b ), or between 0.033 and 0.116 across all parameter combinations studied ( Fig. 3 , Appendix S1: Fig. S1). When 𝑟 P 1 P 2 is between 0.057 and 0.082 with the central parameters, the system can also reach State 3, and it is more likely to do so the higher the initial biomass of P 2 compared to P 1 . Lastly, State 3 is characterized by fast, drastic oscillations across trophic levels. Trophic group and species biomass is moderate on average, but the producers undergo downswings in which producer biomass reaches very low values while available nutrient concentration reaches high values ( Fig. 2c , Fig 4a-b ). This quickly cycling state occurs under higher maximum IGP rates: specifically, 𝑟 P 1 P 2 between 0.057 and 0.097 with the central parameters ( Fig. 2c ) or between 0.045 and 0.267 across all parameter combinations studied ( Fig. 3 , Appendix S1: Fig. S1). Trophic interaction strengths There is a complex relationship between the maximum IGP rate 𝑟 P 1 P 2 and the realized strength of IGP ( eqn. 7 ) and other trophic interactions. The oscillating states are associated with stronger realized IGP while point equilibria are associated with weak realized IGP ( Fig. 5 ). In the static and quickly oscillating states (States 2 and 3), the biomasses of R 1 and P 1 as well as the rates at which other species consume R 1 and P 1 remain low compared to the biomasses and consumption of other species (Appendix S1: Figs. S2 , S3 , S4 ). The distribution of trophic interaction strengths is similar between States 2 and 3, with the major exception that IGP strength is higher in State 3 and very low in State 2 especially as 𝑟 P 1 P 2 increases ( Fig. 5 , Appendix S1: Figs. S3 and S4). In contrast, in State 1 the distribution of trophic interaction strengths is unique. Here, most trophic interactions (rate of consumption of a nutrient or prey by a species per unit biomass of the species) are weaker than in other states, as the biomasses of most species remain low, preventing them from being consumed at a high rate (Appendix S1: Figs. S2 , S3 , S4 ). On the other hand, R 1 maintains a very high biomass, and the consumption rates of R 1 by the herbivores are very high. Crucially, the biomass of P 1 is also much higher than in other states and, accordingly, IGP strength (consumption of P 1 by P 2 ) is higher than in any other state ( Fig. 5 , Appendix S1: Fig. S2). This is true even though the maximum IGP rate 𝑟 P 1 P 2 is at its lowest here, requiring that all else being equal, the consumption of P 1 by P 2 is lowest. Therefore, the high IGP strength in State 1 results from an indirect effect whereby the change in maximum IGP rate redistributes species biomasses, including an increased P 1 : P 2 ratio. Consequently, realized IGP strength (the rate of consumption of P 1 by P 2 per unit biomass of P 2 ) is counterintuitively high under low maximum IGP rates due to the effects of maximum IGP rate on the biomass distribution. Even so, realized IGP strength is generally highest in the cycling states and lowest in the static equilibrium state. Download figure Open in new tab Figure 5: Long-run realized intraguild predation (IGP) strengths resulting from differing maximum IGP rates 𝑟 P 1 P 2 . Parameters match those of Fig. 4 . Realized IGP strength is the consumption rate of the smaller predator by the top predator per unit biomass of the top predator. Causes of cycling Analyses indicate that the strongest consumptive interaction, which differs between State 1 and 3, may induce the oscillations characterizing both of these states. In State 1, herbivory of R 1 by H 1 is the strongest consumptive interaction ( Table 2 ), and only when H 1 is held constant does the food web switch from State 1 to 2. In contrast, in State 3 the strongest consumptive interaction is herbivory of R 2 by H 2 , and oscillations are eliminated only when H 2 is held constant. View this table: View inline View popup Download powerpoint Table 2: Median strength of trophic interactions in each qualitative state. Nutrient uptake strength is measured as nutrient mass taken up per day per biomass of the producer; consumptive interaction strength is measured as biomass consumed per day per biomass of the consumer. Data are shown for the central parameters ( Table 1 ) and 𝑟 P 1 P 2 = 0.045 for State 1, 0.06 for State 2, and 0.08 for State 3. The strongest consumptive interaction is bolded for each state. DISCUSSION We modeled a six-species food web with four trophic levels supported by a limiting nutrient to examine cascading effects of intraguild predation (IGP) strength on the variability of species-level and trophic level biomass. The model food web converges on one of three qualitative states depending on the maximum IGP rate 𝑟 P 1 P 2 (equivalent to the inverse of handling time for the IGP interaction) and other parameters. State 1, characterized by slow oscillations, is the rarest state found only under a small range of low 𝑟 P 1 P 2 values, and is analogous to the “high-production state” in Ceulemans et al. (2019) , so-called due to the relatively high primary production and producer biomass. State 2 is a point equilibrium found under low to moderate 𝑟 P 1 P 2 . Lastly, State 3 is characterized by fast oscillations, is found under moderate to high 𝑟 P 1 P 2 , and is analogous to the “low-production state” in Ceulemans et al. (2019) , so-called due to the relatively low average primary productivity. Biomass oscillations are induced by the strongest consumptive interaction, and are stronger lower in the food web such that the oscillations in producer and available nutrient mass are the most extreme. Due to food web feedback that shifts biomass distributions, IGP strength is counterintuitively highest under very low 𝑟 P 1 P 2 in the slowly oscillating State 1. However, sensitivity analyses reveal that State 1 is rare across parameter combinations ( Fig. 3 , Appendix S1: Fig. S1). Furthermore, in general both higher 𝑟 P 1 P 2 and stronger realized IGP ( eqn. 7 ) are associated with oscillations while lower 𝑟 P 1 P 2 and weaker realized IGP are associated with point equilibria. Therefore, while the relationship between the maximum IGP rate, realized IGP strength, and biomass variability is complex in our model, the results overall support a correlation between stronger IGP in a community and biomass variability for all species. The generally positive association between IGP strength and producer biomass variability fits with our predictions, while the complexity of the association appears to echo the mixed findings of existing theory. Most of the relevant theoretical work suggests IGP is either de-stabilizing or can have various effects on stability depending on assumptions and conditions. However, the effects of the maximum IGP rate parameter in our model are not only complex but also unexpected. Still, the complex relationship between maximum IGP rate and outcomes including biomass and realized IGP strength make sense when accounting for how the food web structure shifts. For example, a very low maximum IGP rate keeps the top predator suppressed to a low biomass, presumably due to food limitation, while allowing the small predator species to reach a very high biomass, presumably due to predation release. So even though the very low maximum IGP rate parameter requires the top predators to consume the small predators at a relatively low rate at a given density of small predators, there are so many small predators that the individual top predators still consume a relatively large amount of small predators and relatively few herbivores, leading to strong realized IGP. The absolute rate of IGP in terms of total biomass consumed per unit time (energy transfer), however, remains moderate as there are few top predators (Appendix S1: Fig. S5). Therefore, realized IGP as we define it (small predator biomass consumed per unit time per unit biomass of top predators) largely reflects the abundance of the small predator relative to the herbivore species, while the maximum IGP rate is a property of the top predator’s behavior and the absolute IGP rate integrates relative and absolute abundances. The relationship between realized IGP strength and the central tendency of trophic level biomasses was also unexpected. Previous work suggested that IGP weakens trophic cascades ( Finke and Denno 2004 , 2005 , Straub et al. 2008 ), such that stronger IGP would lead to lower producer biomass. Instead, we found the opposite ( Figs. 4b , 5 ). This result seems to be related to the balance of the energy pathways, e.g. R 1 → H 1 → P 1 vs. R 2 → H 2 → P 2 . In States 2 and 3, P 1 is suppressed by P 2 to low biomass. However, in State 1 when the maximum IGP rate is especially low, P 2 does not significantly suppress P 1 , so P 1 maintains higher biomass than in other states. This allows P 1 to suppress H 1 to a lower level than in other states, freeing R 1 to reach very high biomass. This is the only state in which the two predators maintain similar biomass, allowing them to more effectively partition prey (herbivore) resources and cause a stronger trophic cascade. Interestingly, the predators’ combined biomass is lowest in this state, as would be expected when (realized) IGP is high, but contrary to general expectations for when a trophic cascade is strongest. Predictions are usually that high IGP will reduce the predators’ total control of their shared prey ( Finke and Denno 2005 , Straub et al. 2008 ). However, in this case, even though the combined biomass of the predators is lower in State 1, their much more balanced biomass appears to allow them to better exploit their shared prey and cause a stronger trophic cascade by affecting both energy pathways and preventing intra-trophic-level biomass compensation. This particular finding highlights the importance of predator species evenness in shaping food web structure in addition to IGP strength. Theoretical studies have found negative ( Pimm and Lawton 1978 , Holt and Polis 1997 , Tanabe and Namba 2005 ), positive ( Křivan 2000 ), and context-dependent ( McCann and Hastings 1997 , Vandermeer 2006 ) relationships between omnivory or IGP and food web stability. There are many differences in assumptions among these models which likely contribute to the various conclusions. For example, unlike the other studies, Křivan (2000) assumed IGP is adaptive, meaning the top predator adjusts its IGP strength to maximize its population growth rate in a form of optimal foraging, which may explain why IGP provides stability in that model. In our model, unlike in most other models of IGP dynamics, we include an alternate herbivore species for a total of three prey species consumed by the top predator. Holt and Huxel (2007) found that adding an alternate prey species to a 3-species IGP food web facilitates coexistence of the predators but also induces unstable dynamics over a wide parameter space, which could help explain why IGP was often associated with oscillations in our model. Besides an additional prey species, we also included additional trophic levels and more complex functional responses for each trophic interaction in the model. Simply put, our model is more complex than the models in the aforementioned studies and could potentially reflect real food webs more accurately, albeit at a sacrifice of generality. Compared to most terrestrial aboveground food webs, our model food web exhibits faster dynamics due to higher growth, consumption, and mortality rates, as well as stronger allometric scaling structuring these rates based on trophic level and organism size. Perhaps in food webs with different characteristics, a different relationship between IGP strength and stability may be found. The time is ripe for more empirical studies guided by IGP theory, which could bridge the gap between the various models and real food web dynamics and provide needed context for further theoretical developments. Empirical tests of the relationship between IGP strength and biomass variability are rare, but a few relevant experiments add to the IGP-stability debate. For example, Fagan (1997) and Holyoak and Sachdev (1998), studying terrestrial arthropod and aquatic protist-bacteria food webs, respectively, concluded that the presence of omnivory was stabilizing. In contrast, other experiments using aquatic protist-bacteria food webs found that IGP increased the variability of bacterial biomass ( Lawler and Morin 1993 ) and predator biomass ( Morin 1999 ). However, all of these experiments only manipulated the presence or density of omnivorous predators, and therefore were not able to separate the effects of IGP from the effects of species identity or competition between the predators (i.e., multiple food web links were added simultaneously, rather than an IGP link only). Chang and Cardinale (2020) overcame this challenge by using mesh to manipulate the degree of IGP without affecting the strength of other species interactions in an aquatic microbial food web, although they did not study biomass stability. Their paired experiment and model found that adding weak IGP to a food web decreases the biomass of the shared prey, but stronger IGP increases it. However, in our model we found a different relationship where realized IGP strength is negatively correlated with the biomass of shared prey (herbivores) within qualitative states. Again, this discrepancy could be due to the larger number of species and nonlinearities in our model which create cascading shifts in biomass throughout the food web as the IGP rate changes. It is also worth noting that the large change in realized IGP strength between the static and oscillating states in our model scarcely affects herbivore biomass, and even the correlation between IGP and herbivore biomass within states is weak. Chang et al. (2020) similarly found that when the predators in their model significantly partitioned their prey resources, as occurs in our model, the effects of IGP strength on predator biomass became weak. Indeed, the more important effects of IGP strength are the effects on biomass variability in our model, hence our focus on variability throughout the paper. Note, we were not able to explore all solutions of our model due to the inability to solve the model analytically. However, the parameters were determined carefully to fit real data, past relevant theoretical work, and key tradeoffs found in real food webs. It is also worth noting that as most parameters are varied away from their central values, the six species coexist within a more restricted range of 𝑟 P 1 P 2 . Perhaps, therefore, the central parameter values are reasonably representative of real food webs which exhibit a wide range of IGP rates. Overall, despite these caveats we are confident in the conclusion that IGP can either generate cycling behavior or maintain stability in model food webs, depending on its strength. We look forward to empirical tests of our findings and for the intuition that this model may bring to understanding the dynamics of aquatic, terrestrial, natural, and managed ecosystems. The results of this study can be applied to biological pest control methods. When multiple predators are used as biological control agents in a single system, releasing or attracting predators which only engage in low rates of IGP may have the effect of enhancing the stability of crop yields. This is particularly likely if the focal species adhere to our model assumptions, such as exhibiting relatively high birth and consumption rates and following allometric scaling as in belowground and planktonic systems. While our results could be applicable to some terrestrial agricultural systems, planktonic systems are also relevant to biological control. Environmental algal technologies, such as cultivating phytoplankton as a way of cleaning wastewater while producing biofuel feedstock, are promising strategies for combating multiple environmental challenges but suffer from the proliferation of zooplankton ( Smith et al. 2010 ). Therefore biological control is promising for improving the economic feasibility of algal cultivation ( Rakowski et al. 2024 ), and this nascent field of algal cultivation biocontrol could benefit from insight into the properties of predators that might improve the stability of algal production. As global change further alters food webs and predation rates, the results of this study imply that increases in IGP rates could run the risk of de-stabilizing populations and key ecosystem functions such as primary production. For example, one study found that warming significantly increased IGP levels among aquatic invertebrates ( Frances and McCauley 2018 ). However, the generality of this relationship across ecosystems remains to be determined, and there may not be reason to expect a universal de-stabilization of ecosystems mediated by changing IGP. Still, the connections explored herein shed light on the intricacies of strong indirect effects in even a relatively simple six-species food web and underline the possibility of strong cascading impacts of changing food web composition and predation rates. Fruitful avenues for future research could include empirically testing whether IGP strength is positively correlated with biomass variability across trophic levels, and improving predictions for when IGP might increase or decrease. AUTHOR CONTRIBUTIONS CJR conceived the idea and built and analyzed the model with guidance from MAL and CEF. CJR drafted the manuscript, and all authors contributed revisions and approved the final version for publication. CONFLICT OF INTEREST STATEMENT The authors report no conflicts of interest. ACKNOWLEDGMENTS We would like to thank Ursula Gaedke and Matthieu Barbier for providing guidance for the model analysis. Christopher Peterson aided in writing the compact form of the model equations. We also thank members of the Farrior and Leibold groups as well as Shalene Jha for providing helpful feedback on the model analysis and earlier versions of the manuscript. The University of Texas at Austin provided financial support. Footnotes OPEN RESEARCH STATEMENT: All code is publicly available on Zenodo at DOI: . https://doi.org/10.5281/zenodo.15733522 REFERENCES 1. ↵ Bond , A. L. , and J. L. Lavers . 2014 . Climate change alters the trophic niche of a declining apex marine predator . Global Change Biology 20 : 2100 – 2107 . OpenUrl CrossRef PubMed Web of Science 2. ↵ Ceulemans , R. , U. Gaedke , T. Klauschies , and C. Guill . 2019 . The effects of functional diversity on biomass production, variability, and resilience of ecosystem functions in a tritrophic system . Scientific Reports 9 : 7541 . OpenUrl CrossRef PubMed 3. ↵ Chang , F. , and B. J. Cardinale . 2020 . Intra-guild predation (IGP) can increase or decrease prey density depending on the strength of IGP . Ecology 101 : e03012 . OpenUrl CrossRef PubMed 4. ↵ Chang , F.-H. , P.-J. Ke , and B. Cardinale . 2020 . Weak intra-guild predation facilitates consumer coexistence but does not guarantee higher consumer density . Ecological Modelling 424 : 109019 . OpenUrl CrossRef 5. ↵ Cucherousset , J. , S. Blanchet , and J. D. Olden . 2012 . Non-native species promote trophic dispersion of food webs . Frontiers in Ecology and the Environment 10 : 406 – 408 . OpenUrl CrossRef 6. ↵ Fagan , W. F . 1997 . Omnivory as a stabilizing feature of natural communities . The American Naturalist 150 : 554 – 567 . OpenUrl CrossRef PubMed Web of Science 7. ↵ Finke , D. L. , and R. F. Denno . 2004 . Predator diversity dampens trophic cascades . Nature 429 : 407 – 410 . OpenUrl CrossRef PubMed Web of Science 8. ↵ Finke , D. L. , and R. F. Denno . 2005 . Predator diversity and the functioning of ecosystems: the role of intraguild predation in dampening trophic cascades . Ecology Letters 8 : 1299 – 1306 . OpenUrl CrossRef Web of Science 9. ↵ Frances , D. N. , and S. J. McCauley . 2018 . Warming drives higher rates of prey consumption and increases rates of intraguild predation . Oecologia 187 : 585 – 596 . OpenUrl CrossRef PubMed 10. ↵ Hammill , E. , E. Johnson , T. B. Atwood , J. Harianto , C. Hinchliffe , P. Calosi , and M. Byrne . 2018 . Ocean acidification alters zooplankton communities and increases top-down pressure of a cubozoan predator . Global Change Biology 24 : e128 – e138 . OpenUrl PubMed 11. ↵ Holt , R. D. , and G. R. Huxel . 2007 . Alternative prey and the dynamics of intraguild predation: theoretical perspectives . Ecology 88 : 2706 – 2712 . OpenUrl CrossRef PubMed Web of Science 12. ↵ Holt , R. D. , and G. A. Polis . 1997 . A theoretical framework for intraguild predation . The American Naturalist 149 : 745 – 764 . OpenUrl CrossRef Web of Science 13. Holyoak , M. , and S. Sachdev . 1998 . Omnivory and the stability of simple food webs . Oecologia 117 : 413 – 419 . OpenUrl CrossRef PubMed 14. ↵ Hu , Z. , J. Yao , X. Chen , X. Gong , Y. Zhang , X. Zhou , H. Guo , and M. Liu . 2022 . Precipitation changes, warming, and N input differentially affect microbial predators in an alpine meadow: Evidence from soil phagotrophic protists . Soil Biology and Biochemistry 165 : 108521 . OpenUrl CrossRef 15. ↵ Jiang , L. , and Z. Pu . 2009 . Different effects of species diversity on temporal stability in single-trophic and multitrophic communities . The American Naturalist 174 : 651 – 659 . OpenUrl CrossRef PubMed Web of Science 16. ↵ Křivan , V . 2000 . Optimal intraguild foraging and population stability . Theoretical Population Biology 58 : 79 – 94 . OpenUrl CrossRef PubMed Web of Science 17. ↵ Lawler , S. P. , and P. J. Morin . 1993 . Food web architecture and population dynamics in laboratory microcosms of protists . The American Naturalist 141 : 675 – 686 . OpenUrl CrossRef PubMed Web of Science 18. ↵ McCann , K. , and A. Hastings . 1997 . Re-evaluating the omnivory-stability relationship in food webs . Proceedings of the Royal Society B: Biological Sciences 264 : 1249 – 1254 . OpenUrl CrossRef Web of Science 19. ↵ McCann , K. , A. Hastings , and G. R. Huxel . 1998 . Weak trophic interactions and the balance of nature . Nature 395 : 794 – 798 . OpenUrl CrossRef Web of Science 20. ↵ Morin , P . 1999 . Productivity, intraguild predation, and population dynamics in experimental food webs . Ecology 80 : 752 – 760 . OpenUrl CrossRef Web of Science 21. ↵ Pimm , S. L. , and J. H. Lawton . 1978 . On feeding on more than one trophic level . Nature 275 : 542 – 544 . OpenUrl CrossRef Web of Science 22. ↵ Polis , G. A. , C. A. Myers , and R. D. Holt . 1989 . The ecology and evolution of intraguild predation: potential competitors that eat each other . Annual Review of Ecology and Systematics 20 : 297 – 330 . OpenUrl CrossRef Web of Science 23. ↵ Post , D. M. , M. E. Conners , and D. S. Goldberg . 2000 . Prey preference by a top predator and the stability of linked food chains . Ecology 81 : 8 – 14 . OpenUrl CrossRef Web of Science 24. ↵ Purvis , A. , J. L. Gittleman , G. Cowlishaw , and G. M. Mace . 2000 . Predicting extinction risk in declining species . Proceedings of the Royal Society of London. Series B: Biological Sciences 267 : 1947 – 1952 . OpenUrl CrossRef PubMed Web of Science 25. ↵ R Core Team . 2024 . R: A language and environment for statistical computing . R Foundation for Statistical Computing , Vienna, Austria . 26. ↵ Rakowski , C. J. , M. A. Leibold , and S. R. Manning . 2024 . Sunfish as zooplankton control agents to improve yields of wastewater-cultivated algae . Journal of Applied Phycology 36 : 2467 – 2476 . OpenUrl CrossRef 27. ↵ Smith , V. H. , B. S. M. Sturm , F. J. deNoyelles , and S. A. Billings . 2010 . The ecology of algal biodiesel production . Trends in Ecology & Evolution 25 : 301 – 309 . OpenUrl CrossRef PubMed 28. ↵ Straub , C. S. , D. L. Finke , and W. E. Snyder . 2008 . Are the conservation of natural enemy biodiversity and biological control compatible goals? Biological Control 45 : 225 – 237 . OpenUrl CrossRef 29. ↵ Tanabe , K. , and T. Namba . 2005 . Omnivory creates chaos in simple food web models . Ecology 86 : 3411 – 3414 . OpenUrl CrossRef Web of Science 30. ↵ Thébault , E. , and M. Loreau . 2005 . Trophic interactions and the relationship between species diversity and ecosystem stability . The American Naturalist 166 : E95 – E114 . OpenUrl CrossRef PubMed Web of Science 31. ↵ Thébault , E. , and M. Loreau . 2006 . The relationship between biodiversity and ecosystem functioning in food webs . Ecological Research 21 : 17 – 25 . OpenUrl CrossRef 32. ↵ Tilman , D. , F. Isbell , and J. M. Cowles . 2014 . Biodiversity and Ecosystem Functioning . Annual Review of Ecology, Evolution, and Systematics 45 : 471 – 493 . OpenUrl CrossRef 33. ↵ Vandermeer , J . 2006 . Omnivory and the stability of food webs . Journal of Theoretical Biology 238 : 497 – 504 . OpenUrl CrossRef PubMed Web of Science 34. ↵ Williams , R. J. , and N. D. Martinez . 2004 . Stabilization of chaotic and non-permanent food-web dynamics . The European Physical Journal B - Condensed Matter 38 : 297 – 303 . OpenUrl View the discussion thread. Back to top Previous Next Posted June 28, 2025. Download PDF Supplementary Material Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Strengthening intraguild predation increases the temporal variability of biomass across all trophic levels in model food webs Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Strengthening intraguild predation increases the temporal variability of biomass across all trophic levels in model food webs Chase J. Rakowski , Mathew A. Leibold , Caroline E. Farrior bioRxiv 2025.06.25.661600; doi: https://doi.org/10.1101/2025.06.25.661600 Share This Article: Copy Citation Tools Strengthening intraguild predation increases the temporal variability of biomass across all trophic levels in model food webs Chase J. Rakowski , Mathew A. Leibold , Caroline E. Farrior bioRxiv 2025.06.25.661600; doi: https://doi.org/10.1101/2025.06.25.661600 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Ecology Subject Areas All Articles Animal Behavior and Cognition (7635) Biochemistry (17697) Bioengineering (13894) Bioinformatics (41951) Biophysics (21455) Cancer Biology (18592) Cell Biology (25507) Clinical Trials (138) Developmental Biology (13380) Ecology (19903) Epidemiology (2067) Evolutionary Biology (24321) Genetics (15610) Genomics (22509) Immunology (17737) Microbiology (40398) Molecular Biology (17182) Neuroscience (88618) Paleontology (667) Pathology (2833) Pharmacology and Toxicology (4825) Physiology (7641) Plant Biology (15158) Scientific Communication and Education (2046) Synthetic Biology (4296) Systems Biology (9825) Zoology (2271)