Non-trophic interactions strengthen the diversity—functioning relationship in an ecological bioenergetic network model

Ecological communities are undeniably diverse, both in terms of the species that compose them as well as the type of interactions that link species to each other. Despite this long recognition of the coexistence of multiple interaction types in nature, little is known about the consequences of this diversity for community functioning. In the ongoing context of global change and increasing species extinction rates, it seems crucial to improve our understanding of the drivers of the relationship between species diversity and ecosystem functioning. Here, using a multispecies dynamical model of ecological communities including various interaction types (e.g. competition for space, predator interference, recruitment facilitation in addition to feeding), we studied the role of the presence and the intensity of these interactions for species diversity, community functioning (biomass and production) and the relationship between diversity and functioning.Taken jointly, the diverse interactions have significant effects on species diversity, whose amplitude and sign depend on the type of interactions involved and their relative abundance. They however consistently increase the slope of the relationship between diversity and functioning, suggesting that species losses might have stronger effects on community functioning than expected when ignoring the diversity of interaction types and focusing on feeding interactions only.


Author summary
The question of how species diversity contributes to the functioning of ecological communities has intrigued ecologists for decades, and is especially relevant in the current context of species extinctions. Ecological communities are not only diverse in terms of the species that compose them but also in terms of the way they interact with each other: for example, species compete for space and for food, eat and facilitate each other. The diversity of ways species interact has rarely been taken into account in the study of ecological communities, although widely acknowledged. Here we show that the diversity of interaction types

