DOI: 10.1038/ncomms12573

OPEN

No complexity–stability relationship in empirical ecosystems Claire Jacquet1,2,3, Charlotte Moritz4,5, Lyne Morissette6, Pierre Legagneux1,2, Franc¸ois Massol7, Philippe Archambault4,8 & Dominique Gravel1,2,9

Understanding the mechanisms responsible for stability and persistence of ecosystems is one of the greatest challenges in ecology. Robert May showed that, contrary to intuition, complex randomly built ecosystems are less likely to be stable than simpler ones. Few attempts have been tried to test May’s prediction empirically, and we still ignore what is the actual complexity–stability relationship in natural ecosystems. Here we perform a stability analysis of 116 quantitative food webs sampled worldwide. We ﬁnd that classic descriptors of complexity (species richness, connectance and interaction strength) are not associated with stability in empirical food webs. Further analysis reveals that a correlation between the effects of predators on prey and those of prey on predators, combined with a high frequency of weak interactions, stabilize food web dynamics relative to the random expectation. We conclude that empirical food webs have several non-random properties contributing to the absence of a complexity–stability relationship.

1 De ´partement de Biologie, Universite´ du Que´bec a` Rimouski, 300 Alle´e des Ursulines, Quebec, Canada G5L 3A1. 2 Quebec Center for Biodiversity Science, Montre´al, Quebec, Canada H3A 1B1. 3 UMR MARBEC, Universite´ de Montpellier, Place Euge`ne Bataillon, F-34095 Montpellier cedex 05, France. 4 Institut des Sciences de la Mer de Rimouski, Universite´ du Que´bec a` Rimouski, 310 Alle´e des Ursulines, Quebec, Canada G5L 3A1. 5 Centre de Recherches Insulaires et Observatoire de l’Environnement, EPHE, PSL Research University, UPVD, CNRS, USR 3278 CRIOBE, F-98729 Moorea, French Polynesia. 6 M-Expertise Marine, 10 rue Luce-Drapeau, Sainte-Luce Quebec, Canada G0K1P0. 7 Unite´ Evolution, Ecologie & Pale´ontologie (EEP), SPICI group, CNRS UMR 8198, Universite´ Lille 1, Baˆtiment SN2, F-59655 Villeneuve d’Ascq cedex, France. 8 Que´bec-Oce´an, De´partement de biologie, Universite´ Laval, Pavillon Alexandre Vachon, 1045, avenue de la me´decine, Quebec, QC, Canada G1V 0A6. 9 De´partement de biologie, Faculte´ des Sciences, Universite´ de Sherbrooke, 2500 Boulevard Universite´, Sherbrooke, Quebec, Canada J1K 2R1. Correspondence and requests for materials should be addressed to C.J. (email: [email protected]).

NATURE COMMUNICATIONS | 7:12573 | DOI: 10.1038/ncomms12573 | www.nature.com/naturecommunications

1

ARTICLE

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12573

he complexity–stability debate1, initiated more than 40 years ago, stems from two apparently conﬂicting observations. On the one hand, complex ecosystems are ubiquitous in nature, as illustrated by diverse tropical forests, coral reefs or intertidal communities. These observations have inspired ecologists to hypothesize that complexity could stabilize ecosystems2,3. On the other hand, theory states that complex random systems are less likely to recover from small perturbations than simpler ones4–6. This prediction was put forth by Robert May7, who studied the relationship between complexity and stability in prandom ecosystems. Ecosystem ﬃﬃﬃﬃﬃﬃ complexity was deﬁned as s SC where S is species richness, C is connectance (the probability that any two species will interact with each other) and s is the s.d. of interaction strength. May7 predicted that a system could be stable only if the criterion pﬃﬃﬃﬃﬃﬃ s SC od was satisﬁed, where d expresses the magnitude of intraspeciﬁc competition. May7 analysed the local stability of randomly generated community matrices. A community matrix is obtained from the linearization around a feasible equilibrium of a system of equations describing the dynamics of the community. The entries of a community matrix quantify the impact of a change in abundance of one species on the dynamics of another species. The real part of the dominant eigenvalue of the community matrix indicates the rate at which a system returns to equilibrium (if negative) or moves away from it (if positive) after small perturbations. It does not guarantee stability following large perturbations (global stability), or that the perturbation will not ﬁrst amplify before vanishing (reactivity)8,9. The stability of a random community matrix can be predicted thanks to the generalization of the circular law10. This theory states that the distribution of the eigenvalues of a S S matrix, whose coefﬁcients are independently sampled from a distribution of mean 0 and variance 1, converges to the uniform distribution in the unit circle in the complex plane, as S-N. The centre of the circle d corresponds to the mean of intraspeciﬁc interaction terms d40 , provided that the variance in intraspeciﬁc Rﬃ is related to interaction terms is not too large11. The radius pﬃﬃﬃﬃﬃ interspeciﬁc interactions and is equal to s SC in random ecosystems, that is May’s complexity measure. Thus, local stability is determined by the combination of two components; one can increase the stability of a system by (i) moving the centre

T

a

b 0

C

+

C=

of the circle to more negative values along the real axis by increasing intraspeciﬁc competition or (ii) decreasing the radius of the circle by reducing the complexity of the system (Fig. 1c). Tang et al12. proposed that another quantity critically affects the stability of more realistic ecosystems, such as predator–prey communities, namely the correlation between coefﬁcients across the diagonal of the community matrix r. They subsequently found that the stability pcriterion for large and random ﬃﬃﬃﬃﬃﬃ community matrices is s SC ð1 þ rÞ Eod where E is the mean of the elements of the community matrix (including zeros). In other words, the correlation between pairs of interactions decreases stability if positive (r40) but increases stability if negative (ro0) with respect to May’s case12. Here we attempt to solve the complexity–stability paradox with a local stability analysis of 116 quantitative food webs sampled worldwide from marine, freshwater and terrestrial habitats. This is the largest data set ever used to test May’s prediction empirically. The complexity–stability relationship has been previously studied with direct observations of energy ﬂows between species, but on a small number of food webs (from one to seven)13–15. Recently, Neutel and Thorne16 reported an absence of complexity–stability relationship in 21 soil food webs, while James et al17. found a weak positive relationship based on 21 food webs from terrestrial and marine habitats. These studies, however, used heterogeneous methodologies, shared several networks, and in several cases, interaction strengths were derived from assumptions rather than from direct observations18,19. The studied food webs were all compiled on the same standard methodology to satisfy the Ecopath modelling framework20. Ecopath is a trophic model, the most widely used tool for ecosystem-based ﬁsheries management, and has also been used to characterize unexploited terrestrial ecosystems21. It relies on a system of linear equations established with the aim of balancing the inﬂows and the outﬂows of each compartment20,22. A large amount of information is included in Ecopath models, such as diet composition, biomass, production and consumption rates of each species, providing an accurate representation of feeding interactions within food webs. Ecopath models provide a unique opportunity to build realistic community matrices with empirical data derived from a standardized protocol. The level of resolution of marine Ecopath models is, however, heterogeneous through food web compartments, with detailed compartments for

c

DCcr × (Q /B)c × erc × Bc Br

–DCcr × (Q /B)c

Imaginary

0

erc × rc ×Br ×Bc

Bc × (P/B)c × DCcr

max

–

–d1

Inefficiency rc ×Br ×Bc

Bc × (Q /B)c × DCcr

C=

+c2,1 +c3,1

–c1,2

–d2

0

+c4,2

–c1,3

0

–d3

+c4,3

0

–c2,4

–c3,4

–d4

R

–d

0

Real

R

Figure 1 | Method summary. (a) Equivalence between Ecopath and Lotka–Volterra models: simpliﬁed diagram of trophic ﬂows between one consumer C and one resource R parameterized with Ecopath model (in blue) and Lotka–Volterra model (in green). B is biomass (t km 2), (P/B)c and (Q/B)c are the production/biomass and consumption/biomass ratios of C respectively (per year), DCcr is the proportion of resource R in the diet of consumer C, erc P=BÞc expresses the efﬁciency of a consumer to convert resource energy into biomass with erc ¼ ððQ=B Þ . (b) Community matrix construction: derivation of c

community matrix elements for the simpliﬁed food web presented in diagram A, and an example of community matrix structure observed in real food webs. (c) Measure of stability: the eigenvalues of a large community matrix are uniformly distributed on a circle on the complex plane (axes cross at the origin). On the real axis, the dominant eigenvalue Re(lmax) ¼ R d, where the centre of the circle d is equal to the mean of intraspeciﬁc interaction terms d40 , pﬃﬃﬃﬃﬃﬃ and the radius R is related to interspeciﬁc interaction terms (that is, off-diagonal elements of the community matrix) and is equal to s SC in random pﬃﬃﬃﬃﬃﬃ matrices. For predator-prey communities, the eigenvalues are uniformly distributed on an ellipse with a horizontal half axis R ¼ ð1 þ rÞs SC (ref. 11). 2

NATURE COMMUNICATIONS | 7:12573 | DOI: 10.1038/ncomms12573 | www.nature.com/naturecommunications

ARTICLE

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12573

collected ﬁshes and more aggregated compartments for plankton and invertebrates. We translated parameters of the Ecopath models into interaction coefﬁcients of the Lotka–Volterra interaction model following the same approach as De Ruiter et al13. Interaction coefﬁcients from all pairwise interactions of a food web make the interaction matrix A ¼ [aij]. Because of the equilibrium assumption of Ecopath models, a community matrix C can be constructed for each food web by multiplying the interaction matrix A with species biomass (Fig. 1b). We measured food web stability using the real part of the dominant eigenvalue of the community matrix C to be directly comparable to May’s approach. The diagonal elements of the community matrices were set to 0 to focus on the effect of interspeciﬁc interactions on stability. Note that Re(lmax) will be positive, since R40 (Fig. 1c). This method is comparable to other studies that calculated stability by assessing the level of intraspeciﬁc interaction needed for all eigenvalues in a community matrix to have negative real parts14,16,19. We show that complexity is not related to stability in empirical ecosystems. We ﬁnd that the intrinsic energetic organization of food webs creates a high frequency of weak interactions and a correlation between pairs of interactions. These non-random properties are highly stabilizing and contribute to the absence of a complexity–stability relationship. Results Complexity–stability relationship in empirical ecosystems. We ﬁrst investigated the relationship between stability and classic Marine

descriptors of ecosystem complexity23, that is, species richness S, connectance C and s.d. of interaction strengths s. We observed no relationship between food web stability and species richness, neither with connectance nor with s.d. of interaction strength (Fig. 2). Further analyses revealed that this result was robust to the variability of sampling intensity among the 116 food webs and to uncertainty related to Ecopath parameter estimates (Methods section, Supplementary Figs 1,2 and 3). We neither found signiﬁcant complexity–stability relationship using the stability criterion derived by Tang et al12. that integrates correlation between pairs of interactions and mean of interaction strengths (Supplementary Fig. 4). The absence of a complexity–stability relationship in empirical food webs demonstrates that the random matrices studied by May7 to derive stability criteria deviate signiﬁcantly from empirical systems. As May7 stated in the re-edition of his book, his theory provides the baseline against which we should compare empirical systems and ﬁnd the non-random features stabilizing them. We therefore investigated further the mechanisms preventing the negative relationship between complexity and stability to occur. Correlationpﬃﬃﬃﬃﬃ between complexity parameters. May’s stability ﬃ criterion s SCod indirectly predicts that for complex systems to persist, interaction strength should be weaker in species-rich and highly connected systems7,24. In other words, complex ecosystems could persist in nature provided that S, C and s are not independent. Inﬃ the same way, the inequality derivedpﬃﬃﬃﬃﬃ byﬃ pﬃﬃﬃﬃﬃ predicts that r and s SC Tang et al.12, s SC ð1 þ rÞ Eod, should be correlated in feasible ecosystems. We therefore

Freshwater

a

Terrestrial

b

Stability Re (λmax)

6 5 4 3 2 1 0 10

20 30 40 Species richness S

c

50

0.2

0.3

0.4 0.5 0.6 Connectance C

0.7

0.8

d

Stability Re (λmax)

6 5 4 3 2 1 0 0

20

40 60 80 100 Standard deviation of IS

0

50

100 150 Complexity SC

200

Figure 2 | Food web stability related to complexity parameters in 116 food webs. (a) Number of species S (linear regression: P ¼ 0.97, R2o10 5), (b) Connectance C ¼ (L/S2) where L is the number of links (P ¼ 0.98, R2o10 6), (c) Standard deviation of interaction strengths s (P ¼ 0.1, R2 ¼ 0.02), pﬃﬃﬃﬃﬃﬃ (d) May’s complexity measure s SC (P ¼ 0.02, R2 ¼ 0.04). Stability is measured as Re(lmax) for marine (blue), freshwater (green) and terrestrial ecosystems (orange). Food webs with eigenvalues close to zero are the most stable. All quantities are dimensionless. NATURE COMMUNICATIONS | 7:12573 | DOI: 10.1038/ncomms12573 | www.nature.com/naturecommunications

3

ARTICLE

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12573

hypothesized that, contrary to randomly built ecosystems, parameters describing complexity are not independent in nature. We found that the s.d. of interaction strength s across the 116 food webs was negativelypcorrelated to the product of ﬃﬃﬃﬃﬃﬃ species richness and connectance SC (Fig. 3a) and contrary to expectations, pﬃﬃﬃﬃﬃﬃwe observed a slightly positive correlation between pﬃﬃﬃﬃﬃﬃ r and s SC (Fig. 3b). The correlation between s and SC decreased the overall complexity and conﬁrmed the existence of feasibility constraints on communities. However, we still observed higher values of s than predicted by May’s stability criterion and this observation did not explain the absence of complexity– stability relationship in empirical systems. Non-random properties of empirical community matrices. Random matrix theory supposes that interaction coefﬁcients are independent and identically distributed in the community matrix. Marine

a

Freshwater

Terrestrial

120 100

80 60 40 20 0 2.0

2.5

3.0

3.5

4.0

4.5

SC

Pairwise correlation

b

0.0 −0.2 −0.4 −0.6 −0.8 0

50

100

150

200

SC

Figure 3 | Correlation between complexity parameters in real food webs. (a) s is the s.d. of interaction strengths, S the number of species and C the pﬃﬃﬃﬃﬃﬃ connectance. The product SC was negatively correlated to s (Spearman’s rank correlation (Po10 13, r ¼ 0.64). (b) r¼ corrðcij ; cji Þ is the pﬃﬃﬃﬃﬃﬃ correlation between pairwise interactions. The product s SC was positively correlated to r (P ¼ 0.02, r ¼ 0.22).

However, many studies on the complexity-stability relationship suggest that real ecosystems have non-random structural properties promoting their stability despite their complexity23. We focused on four non-random properties observed in our empirical community matrices, and then investigated their contribution to stability with randomization tests. (i) Pyramidal structure of interaction strength13,17,18: we found that interaction strength was related to trophic level, the occurrence of strong interactions being more likely at low trophic levels. Species biomass distribution affected the mean and variance of row i, since cji ¼ aij Bj. Consequently, rows had different means and variance, a feature we call row structure. We hypothesized that food webs without this row structure are less stable than real food webs. (ii) Interaction strength topology14,15,18: trophic structure determines the position and the direction of interaction strength (that is, ‘who eats whom’), and creates a non-random topology of interaction strengths. We hypothesized that food webs with a random topological structure are less stable than real food webs. (iii) Correlation between pairs of predator-prey interactions12,25,26: we found a correlation between pairs of interaction strengths cij and cji in community matrices, since cji ¼ ( cij eij Bj)/Bi (Fig. 1b). We therefore hypothesized that food webs with an empirical topological structure, but with a null correlation between pairs of interaction strengths, are less stable than real food webs. (iv) Interaction strength frequency distribution: in agreement with previous studies13,27–29, we observed a leptokurtic distribution of interaction strengths (high proportion of weak interactions). Consequently, we hypothesized that food webs with a highly peaked and long tailed distribution of interaction strengths are more stable than ﬂatter distributions, such as the normal distribution. Randomization tests. We performed eight randomization tests to remove one or several properties of natural food webs and computed stability of the permuted community matrices (called H1–H8 at Table 1, see Methods section for details). We used this method to determine whether these properties had a signiﬁcant effect on the distribution of eigenvalues across the 116 food webs, and their impact on the complexity–stability relationship. Randomization tests removed some non-random features of empirically built community matrices, generating matrices more similar to the ones expected under the random matrix theory, in which elements are drawn from a standardized distribution. The distribution of eigenvalues of the permuted food webs was compared to stability of the original food webs. We found that each of the four structural properties enhanced food web stability (Fig. 4a). The removal of the empirical distribution of interaction

Table 1 | Properties conserved by each randomization test (indicated by a ü). Hypothesis H1 H2 H3 H4 H5 H6 H7 H8

Row structure

Topology ü ü ü ü

Pairwise correlation ü ü ü ü

Frequency distribution ü ü ü ü

The column ‘row structure’ speciﬁes whether the pyramidal structure of interaction strength is conserved or not in randomization tests H1–H8. Similarly, ‘topology’ corresponds to interaction strength topology (‘who eats whom’), ‘pairwise correlation’ corresponds to the correlation between pairs of predator–prey interactions and ‘frequency distribution’ corresponds to the leptokurtic distribution of interaction strengths (high proportion of weak interactions).

4

NATURE COMMUNICATIONS | 7:12573 | DOI: 10.1038/ncomms12573 | www.nature.com/naturecommunications

ARTICLE

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12573

strengths (with many weak interactions, H4) had the strongest impact on stability, followed by the removal of correlation between pairs of interactions (H3). Note that in all the randomization tests, the pyramidal structure of interaction strength was removed. Stability decreased when only this property was removed, keeping empirical topology, pairwise correlation and interaction strength distribution (H1). The randomization of interaction strength topology (H2) was also destabilizing, but to a very lesser extent compared with others non-random properties (Fig. 4a). Randomization tests resulted in some cases in a negative complexity–stability relationship, although weaker than one should expect from the random matrix theory. Even if randomized matrices conserved the same S, C and s2 as original ones (and thus their correlation, presented in Fig. 3a), we found a negative complexity–stability relationship when we normalized interaction strength distribution (H4, linear regression: Po10 16, R2 ¼ 0.64) and removed correlation between pairs of interactions (H3, linear regression: P ¼ 10 7, R2 ¼ 0.2, Fig. 4c). The removal of the pyramidal structure of interaction strengths and the topology found in empirical ecosystems did not affect the

relationship between complexity and stability (linear regressions, H1: P ¼ 0.38, R2 ¼ 0.002, H2: P ¼ 0.2, R2 ¼ 0.006, Fig. 4c). All food web properties contributed to stability, but clearly, the leptokurtic distribution of interaction strength had the strongest impact on the complexity–stability relationship. We found a signiﬁcant negative relationship between stability and complexity when we removed this property (H4, Fig. 4c). Conversely, when we only kept the empirical distribution of interaction strengths (H5), the slope of the complexity–stability relationship was signiﬁcantly ﬂatter than in the random case (H8). Topology of interaction strengths (H7) or pairwise correlation (H6) alone did not signiﬁcantly affect the complexity–stability relationship (Fig. 4d). We conclude that May’s stability criterion does not apply to empirical ecosystems because of their structure, which has several stabilizing non-random properties. First, the high frequency of weak interactions balanced the destabilizing effect of complexity (H4). Interestingly, we observed a strong positive correlation between the kurtosis k (index of the peakedness of the interaction strength distribution) and species richness in real food webs (Supplementary Fig. 5). Thus the probability of having many

H5 Interaction strength distribution only H6 Pairwise correlation only H7 Topology only H8 Random food web

Real food web H1 Random row structure H2 Random topology H3 Random pairwise correlation H4 Random interaction strength distribution

b

Density

a 4

4

3

3

2

2

1

1

0

0 0.0

0.5 1.0 1.5 Stability log10 (Re(λmax))

2.0

c

0.0

0.5 1.0 1.5 Stability log10 (Re(λmax))

2.0

d 50

Stability Re(λmax)

80 60

30

40 20

10 0

0 0

50

100 150 Complexity SC

200

0

50

100 150 Complexity SC

200

Figure 4 | Complexity–stability relationship in empirical and permuted food webs. Frequency distributions of eigenvalues of real and permuted food webs : (a) permutation tests H1 to H4, (b) permutation tests H5 to H8. Eigenvalues are on a logarithmic scale and dimensionless. Permutation tests were carried out 1,000 times for each food web. Eigenvalue distributions were smoothed using a kernel density estimate of 0.28. (c) Stability of real and pﬃﬃﬃﬃﬃﬃ permuted food webs related to complexity (permutation tests H1–H4). Stability is measured as Re(lmax) and s SC corresponds to complexity. Statistics of the linear regression between complexity and stability: real food webs (slope ¼ 0.005, P ¼ 0.02, R2 ¼ 0.04), H1: random row structure (slope ¼ 0.003, P ¼ 0.38, R2 ¼ 0.002), H2: random topology (slope ¼ 0.006, P ¼ 0.2, R2 ¼ 0.006), H3: random pairwise correlation (slope ¼ 0.06, P ¼ 10 7, R2 ¼ 0.2), H4: random interaction strength distribution (slope ¼ 0.24, Po10 16, R2 ¼ 0.65). (d) Stability of permuted food webs related to complexity (permutation tests H5–H8). Statistics of the linear regression between complexity and stability: H5: empirical distribution of interaction strengths only (slope ¼ 0.06, P ¼ 10 9, R2 ¼ 0.25). H6: pairwise correlation only (slope ¼ 0.3, Po10 16, R2 ¼ 0.63). H7: topology only (slope ¼ 0.32, Po10 16, R2 ¼ 0.6). H8: random food webs (slope ¼ 0.33, Po10 16, R2 ¼ 0.66). NATURE COMMUNICATIONS | 7:12573 | DOI: 10.1038/ncomms12573 | www.nature.com/naturecommunications

5

ARTICLE

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12573

weak interactions increased with species richness. The negative correlation between pairs of interaction strengths cij and cji is also a strong stabilizing property of empirical community matrices (H3). Finally, the non-random topology of interaction strengths (H2) was also stabilizing, as suggested by previous studies13,14,16,18. Discussion The relevance of local stability analysis to study real ecosystems may be questioned. More general and realistic deﬁnitions of stability have been introduced during the complexity–stability debate, such as persistence, variability, resilience or resistance30. Indeed, local stability analysis only tests the impact of small perturbations on ecological dynamics, and may not apply to large and/or cumulative perturbations typical of most empirical studies. It neither considers the covariance among species and thus the stability of the aggregated properties of the community31. However, it allows the use of analytically tractable community matrices, and thus the investigation of May’s complexity–stability relationship on real ecosystems. Our study yields new insight on the complexity–stability debate. Random matrix theory cannot predict the stability of real ecosystems because interaction strengths are not independent and identically distributed in empirically derived community matrices. Trophic structure creates a negative correlation between pairs of interactions and a non-random distribution of interaction strengths, with many weak interactions and few strong ones at the bottom of the food webs. The likely explanation for the strong effect of the leptokurtic distribution of interaction strengths is the size of the community matrices we investigated. Random matrix theory is performed in the limit of inﬁnitely large matrices and all distributions are expected to converge in systems of several hundreds of species11. The community matrices we investigated had between 6 and 54 species. A detailed investigation of some community matrices revealed that small modules (two to ﬁve species) were often responsible for extreme eigenvalues. These modules could drive strong negative or positive feedbacks16 and thus dominate the dynamics of the entire system. Random matrix theory could provide a sufﬁcient approximation for large ecosystems, but needs to be reﬁned for smaller and realistic food webs such as the ones we investigated. The study of small community matrices might require a different theoretical framework. For instance, Neutel and Thorne16 showed that the stability of a dynamical system could be predicted from the analysis of feedback loops. However, this approach requires knowledge of all of the elements of the community matrix and does not provide a statement about the expected relationship between S, C, s and the occurrence of feedback loops. Such a theory would be needed to make quantitative predictions about the stability of a system with estimates of only few state variables. Our food web dataset provided a great opportunity to study the effect of interspeciﬁc interactions on the relationship between complexity and stability and to demonstrate the existence of a negative correlation between S, C and s in empirical ecosystems. We had, however, no information about the strength of intraspeciﬁc interactions, which is a strong stabilizing mechanism. Our analysis thus focused on the radius of the distribution of eigenvalues in the complex plane, ignoring the location of the centre (Fig. 1c). It is possible that the absence of relationship between complexity and stability results from a positive correlation between thepﬃﬃﬃﬃﬃ strength of intraspeciﬁc ﬃ interactions d and complexity s SC . Here we hypothesized that the food webs we studied were mainly top-down controlled, and that the strength of intraspeciﬁc interactions was negligible in comparison to interspeciﬁc interactions. Nonetheless, we 6

evaluated the sensitivity of our ﬁndings to the addition of intraspeciﬁc interaction terms proportional to species equilibrium biomass, since cii ¼ aii Bi. In agreement with random matrix theory and previous studies11,14,19, the addition of intraspeciﬁc interactions was stabilizing, but had no effect on the correlation between complexity and stability (Methods section, Supplementary Fig. 6). Our results emphasize that further empirical investigations should better consider the relationship between ecosystem complexity and density dependence. The analysis of empirically derived community matrices, combined with the observation of a complexity–stability relationship when their non-random structural properties were removed, demonstrates that the properties captured by Ecopath models contribute to the stability of complex food webs. Further empirical investigations are necessary to better approach real ecosystems and to study the stabilizing effect of the properties ignored or poorly described in Ecopath models, such as species age structure, energy ﬂows from detrital pool or external inputs. We showed that complexity is not related to stability in empirical ecosystems, a question that has stimulated ecological research for four decades. We found that the intrinsic energetic organization of food webs is highly stabilizing and allows complex ecosystems to recover from perturbations. Coexistence also constrains the feasibility of ecosystems, imposing a non-random structure of interactions and a correlation between S, C and s that decreases the overall complexity24. The non-random structure of food webs occurs from the successive addition of consumers having an increasingly large diet, which causes a growing frequency of weak interactions. The complexity–stability debate has contributed to the development of productive research that have pointed out the key role of the structural properties of real ecosystems. Methods Ecopath modelling framework. We compiled 116 Ecopath food web models from published studies (Supplementary Table 1). Ecopath provides a quantitative overview of how species interact in a food web. Species sharing the same prey and predators and having similar physiological characteristics are aggregated in trophic species. The dynamics of each species i is described by the difference between biomass production and biomass losses due to harvesting, predation or other unspeciﬁed sources. It can be expressed as: i Xh dBi ¼ Bi ðP=BÞi Yi Bj ðQ=BÞj DCji M0i Bi ð1Þ dt j where Bi (t km 2) and (P/B)i (per year) are biomass and production/biomass ratio of species i, respectively, Yi (t km 2 per year) corresponds to ﬁshery yields, (Q/B)j (per year) is consumption/biomass ratio of predator j and DCji is the proportion of species i in the diet of predator j. Other mortality sources, M0i (per year), can be expressed as (1 EEi) (P/B)i, where EEi is called the ecotrophic efﬁciency of i, corresponding to the fraction of the production that is used in the food web. The Ecopath model assumes mass-balance, meaning that all species biomass are at equilibrium (dBi/dt ¼ 0). Input parameters (that is, biomass, production/biomass and consumption/ biomass ratios, ﬁshery yields, and diet composition) can have different origins: ﬁeld sampling (for example, trawl survey), derived from similar Ecopath models, or known empirical relationships. Ecopath with Ecosim software includes routines that estimate missing parameters based on the mass-balance hypothesis and the generalized inverse method for a system of n linear equations and m unknowns (see Christensen et al.22, p. 12–15). In general, the biomasses, production/biomass and consumption/biomass ratios are entered for all groups to estimate ecotrophic efﬁciency, which is difﬁcult to measure in the ﬁeld32. The Ecoranger module, also included in Ecopath with Ecosim software, can be used to explore the effect of uncertainty in input data on estimated parameters. This module calculates probability distributions of output parameters based on the conﬁdence intervals of input parameters speciﬁed by the users32. Full details of the Ecopath modelling approach and the Ecopath with Ecosim software can be obtained from www.ecopath.org. Parameterization of Lotka–Volterra interaction coefﬁcients. We used the method from De Ruiter et al13. to derive the community matrices from Ecopath models (Fig. 1a): assuming direct dependence of feeding rates on predator

NATURE COMMUNICATIONS | 7:12573 | DOI: 10.1038/ncomms12573 | www.nature.com/naturecommunications

ARTICLE

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12573

population density, we calculated the per capita effect of predator j on the growth rate of prey i as aij ¼ ((Q/B)j DCji)/Bi. Effects of prey on their predator are deﬁned as predator growth resulting from this predation. Consequently, effect of the prey i on the predator j is related to effect of the predator on the prey according to: aji ¼ eij aij, where B is biomass, DCji is the proportion of species i in the diet of predator j, eij is the efﬁciency with which j converts food into biomass, from ðP=BÞ

feeding on i: eij ¼ ðQ=BÞj and (P/B)j and (Q/B)j are predator production/biomass j

and consumption/biomass ratios respectively. We obtained the following Lotka–Volterra interaction equation: XS dBi ¼ Bi bi þ aij Bj aii Bi ð2Þ j¼1 dt where bi is the intrinsic growth rate (that is, the intrinsic rate of increase for autotrophs, and natural mortality and losses for heterotrophs), Bi and Bj are, respectively, biomass of species i and j, interaction strength aij corresponds to the per capita effect of species j on the growth rate of species i and aii represents the per capita self limitation of species i. Assuming mass-balance, we obtain the following P expression for intrinsic growth rate: bi ¼ Sj¼1 aij Bj þ aii Bi : Correlation between pairs of interactions. Pairwise correlation was calculated Eðcij cji Þ Eðcij Þ2 where E(cij) is using the formula from Tang et al.12: r¼corr cij ; cji ¼ V the mean of the off-diagonal elements cij of the community matrix, their variance is V and E(cijcji) is the mean of the products of the pairs cijcji. Randomization tests. Reported dominant eigenvalues of randomized food webs corresponded to the mean of 1,000 replicates. All permutation tests conserved S, C and s. To randomize the pyramidal structure of interaction strengths (H1), we swapped pairs of predator-prey interactions (the pair cij/ þ cji was replaced by the pair ckl/ þ clk). This permutation only changed row structure (mean and variance) and did not change topology, frequency distribution of interaction strengths nor correlation between pairs of interactions. To randomize interaction strength topology (H2), we swapped the element of the community matrix cij with the element cji. This permutation only removed food web topology and did not change the frequency distribution of interaction strengths or pairwise correlation. To remove correlation between pairs of predator-prey interactions (H3), we permuted off-diagonal elements of the community matrix, keeping the topological structure and the frequency distribution of interaction strengths. Positive and negative terms were permuted separately to keep identical averages of positive and negative interactions. For randomization of interaction strength distribution (H4), we created a random community matrix in which off-diagonal elements were picked from a bivariate normal distribution N2(m, S) where the mean vector m is composed of the mean of positive (m þ ) and the mean of negative (m ) terms, and S is the covariance matrix between positive and negative terms of the original community matrix. Original pairs of positive/negative terms were replaced by positive/negative terms from the bivariate normal distribution. For large random community matrices, the correlation between pairwise interactions is expected to be identical to the original community matrix. Randomization test H5 only kept frequency distribution of interaction strengths. This test is a combination of permutation H2, that randomizes the topology of interaction strengths, and permutation H3, that removes pairwise correlation. Randomization test H6 only kept pairwise correlation; this test is a combination of permutation H2 and randomization H4, that creates a random community matrix in which off-diagonal elements are picked from a bivariate normal distribution. Randomization test H7 only kept the topology of interaction strengths, which is a combination of tests H3 and H4. Randomization test H8 created community matrices in which elements were identically and independently distributed, that is food webs with a random topology, a normal distribution of interaction strengths and no correlation between pairs of interactions. This test corresponds to test H2, that randomizes the topology of the community matrix, combined to a randomization that creates a community matrix in which positive and negative off-diagonal elements are picked from a normal distribution N(m þ ,s2 þ ) and N(m ,s2 ), where m þ and m are the mean and s2 þ and s2 the variance of positive/negative elements of the original community matrix. Parameter uncertainty. We investigated the impact of parameter uncertainty on our ﬁndings. In the section ‘Interspeciﬁc interaction terms of the community matrix’, we evaluated the sensitivity of our results to variability in interspeciﬁc interaction terms. The parameters used to build empirical community matrices come from Ecopath data and each of them bears some uncertainty22. Consequently, we tested whether the introduction of variability in input parameters could bias the complexity–stability relationship. In the section ‘Intraspeciﬁc interaction terms of the community matrix’, we evaluated the robustness of our results to the addition of density dependence. Because Ecopath models depict exclusively trophic interactions between species, we had no empirical information about the strength of intraspeciﬁc interactions and we decided not to model density dependence in the Lotka–Volterra model. Our method is comparable to other studies that calculated stability by assessing the level of intraspeciﬁc interaction needed for all eigenvalues in a community matrix to

have negative real parts (diagonal dominance)6,16,19. These studies assumed that all diagonal elements cii of the community matrix are the same. However, to obtain the community matrix, the interaction matrix A is multiplied by species biomass, which means that diagonal elements are non-constant: cii ¼ aii Bi. We therefore evaluated the robustness of our results to the addition of diagonal elements structured by species biomass. Finally, in the section ‘Food web resolution’, we assessed the impact of food web resolution level on the complexity–stability relationship. Ecopath model is mainly used for ecosystem-based ﬁsheries management and the level of resolution of several food webs is not homogeneous through all ecological compartments. Harvested ﬁshes are generally resolved at the species level, while species at the bottom of the food web, such as plankton and invertebrates, are highly aggregated. We therefore analysed the complexity–stability relationship on a subset of the best resolved Ecopath food webs. Overall, we found the same qualitative results than our main study. We conclude that our ﬁndings are robust to (i) input parameter uncertainty, (ii) addition of non-zero diagonal elements in community matrices and (iii) differences in food web resolution level. Interspeciﬁc interaction terms of the community matrix. We ran sensitivity analyses to determine how uncertainties in parameter estimates could affect the results of the study. For each input parameter, we tested if uncertainty biases (that is, overestimates or underestimates) food web stability and the variables determining complexity, and if our ﬁndings are qualitatively affected by this bias. We used the following parameters from Ecopath data to determine the interspeciﬁc terms of a community matrix: (i) biomass B, (ii) consumption/biomass ratio (Q/B), (iii) production/biomass ratio (P/B) and (iv) diet composition DC. Uncertainty in these parameters could inﬂuence our results through the dominant eigenvalue Re(lmax), through the standard of interaction strength s pﬃﬃﬃﬃﬃdeviation ﬃ (related to May’s complexity criterion s SC p),ﬃﬃﬃﬃﬃor ﬃ through the pairwise correlation r (related to Tang’s complexity criterion s SC ð1 þ rÞ E). We used a resampling procedure to evaluate the sensitivity of our results to parameter uncertainty. For each of the 116 Ecopath models, we proceeded as follows: we resampled each parameter, B, (Q/B) and (P/B), 1,000 times from a normal distribution N(m,s) with m ¼ Xi, s ¼ Xi/10 (corresponding to a CV ¼ 10%) and Xi is the reported value of parameter X for species i. We chose a CV of 10% because higher values could lead to negative P/B. We built a matrix of diet composition in which predators have no prey preferences (that is, they are opportunistic feeders, attacking prey in proportion to their availability). The proportion of prey i in the diet of predator j corresponds to the ratio between biomass of i and total biomass of all j’s prey species. (i) Biomass: for each of the 116 Ecopath models, 1,000 community matrices were built from the resampled values of B. Diet composition, production/biomass and consumption/biomass ratios were kept constant and corresponded to the values reported in Ecopath data. (ii) Consumption/biomass ratio: for each of the 116 Ecopath models, 1,000 community matrices were built from the resampled values of Q/B. Biomass, diet composition and production/biomass ratio corresponded to the values reported in Ecopath data. (iii) Production/biomass ratio: for each of the 116 Ecopath models, 1,000 community matrices were built from the resampled values of P/B. Biomass, diet composition and consumption/biomass ratio corresponded to the values reported in Ecopath data. (iv) Diet composition: for each of the 116 Ecopath models, a community matrix was built using the matrix of diet composition in which predators have no prey preferences. Biomass, production/biomass and consumption/biomass ratios corresponded to the values reported in Ecopath data. We calculated the dominant eigenvalue Re(lmax), the standard deviation of interaction strengths s and the correlation between pairwise interactions r of these community matrices and compared their values to the ones found in the original community matrices (Supplementary Fig. 2). We found that uncertainty in the estimation of B, P/B and Q/B had no effect on food web stability or the variables determining complexity. The absence of diet preferences was destabilizing and also decreased s.d. of interaction strength. However, we found that the deviation between the original leading eigenvalues and the ones obtained after the addition of variability in input parameters was not correlated to complexity (Supplementary Fig. 3). The complexity–stability relationship would have been biased if, for instance, the addition of variability in the biomass estimates would have a more profound impact on the leading eigenvalue of highly complex webs than the one of simpler food webs. In agreement with Barabas et al33. these results demonstrate that our ﬁndings are robust to the addition of variability in interspeciﬁc interaction terms of the community matrices. Intraspeciﬁc interaction terms of the community matrix. The diagonal elements cii of the community matrix express the strength of density dependence, which is highly stabilizing as it moves the dominant eigenvalue to more negative values: pﬃﬃﬃﬃﬃﬃ Re(lmax) ¼ R cii , where the radius of the unit circle R corresponds to s SC in random matrices and cii is the mean of diagonal elements (Introduction section and Fig. 1c). Our aim was not to assess the local stability of empirical food webs but to investigate the relationship between stability and complexity using realistic

NATURE COMMUNICATIONS | 7:12573 | DOI: 10.1038/ncomms12573 | www.nature.com/naturecommunications

7

ARTICLE

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12573

community matrices on a large gradient of complexity. In the main text, we therefore set all the diagonal elements to 0 to focus on the effect of interspeciﬁc interactions on stability. Our method is comparable to other studies that calculated stability by assessing the level of intraspeciﬁc interaction needed for all eigenvalues in a community matrix to have negative real parts (diagonal dominance)6,16,19. These studies assumed that all diagonal elements cii of the community matrix are the same. However, to obtain the community matrix, the interaction matrix A is multiplied by species biomass and the diagonal of an empirically derived community matrix should be structured by species biomass, as cii ¼ aii Bi*. Here we investigated how a non-constant diagonal, in which elements cii are proportional to species biomass, affects the dominant eigenvalues reported in our analysis. We compared the dominant eigenvalues of community matrices with aii ¼ 0 (original food webs) and aii ¼ 0.01 or 0.1 (Supplementary Fig. 6). We found that the addition of the intraspeciﬁc interaction terms was stabilizing, but had no effect on the absence of complexity–stability relationship. Food web resolution. Ecopath is mainly used for ecosystem-based ﬁsheries management. Consequently, the structure of food webs parameterized with Ecopath is often biased, with detailed compartments for harvested ﬁshes and more aggregated compartments for species at the bottom of the food web (that is, plankton and invertebrates). Food web resolution inﬂuences the estimation of species richness, connectance and interaction strength. To assess the robustness of our analysis, we investigated the complexity–stability relationship on a subset of the best resolved Ecopath models. We measured the amount of aggregation of each model, based on the criterion that groups with taxonomic name were more resolved than groups with trophic function names. We deﬁned four resolution levels and qualiﬁed one level for each species with the following indices: taxonomic species (that is, greenland turbot, index ¼ 1), family/class (that is, whales, gadoids; index ¼ 0.7), trophic function (that is, small demersal ﬁsh; index ¼ 0.4) and general name (that is, benthos, ﬁsh; index ¼ 0.1). Resolution indices RI of Ecopath models correspond to the mean resolution index of species within each food web and are listed in Supplementary Table 1. We studied the complexity–stability relationship on a subset of the 37 best resolved models (with RIZ0.7) and found results similar to the overall analysis (Supplementary Fig. 1). Data availability. The data that support the ﬁndings of this study are available from the corresponding author on request.

References 1. McCann, K. S. The diversity-stability debate. Nature 405, 228–233 (2000). 2. MacArthur, R. H. Fluctuations of animal populations and a measure of community stability. Ecology 36, 533–536 (1955). 3. Paine, R. Food web complexity and species diversity. Am. Nat. 100, 65–75 (1966). 4. Levins, R. Evolution in Changing Environments: Some Theoretical Explorations (Princeton University Press, 1968). 5. Gardner, M. R. & Ashby, W. R. Connectance of large dynamic (cybernetic) systems: critical values for stability. Nature 228, 784 (1970). 6. May, R. Will a large complex system be stable? Nature 238, 413–414 (1972). 7. May, R. Stability and Complexity in Model Ecosystems (Princeton University Press, 2001). 8. Neubert, M. G. & Caswell, H. Alternatives to resilience for measuring the responses of ecological systems to perturbations. Ecology 78, 653–665 (1997). 9. Tang, S. & Allesina, S. Reactivity and stability of large ecosystems. Front. Ecol. Evol. 2, 1–8 (2014). 10. Tao, T., Vu, V. & Krishnapur, M. Random matrices: Universality of ESDs and the circular law. Ann. Probab. 38, 2023–2065 (2010). 11. Allesina, S. & Tang, S. Stability criteria for complex ecosystems. Nature 483, 205–208 (2012). 12. Tang, S., Pawar, S. & Allesina, S. Correlation between interaction strengths drives stability in large ecological networks. Ecol. Lett. 17, 1094–1100 (2014). 13. De Ruiter, P. C., Neutel, A.-M. & Moore, J. C. Energetics, patterns of interaction strengths, and stability in real ecosystems. Science 269, 1257–1260 (1995). 14. Neutel, A.-M. et al. Reconciling complexity with stability in naturally assembling food webs. Nature 449, 599–602 (2007). 15. Emmerson, M. & Raffaelli, D. Predator-prey body size, interaction strength and the stability of a real food web. J. Anim. Ecol. 73, 399–409 (2004). 16. Neutel, A.-M. & Thorne, M. A. S. Interaction strengths in balanced carbon cycles and the absence of a relation between ecosystem complexity and stability. Ecol. Lett. 17, 651–661 (2014). 17. James, A. et al. Constructing random matrices to represent real ecosystems. Am. Nat. 185, 680–692 (2015). 18. Yodzis, P. The stability of real ecosystems. Nature 289, 674–676 (1981). 19. Neutel, A.-M., Heesterbeek, J. a. P. & De Ruiter, P. C. Stability in real food webs: weak links in long loops. Science 296, 1120–1123 (2002).

8

20. Christensen, V. ECOPATH II - a software for balancing steady-state ecosystem models and calculating network characteristics*. Ecol. Model. 61, 169–185 (1992). 21. Legagneux, P. et al. Arctic ecosystem structure and functioning shaped by climate and herbivore body size. Nat. Clim. Change 4, 379–383 (2014). 22. Christensen, V., Walters, C. J. & Pauly, D. Ecopath with Ecosim: a User’s Guide (Fisheries Centre, University of British Columbia, 2000). 23. Dunne, J. a. in Ecological Networks: Linking Structure to Dynamics in Food Webs (eds Pascual, M. & Dunne, J. A.) 27–86 (Oxford University Press, 2006). 24. Bastolla, U., La¨ssig, M., Manrubia, S. C. & Valleriani, a. Biodiversity in model ecosystems, I: coexistence conditions for competing species. J. Theor. Biol. 235, 521–530 (2005). 25. Montoya, J. M., Woodward, G., Emmerson, M. & Sole´, R. V. Press perturbations and indirect effects in real food webs. Ecology 90, 2426–2433 (2009). 26. De Angelis, D. L. Stability and connectance in food web models. Ecology 56, 238–243 (1975). 27. Paine, R. Food-web analysis through ﬁeld measurement of per capita interaction strength. Nature 355, 73–75 (1992). 28. McCann, K. S. & Hastings, A. Weak trophic interactions and the balance of nature. Nature 247, 337–345 (1998). 29. Berlow, E. L. Strong effects of weak interactions in ecological communities. Nature 398, 330–334 (1999). 30. Grimm, V. Babel, or the ecological stability discussions : an inventory and analysis of terminology and a guide for avoiding confusion. Oecologia 109, 323–334 (1997). 31. de Mazancourt, C. et al. Predicting ecosystem stability from community composition and biodiversity. Ecol. Lett. 16, 617–625 (2013). 32. Christensen, V. & Walters, C. J. Ecopath with Ecosim: methods, capabilities and limitations. Ecol. Model. 172, 109–139 (2004). 33. Barabas, G. & Allesina, S. Predicting global community properties from uncertain estimates of interaction strengths. J. R. Soc. Interface 12, 20150218 (2015).

Acknowledgements We thank students from the Chair in Biogeography and Metacommunity Ecology at UQAR for helpful comments on earlier versions of the manuscript. C.J. was supported by a grant from the Ministry of Higher Education and Research of France. C.M.’s post-doctoral fellowship grants were provided by Ressources Aquatiques Que´bec and the NSERC through the Canadian Fisheries Research Network and the FQRNT. Financial support was provided by the Canada Research Chair program and a NSERC Discovery grant to D.G.

Author contributions C.J., C.M., L.M., P.L., F.M., P.A. and D.G. designed research; C.J., C.M., P.L. and D.G. conducted research. C.J., P.L., F.M. and D.G. contributed to the analytical tools. C.J. and D.G. wrote the paper and C.J., C.M., L.M., P.L., F.M., P.A. and D.G. edited the paper.

Additional information Supplementary Information accompanies this paper at http://www.nature.com/ naturecommunications Competing ﬁnancial interests: The authors declare no competing ﬁnancial interests. Reprints and permission information is available online at http://npg.nature.com/ reprintsandpermissions/ How to cite this article: Jacquet, C. et al. No complexity–stability relationship in empirical ecosystems. Nat. Commun. 7:12573 doi: 10.1038/ncomms12573 (2016). This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

r The Author(s) 2016

NATURE COMMUNICATIONS | 7:12573 | DOI: 10.1038/ncomms12573 | www.nature.com/naturecommunications