The structure and robustness of ecological networks with two interaction types

Until recently, most ecological network analyses investigating the effects of species’ declines and extinctions have focused on a single type of interaction (e.g. feeding). In nature, however, diverse interactions co-occur, each of them forming a layer of a ‘multilayer’ network. Data including information on multiple interaction types has recently started to emerge, giving us the opportunity to have a first glance at possible commonalities in the structure of these networks. We studied the structural features of 44 tripartite ecological networks from the literature, each composed of two layers of interactions (e.g. herbivory and pollination), and investigated their robustness to species losses. Considering two interactions simultaneously, we found that the robustness of the whole community is a combination of the robustness of the two ecological networks composing it. The way in which the layers of interactions are connected to each other affects the interdependence of their robustness. In many networks, this interdependence is low, suggesting that restoration efforts would not automatically propagate through the whole community. Our results highlight the importance of considering multiple interactions simultaneously to better gauge the robustness of ecological communities to species loss and to more reliably identify key species that are important for the persistence of ecological communities.


Fig A:
Map showing the location of the 6 studies we used to build the database of networks.The color and shape of the points represents the different types of networks, and the text the name of the site, the country, and the surname of the leading author of the study.The map has been done by the authors using package "rnaturalearth" [2] in R [3] .This is also true when looking at the average degree heterogeneity by species set (see Fig CA and C).This high degree heterogeneity means that empirical tripartite networks, contrary to the expectation in random networks where all nodes have a similar connectivity, have few nodes with a high connectivity (generalists) and many nodes with a lower connectivity (specialists), as shown in Fig D .This 'long tailed' degree distribution is a well established feature of plant-pollinator [9,10], plant-seed disperser [9] and, to a lesser extent, food webs [11,12] that we recover in the tripartite networks in our data-set as well.Finally, looking only at the degree heterogeneity of the nodes in the shared set, we find that in 50% of networks the degree heterogeneity of the nodes in shared set in the tripartite network is higher than the degree heterogeneity of the subset of shared set nodes in the two bipartite networks (Table B), and in all cases the degree heterogeneity of the shared set nodes in the tripartite is at least larger than that of one of the two subsets (Table B).This means that the two layers are connected in such a way that the degree heterogeneity of the shared set nodes is increased.Horizontal grey lines mark the limits of the confidence interval of 1.96 (95%) and 2.33 (98%) for rejecting the null hypothesis.for the three different sets of species, from left to right: insect herbivores, insect pollinators, plants and all nodes combined (i.e. the merged network).The continuous black line is the empirical cumulative degree distribution, with the best power law distribution fit in a dashed blue line.The continuous red line is the cumulative degree distribution in null model "1" ensemble (100 randomizations), and the dashed red line the best exponential fit to that distribution.The insets show the degree distribution of the empirical network (in grey) and that of 100 randomizations (in red).
Apart from this heterogeneity in the species degree, in empirical networks the species are not randomly connected, rather the degree of a node can depend on the degree of its neighbours.This is known as 'mixing patterns' or degree-degree correlations, and is measured using correlation coefficients between the degree of a node and the degree of its neighbours.Pearson's correlation coefficient between the degree of a node and the average degree of the node's neighbours (r) was found to be negative in all empirical networks, indicating that they are all disassortative, i.e. that species with many connections tend to be linked to species with few connections (see Fig BB).Since degree heterogeneity is known to promote degree disassortativity [13], we also studied if the empirical networks were more disassortative than expected given their degree distribution by comparing the observed values to those of the constant "K" ensemble of null model 4. In line with previous findings in ecological plant-mutualistic networks [13], we find that most of the observed tripartite networks are not significantly more disassortative than the null expectation (with the exception of MM networks, which exhibit a higher disassortativity than expected), as the value of r in the empirical networks is consistent with those of the "K" ensemble average (null model 4) with a confidence con 95%(98%) in 86%(90%) of the networks (

see Fig BD and Table A above).
When we look at the interaction layer scale (i.e. at the two unilayer networks composing the tripartite network), we find that all pollination, dispersion, ant-mutualism and a majority (80%) of herbivory networks are disassortative, while most of parasitoids networks (80%) are assortative (i.e. have r > 0), although most herbivory and parasitism networks exhibit weak degree-degree correlations (values of r close to 0), as shown in Fig CB .However, as with the tripartite networks, when we compare with the randomizations in null model 4 they tend to not be significantly more assorative/disassortative than the null expectation (Fig CD ), with 75%(91%) or interaction layers having a value of r compatible with that of the random ensemble with a 95%(98%) confidence interval.This means that most of the assortativity/disassortativity found in the empirical networks stems from their low/high degree heterogeneity.B: Degree heterogeneity of the shared set species in the tripartite networks in the data-set.Different columns show the name of the layer, the number of species in the layer (N), the degree heterogeneity of all shared set species in the tripartite network (LS_HD), and also for the subset present in each of the two bipartite networks that compose the tripartite.Networks marked with a star indicates those where the degree heterogeneity of the shared set species in the merged (i.e.tripartite) network is larger than that of the subset of shared set species included in each of the two bipartite networks composing the tripartite.

