Testing for top‐down cascading effects in a biomass‐driven ecological network of soil invertebrates

Abstract To investigate the structural changes of a food‐web architecture, we considered real data coming from a soil food web in one abandoned pasture with former low‐pressure agriculture management and we reproduced the corresponding ecological network within a multi‐agent fully programmable modeling environment in order to simulate dynamically the cascading effects due to the removal of entire functional guilds. We performed several simulations differing from each other for the functional implications. At the first trophic level, we simulated a removal of the prey, that is, herbivores and microbivores, while at the second trophic level, we simulated a removal of the predators, that is, omnivores and carnivores. The five main guilds were removed either separately or in combination. The alteration in the food‐web architecture induced by the removal of entire functional guilds was the highest when the entire second trophic level was removed, while the removal of all microbivores caused an alteration in the food‐web structure of less than 5% of the total changes due to the removal of opportunistic and predatory species. Omnivores alone account for the highest shifts in time of the numerical abundances of the remaining species, providing computational evidence of the importance of the degree of omnivory in the stabilization of soil biota.


| 7063
CONTI eT al. correlation between aboveground farming intensity and the belowground functional diversity is known (Mulder, De Zwart, Van Wijnen, Schouten, & Breure, 2003;Mulder, Cohen, Setälä, Bloem & Breure, 2005) and makes some assumptions on any constant biomass distribution restrictive. Hunt et al. (1987) assumed that biomass inputs exactly balance biomass outputs at all times in each soil compartment of the detrital food web. But this steady-state assumption not always holds, as the biomass decreases of microorganisms, microfauna and mesofauna indicate negative effect on the soil buffer capacity (Hunt & Wall, 2002;Wall, Nielsen & Six, 2015) due to an increased pressure, or even soil exploitation, caused by intensive management practice Mulder, Den Hollander, & Hendriks, 2008).
There is a consensus regarding the evidence that any loss of keystone species severely disrupts ecosystem functioning (IPBES, 2019), but the functional impacts on the food-web architecture of such decimations are almost unknown. Functional guilds are keystone units and can be defined as trait-driven groups of species with key roles in community architecture and therefore ecosystem functioning (Mouquet, Gravel, Massol, & Calcagno, 2013;Power et al., 1996). Recently, Brose et al. (2019) pointed out the importance of preypredator relationships, demonstrating that metabolic rate and functional traits of predatory species are more important in determining the interactions between the weight of the prey and the weight of the predator. Despite this recent meta-analysis, it seems often difficult to extrapolate such outcomes to a much wider context, like behavioral ecology, ecosystem services, or ecosystem functioning (QUINTESSENCE, 2016). Hence, such an extrapolation has to rely upon undiluted mathematical evidences (Bourne, Brenner, & Eisen, 2005;Cohen, 2004).
Mathematically, food webs can be seen as complex ecological networks (Pascual & Dunne, 2005), that is, as graphs whose nodes represent the species present in the ecosystem and whose edges, or links, symbolize the prey-predator relationships between the various species. In particular, a food web is represented by a directed graph, because the connections are expressed by oriented links, which connect the prey/resource to the predator/consumer. In other words, the direction of a given link follows the flow of energy that, through predation, passes from the prey to its predator. Each node/ species will be therefore characterized by a degree that defines the total number of (ingoing and outgoing) links connecting it to the other nodes. These links reflect the architecture of any food web and characterize the energy flux across trophic levels and therefore ecosystem functioning as a whole.
Big data at all biological scales became a central feature of research and discovery in the life sciences (Bourne et al., 2005;QUINTESSENCE, 2016). Our aim is to propose here a simple method to quantify how disproportionate the impact of less functional diversity can be and to illustrate its application with one real soil food web.
In particular, we rebuilt the considered ecological network (Mulder & Elser, 2009) within a fully programmable multi-agent environment (Wilensky, 1999) in an attempt to figure out the cascading effects by weakening the functional diversity according to the of soil invertebrates through a simulated removal of entire functional guilds.

