Impact of biodiversity loss on production in complex marine food webs mitigated by prey-release

Public concern over biodiversity loss is often rationalized as a threat to ecosystem functioning, but biodiversity-ecosystem functioning (BEF) relations are hard to empirically quantify at large scales. We use a realistic marine food-web model, resolving species over five trophic levels, to study how total fish production changes with species richness. This complex model predicts that BEF relations, on average, follow simple Michaelis–Menten curves when species are randomly deleted. These are shaped mainly by release of fish from predation, rather than the release from competition expected from simpler communities. Ordering species deletions by decreasing body mass or trophic level, representing ‘fishing down the food web’, accentuates prey-release effects and results in unimodal relationships. In contrast, simultaneous unselective harvesting diminishes these effects and produces an almost linear BEF relation, with maximum multispecies fisheries yield at ≈40% of initial species richness. These findings have important implications for the valuation of marine biodiversity.

). The function derived from a fitted truncated Pareto distribution is also shown.   The five model food webs were produced from the third five of 20 runs of the PDMM assembly algorithm, whereas the empirically derived values pertain to temperate shelf communities in the Northeast Atlantic. S denotes "species richness", M mat denotes "Maturation body mass" and B a denotes "Biomass density" (Biomass, B, divided by system area, a). In calculating the slopes of the diversity spectra, lower bounds of 1 kg were used except for Model food web 15, for which a lower bound of 3 kg was used. The five model food webs were produced from the last five of 20 runs of the PDMM assembly algorithm, whereas the empirically derived values pertain to temperate shelf communities in the Northeast Atlantic. S denotes "species richness", M mat denotes "Maturation body mass" and B a denotes "Biomass density" (Biomass, B, divided by system area, a). In calculating the slopes of the diversity spectra, lower bounds of 1 kg were used except for Model food webs 16 and 20, for which lower bounds of 35 kg and 4 kg were used, respectively.

Further details on the Population-Dynamical Matching Model (PDMM)
Overview of PDMM: Model elements. The PDMM was originally designed to extend the understanding of food-web topology represented by the Matching Model 22 , by adding an explicit representation of interspecific size-structure and biomass dynamics through a system of differential equations 23 . Food-web topology and the strengths of trophic (feeding) interactions in the PDMM are determined by the body sizes of consumer and resource species together with abstract trophic traits of consumers and resources 24 . These emerge from a community assembly process consisting of a cycle of species introductions and simulation of food-web dynamics to a new equilibrium 23 . In subsequent studies, the PDMM was developed and used to represent temperate shelf communities in the Northeast Atlantic [25][26][27] . The model has been shown to reproduce structural 26 and dynamic 25,27 properties of temperate shelf food webs. For the present study we have incorporated additional minor modifications to the PDMM to improve control over species richness in mature communities during community assembly and to increase its performance in reaching equilibrium dynamics. In every other aspect, it is identical to the PDMM implementation in Fung et al. 26 In the remainder of this subsection and the next two subsections ("Overview of PDMM: Community assembly" and "Overview of PDMM: Population dynamics and foodweb structure"), the structure of the PDMM is summarised. The next three subsections ("Detailed description of PDMM: Community assembly", "Detailed Description of PDMM: Population dynamics of producer species" and "Detailed Description of PDMM: Population dynamics of consumer species") give a complete model description, highlighting changes made in comparison with Fung et al. 26 The parameter values used for the PDMM implementation in the present study are also largely the same as those in Fung et al. 26 , except for changes arising from the minor modifications made. These changes in parameterisation of the PDMM are described in the last subsection of this section ("Model parameterisation").
Overview of PDMM: Community assembly. The PDMM uses a stochastic, iterative assembly algorithm to construct dynamic model food webs with many species and realistic interspecific body size and trophic link structures 23,26 . One model food web is produced from each run of the algorithm. The algorithm starts with an empty model food web (with no species) and then repeatedly adds sets of new species to this web, representing invasions 23,26 .
Apart from the first set of species, which have body sizes that take fixed default values and abstract trophic traits with values chosen from random distributions, new species have trait values that are determined by randomly changing those of a species already found in the food web. This has the effect of implicitly emulating the effect of phylogenetic structure in a hypothetical pool of species from which species invade, as well as evolution of species in this pool 28 . The resulting phylogenetic correlations of the trait values of species within the food web have been found to be important determinants of food-web structure 22,[29][30][31] . In between each set of species additions, food-web biomass dynamics are simulated to represent the dynamic response of the food web to the species additions 23,26 . During simulation of dynamics, extinct species are removed from the food web, representing a natural selection filter for the species subset that can coexist dynamically. The resulting food web emerges from a simulation of natural dynamic evolutionary and population processes, rather than being pre-specified. This gradual assembly of model food webs using a developmental approach (sensu Taylor 32 ) enables the construction of food webs with up to thousands of dynamically coexisting species 26 .