Introduction
Despite the wide recognition of the coexistence of multiple interaction types linking species in nature [1][2][3], research on ecological networks has been massively dominated by studies on a single interaction at a time (e.g. trophic, competitive or mutualistic; e.g. [4][5][6]). The implications of the diversity of interactions for ecological community dynamics and resilience remains therefore largely unknown, despite a recent growing interest in the ecological literature [7][8][9][10]. Among interaction types, feeding has massively dominated the literature [2], leading to the analysis of the structural properties of food webs on data sets and to the use of modeling to investigate the functional consequences of these structures (e.g. [4,[11][12][13][14][15][16]). Early on, Arditi and colleagues [17] proposed to integrate non-trophic interactions in such dynamical models as modifications of trophic interactions (so-called 'rheagogies'). Building on that idea, Goudard and Loreau [18] investigated the effect of rheagogies on the relationship between biodiversity and ecosystem functioning (BEF) in a tri-trophic model. They showed that ecosystem biomass and production depended not only on species richness but also on the connectance and magnitude of the non-trophic interactions.
Several studies have investigated the role of incorporating specific interactions in food webs. For example, incorporating interspecific facilitation in a resource-consumer model allowed species coexistence in communities of plants consuming a single resource [19]. This increase in species diversity also happens in ecological communities with higher trophic levels including both trophic and facilitative interactions [3]. In the same model, intra-and interspecific predator interference increased species coexistence as well in multi-trophic webs, although to a lesser extent than facilitation among plants [3].
More generally, the joint effect of several interaction types is expected to affect community functioning and stability. Extending May's work, Allesina and Tang [20] showed that communities including a mixture of mutualistic and competitive interactions with equal probability were less likely to be stable than random ones (i.e. where interactions between species are randomly chosen), themselves being less stable than predator-prey communities (i.e. in which interactions come in pairs of opposite sign). Using a similar approach, Suweis and colleagues [21] explored the effect of mixing mutualistic and predator-prey interactions on stability, and showed that, without making any further hypothesis, increasing the proportion of mutualistic interactions tend to destabilize the community. Conversely, in a spatially explicit model including both mutualism and antagonism, Lurgi et al. [9] found that increasing the proportion of mutualism increased the stability of the communities. Addressing the relationship between structure and stability, Sauve et al. [8] showed that the role of nestedness and modularity-structural properties that were shown to promote stability in their single interaction types networks (more specifically in mutualistic networks for nestedness and in antagonistic networks for modularity)-was weakened in networks combining mutualistic and antagonistic interactions. Note that this result contrasts with Allesina and Tang [20]'s result on community matrices who showed that, for mutualistic interactions, nested matrices were less likely to be stable than unstructured matrices.
Combining dynamical models with an empirical network analysis including all known non-trophic interactions between the species of intertidal communities in central Chile [22], Kéfi et al. [10] found that the specific ways in which the different layers of interactions are structured in the data increased community biomass, species persistence and tend to improve community resilience to species extinction compared to randomized counter-parts. More recently, García-Callejas et al. [23] used a dynamical model to investigate the effect of the relative frequency of different interaction types on species persistence and showed that persistence was more likely in species-poor communities if positive interactions were present, while this role of positive interactions was less important in species-rich communities.
Altogether, these studies suggest that the joint effect of several interaction types could alter fundamental properties of ecological systems-such as species coexistence, production and community stability-with however a clear lack of consensus on how. So far, most studies have addressed these questions with specific subsets of non-trophic interactions [3,8,18,19], in small species modules [24,25], in networks with limited numbers of trophic levels [19] or with unrealistic trophic structure [18]. Only a few studies have extended these approaches to complex networks of interactions with a diversity of interaction types (see e.g. [9,23,26]). We therefore still lack a clear view on the overall role of the diversity of interaction types per se for species diversity and community functioning, and especially how they may affect the relationship between diversity and functioning.
In the 90ies, because of the raising awareness of the increase in species extinction rates, the long-lasting interest on the origin and maintenance of species diversity shifted toward the study of the consequences of biodiversity, and especially of its loss, for ecosystem functioning [27]. This became an entire sub-field of ecology referred to as 'Biodiversity and Ecosystem Functioning' (so-called BEF) and lead to decades of experimental and theoretical research investigating how diversity affects functioning (see [28][29][30][31][32][33] for reviews). Results of experimental studies suggests that more diverse communities generally produce more biomass than less diverse ones [34,35]. Theoretically, the question has been addressed as well; models have long focused on plant communities (i.e. a single trophic level) (e.g. [36]), but have more recently started to expand these investigations to more complex, realistic communities (e.g. [37][38][39]). Until now, as far as we know, studies had not specifically investigated the role of the diversity of interactions types on the shape of the BEF.
Here, using a bioenergetics resource-consumer model in which broad categories of nontrophic interactions were introduced [3], we systematically investigated the functioning of 'multiplex' ecological networks, i.e. how multiple interactions (their abundance and intensity) affect species coexistence, community functioning (biomass and production), and the relationship between diversity and functioning. Our model includes, in addition to the consumerresource interactions, competition for space among sessile species, predator interference, refuge provisioning, recruitment facilitation as well as effects that increase or decrease mortality.

The dynamical model
The trophic model. We used an allometric-scaling dynamic food web model [12,40]. These models have been used extensively to explore the dynamics and stability of complex ecological networks [12,13,41]. The food web model consisted of plants (primary producers, at the base of the network, which consume nutrients not explicitly modeled here), and consumers (animals which eat plants and other consumers). The number of species (plants and consumers) and the structure of the web (who eats whom) were initially determined based on the niche model [42] (for details see section 'Simulated networks' in the 'Numerical simulations' part). We mapped dynamical equations to that food web skeleton. The change in species i's biomass density B i (in ½mass� ½area� ) was described by an ordinary differential equation of the general form: where the first term describes plant growth; the second term describes the biomass gained by the consumption of other species j; the third term describes mortality due to predation, summed over all consumers k of species i; the fourth term represents the metabolic demands of species i; the last term is the natural mortality of species i. More precisely: • r i is the intrinsic growth rate of primary producers (in [time] −1 ; r i is positive for primary producers and null for other species); • G i is the growth term described in Eq (2) below; • � o j is a conversion efficiency (dimensionless) which determines how much biomass eaten of resource j is converted into biomass of consumer i; • F ij is the functional response, i.e. the rate at which consumer i feeds on resource j (see Eq (3) below; in [time] −1 ); • x i is the metabolic demand of consumer species i (in [time] −1 ); note that for basal species, metabolic demand is already taken into account in the intrinsic growth rate r i ; • d i the natural mortality rate (in [time] −1 ).

Plant growth.
We assumed a logistic growth for basal species: with K i the carrying capacity of the environment for species i (in ½mass� ½area� ). Functional response. We used a multi-prey Holling-type functional response. The feeding rate of species i on species j is expressed as: where: • w i is the relative consumption rate of predator i on its prey, which accounts for the fact that a consumer has to split its consumption between its different resources (dimensionless); • a ij is the capture coefficient in ½area� ½time� � ½area� q ½mass� q (the attack rate here is a ij B q j which has the unit ½area� ½time� ).
• 1+q is the Hill-exponent, where the Hill-coefficient q makes the functional response vary gradually from a type II (q = 0) to a type III (q = 1) [43] (dimensionless); • h ij is the handling time in ½time� ½mass� .
• Note that m i , the body mass of species i, is needed here for model and unit consistency [39].
Introducing non-trophic interactions. We introduced in this model a number of nontrophic effects found to be frequent ones mentioned in the literature [3,10,22]. We made the relevant parameters of the trophic model become a function of the density of the species source of the effect [3,18]. As a first approximation, we assume all such dependencies to have a similar, linear shape.
Competition for space. We add a space-dependent term, g i , which affects species' net growth rate: where g i -the competition for space term-is evaluated as in Eq (5), l refers to all the species that potentially compete for space with each other, which need to be sessile (excluding intraspecific competition, i.e. l different from i), c 0 is the overall intensity of competition for space and γ il is the strength of competition exerted by species l on species i, which is assumed to increase with the amount of space occupied by each individual of species l (see upcoming subsection on 'Parameter values used' in the 'Numerical simulations' part). Note that the element γ il is zero if either i or l is non-sessile and even if both species are sessile it is non-zero only with a certain probability (see 'Simulated networks' in 'Numerical simulations' part). This makes competition for space asymmetric, as some species can have a large negative effect on others but are not negatively affected themselves. Also note that because plants have a higher probability of being sessile than other species in the web (see upcoming subsection on 'Simulated networks' in the 'Numerical simulations' part) competition for space occurs more frequently among plants than among other species of the network. Finally, note that if γ il is null for all l, Eq (4) is identical to Eq (1). Competition for space was assumed to only operate if the net growth rate of the target species (the first term in between brackets in Eq (4)) is positive, i.e. if ðr i G i þ P j2prey � o j F ij À x i Þ > 0. 0therwise, g i is set to 1.
Predator interference. We introduced predator interference in the feeding rate as follows: where s denotes all the predators of prey j and δ si is the strength of interference competition between predators s and i (Beddington-DeAngelis type, [12,44]). Again, even if two predators share a prey, this term is non-zero only with a certain probability, and its values is assumed to depend on the differences of the body mass of the two predators (see upcoming subsection on 'Parameter values used' in the 'Numerical simulations' part). The constant i 0 is the overall intensity of predator interference. Effects on mortality. A number of negative interactions lead to a decrease in the survival of the target species i (e.g. whiplash). Some species might also increase the survival of target species (e.g. improvement of local environmental conditions). We summarized these two types of effects on the mortality rate, d i as follows: with n 0 and p 0 the overall intensities of negative and facilitative effects on mortality (i.e. resp. increase and decrease in mortality), and n and p the interaction matrices containing zeros and ones.
Refuge provisioning from predators. Refuge provisioning can happen in different ways: a species can protect another from abiotic stress (e.g. decreasing its mortality-see previous example) but a species can also protect another from its predator (e.g. affecting the attack rate of the predator). A refuge provision from species k to species j from its predators can be modeled as follows. For all the predator species i of j: where a ij new tends to 0 in the presence of facilitators, r 0 is the intensity of the refuge effect and ϕ is the interaction matrix containing zeros and ones. Effects on recruitment. Species may increase (e.g. habitat amelioration) the recruitment of new plants in the community. We therefore created the term e i which is multiplied by the growth rate of the species: with e 0 the overall intensity of facilitative effects on recruitment, and η the interaction matrices containing zeros and ones. Note that this term only applies to plants.

Numerical simulations
Simulated networks. Trophic networks were generated with the niche model [42] starting with 100 species including a fixed number of 20 plants (primary producers) and a connectance of 0.06 (i.e. about 600 trophic links per network) [45]. Cannibalism was allowed but trophic networks containing cycles were discarded.
We imposed a fixed number of 33 sessile species in each network: for each species, we uniformly drew a 'mobility' trait with probability 0.2 for plants and 0.8 for other species, and we repeated the procedure until the network contained 33 sessile species. As a consequence of this choice of probability, plants have a much higher chance of being sessile in our networks than other species in the web.
The location of the non-trophic links was chosen randomly in the trophic web, but following a number of basic rules inspired from the Chilean data set [10,22]. Competition for space was drawn between two sessile species, interference between two mobile predators that have at least one prey in common, refuge provisioning from a sessile species to a prey (i.e. a species that has at least one consumer), and recruitment facilitation from any species to a plant. Effects on mortality were drawn between any pair of species (i.e. no constraint). Non-trophic links were drawn only between different species because we focused on the role of inter-specific interactions. Note however that our simulations all include intra-specific interference with i 0−intra = 0.8, as in [39], without which the species diversity at steady state is much lower.
With the previous rules satisfied, the interaction probability (i.e. the probability for a nonzero element of the corresponding non-trophic interaction matrix) between two species was respectively set to 0.098, 0.15, 0.01, 0.01, 0.033 and 0.063 for competition, interference, increase and decrease in mortality, refuge provisioning and recruitment facilitation. These settings allowed to get, on average, 100 non-trophic links for any of the six types of non-trophic interactions (because of the imposed rules, the probabilities need to be different for each interaction type). This means that a simulation started with about 600 trophic and 100 non-trophic links of a given non-trophic interaction type.
Simulations setup. In the first part (Fig 1; see also S1 and S2 Figs), simulations were first run with trophic links only, and 100 trophic networks were selected in which no disconnected plant (i.e. plants that have no consumer) was present at the end of the dynamics (as in [39]). For each of these trophic networks, we drew 50 non-trophic networks of each non-trophic interaction type (to study the effect of each type of non-trophic interaction individually). Again, we only kept the networks in which no disconnected plant was present at the end of the dynamics with the non-trophic links. For each non-trophic interaction type, we defined a range of non-trophic intensity values as follows. We started from a minimum non-trophic intensity value, and we linearly increased this value to a maximum. The minimum (resp. maximum) values were selected so that they correspond to a 2.5% (resp. 10%) variation in species diversity at the end of the simulation in the case with compared to without NTIs. This was the case for all NTIs except for positive effects on mortality for which reaching an effect of 10% change in diversity was not possible even for very high values of p 0 (see colored lines in Fig 1). With this procedure, we put all non-trophic interactions on equal footing, which allowed comparing their effect on outcome variables. This procedure allowed us to compute the slope of the effect of the non-trophic interaction intensity on final species diversity using a linear regression, which is an indicator of the strength of the non-trophic interaction. The slope of this regression is displayed on top of each panel of Fig 1. In the second part (Figs 2 and 3; see also S4 Fig), decrease in mortality (also referred to as 'positive effects on mortality' in the text) were discarded because of their lack of significant effect on species diversity. For the other five non-trophic interaction types, we selected 1000 simulated trophic networks (again with no disconnected plants), and for each of these networks, we repeated the following procedure 100 times: we drew five non-trophic networks (one per non-trophic interaction type), with one quarter of the previous interaction probabilities for the four negative non-trophic interaction types. This way, we ensure that we have added an equal number of positive and negative links in the networks (i.e. there are as many positive as negative links in the networks studied). We uniformly drew the non-trophic interaction intensities in the same range of values as in the first part (see above). Hence, in this set of simulations, we started with about 600 trophic links, 100 positive (facilitation for recruitment) and 100 negative non-trophic links (about 25 for each of the four types: competition for space, interference, refuge and increase in mortality). In the sensitivity analysis (S5 Fig), we reproduced the same procedure for each of the parameter values investigated (namely the Hill coefficient, expo and capture coefficient a 0 -see after).
In the same vein, we did another set of simulations in which the intensities of the non-trophic interactions were fixed (to the values leading to 10% variation in species diversity; i 0 = 3, r 0 = 1.75, c 0 = 0.012, e 0 = 1.8 and n 0 = 3), but we varied the number of links of each of the non-trophic interaction type (S3 Fig). We simulated 100 trophic networks with no disconnected plants and, for each trophic network, we drew 100 non-trophic interaction networks for each non-trophic interaction type of varying size: by modulating the previously mentioned Simulation runs. We used the GNU Scientific Library (https://www.gnu.org/software/ gsl/) solver with the embedded Runge-Kutta-Fehlberg (4,5) method. For each network, numerical simulations were run until steady state was reached (we set a maximum time t = 5000 which we observed to be sufficient). During the dynamics, we set a species to extinction when its biomass was very small (<1e-6). To fairly compare results with and without non-trophic interactions, we used the same initial conditions in both cases. At steady state, we measured: diversity, i.e. the number of surviving species (biomass> = 1e-6), total biomass (sum of the biomass of all surviving species at steady state) and total production (the sum of the intrinsic growth of basal species and food uptake minus the respiration of consumers, over all surviving species, i.e. first term in Eq (4)). We also evaluated the normalized ratio of each of along the x-axis were chosen such that it lead to a ±2.5% (green dashed line; minimum value of the parameter range) to ±10% (blue line; maximum value of the parameter range) change in species diversity, relatively to the case without NTI. Note that the y-axes of the different panels differ. The NTIs were categorized into 'positive' (i.e. beneficial, recruitment facilitation) vs 'negative' (i.e. detrimental; interference, refuge provisioning, competition for space and increase in mortality) based on their effect on diversity.
https://doi.org/10.1371/journal.pcbi.1007269.g001 these metrics in the case with and without non-trophic interactions. For instance, we call diversity ratio (in particular in the legends of the figures) the difference between the diversity with and without non-trophic interactions divided by the diversity without non-trophic interactions.
Parameter values used.
• � 0 j was set to 0.45 if the resource j is a plant and to 0.85 otherwise; • γ ij , the effect strength of competition for space of species j on species i, is assumed to depend on body mass such that g ij ¼ g 0 m 2 3 j , with γ 0 = 1.
• δ ij is the term of interference between predator i and predator j; we assume that the more similar the body masses of the two predators, the stronger the interference between them, such that d ij ¼ • the Hill exponent was tested with q = 0.3 and 0.7 (in addition to the value 0.5 used for the main simulation) • we tested a 0 ¼ a 0 sesscons ¼ a 0 sessres ¼ 10 (in this case, we had to set h 0 = 0.1 to avoid massive extinctions) and 250 (in addition to the value 50 used for the main simulation) • the parameter expo was set to 25 and 75 (in addition to the value 50 used for the main simulation) We focused on these three parameters because they are linked with both the largest variation in real ecological systems and the largest measurement uncertainties. Also, they are known to strongly affect biomass flow patterns and dynamical stability in food webs. Others parameters of the bioenergetics model are either fixed by defining the scales for time and biomass density (e.g. r 0 ) or do not show much variation in natural systems (e.g. the assimilation efficiencies).
The ranges over which the parameters were varied in the sensitivity analysis were chosen as broad as possible, but with the constraint of still enabling enough species survival in the food webs (i.e. in the networks without NTIs). This restriction is useful as the aim of this study is not to explore conditions for persistence in food web models, but to focus on the effect NTIs have on diversity, ecosystem functioning, and the relationship between both.

Results
In what follows, we use NTI(s) to refer to non-trophic interaction(s).

Effect of the presence and intensity of each NTI
We ran community dynamics with or without NTIs, and evaluated the relative difference in community characteristics at steady state obtained in the presence compared to in the absence of each NTI. This allowed comparing the effects of the different NTI types (for a range of interaction intensities; see Methods). At steady state, the measured characteristics of the communities were: species diversity (the number of species which survived at steady-state, i.e. whose biomass was above a threshold level), total biomass (the sum of the biomass of all surviving species at steady state) and total production (the sum of the intrinsic growth of basal species and food uptake minus respiration of consumers, over all surviving species; first term in Eq (4) in Methods, 'The dynamical model').
The following NTIs were introduced, one at a time, in the consumer-resource model: i) predator interference, which can occur between two predators which share at least one prey, ii) refuge provisioning which can happen if a species protects another from its predator (e.g. affecting the attack rate of the predator), iii) competition for space which occurs predominantly between sessile species, iv) recruitment facilitation which happens when some species increase the recruitment of new plants in the community (e.g. by habitat amelioration), v) increases in mortality when some species decrease the survival of others (e.g. because of whiplash) and vi) decreases in mortality when some species increase the survival of others (e.g. by improving the local environmental conditions).
We found that interference had a negative effect on diversity and community production and a weak (negative) effect on biomass (1st row of Fig 1). Through time, interference decreased the consumption of some of the predators; this initially favored some of the basal and intermediate species (that were less consumed), and eventually lead to the extinction of some of the intermediate and top predators, as well as to a decrease in their total biomass (1st rows in S1 and S2 Figs). Primary producers, some of which were relieved from consumers, exhibited a slight gain in biomass.
Refuge provisioning had similar overall effects, but with a larger effect on biomass than the one of interference (2nd row of Fig 1). In this case, species benefiting from refuges remained in the system but were less accessible resources. This lead to a loss of biomass and subsequent extinctions of some consumers (which could not access their prey; especially top predators), while their resources remained under protection and gained a bit of biomass (except for those whose protector went extinct) (see S1 and S2 Figs, 2nd row).
Competition for space had a strong negative effect on all variables (which affected all trophic levels), while recruitment facilitation had a positive effect on all community characteristics (but affected only consumer and predator species; 3rd and 4th rows of Fig 1 and S2 Fig). Through time, these effects tend to first affect the basal species, then the intermediate and eventually the top predators (see S1 Fig).
Modifications of mortality rates produced very weak effects overall. Increasing mortality had a negative effect on diversity and biomass and no effect on production (5th row of Fig 1).
Decreasing mortality had a very weak positive effect on diversity and biomass and no effect on production (6th row of Fig 1).In what follows, we did not consider decreases in mortality (also referred to as 'positive effects on mortality') which was the weakest of all NTIs considered overall and focused instead on the five remaining NTIs, namely interference among predators, refuge provisioning, competition for space, recruitment facilitation, and increase in mortality (also referred to as 'negative effects on mortality').
Overall, the most influential NTIs among the ones studied were competition for space and recruitment facilitation, in terms of both diversity and functioning (see slopes linking the parameter values to see the extent of the effects in Fig 1). The effect of competition for space on diversity was two orders of magnitude larger than those of all the other NTIs (namely interference, refuge provisioning, recruitment facilitation and increase in mortality). Regarding biomass, the effect of competition for space was two orders of magnitude larger than the one of recruitment facilitation, which was itself an order of magnitude larger than the one of all the other NTIs (namely interference, refuge provisioning, and increase in mortality).
For all NTIs except for competition for space, effects seemed to be stronger on intermediate and top trophic levels at steady state (see S2 Fig). Regarding species diversity, this was partly due to the fact that plant species already all persisted with trophic interactions only, andbesides competition for space (and to a much lesser extent increases in mortality)-the other NTIs were not able to lead to plant extinctions, because their effects either corresponded to a decrease in plants consumption (interference, refuge) or to a positive effect on plants (recruitment facilitation, decrease in mortality). Therefore, the NTIs studied here had very little opportunities for affecting plant species diversity. Conversely, the NTIs studied had more leverage on intermediate and higher trophic levels where species did not all persist in webs with trophic interactions only, and where they could therefore either increase or decrease species diversity. Regarding biomass, effects seemed to first affect basal species but then climb up the food web to eventually affect the top predators more strongly (see S1 Fig).
This first set of simulations helped us categorize the NTIs studied into 'positive' (i.e. beneficial; recruitment facilitation) vs 'negative' (i.e. detrimental; interspecific predator interference, refuge provisioning, competition for space and increase in mortality) based on their effect on diversity.