| Data sampling
We used a reference data set of empirically observed soil invertebrates in the sandy soils of one Dutch pasture with former low-pressure management (Mulder & Elser, 2009). Three replicate samples of about 5 m 2 from the upper 10 cm of soil for the fauna were taken.
Bulk samples of 50 soil cores (diameter 2.3 cm) were used to extract the microfauna, and two soil cores (diameter 5.8 cm) were used to extract the mesofauna.
Extraction of free-living nematodes was performed within one week of core sampling using Oostenbrink funnels, and all the elutriated nematodes were collected; ecto-and endoparasitic nematodes were recovered with centrifugal flotation. All nematode individuals were counted, and ~150 randomly chosen specimens were identified and measured under a light microscope (Mulder & Vonk, 2011). Enchytraeid worms (Oligochaeta: Enchytraeidae) were sampled by wet extraction and microarthropods (Acarina and Collembola) by dry extraction (Cohen & Mulder, 2014). In both sampling protocols, the heat was increased gradually with incandescent bulbs, and the invertebrates escaped by moving downward. For enchytraeids and microarthropods, the abundances for 1 m 2 × 10 cm depth were derived from the surface and the bulk density of the soil samples .
All these organisms live in a dark and intricate world, interacting in a detrital food web and in close contact with the soil. Each organism has its own trait-driven function in soil biota, giving the soil its exclusive properties, and an interaction matrix was created based on an inventory of multitrophic interactions of soil food webs that provide all links consistent with literature-derived guilds (Mulder & Elser, 2009 : their table   S2). Despite the observation that energetic equivalence rule is rarely supported within local communities (Morlon et al., 2009), our reference soil food web was particularly stable according to the Eltonian rule (Elton, 1927). The lumped dry weight of all the sampled invertebrates of the first trophic level was exactly 10.28 times the lumped dry weight of all the sampled invertebrates of the second trophic level, as expected according to the energetic equivalence rule (Mulder & Elser, 2009).