Overview of PDMM: Population dynamics and food-web structure. Food-web dynamics
in the PDMM are specified by linked differential equations that describe gains and losses in the biomass of each species' population due to trophic interactions, metabolism and natural mortality. Equations for producers differ from those for consumers, to reflect their differences in energy acquisition and growth. Specifically, if the S species in a PDMM food web are indexed such that the S P producer species are allocated indices of 1 to S P and the S  S P consumer species are allocated higher indices, then the rate of change of the biomass of a producer species population i, 1  i  S P , is given by: where B i denotes the biomass of a producer species i, the coefficient G i denotes the realised net growth rate of i (discounting losses due to metabolism), and f ij is the functional response for consumer species j feeding on i. Producer species i grows at a maximum net growth rate  max, i in the absence of other producer species and consumption by consumer species. This maximum net growth rate is reduced in the presence of other producer species, due to competition for limiting resources such as light and nutrients and also transmission of diseases (by pathogens such as viruses). The exact form of G i used for the PDMM version in this study differs from that for the PDMM version used in Fung et al. 26 and is described in detail in the subsection "Detailed description of PDMM: Population dynamics of producer species", together with the underlying rationale.  max, i scales with the maturation body mass of i, one of the traits modelled for each species, according to an empirically-derived allometric equation 6 .
The functional response f ij is a non-linear Holling type II response extended to include a mechanistic theory of prey-switching among multiple resource species 33 . Functional responses are parameterised by the trophic interaction strengths, c kj , between pairs of species 33 . In the PDMM, these parameters depend on the traits of the interacting species.
Specifically, the value of c kj is a product of two factors. The first factor describes the dependence of trophic interaction strengths on the ratio between consumer and resource body masses 8,9,34,35 . The second factor describes the dependence of interaction strengths on the abstract trophic traits of the interacting species. Abstract trophic traits of a species are given by numerical vectors of "foraging" traits and "vulnerability" traits, each of length D, and trophic interactions become stronger when the foraging traits of the consumer j are closer to the vulnerability traits of the resource k [22,24]. Thus, consumption rates are not exclusively determined by the size preference of consumers (preferred ratio between maturation body mass of a consumer species to that of a consumed species), implying that the trophic level of a consumer species does not necessarily increase with its body size (maturation body mass), as observed 8  The rate of change of the biomass of a consumer species population i, S P 1  i  S , is given by: where  is the consumer assimilation efficiency and l c, i is the metabolic loss rate of consumer species i. Parameters of the functional response and metabolism are chosen such that populations of consumer species grow at a maximum net rate of r max, i when resources are not limiting and predation is absent. Both r max, i and l c, i scale with the maturation body mass of i according to empirically-derived allometric equations 5,36 . In equation (1.2), the first term on the right-hand side represents the rate at which consumer species i produces biomass following consumption, and the sum of these terms for all model fish species is used as the measure of production in this study. It is noted that the terms represent the biomass assimilated by fish species and therefore available for human harvesting, which is more directly relevant in the ecosystem service context considered in this study. In other contexts, the biomass ingested may be more appropriate. However, the biomass ingested is simply the biomass assimilated divided by  . Thus, the results derived in this study for biomass assimilated would also hold if the biomass ingested was used instead, up to a constant scaling factor.  If there are no species in the existing food web of the same type as species k, then the trophic trait values of the latter are assigned default values. Following Fung et al. 26 , M mat,k is set to a default value of 10 -10 kg for producers and a default value of 10 -7 kg for consumers. In addition, the trait vectors V k and F k are sampled from the corresponding hyperspheres assuming even probability distributions 26 .
Sets of trophic trait values for each new species k are sampled as described above until a set is found giving a positive growth rate (invasion fitness 36 ) in the invaded model food web.
This leads to a set of new species that can all invade the existing community. These are added to the food web and food-web dynamics with the new species are simulated until an equilibrium is reached or 500 yrs have elapsed. Food-web dynamics are considered to have reached an equilibrium when the biomass dynamics of all species in the food web have reached an equilibrium, according to criteria as described in Fung et al. 26 The waiting time of 500 yrs is longer than the 200 yrs used in Fung et al. 26 and is used to increase the probability that the dynamic response of a model food web to a set of species additions is fully captured.
An equilibrium was reached during the waiting time of 500 yrs for 92% of approximately 1,100 iterations of the assembly algorithm used in the present study. interactions among producer species were controlled by an additional set of abstract traits. In the present PDMM implementation a logistic model for G i is used instead that allows finer control of producer species richness, following an approach used frequently in food-web models [37][38][39] . Specifically, producer species i grows logistically according to where  max, i is the maximum net growth rate of i,  ij measures the negative effect of producer species j on the growth rate of i due to competition for limiting resources or transmission of diseases, and K j is the carrying capacity of j.  max, i is given by an allometric equation relating it to the maturation body mass of i, following an empirical study 6 : The value of  ij is set to one when i  j , reflecting intraspecific competition that is always present, and is otherwise set to a random value that may or may not be zero, reflecting random interspecific antagonistic interactions. Mathematically,  [17]. The carrying capacity of species j, K j , is derived from parameters for which empirical data are available, following the approach of previous implementations 23, 26 . In a monoculture of producer species j, the maximum net rate of biomass production in the model is  max, j K j . This is set equal to the empirical maximum net rate of biomass production attained by a producer species in monoculture or when dominating a community, which is taken to be a constant based on the observed size-independence in real systems 40 . Thus,  max, j K j  P N,max a , where P N,max is the empirical maximum net rate of biomass production attained by a producer species in monoculture per unit planar area and a is the empirical planar area of the system modelled. This gives, (1.9) The loss rate of producer species i due to predation by consumer species j is determined by the functional response f ij , details of which are presented in the subsection "Detailed description of PDMM: Population dynamics of consumer species" below. This has the same form as in Fung et al. 26 and is given by which is a Holling type II response extended to include a mechanistic theory of preyswitching 33 . In this functional response, b j is the attack rate of consumer species j, T j is the handling time of consumer species j, s ik is the similarity between two prey species i and k with respect to prey-switching, and c ij is the trophic interaction strength between resource species i and consumer species j. b j scales allometrically with the maturation body mass of the consumer j [7,11,12], M mat, j , and can thus be written as Following Fung et al. 26 , T j can be expressed as where  is the consumer assimilation efficiency, r max, j is the maximum net growth rate of consumer species j, and l c, j is the metabolic loss rate of consumer species j. Both r max, j and l c, j scale allometrically with M mat, j (subsection "Overview of PDMM: Population dynamics and food-web structure"): (1.14) The coefficient s ik is modelled as a decreasing function of the Euclidean distance between the vulnerability trait vectors of the two resource species i and k [26]. Specifically, and where w c controls how strongly c ij depends on the proximity of the foraging and vulnerability trait values, and is called the consumer niche width;  c determines the rate of decrease of c ij when the consumer-resource maturation body mass ratio M mat, j M mat,i increases above the preferred ratio  c ; and  c determines the rate of decrease of c ij when M mat, j M mat,i decreases below  c . From the subsection "Model parameterisation" below, the empirically derived value of  c  10 2 [8,9]. In addition, following the logic of Fung et al. 26 ,  c and  c were chosen to be 0.07 and 0.25 respectively, such that c ij decreases slowly when the resource species i decreases below the preferred maturation body mass of consumer species j and decreases quicker when it increases above the preferred mass. These values ensure that consumer species j feeds on resource species with a wide range of body masses, which is typical for marine organisms that change their body masses over many orders of magnitude as they grow from larvae to adults (Rossberg 35 , Section IX.C). For example, cod (Gadus morhua) feed on resource species with a range of maturation body masses spanning more than 7 orders of magnitude, representing small plankton to adult fish 1,3,18,41 .
It is noted that the trophic interaction coefficient c ij depends not only on maturation body mass, but also on the foraging traits of consumer species j and the vulnerability traits of resource species i. Therefore, species j does not necessarily feed most on resource species with the preferred body mass, since consumer preferences are also determined by other traits. For example, Supplementary Fig. 1 shows the realised consumer-resource maturation body mass ratio for each fish species in one of the PDMM food webs used, on a logarithmic scale (Model food web 1 in Supplementary Table 3). This ratio is calculated as the weighted  rather than that of the Celtic Sea. However, following the mathematical argument by Fung et al. 26 , changing a has no major effect on model dynamics. C lc is now set to 0.259 kg 1/4 yr -1 rather than 0.298 kg 1/4 yr -1 . The new value is the average of two empirical values from Brown et al. 5 for the metabolic loss rate of consumer species, the first derived using empirical data for fish and the other derived using empirical data for invertebrates, including aquatic types. The value of  c is set to 10 2 instead of 10 3 . Both values are consistent with the empirical range of 10 2 -10 3 found using consumer-resource body mass ratios for a demersal fish community 8 and a demersal fish and invertebrate community 9 , both in the North Sea.
There are no empirical data to directly constrain  and p  , the remaining two new parameters in the PDMM implementation used here. After all parameters for which empirical data is available had been fixed as described above, the values of these two parameters were varied heuristically together with those of 14 other parameters lacking empirical data (Supplementary Table 1), until model food webs could reasonably reproduce values of key properties of temperate shelf communities in the Northeast Atlantic (see section "Full details of validation of PDMM food webs" for these properties and corresponding empirical ranges). This heuristic approach follows that in Rossberg et al. 23 and Fung et al. 26 ; it is necessary because assembly of each model community is computationally expensive, lasting from days to weeks even with the use of multi-threading over eight processors on a computer cluster. The method gave   0.4 and p   0.0025, and values of the other parameters that are the same as in Fung et al. 26 , except for