Combined effects of the NTIs on species diversity
Next we mixed the five remaining NTIs together, with NTI intensities picked at random within predefined ranges to study the joint effect of the NTIs considered. These predifined ranges were chosen so that each NTI increases or decreases the diversity of the system by 2.5% to 10% compared to the case without NTI (see Methods and Fig 1). Pre-defining these ranges for each of the NTI taken individually allows to put all NTIs on comparable grounds.
Not unexpectedly, the effect of the presence of the NTIs depended on the relative number of links of the different NTIs and on their intensities. When all interaction types were together with an equal proportion of positive and negative NTIs, networks with NTIs tended to have a smaller species diversity than networks without NTIs (Fig 2A). In other words, NTIs lead to extinctions of species compared to simulations run with feeding interactions alone. There were also quite a few number of cases where the net effect on diversity was null.
There was nonetheless a large fraction of cases where NTIs tended to enhance species diversity; these were clearly cases where beneficial NTIs were present and strong (orange and red areas on Fig 2B). It was noteworthy that the NTI values were all chosen at random for each of the simulations, so all combinations of intensity values were possible and present across simulations, but our results showed that positive effects of NTIs on diversity always happened when the beneficial NTI (recruitment facilitation) was strong while the detrimental NTIs were weaker. Now fixing the intensities of all NTI links to their maximum value (corresponding to a 10% effect on diversity ratio; see Fig 1) and focusing on their relative abundance, we found that a greater number of recruitment facilitation links tend to favor positive effects on diversity while increasing the number of interference, refuge or competitive links pushed toward negative effects on diversity (S3 Fig).