Interdependence
We study how the Robustness of the two animal species sets are correlated in the tripartite networks of our database when plants were randomly driven to extinction.A positive correlation means that the same plant species that are important for the survival of the two sets of animal species (i.e.plants important for pollinators are also important for herbivores).However, if the relevant plants are different, deleting a given plant will cause a large loss of animal species in one set, but not have much effect in the other set, rendering no interdependence.
We find that, in general, the correlation between the Robustness, henceforth called interdependence (I) was either positive or null in all the networks.However, some MA and MM networks (∼20% of all networks in the database) show a negative correlation (albeit very weak in most cases).The value of I found in AA networks (Fig EA) is on average significantly higher from that found in MM and MA networks, mirroring previous results where the correlations between species sets that where linked by feeding interactions were found to be stronger [7].
Comparison with the four null models indicates that, for AA networks, the positive interdependence may be a result of the process in which secondary extinctions take place (i.e.cascading extinctions), since the empirical value is not significantly different from the null expectation in most cases (see Fig E B to E).As for MM and MA networks, comparison with all null models shows that the structure of empirical networks (in particular the degree heterogeneity of the species, Fig E B and C) makes the two layers more decoupled than expected by chance, since the null models predict in fact a negative Interdependence.In the two null models where degree heterogeneity is conserved the empirical interdependence is not different from the null expectation in any case (Fig E D and E).

Effect of structural features on interdependence
The different ways in which the two interaction layers are connected is reflected also in how the structural features of the connector nodes are able to explain the interdependence.In MM and MA networks, these structural features play an important role: H C and P R c alone account for ∼55% and ∼20% of the variance in interdependence, respectively (see Fig FE and F).In AA networks, in contrast, the structural metric that best explains the variance in interdependence is C, but only ∼16% of it (see Fig FD).In MM and MA networks, a multiple linear regression explains 64% of the interdependence as a function of H C (how frequently hubs are connector nodes) and P R c (the average participation ratio of the connector nodes) and degree heterogeneity, as shown in Table C.No structural metric explains more than 16% of the variance in interdependence in AA networks.
Regarding the presence of (weak) anti-correlations found in MA and MA networks, one could think that their origin lies in the fact that not all the targeted plant species are shared between the two layers in MA and MM networks.In those cases, secondary extinctions can happen in only one of the two interaction layers, while the other remains intact, introducing a  In fact, we find that the quantity that best correlates with I in MM and MA networks (r ≈ 0.75) is the tendency of the 20% of the most connected nodes in the shared set to be connector nodes (FigFE).That is, the fact that the generalist species in the shared set tend to be involved in the two interaction layers or not makes the difference in terms of interdependence in MA and MM networks.If the generalist plants tend to be involved in both interactions, then the interdependence is larger.If, on the contrary, the generalist plants are involved in only one of the two interaction layers, then we see no correlation, as in MM networks, or even an anti-correlation.In AA networks, on the contrary, the structural metric that best correlates with I (albeit weakly) is the proportion of connector nodes in the shared set (C). Since, as said before, for a parasitoid to undergo extinction the herbivore host has to go extinct first, it is not a surprise that the higher the number of connectors, the higher the interdependence of the two interaction layers.