Full details of validation of PDMM food webs
The 20 complex model food webs used in the present study were generated using 20 independent runs of the PDMM assembly algorithm. In each run, the algorithm was stopped consumed by a fish species, counting all resource species that contribute >1% to the diet of a fish species 43 . For (vi), the trophic level for a fish species was calculated as 1 plus the weighted average trophic level of the species that it consumes, with weights given by the biomass intake rates normalised by the total biomass intake rate, i.e. proportional contributions to diets. For the present study, these six properties were supplemented by four more properties: (vii) the diet-partitioning exponent, which is estimated as the slope of a linear regression of the logarithm of the dietary diversity of fish species against the logit of the proportional diet contribution above which a consumed species is counted as a resource 16,17 ; (viii) the slope of the diversity spectrum, which is the slope of a fitted probability density function of asymptotic body masses of fish species on a log 10 -log 10 scale 19 ; (ix) the slope of the biomass size-spectrum, which is the slope of a linear regression of biomass against individual body mass for fish species on a log 10 -log 10 scale 20 ; and (x) the range of biomass densities for fish species. Use of (vii), (viii) and (ix) allowed more rigorous quantitative assessments of trophic link structure, size-structure among species and sizestructure among individuals, respectively. Use of (x) allowed a quantitative assessment of fish abundances. Together, the 10 food-web properties examined are major determinants of food-web structure and dynamics.
In calculating the slope of the diversity spectrum for a model food web, the method detailed in Section S10.2 of the Online Appendix of Reuman et al. 19 19 The slope of the diversity spectrum is derived as the value of the shape parameter for the truncated Pareto distribution minus one, following calculations in Reuman et al. 19 (Section S10.2 of their Online Appendix).
To calculate the slope of the biomass size-spectrum of a model food web, the graph of biomass against individual body mass for the model fish species was drawn first, on a log 10log 10 scale. Since the PDMM does not resolve the intraspecific size-structure of model species, we assumed that each model fish species at equilibrium has an intraspecific sizestructure of the form derived by Rossberg 35  Fish species richness of nine model food webs were within the empirical range, with richness for the remaining 11 webs being within 23% of the empirical range (Table 1; Supplementary Tables 3-6). Model phytoplankton species richness consistently exceeded the highest empirical estimate, by 51-74% (Table 1; Supplementary Tables 3-6). The dietary diversity of fish species 16 quantifies the trophic link density among fish species, and thus implicitly quantifies their trophic niche widths. Model values for this property fell within the empirical range (Tables 1; Supplementary Tables 3-6). In addition, model values of the diet partitioning exponent 16 , which reflects the distribution of trophic link strengths, also fell within the corresponding empirical range (Table 1; Supplementary Tables 3-6).
Supplementary Fig. 3 graphically depicts the relationship between dietary diversity and diet threshold for 10 of the 20 model food webs. In every web, the maturation body masses of all phytoplankton and fish species fell within the empirical ranges (Table 1; Supplementary   Tables 3-6). Furthermore, >86.2% of fish species in any one food web had a trophic level within the empirical range, whereas >96.4% of fish species had a biomass density within the empirical range (Table 1; Supplementary Tables 3-6). Moreover, in every web, the maturation body masses, trophic levels and biomass densities described cover 54.6-100% of the corresponding four empirical ranges (using a log 10 scale for body masses and biomass densities), except the biomass densities for one web, which cover 41.3% of the corresponding empirical range (Table 1; Supplementary Tables 3-6).
To further characterise the distribution of species among different size classes, we calculated the slope of the diversity spectrum for each model food web 19 . Using a lower asymptotic body mass bound of 1 kg [19] (0 on a log 10 scale), slopes for 11 food webs fell within the empirical range, whereas slopes for five food webs were within 9% of the empirical range (Table 1; Supplementary Tables 3-6). Slopes for the remaining four webs (0.0000976-0.654) fell outside the empirical range, but fell within the range when the lower asymptotic body mass bound was increased to 3-35 kg (0.477-1.54 on a log 10 scale; Table 1; Supplementary   Tables 3-6). All fitted functions (considering a 3-35 kg lower asymptotic body mass bound for the four webs mentioned) gave a good fit to the underlying data from the model food webs, with R 2  0.92 ( Supplementary Fig. 4 shows the diversity spectrum for one of the 20 PDMM food webs together with the fitted spectrum). We also estimated the slope of the biomass size-spectrum for each model food web 20 . Considering the body mass range of 0.1-10 kg [20] (-1 to 1 on a log 10 scale), slopes for 16 food webs fell within the expected range predicted by applying macroecological theory to empirical data, with slopes for the remaining four webs being more negative by 9-114% (Table 1; Supplementary Tables 3-6; Supplementary Fig. 5 shows the biomass size-spectra for 10 of the PDMM food webs). The linear regressions were generally a good fit to the underlying data from the model food webs, with R 2  0.67 for 12 webs; however, they gave poor fits for eight webs, with R 2  0.14. It is noted that the biomass size-spectra for the model webs all exhibit a declining or approximately linear trend, such that total biomass of large individuals tends to be no greater than that of smaller individuals. This does not contradict the slight increasing trends observed for species population biomass with maturation body mass, since there are typically fewer species at large maturation body masses ( Supplementary Fig. 2).