Combined effects of NTIs on the Biodiversity-Ecosystem functioning relationship
How did these effects on species diversity translate into community functioning? Using the previous simulations where NTI intensities were picked at random, we found that both in food webs (i.e. in ecological networks without NTIs) and in ecological networks with NTIs, the relationship between species diversity and biomass at steady state was positive (Fig 3A and 3B; this was also the case for production: see S4 Fig)Strikingly, in presence of NTIs, the relationship was significantly stronger than in their absence (ANCOVA p-value<1e-16; comparing slopes in Fig 3A and 3B). We checked that this result was robust to changes in the value of major parameters of the model (namely the Hill exponent which determines the shape of the functional response, the parameter expo which determines how species body mass depends on their trophic level, and the capture coefficient a 0 of consumers; S5 Fig).
Plotting the biomass ratio (i.e. the variation in biomass with NTIs compared to without NTIs) as a function of the species diversity ratio suggested that when NTIs contributed to a gain in species, this generally translated into a gain in biomass as well (Fig 3C). Actually, networks with NTIs tended to gain biomass (compared to networks without NTIs) even when there was no gain (or even a weak loss) in diversity (see boxes at diversity ratios of -0.1 and 0 in Fig 3C).When there was a small loss of diversity in presence of NTI (-0.1-0), the remaining species took advantage of these extinctions and gained biomass. When there was a gain in species diversity compared to the case without NTIs, this was often happening because of the presence of beneficial NTIs (Fig 2B), and those beneficial links lead to a considerable increase in biomass as well. There was, however, a large variability around these trends due to the fact that each simulation corresponded to a different combination of NTI intensities.