Robustness
All types of tripartite networks seem to react in a similar way (with Robustness in the same range of values) to the three extinction scenarios explored, although some differences between network types were significant (see Fig GA to C).In line with previous studies on networks with only one interaction type, all ecological tripartite networks were most fragile when plants were selectively attacked targeting the most connected plants first (DD scenario), and the least fragile when plants were attacked selecting the specialists plants first (ID scenario), as previously reported in mutualistic plant-pollinator, food-webs and host-parasite networks (albeit fish-parasites, not herbivore-parasites) [14,10,15,16,17].In the random extinction scenario (RND scenario), Robustness showed intermediate values, with AA networks being the more fragile of the three different type of networks and MM networks the more robust (see Fig GA).We did a multiple linear regression of the Robustness (R) in the RND extinction scenario as a function of the different structural features (Table D)

Effect of basic structural features in Robustness
Degree heterogeneity is associated with an increase in the Robustness of all tripartite networks in the RND and ID scenarios (r RN D = −0.62,r ID = −0.75),but increasing it in MA and MM networks in the DD scenario (r DD = 0.73).The positive effect of degree heterogeneity on robustness was already reported for bipartite mutualistic networks in [18] (through nestedness, but it was also shown that nestedness is consequence of degree heterogeneity [13]).These results are also in line with previous studies showing that ecological networks are quite robust to the extinction of the most specialist species, but quite fragile if the most generalized species were the ones going extinct [11,14,10].This robustness to random extinctions has largely been explained by their heterogeneous structure, in food webs [11] and in plant-mutualistic networks [10], and is also present in the tripartite networks in our data set.

Relevance of the empirical structure
To gauge the relevance of the empirical structure in determining Robustness in the RND extinction scenario, we tested how the Robustness of the empirical networks compared to that of their different randomizations.AA empirical networks where consistently more robust than their randomizations when the original degree distribution was lost (i.e., in all null models except "4", see Fig JB to D). Empirical MA and MM networks, on the other hand, were less robust than their randomizations in ensemble "1" (and MM networks also in ensemble "2").This indicates two things.First, that the degree distribution is a mayor driver of Robustness, since as long as it is kept constant (null model "4") the Robustness of all types of empirical networks is not significantly different from the random expectation.Second, that the degree distribution of empirical AA networks is such that makes them more robust to the loss of plant species.
We also compared the robustness of the MA and MM tripartite networks with that of the two bipartite networks that compose the tripartite to see if considering multiple interactions at the same time enhanced or hindered the robustness of the whole community.This analysis was only possible in MA and MM networks, since in AA networks plants appear only in one of the two bipartite networks -plant-herbivores-making it impossible to quantify Robustness in the other bipartite network -herbivore-parasite-We find that, in all cases, considering the merged network did not increased the robustness to plant extinction, a result in line with previous findings in multiple-mutualistic ecological networks [8].Rather, the robustness of the whole community (R) was in between that of the two bipartite networks composing the tripartite network (Fig KD ).