Interpretation of model validation results.
Food-web structure can be characterised by the number, biomasses, body sizes and trophic link strengths of constituent species 45,46 .
Together, these structural elements are fundamental determinants of energy flows between species, so reproducing these elements is expected to capture essential aspects of food-web dynamics and stability. We compared values of 10 properties quantifying these elements for the 20 complex model food webs used with empirically derived values from temperate shelf communities. There was largely good agreement, allowing the model webs to be taken as realistically representing the food-web structure of temperate shelf communities. Although model phytoplankton species richness exceeded the highest empirical estimate by 51-74% (Table 1; Supplementary Tables 3-6), empirical counts are likely to be underestimates because of incomplete sampling and identification 47 . The slopes of the model biomass sizespectra were largely within the range derived by applying macroecological theory to empirical data 20 , but the corresponding linear regressions used to calculate the slopes often gave poor fits to the model spectra. This could be because of the presence of food-web effects independent of body size in the PDMM 23 , which introduce extra variation not captured in macroecological theory 20 .
The realistic structure of the model food webs is expected to give realistic production values.
Indeed, Heath 48 estimated that fish in the North Sea consume 13.9-19.5 gC m -2 yr -1 over the period 1973-1999. Using conversion factors of 1 gC = 1/0.4 g dry weight 7 and 1 g dry weight = 1/0.32 g wet weight 49 , and assuming an assimilation efficiency of 0.6 [10,50], total fish production in the North Sea was 65.2-91.4 g m -2 yr -1 . This is in good agreement with the mean value of 96.6 g m -2 yr -1 for the model webs, especially considering that the empirical range could be an underestimate due to historical fishing 20 .