| Analysis of the food web: a modified Lotka-Volterra model
In contrast to previous studies on simulated species extinction (Ives & Cardinale, 2004), we want to disentangle the cascading effects of the removal of selected functional groups (here, entire guilds of soil invertebrates identified at genus or family). Using the data discussed in the previous paragraph, a real food web was built (Table S1) within NetLogo (Wilensky, 1999). For these 62 taxa (hereafter, just species), we know the abundance X i , the body mass M i , the biomass B i (given by B i = X i M i ), and the value of the growth rate r i in condition of interaction, with i = 1, 2, ..., n. All identified soil invertebrate "species" fell into five main guilds, and the independent trophic links among guilds (from any possible prey to its consumer) were inferred from published literature. The complete inventory of multitrophic interactions consistent with literature-derived guilds is fully downloadable from Mulder and Elser (2009).
In Figure S1, we show this ecological network, with the numbered nodes/species placed in a circular layout, where each group of species is distinguished by a different color-as explained in detail in Table S1.
The size of each node in online Figure S1 is proportional to the base-10 logarithm of the abundance of the corresponding species (Log X i ).
As multitrophic interactions between basal consumers and allochthonous resources are donor-controlled, that is, according to Polis, Anderson, and Holt (1997) "consumers benefit from but do not affect resource renewal rate," we postulated constant allochthonous resources. Our lemma is therefore that the faunal populations have an unlimited resource supply of bacteria, fungi, and roots. This is because we are focusing only on the prey-predator interactions among invertebrates and not on all the resource-consumer interactions occurring in soil biota. Hence, we keep in our simulations microbial and plant biomasses constantly available for grazing by basal (specialized) species and nonbasal (omnivore) species (cf. Hunt et al., 1987;Polis et al., 1997). Then, our assumptions will be: • In the absence of predators, the population of the prey would grow proportionally to its size.
• In the absence of prey, the population of the predator would decline proportionally to its size, meaning extinction of that population.
• When both predator and prey are present, the interspecific effect of the predation has to be represented as a decrease in the population of the prey and an increase in the population of the predator (Supporting information).
When in an ecosystem there are more than one prey and one predator species, assuming a fixed amount of foraging effort, the Lotka-Volterra equations can be extended as in Hodzic, Selman, and Hadzikadic (2016): where i, j = 1, 2, ..., n, being n the total number of species in the ecosystem. Coefficient A i is the intrinsic growth rate of the i-th species, while coefficients A ij in the summation take into account the effect of the predation. We would like that the term within parenthesis can be interpreted as a sort of generalized growth rate for species X i , which considers the interaction of X i with the other species. We can therefore build a community matrix A, relative to the ecosystem, whose elements A ij weigh the effect of predation between pairs of species: Hence, A ij ≠ 0 when species i and j are linked by predation phenomena (in particular, A ij > 0 when species i preys species j, and A ij < 0 when species j preys species i), while A ij = 0 when species i and j do not have any predation connection. Nonzero diagonal elements A ii represent the phenomenon of cannibalism when an individual of species i preys on another individual of the same species.
Our starting point is the logistic equation (Supporting information), which takes into account the intraspecific regulation due to the presence of limited resources, applied here to the abundance X i of each of the 62 species as in Kondoh (2005): where r i and K i are, respectively, the growth rate and the carrying capacity of the i-th species. In order to adapt these equations to real data, we considered that, in addition to the intraspecific interactions, also interspecific interactions occur among species. For this reason, the growth rate r i must also consider the effect due to predation.
According to the prescription of Equation (1), we assume that, in the presence of interaction, r i is expressed as follows: where A ij are the elements of the community matrix A, and r i0 is the growth rate of the i-th species in the absence of interaction. By combining the Equations (2) and (3), one finally gets: As we have already seen, in general the elements A ij of the community matrix weigh the food interaction among pairs of species (A ij = 0 when species i and j do not have any predation connection). In this study, we decided to link these weights to the biomasses of the prey species. Such a biomass-driven perspective, in fact, focuses on groups of species and on the lumped biomass values of the populations that constitute that functional group (Moore & De Ruiter, 2012).
Let us consider, for example, a generic node/species i with degree 2 which is, simultaneously, a predator for a given node/species m (therefore, an ingoing link will exist from node m toward node i) and a prey for another node/species l (in this case, a directed outgoing link will exist from node i toward node l). In this case, we will define the coefficient A im as the following ratio: where at the numerator there is the biomass B m of the prey species m and at the denominator the summation of the biomasses of all the prey species of node i (included B m ): The coefficient is positive because the flux of energy goes from m to i; therefore, after an encounter, species i will have an increase in abundance. Similarly, we defined the coefficient A il as follows: where at the numerator there is the biomass B i of the prey species i and at the denominator the summation of the biomasses of all the prey species of node l (included B i ). In this case, the coefficient is negative because the flux of energy goes from i to l; consequently, after an encounter, species i will have a decrease in abundance. The rationale behind these definitions-which, of course, can be applied to node/ species with any degree-is that, when a predator has a diet based on a few prey species, he will consume a greater quantity of each of them depending on their single biomass in relation to the total biomass of its prey species. If, on the contrary, the same predator preys many species, it will consume a smaller quantity of each of them in relation to their biomass compared to the total.
In order to proceed with the calculation of the abundances through Equation (4), we should know the term r i0 which can be inferred from real data by making the following assumption. At the time of sampling, the system was-according to Eltonian rules-in a state of stability in which the abundance of each species, in the presence of interactions with all the other ones, had reached its carrying capacity; thus, X i = K i . Therefore, according to the Equation (3), it is possible to obtain the value of the net growth rate without interaction, r i0 , starting from that one measured when there is interaction, r i , considering all the species in their stationary state K i : The parameter α is a coupling coefficient that can be considered as a measure of the interaction strength of a given species within the rest of the food web, and it is chosen so that the carnivores and part of the omnivores (species 57, 58, 59, 60 and 61 of the Table S1) have a negative r i0 . In fact, in the absence of interaction between species (thus without possibility of predation), carnivores must have a negative growth rate. The same happens for omnivores whose diet is composed of animals rather than plants.
Considering our lemma (we kept bacteria, fungi, and roots constantly available to soil invertebrates) and under the further plausible assumption that, in the absence of interaction among species, the carrying capacity of each species, say K i0 , would be greater than the same one in the presence of interaction, that is, K i , we postulate that: Notice that, despite this prescription, in the absence of interaction, species with r i0 < 0 will tend to extinction.
Summarizing, we can effectively rewrite our dynamical equations as follows: with β = 0 in the case in which there is no interaction and β = 1 in the case in which there is interaction, of strength α, between the species.
Equation (9), applied to each node of our ecological network, allows us to simulate the dynamical evolution of the system in several representative scenarios, where different kinds of perturbations will be realized in order to study the reaction of the species. Notice that, despite Equation (8), in the absence of interaction(s), species with r i0 < 0 will tend to extinction as expected. All the simulations were done by choosing the initial abundance of the species in the interval: so that they cannot exceed their carrying capacity. For each scenario, starting from the initial conditions (10), at each time step the populations X i (t) of all species are updated by numerically integrating equation (9) until the system has reached a condition of stability. Notice that a variation of X i (t) for a given species implies a variation of its biomass In turn, this can induce cascading effects on some of the other species, according to the following rules: When the abundances of carnivores or omnivores start to decrease in time, it is reasonable to assume that a corresponding increase in the carrying capacity K i of their prey will occur. The latter can be obtained starting from the potential variation of the biomass of the i-th prey at time t, that is calculated as the product of its own current biomass B i (t) times the ratio between the total biomass of dead predators and the total biomass of all their prey, that is: This quantity can be translated into a consequent potential increase in the prey's abundance, allowed by the decrease of its predators, and this increase can therefore be added to its carrying capacity, so that: for herbivores, fungivores, and bacterivores (first trophic level), and tot_biomass_dead_predators tot_biomass_prey for carnivores and omnivores (second trophic level). The multiplicative factor 10 in Equation (13) was inserted according to the Eltonian rule across adjacent trophic levels (Elton, 1927).
The new value for K i t + 1 will be inserted in the Equation (9), thus influencing the further dynamical evolution of the system. The same rule does not need to be applied if the decrease concerns the abundances of herbivores, fungivores, and bacterivores since these grazing species are according to the lemma only prey; therefore, they cannot induce variations in the carrying capacity of other species.
In order to quantify structural changes and to compare one single simulation to others, that is, the results of the simulations carried out by removing either guild, we introduce an Alteration Index (AI), defined for each ecosystem as follows: where X sk and X fk are, respectively, the abundance of species k-th calculated after 100 time steps, that is, in the steady state, and the abundance of the same species calculated at the end of the simulation. In other words, AI considers the sum of the absolute variations in abundance that the species undergo due to the forced removal of some other species, normalized with respect to their abundance in the steady state. Our AI is of particular interest as the stability of depaupered food webs remaining after deleting functional groups has not been examined systematically by Hunt and Wall (2002). These authors state that such a new food-web stability may be critically important to ecosystem function, and therefore, we see AI as a measure of the alteration of the ecosystem due to the introduced perturbation. Note that here the summation includes only the guilds for which forced removal does not occur, since we are interested in quantifying only the direct effects of the perturbation (without including the perturbation itself).