Using more advances extinction algorithms
Through the main text, we employ the classical co-extinction algorithm [10] to measure robustness to plant losses.Although this approach lacks realism (since there are no ecological dynamics behind it, nor any information about species preferences), it has proven useful in understanding the threat that biodiversity loss poses to ecosystem services and functioning [12,10,7].
To verify that the classical approach provides a lower bound on the damage caused on the ecological community (i.e. is the most conservative approach), we implemented the stochastic version of the co-extinction algorithm in mutualistic-antagonistic networks [19].In this stochastic version, a species i will undergo a secondary extinction following the extinction of plant j with a probability P ij = R i .dij , where d ij is the dependency of species i on j (interaction weight), and R i represents the intrinsic demographic dependence of species i on mutualism (for example, it may reflect the intrinsic dependence of a pollinator on nectar or on other resources, so that a lower value of R i would mean that that animal species does not rely completely on mutualism to survive, or that a plant can also rely on self pollination to survive).We simulated extinction cascades using this algorithm in the following manner: We erased plants following the extinction sequence.At each plant extinction, their animal partners had a probability P i j of undergo extinction.In our case, to simulate obligate mutualism we fixed R i = 1 for all animal species, but we continued using a bottom up control scenario by fixing R 0 for plants (i.e.plants could only go extinct in their turn on the extinction sequence, even if they had no mutualistic partners).
Below we show the comparison between our results in the main text (using the classical approach) and those obtained with the stochastic model on weighted networks, when available.When weights were not provided, we considered all equal weights of 1, which stills allows to implement the stochastic version of the co-extinction model.Fig L shows that the values of robustness (left) and Interdependence (right) are similar independently of the algorithm used, in particular for interdependence they are almost equivalent.Note that in the left figure, the values in the horizontal axis (robustness) are always higher than their stochastic weighted counterparts in the vertical axis, which corroborates that the classical approach at quantifying robustness is the most conservative, providing a clear lower bound that does not depend on any further parameter.
Fig M shows that the result of robustness of the merged network being a composition of the robustness of larger and smallest interaction layers also holds also in the case of using the stochastic weighted version.
In light of these results, the classical approach provides a lower bound on the damage caused to the ecological community.Given the fact that is also a more parsimonious model, without any parameter, we decided to use the classical approach in the main text.