Discussion
Using a bioenergetic model in which six types of NTIs were incorporated, we found that these NTIs in isolation and jointly affected significantly species diversity and community functioning (biomass and production), consistently with previous studies addressing the role of the diversity of interaction types in module or network contexts [3,10,[17][18][19][20]23].
Overall, when taken together and with a balanced number of beneficial and detrimental interactions (as defined by their individual effects on diversity), the presence of NTIs tended to have a slightly negative effect on species diversity. This is in agreement with Goudard and Loreau [18] who studied NTIs that are modifications of feeding links; with equal numbers of positive and negative effects, they found a decrease in the total number of species when NTIs are incorporated. In our case, we did not expect this result since the range of NTI intensities spanned was selected such that each interaction type had equivalent effects on diversity when taken individually. Yet, despite controlling for both the number and the intensities of NTIs, the joint effect of the NTIs, when simultaneously incorporated in the model, tended to be negative for the species diversity of the resulting communities. It is noteworthy that the only beneficial NTI studied in this model, recruitment facilitation, operates on a single trophic level, namely the plants. When considered alone (i.e. without any other NTI), we showed that the positive effect of recruitment facilitation on diversity of the entire networks was entirely due to indirect effects on consumer species (trophic level > 2, S2 Fig). Thus, the negative joint effect of NTIs on diversity suggests that recruitment facilitation has less leverage on diversity than the other NTIs, and more specifically, that the fact that its positive effect on diversity comes about only by indirectly affecting species on higher trophic levels, does not allow it to compensate for the more direct negative diversity effects of some of the detrimental NTIs.
Surprisingly, we found interference between predators to have a negative impact on diversity, which contrasts with other studies reporting stabilizing effects on population dynamics and positive effects on diversity [12,49,50] when interference is included via a Beddington-DeAngelis functional response [51,52]. However, these studies included interference either as a purely intra-specific effect or assumed that intra-specific interference was stronger than inter-specific interference. Here, intra-specific interference was present in all simulations at a fixed value, and our aim was to study the effect of changes in the intensity of inter-specific interference. In this sense, our result that interference reduces diversity is reflecting classic results from competition theory, namely that competition is destabilizing if it is stronger between species than within species [53]. As interference was only included for predators that are already competing for at least one common prey species, the decline in diversity can be attributed to an increased effect of the competitive exclusion principle [54]. It has to be noted, however, that the complex network structure of trophic and non-trophic interactions provides a plethora of niches, which reduces the direct applicability of this principle [13].
Interestingly, we found that NTIs affected the relationship between diversity and functioning, and this result seems to be robust to changes in the value of key parameters of our model. Despite the array of possible effects from the NTIs, individually and combined, the relationship between species diversity and biomass was found to be significantly stronger in networks with NTIs than in networks without them. Again, this was not necessarily expected since simulations were run with all NTIs together, whose intensity values were picked at random in a range such that their effects on diversity was controlled; we could therefore have expected, e.g. compensatory or negative effects since most of the NTIs studied here tended to have negative effects on diversity (Fig 1). The effect of NTIs on the slope of the diversity-biomass relationship means that when species-rich networks gain even more species, that goes with a disproportionately higher gain in biomass in the presence than in the absence of NTIs. This also means that, conversely, species-poor communities lose more biomass with additional species loss with than without NTI. This is due to the fact that species-rich communities are communities in which beneficial interactions are present and strong, while species-poor communities are communities where detrimental NTIs operate. This result is interesting in that it suggests that species loss may have stronger consequences on community functioning than expected if ignoring non-feeding interactions. Further work is needed to see if this result extends to other models as well as to real ecosystems.
Of course, our study presents a number of limitations. The strongest NTIs studied here, namely competition for space and facilitation for recruitment, both affect mainly plants in our model. It is therefore unclear whether these NTI types appear to exert stronger effects because plants are the affected species. This could be a topic of further investigations.
Moreover, we have focused on a selection of six NTIs that are the major ones found to occur in the Chilean web [10,22], but other interaction types not present in this data set are known to be frequent and important in nature. Examples include parasitism, effects on resource availability, plant dispersal or animal movement. Further work could introduce these other interactions in a single framework.
We also assumed here the intensity of the interactions between pairs of species to be constant through time. Some interactions may however be context-dependent (beyond the biomass or abundance of other species). For example, adaptive inducible defenses are a form of phenotypic plasticity that affects the strength of predator-prey interactions [55]. Changes in morphological, behavioral, or life-historical traits in response to chemical, mechanical or visual signals from predators have been reported in the literature for a number of organisms [56]. These responses can moreover occur with a lag, given that the expression of defenses may involve considerable time, relative to the organisms life-cycle [57]. The intensity, and even the type of interactions, could also change with e.g. changes in abiotic factors such as climate.
We have no information regarding the relative importance or intensities of the different interaction types. We proposed a way of putting all interaction types on equal footing regarding their effect on species diversity. This is however a debatable choice-we could for example have chosen to make NTIs comparable regarding their effect on biomass. In nature, it is likely that interaction intensities are not equivalent and that some of them are much stronger than others. Making progress along these lines requires experimental work aiming at quantifying different interactions types, which involves a number of challenges [1].
In this study, NTIs were plugged in the food web randomly although with a number of constraints based on our knowledge of the Chilean web [10,22]. We did not explicitly investigate the role of the structure of the NTI network, despite the fact that previous studies have suggested that it might play an important role [10,58]. This remains a difficult task since it can be necessary to take into account the dependency between the structure of the different layers (e.g. NTI types, as observed in [10]); however this is a promising avenue of future research.
Previous studies have used the community matrix approach focused on net effects between species to investigate the role of the diversity of interaction types [7,20]. This approach has a number of advantages, including the fact that it allows analytical predictions and generalizations. Here, we chose to focus on a more mechanistic approach, starting from the mechanism of the NTIs without assuming their net effect. For example, we had initially assumed that refuge provisioning would be a beneficial NTI, meaning that it would have a positive effect on species diversity. We however found the opposite in the model simulations-indeed, refuge provisioning protects prey from their consumers but it also deprives consumers from their resource and it seems that this latter effect has stronger consequences at the community scale. In a dynamical model, Gross [19] showed that interspecific facilitation among plants allowed the maintenance of species diversity despite the fact that the net effects measured among plants remained negative. These results highlight that insights gained from the analysis of few-species systems cannot be easily translated into the dynamics of complex communities. Focusing on the net effects between species may conceal important coexistence mechanisms when species simultaneously engage in both detrimental and beneficial interactions and stresses the importance of working with mechanistic models to better understand the consequences of NTIs for community diversity and functioning. Nonetheless, a more mechanistic approach implies, for each of the NTI identified, to model it in one specific way, matching our knowledge of how these interactions operate in nature (in our case here, having the Chilean web in mind [10,22]). For example, regarding refuge provisioning, alternative ways of how it could affect trophic interactions are certainly conceivable. It could increase the Hill coefficient of all predators of the protected prey, e.g. to mimic the fact that prey at low density would become better protected, but if prey density is too high, the predators would still see them. This could have a different effect on diversity than found here.
Our study is a step toward getting a better understanding of the dynamics of multiplex ecological networks (i.e. including several interaction types among a set of species), and more precisely of the role of NTIs on community functioning. Our model results suggest that, when simultaneously included, and assembled according to simple rules reflecting observations in nature, NTIs tend to mechanically strengthen the BEF, making the dependency between the number of species present in the community and the functioning of this community (in terms of biomass or production) stronger. This result has important consequences for predicting the consequences of species loss on community functioning.
Supporting information S1 Fig. Relative change in species diversity (left column) and biomass (right column) per trophic level through time for 500 simulations for the main NTI types (5 first rows) and the NTI all together (6th row). Each simulation starts with 600 trophic links and 100 nontrophic links. Note that on the x-axis, time is on a log-scale. Green: TL = 1, yellow: TL between 2 and 3, red: TL> 3. TL refers to the prey-averaged trophic level measured as one plus the mean trophic level of all the species resources, where the trophic level of a resource is the chain length from the resource to a basal species [59]. Species whose TL is 1 are primary producers. (PDF)