| COMPUTATIONAL RE SULTS
Let us now discuss in detail the simulations performed in nine different scenarios: removing separately the five main guilds or removing  (2002)).
In the simulation shown in Figure 1a, we set β = 0 in order to test the behavior of the system in the (unrealistic) scenario in which there were no interactions among species. As expected, species with r i0 < 0 (i.e., the carnivorous and omnivorous species F I G U R E 2 Graphic representation of the simulation results in the case all herbivores (a), bacterivores (b), or fungivores (c) are removed from the ecosystem. Left, the temporal evolutions of species abundances and right, the representations of the food web at the end of these three simulations [for the upper two panels, the first 500 time steps are shown; in the lower panel, fluctuations occur after 700 time steps] 57, 58, 59, 60, and 61) become extinct. The other species, instead, rapidly reach their own carrying capacity. On the other hand, in the simulation shown in Figure 1b, we set β = 1, that is, we consider the standard scenario in which there is interaction between species, but the system is not disturbed. In this case, all species reach their carrying capacity and the system goes rapidly in balance. This is the condition with which we have to compare all the scenarios, in which we always set β = 1, but the soil system is perturbed by the forced loss of entire guilds.
In summary, we will perform nine simulations differing one from each other concerning the functional guilds which were removed from the food web. We removed either one or more guilds (even an entire trophic level in the case of the removal of all omnivores and carnivores) to forecast an evolution of such an artificially depaupered food web.
• Removal of either all the herbivorous species together or all the bacterivorous species together. The soil system does not seem to be affected by any of these disturbances, and all other species settle at their carrying capacity (Figure 2a and b), and the alterations of the food-web architectures, as computed by Equation (15) terivorous enchytraeids (Mulder et al., 2008), and therefore, we would have expected stronger cascading effects.
• Removal of all the fungivorous species. In this case, being the fungivorous Tylenchidae (node 15) the most abundant species of the entire system, its removal causes a slight decline of the second most abundant species of the system, the herbivorous Helicotylenchus (node 2), which decreases from 10 5.9 to 10 5.82 . A possible expla-  Table S1) from 10 4.6 to 10 4.56 is observed (Figure 2c), as resources do not seem to be sufficient. Anatonchus is well known to be a stress-resistant nematode, capable to survive and dominate in hostile environments (Fiscus & Neher, 2002;Neher, Wu, Barbercheck, & Anas, 2005).
• Removal of all the herbivorous and fungivorous species together. The abundance of the carnivorous nematode Anatonchus (node 43 of Table S1) decreases from 10 4.6 to 10 4.48 and, as a consequence, a slight increase is also observed in some species of bacterivorous nematodes and enchytraeid worms (Figure 3a), and the alteration is correspondingly low (AI = 1.82). Being Anatonchus the most abundant predatory nematode, the lack of herbivorous and fungivorous prey not only affects this predatory species but also enhances the occurrence of bacterivorous prey at a lower trophic level.
• Removal of all the herbivorous and bacterivorous species together.
The carnivorous mite Dendrolaelaps (node 47 of Table S1) and the truly omnivorous mites (nodes 57-61) become extinct. Other carnivores (43-46) succeed in adapting their population size to the much fewer resources available. There is also a slight increase in fungivores and in the remaining omnivores (Figure 3b). At the end, this simulation is resulting in a strong re-assemblage of the food web with an AI equal to 11.53 and the fluctuations in the numerical abundances of so many species were rather unexpected given the lack of response of the food web to the removal of herbivores alone (Figure 2a) or bacterivores alone (Figure 2b). In short, the sum of the effects due to the removal of both guilds was highly relevant even if the separate removal of these two guilds did not influence our web.
• Removal of all the fungivorous and bacterivorous species together.
Carnivores 44, 45, and 46 become extinct, but even the abundance of second most abundant species, the herbivorous Helicotylenchus (node 2), decreases. There is still a slight increase in herbivorous and omnivorous species ( Figure 3c); in fact, AI reaches 12.97. Only the carnivores 43 and 47 seem to succeed in adapting their number to the fewer resources available. As aforementioned, Anatonchus is well known as a stress-resistant nematode, and therefore, its successful performance is expected. The same holds for the mite Dendrolaelaps, because as generalist its trophic niche will be much larger than that of other more specialized carnivores. As in the previous figure, the effect due to the removal of two guilds at the same time was really relevant, although it should be mentioned that these effects are likely due to the removal of fungivores (Figure 2c) more than the removal of bacterivores ( Figure 2b). In any case, the trophic effects due to the combined removal of microbivores in Figure 3c were opposite for the first trophic level (mostly increase in herbivorous species) and the second trophic level (mostly decrease in carnivorous species).
• Removal of either all the omnivorous species or all the carnivorous species. Omnivores are known to have the capability to modify their feeding behavior based on the resource profitability, either by switching prey species (Holt & Polis, 1997)  alteration (AI = 195.75), the highest alteration due to the deletion of a single guild. In contrast, Figure 4b shows that the removal of only the carnivorous species causes only a disproportionally small impact on the food-web architecture (AI = 29.41).
• Removal of all the invertebrates at the second trophic level (omnivorous and carnivorous species). Omnivores and carnivores ate the producers in proportion to their biomass, with omnivores sharing no evident random preference. As expected, an increase in the abundance of all other species at the first trophic level is observed, showing the response of the system to the lack of top-down control ( Figure 5). This scenario is the one with the highest alteration value despite it does not lead to the secondary extinction of any species because the entire upper trophic level was removed. In fact, the removal of omnivores and carnivores causes a significant increase in their prey according to the Equation (13) (2000) and Kratina et al. (2012) that a major effect of omnivory is to lessen the risk of species extinctions following the loss of a herbivore, it is surprising that so many omnivores got extinct. However, missing top-down control by secondary consumers in general, and by omnivores in particular, independently increased the primary consumers and confirmed a recent meta-analysis (Mancinelli & Mulder, 2015). These authors concluded in fact that omnivores, either alone or with predators, exerted a much stronger negative effect than solely predators in terrestrial systems. Seen that the species loss can be assessed when guilds disappear, this method visualizes alterations in the food-web architecture.

ACK N OWLED G M ENTS
The fieldwork, laboratory sampling, and survey were supported by the RIVM Directorate (Dutch grant EIA-S607001) to Christian Mulder. All authors are very grateful to the two anonymous referees for their extensive and highly constructive comments.

CO N FLI C T O F I NTE R E S T
The authors declare no competing interest.