Partitioning variation in production
By performing random deletion experiments with 20 PDMM food webs, the relationship between (normalised) fish species richness and mean total fish biomass production was derived in the main text, together with the standard errors in this production (Fig. 1). The standard error at a particular richness value measures the variation in the mean production, not in the production values. Thus, in Supplementary Fig. 10, we show the mean relationship together with the standard deviations in production. This variation in production exhibits an increasing trend with richness. However, the coefficient of variation remains largely similar across the entire range of richness values ( Supplementary Fig. 8).
In order to obtain a deeper insight into the sources of variation, the variance in production was partitioned according to variance between food webs and within food webs. This was done by considering each fish species richness value in turn, calculating the sum of squared differences between all production values and the mean production for the species richness value considered, and then partitioning this sum into contributions due to differences between and within webs. More precisely, at a particular richness value of    x , denote the number of corresponding production values produced from the 10 random deletion experiments for PDMM food web i (1  i  20 ) by n x,i and the production values by P x, ij (1  j  n x,i ). Then the sum of squared differences is where P x,  is the mean over all i and j. The contribution to this sum due to variation in production between food webs is calculated as where P x, i is the mean over all j for given i. The contribution to the sum due to variation in production within food webs (arising from use of 10 random deletion experiments for each web) is calculated as Supplementary Fig. 9 shows that SS x,b dominates SS x,w for a large range of  (approximately   0.1) and tends to increase with  . This shows that increasing variation in production is due predominantly to greater differences between food webs rather than differences within food webs.