Plant importance
The histogram of importance of plants indicates that the ranking is particularly sensitive in the case of the more important plants, both in the particular case that we use as an example in the main text (

Null models and Z-scores
To assess the importance of network structure in determining a certain network feature, we compared measurements of that feature performed on empirical networks with measurements performed on randomized versions of those networks keeping some properties fixed.We use four different null-models (see Fig 5 in main text), going from the least restrictive to the most restrictive (left to right in the figure), they are: • "1", or "NL": the number of nodes in each of the three species sets (N a , N b , and N l ), as well as the number of links in each interaction layer (L α , L β ) are kept constant.The randomization is performed in both interaction layers (α and β) separately, by linking L α (L β ) pairs of nodes from the two different sets involved in the interaction: a and l (b and l), ensuring that each of the nodes in the original network receives at least one link.Networks generated in the "constant NL" ensemble have the same size and connectance as the empirical networks, but erase any other structure present in the original interactions.In particular, the degree distribution of each species set is lost as well as the degree-degree correlations between nodes of species sets.When a metric is similar in the empirical networks and in the constant NL ensemble, one knows that the number of species and the connectance are the relevant features determining that metric.• "2", or "NL2": the number of nodes in each species set (N a ,N b , N l ) is kept constant, as well as the number of interactions in each interaction layer (L α ,L β ), and the degree of all nodes except those in the shared set.This null model allows the study of the effect of the degree distribution of the shared set species has in a given metric.When a metric is similar in the empirical network and in this random ensemble, one knows that the size of the sets and connectance are the the relevant features determining that metric (i.e. the degree distribution of the shared set species does not play a relevant role).• "3", or "K2": the number of nodes in each species set (N a ,N b , N l ) is kept constant, as well as the degree distribution of the species in both interaction layers.However, the nodes in the shared set do not keep their total degree (i.e.we beak the correlation between the degree of the shared set nodes in each interaction layer).This null model allows the study of the effect the correlation between the degree of connector nodes in each interaction layer has in a given metric.When a metric is similar in the empirical network and in this ensemble one knows that the correlation in the degrees of the connector species in different interaction layers is not relevant in determining that given metric.• "4", or "K": the number of nodes in each species set (N a ,N b , N l ), the total number of links in each interaction layer (L α ,L β ), and each node's degree (the number of links of a given node) are kept constant.However, the identity of neighbours are reshuffled within a layer.The reshuffling of the empirical networks is performed using the 'curveball algorithm', that ensures an unbiased sampling of random matrices [20].Networks generated in the "constant K" ensemble have exactly the same degree distribution as the empirical networks, but erase any other structure present in the original interactions, in particular the degree-degree correlation.When a metric is significantly different in the empirical network and in this null model, one knows that it is due to other structural features than degree distribution.When a metric is similar in the empirical networks and in the constant K ensemble, one knows that size and degree distribution are the relevant features determining that metric.
To compare the value of a given metric in an empirical network with that obtained in the random ensemble, we measure the Z-Score of that metric.The Z-Score quantifies the number of standard deviations by which the value of a raw score (i.e., the value measured in the empirical network) is above or below the mean value of what is being measured in the random ensemble, and is defined as: where x is the average x measured in the random ensemble and σ x its standard deviation.Positive/negative values mean that the measure in the empirical network is higher/lower than in the random-ensemble.The relevant Z score values when using a 95% confidence level are -1.96 and +1.96.The p-value associated with a 95% confidence level is 0.05.If the Z-Score is between -1.96 and +1.96, the p-value will be larger than 0.05, and one cannot reject the null hypothesis; Conversely, for a measure to be outside the 98% confidence interval (p<0.02) the absolute value of the Z-score has to go up to 2.33.
In the case where one wants to compare two distributions, for example when comparing the Robustness in the empirical network and the random ensemble (note that the Robustness of the empirical networks is a distribution of values rather than a number) we use the z-test.In this case, we compute the Z-score as: where x1 is the mean value of sample one (empirical network), x2 is the mean value of sample two (randomized ensemble), σ x1 is the standard deviation of sample one divided by the square root of the number of data points, and σ x2 is the standard deviation of sample two divided by the square root of the number of data points. .

D
Fig B:Boxplots of the basic structural features in the empirical tripartite networks of our data-set grouped (and color coded) by the sign of the interactions involved.A: Degree Heterogeneity (σ k / < k >) B: Degreedegree correlations (r).C: Z-score of the degree heterogeneity (σ k / < k >) in null model "1".D: Z-score of the degree-degrees correlations (r) in null model "4" in the tripartite networks.In B and D horizontal grey lines mark the limits of the confidence interval of 1.96 (95%) and 2.33 (98%).Differences between groups are measured by independent t-test (*** p < 0.001, ** p < 0.01, 'ns' not significant).

D
Fig C: A: Degree heterogeneity of the empirical tripartite ecological networks (σ k / < k >) at the species set level B: Degreedegree correlations (r) at the interaction layer scale.C: Boxplot of the Z-score of the degree heterogeneity (σ k / < k >) by species set in null model 1.D: Boxplot of the Z-score of the degree-degree correlations by interaction layer in null model 4.Horizontal grey lines mark the limits of the confidence interval of 1.96 (95%) and 2.33 (98%) for rejecting the null hypothesis.
Fig D:Cumulative degree distribution in a multipartite 'mutualistic-antagonistic' network of pollination-herbivory compared with its randomization in null model 1.Figures show the degree distribution for the three different sets of species, from left to right: insect herbivores, insect pollinators, plants and all nodes combined (i.e. the merged network).The continuous black line is the empirical cumulative degree distribution, with the best power law distribution fit in a dashed blue line.The continuous red line is the cumulative degree distribution in null model "1" ensemble (100 randomizations), and the dashed red line the best exponential fit to that distribution.The insets show the degree distribution of the empirical network (in grey) and that of 100 randomizations (in red).
Fig D:Cumulative degree distribution in a multipartite 'mutualistic-antagonistic' network of pollination-herbivory compared with its randomization in null model 1.Figures show the degree distribution for the three different sets of species, from left to right: insect herbivores, insect pollinators, plants and all nodes combined (i.e. the merged network).The continuous black line is the empirical cumulative degree distribution, with the best power law distribution fit in a dashed blue line.The continuous red line is the cumulative degree distribution in null model "1" ensemble (100 randomizations), and the dashed red line the best exponential fit to that distribution.The insets show the degree distribution of the empirical network (in grey) and that of 100 randomizations (in red).

E
Fig E: A: Interdependence (I) in the tripartite networks in our data set when plants are randomly driven to extinction, grouped and color coded by the signs of their interaction layers.B to E: Z-score of interdependence (I) in the different null models.The horizontal lines represent the Z-score values associated with a 95% and 98% confidence interval, z=1.96 and z=2.33 respectively.

F
Fig F: Correlation of interdependence (I) vs the structural features studied: A) degree heterogeneity (σ k / < k >), B) degree heterogeneity of the shared set species (σ k / < k > LS ), C) degree-degree correlations (r), D) proportion of connector nodes inside the shared set (C), E) Proportion of shared set hubs that are connectors (H C ), and F) average participation ratio of the connector nodes (P R C ) on the tripartite networks in our data set.The colors of the points represent the different types of tripartite networks according to the sign of their interaction layers (see legend).The values in the lower right side are the Pearson correlation coefficients considering only MM and MA networks (blue), only AA networks (violet).