S2 Fig. Relative change in each of the measured variables (3 columns) as a function of the intensity of the NTI by trophic level, for each NTI type (6 rows).
Note that each NTI has its own y-axis. Green: TL = 1, yellow: TL between 2 and 3, red: TL> 3. TL refers to the prey-averaged trophic level measured as one plus the mean trophic level of all the species resources, where the trophic level of a resource is the chain length from the resource to a basal species [59]. Species whose TL is 1 are primary producers. Note that the y-axis differ for the different NTIs. (PDF) S3 Fig. Frequency of simulations leading to a given diversity ratio, i.e. variation in species diversity in the case with compared to without NTI links. (to match the term defined in the Methods). The intensity of each NTI type is now fixed and we vary the relative number of links of different NTI types when put together. The left column correspond to situations where there are twice more detrimental than beneficial links, the middle column shows results with equal number of beneficial and detrimental links, and the right column correspond to cases where there are twice more beneficial than detrimental links. The panels correspond to different combinations of NTIs: I for interference, E for recruitment facilitation, N for negative effects on mortality, R for refuge, C for competition. At fixed intensity and with equal number of links, detrimental links tend to take over (slightly). There are configurations of relative abundance of the different types of NTIs in which positive and negative effects can balance each other. (PDF)  Fig 3A and 3B to changes in the value of major parameters of the model, namely the Hill coefficient, which determines the shape of the functional response, the parameter expo which determines how species body masses depend on their trophic level, and the capture coefficient a 0 of the consumers. See Methods, part 'The dynamical model', for more details about the parameters and how they contribute to the dynamical equations. Independently of the combinations of parameter values found, the slope of the BEF is stronger in the presence than in the absence of NTIs. (PDF)