Mathematical analyses of BEF relations
Mean-field theory. In a Lotka-Volterra model describing the food-web dynamics of S fish species, the rate of change of the biomass density of species i in kg m -2 ,  i , is given by: where r i is the intrinsic population growth rate of species i in yr -1 and  ij is the interaction coefficient quantifying the effect of species j on i in kg -1 m 2 yr -1 . r i represents gains in the biomass of fish species i because of consumption of non-fish species minus losses due to non-predation processes, largely metabolism. We first examine the general case where the interspecific interaction coefficients can be different from each other, such that  ij does not have to be the same as  ij . This general case encompasses food webs with their asymmetric interactions. In the subsection "Implications of competition symmetry" below, we will examine the special case where the interspecific interaction coefficients are symmetric, i.e.
 ij   ji , in order to provide a contrast to results for the general case that we derive in this subsection. The symmetric case pertains to species that partake predominantly in competitive rather than trophic interactions.
In the general case, we apply a mean-field approximation, under which it is assumed that interspecific interaction strengths are largely independent, such that Here, the angled brackets denote expectation values for random choices of quantities with unsummed indices inside the brackets 51,52 . When evaluating the expected interaction coefficient,  ij , there is a need to distinguish between the expected contribution from off-diagonal terms,  ij i j , and the expected contribution from diagonal terms,  ii , because the latter tends to be systematically larger than the former. This gives Using the mean-field approximation, the equilibrium conditions for the dynamic system specified by equation 4.1 can be evaluated to give [51,52]. The mean total biomass density,  , can then be obtained as   S  i . It is assumed that the mean total biomass production per unit area, P , is equal to  multiplied by a constant C with units of yr -1 . Under this assumption, the mean-field approximation Denoting the number of fish species in the pristine (unfished) state as S pris , the equation for P can be rearranged to give a Michaelis-Menten (MM) function: where   S S pris is the normalised fish species richness, A  C r i  ij i j and Thus, (4.8) In random-deletion simulations using 20 PDMM food webs (10 replicates per web), the average s across all replicates and food webs is 0.26, which is small. Thus, in equation 4.8, s can be taken to be 0 as a simplifying approximation, which gives The first term on the right of P equals  P S, which is the expected immediate loss in P due to random deletion of a fish species. Therefore, the second term on the right must represent the expected change in P due to dynamic responses following random deletion of a fish species. This can be verified by observing that the second term is equal to represents the expected change in production of a fish species resulting from random deletion of another fish species. Thus,  d P S   dS multiplied by S represents the expected change in the total production of all remaining fish species following the random deletion, which is the expected change in P due to dynamic responses following the deletion. In Fig. 2a in the main text, the two terms of equation 4.9 are plotted as functions of normalised fish species richness,  . To achieve this, P is rewritten as (4.10) The second term, corresponding to the effect of dynamic responses on P, is consistently higher by small amounts than simulation results using the 20 PDMM food webs. This discrepancy is a consequence of the approximation of disregarding secondary extinctions ( s  0) in the derivation of the second term. When including this effect in the derivation (using equation 4.8), the corresponding correction to the second term is   increases from one equilibrium state to another during a modelled community assembly process, provided that interactions are symmetric.
Assume the community has reached a locally asymptotically stable equilibrium with S species. How will Q change if a subset of the S species in the community is selected and deleted from the community, and populations are allowed to relax until a new equilibrium is reached? To answer this question, call the original S-species equilibrium community the high-diversity community. Denote as the low-diversity community the community at the new equilibrium after species deletion -this community is obtained from the high-diversity community by deleting the subset of selected species as well as any species that reach zero abundance during dynamics to the new equilibrium (secondary extinctions). Because ij  is symmetric, the local asymptotic stability of the high-diversity community implies its global stability. * Thus, if one re-inserts all species deleted from the high-diversity community into the low-diversity community at low abundances, the original high-diversity community will re-emerge. By the considerations above, the value Q will increase in this process. This implies directly that Q is lower for the low-diversity community than it is for the high-  as a measure of total community production. This will generally be different from the measure of total system production used in the main text. The Lotka-Volterra model is simply too abstract to isolate gross production among the various biological contributions to population dynamics. On the other hand, the approximation that production is proportional to biomass for a given species (constant production/biomass ratio or turnover rate) is frequently applied in the ecological * Local stability of the system at an inner ( 0   . From Sylvester's law of inertia 58 it then follows that the matrix α  has the same number of negative and zero eigenvalues as J , and no positive eigenvalues. Now, consider first the case that zero is an eigenvalue of α , and let u be a vector in the null space of α . It is then not difficult to verify that with any equilibrium i  of the Lotka-Volterra system and any real number  , the point Despite these caveats, however, we have shown that interaction symmetry is essential for obtaining a steady increase of production with each species added to a community and a steady decrease of production with each species deleted, where in both cases secondary extinctions might occur as a result of the manipulations. For strongly asymmetric interactions, this cannot be expected. However, as detailed in the subsection "Mean-field theory" above, in the specific case of random additions and deletions of species, these trends can still be expected to be seen on average. Thus, the key difference between the general case where interactions can be asymmetric and the specific case where interactions are symmetric is that in the former case, production increases on average with increasing richness, whereas in the latter case, production always increases with increasing richness.