F
Fig H: Effect of degree heterogeneity (σ k / < k >) and degree-degree correlations (r) on the Robustness (R) of the tripartite networks in our data set.The colors of the point represent the different types of tripartite networks according to the sign of their interaction layers (violet AA, green MA, and blue MM).The values in the upper right side are the Pearson correlation coefficients considering only MM and MA networks (blue), and only AA networks (violet).The lines represents the best linear regression for each of the correlations.

I
Fig I: Effect of structural features of the shared set on the Robustness (R) of the tripartite networks in our data set.The color of the points represents the different types of tripartite networks according to the sign of their interaction layers (violet AA, green MA, and blue MM).The values on the side are the Pearson correlation coefficients considering only MM and MA networks (blue), and only AA networks (violet).The lines represent the best linear regression for each of the correlations.

E
Fig J: A: Robustness (R) of the tripartite networks in our data set when plants are randomly driven to extinction, grouped and color coded by the sign of the interaction layers.B to E: Z-score of R in the different null models.The horizontal lines represent the Z-score values associated with a 95% and 98% confidence interval, z=1.96 and z=2.33 respectively.

Fig K :
Fig K:Upper panel: correlation of the Robustness of the whole merged network (R) with A) the Robustness of the larger interaction layer (R L ), B) the Robustness of the smaller interaction layer (R S ), and C) with the best estimation of Robustness (R (est) ) as a composition of the Robustness of the two interaction layers.Lower panel: Comparison of the Robustness of the tripartite networks (in blue) with that of the two interaction layers composing them (R L in green and R S in red) and with the best estimated Robustness as a composition of the Robustness of the two interaction layers (formula shown).

FigC
Fig L: Comparative of results between the stochastic co-extinction algorithm in weighted networks (when available) with the results of the classic co-extinction algorithm used in the main text, for Robustness (left) and Interdependence (right).
Fig N: Upper row: Distribution of plant importance values in the network used as example in Fig. 4. The histograms depict the amount of plants with a given importance in the whole community (left), with a given importance for pollinators (center) and for herbivores (rigth).Lower row: Ranking of importance (based on importance) vs importance.Inset values show the spearman rank correlation coefficient (ρ) and pearson correlation coefficient (r).
Fig O: Distribution of plant importance values considering all the communities together.The histograms depict the amount of plants with a given importance in the whole community (left), with a given importance for animals connected through positive interactions -pollination or seed dispersion-(center) and for animals connected through negative interactions -herbivory-(rigth).

Table C :
Multiple linear regression of interdependence (I) vs structural features