Stock-recruitment relations controlled by feeding interactions alone

Empirical stock recruitment (S-R) relationships exhibit clear evidence of density dependence. Since the early works of Ricker (1954), Schaefer (1954), Beverton and Holt (1957) and others, characterizing and interpreting these relationships, much progress has been made in linking this density-dependence to ecological mechanisms. Reviewing our current understanding, Houde (2008, 2009) stresses the distinction between a variety of oceanographic, climatic, and ecological mechanisms controlling recruitment and manifestly density-dependent mechanisms contributing to population regulation. While all controlling mechanisms combine to generate the observed large variability in recruitment even at constant spawning-stock biomass (SSB), only regulating mechanisms can explain the observed non-linear dependence of recruitment R on SSB. Among possible regulating mechanisms that have been considered are competition for refuge or territory, diseases or parasitism at various life stages, larval competition for oxygen, and densitydependent trophic (feeding) interactions, i.e., food intake and predation mortality, with cannibalism as a special case (Houde 2009). The goal of the present study is to test the plausibility of the hypothesis that S-R relationships are primarily shaped by the density dependencies of trophic interactions. In addition to its relevance for the interpretation of S-R data, this hypothesis also has a bearing on questions such as to what extent trophic interactions are relevant for management models or what the mechanisms are that determine the carrying capacity of individual stocks and entire communities (Batchelder and Kim 2008). Management aiming at maximum sustainable yield (MSY), in turn, hinges on a good understanding of density dependencies and carrying capacity. This is clear already from the textbook analysis of MSY for the logistic model (Conroy and Carroll 2009), where the optimal population biomass comes out as exactly half the carrying capacity.


Introduction
Empirical stock recruitment (S-R) relationships exhibit clear evidence of density dependence.Since the early works of Ricker (1954), Schaefer (1954), Beverton and Holt (1957) and others, characterizing and interpreting these relationships, much progress has been made in linking this density-dependence to ecological mechanisms.Reviewing our current understanding, Houde (2008Houde ( , 2009) ) stresses the distinction between a variety of oceanographic, climatic, and ecological mechanisms controlling recruitment and manifestly density-dependent mechanisms contributing to population regulation.While all controlling mechanisms combine to generate the observed large variability in recruitment even at constant spawning-stock biomass (SSB), only regulating mechanisms can explain the observed non-linear dependence of recruitment R on SSB.Among possible regulating mechanisms that have been considered are competition for refuge or territory, diseases or parasitism at various life stages, larval competition for oxygen, and densitydependent trophic (feeding) interactions, i.e., food intake and predation mortality, with cannibalism as a special case (Houde 2009).
The goal of the present study is to test the plausibility of the hypothesis that S-R relationships are primarily shaped by the density dependencies of trophic interactions.In addition to its relevance for the interpretation of S-R data, this hypothesis also has a bearing on questions such as to what extent trophic interactions are relevant for management models or what the mechanisms are that determine the carrying capacity of individual stocks and entire communities (Batchelder and Kim 2008).Management aiming at maximum sustainable yield (MSY), in turn, hinges on a good understanding of density dependencies and carrying capacity.This is clear already from the textbook analysis of MSY for the logistic model (Conroy and Carroll 2009), where the optimal population biomass comes out as exactly half the carrying capacity.
Our approach to this set of interrelated problems is to construct a species-rich fish-community model in which trophic interactions are the only source of density dependence, and to study in detail the resulting S-R relationships, relating them to observations.The model we study is a variant of the Fish Community Size Resolved Model (FCSRM) (Hartvig et al. 2011;Andersen and Pedersen 2010;Houle et al. 2012), developed based on theory by Andersen and Beyer (2006).In this model, individual fish interact exclusively through predator-prey relationships.Just as in the original model formulation by Hartvig et al. (2011), individual reproductive output is a linear function of food intake (with negative intercept), that is, S-R relationships are not a priori built into the model.Density dependencies in the model arise only from the dependence of growth and fecundity on food availability and the dependence of mortality on the abundance of predators.Individuals are characterized by their momentary body mass and their species identity, with the latter determining prey, predators and asymptotic body mass.Intentionally, the model describes fish physiology and ecology in a highly simplified form, so simplifying the interpretation of model outputs.Yet, the model is realistic enough to choose all model parameters on empirical grounds, and the resulting model communities resemble communities of interacting fishes in many aspects (Andersen et al. 2008;Hartvig et al. 2011;Rossberg 2012).
The simulated S-R relationships for each model species were obtained by adding a constant µ a to the mortality of that species (leaving that of other species unchanged) and recording SSB and recruitment R as µ a gradually increases.The resulting pressure-response relationships are similar to those produced by fishing-where µ a would correspond to fishing mortalitybut avoid complications such as those resulting from mixed (multi-species) fisheries or gear selectivity by size.The simulated S-R data were then evaluated with respect to visual patterns, the best-fitting functional relationships, heteroscedacity of variations around this curve, and their steepness (Mace and Doonan 1988), revealing surprisingly good agreement with observations.

Simulation model
Most aspects of the FCSRM have been described and motivated in great detail by previous authors (Andersen and Pedersen 2010;Hartvig et al. 2011;Houle et al. 2012).Here, the model structure is outlined just as much as required to interpret simulation results.Only two innovations that we introduce for the present study are described in more detail.
In view of the tremendous numbers of offspring produced by most fishes, individual-based models that track the full life cycle of each individual can be computationally very expensive.Instead, the FCSRM keeps track of the body-size distribution of individuals in each population.This distribution is represented by a time-dependent histogram with bins equally spaced along the logarithmic body-mass axis.Bin width is approximately log 1.1, corresponding to a factor of 1.1 in body mass.
The modelled fish feed either on other fish or on a community of smaller non-fish species (e.g.zooplankton), summarily described by a resource size spectrum.The size-resolved population dynamics of the resource size spectrum is described by a simple semi-chemostatic growth model without any coupling between size classes.
Food intake of individual fish is determined by a bell-shaped predator-prey size ratio (PPMR) window (here in a modified form following Houle et al. (2012)), a trophic link strength θ ij , depending on the species identities of prey i and predator j, and a Type II functional response to prey availability.Assimilated food is first discounted by density-independent respiration, and the remainder is invested into growth, g j (m), and reproduction, g r,j (m), both of which depend on individual body mass m and species identity j.The proportion g r,j (m)/[g j (m)+ g r,j (m)] of assimilated intake invested into reproduction is controlled by a maturation selection function, which is described as a product of a maturity ogive and a functional form chosen such that adults grow following approximate von-Bertalanffy trajectories (Hartvig et al. 2011).Different species have different asymptotic body masses M j .Maturation body mass, defined as the inflection point of the maturity ogive on the log m axis, or the body mass at 50% maturity, is set to be m * j = M i /4 (Beverton 1992).The size of the offspring produced, on the other hand, is chosen following Hartvig et al. (2011) to be the same value m 0 = 0.5 mg for all species (Andersen et al. 2008).
Natural mortality in the FCSRM has three components: (1) food intake by all individuals is matched by corresponding mortality of their prey; (2) if food is insufficient to compensate for respiration losses, starvation mortality sets in; and (3) a small additional background mortality accounts for death of predominantly large adults for reasons not explicitly described by the model.
Standard allometric scaling laws are employed for the dependencies of maximal food-intake rate, respiration rate, and background mortality rate on the body mass of individuals.Similarly, an allometric law controls the scaling of search/attack rates with body mass, as a third factor together with PPMR window and link-strength matrix.
Simulations are initialized with S ini species with asymptotic body masses sampled at random from an even distribution between 0.1 g and 100 kg on a logarithmic axis.Initial population sizes and structures are chosen according to generic theoretical considerations (Andersen et al. 2008), and simulations run until a quasi-steady state is reached.Populations reaching very low biomasses are removed as extirpated.The quasi-steady state finally reached tends to exhibit complex oscillations, either regular or chaotic (Hartvig et al. 2011).
The FCSRM in the form originally proposed by Hartvig et al. (2011) admits coexistence of only a limited number of species (at most 9, on average 5).This is likely a result of strong interspecific competition.To overcome this limitation, the model was modified in two crucial ways.First, rather than sampling the entries of the link-strength matrix θ ij from an even distribution between 0 and 1 (Hartvig et al. 2011), we sampled them as with σ = 2 √ 2 ln S ini and ξ ij denoting S ini × S ini independent standard normal variables.This log-normal distribution of link strengths is suggested by the observed statistics of fish diet partitioning (Rossberg et al. 2011, appendices).The term −σ 2 /2 in Eq. [2] guarantees that the community mean of θ ij is close to 1.With link strengths chosen as in Eq. [2], diets are more narrow and competition is reduced as compared to the model variant by Hartvig et al. (2011).
The second modification we introduced carries this idea over to feeding on the resource size spectrum.Assuming that different species are specialized to feed on different parts of the resource spectrum, each species is assigned its "own" resource spectrum, which dynamically responds to predation pressure from this species alone.
With these two modifications, we find that arbitrary large numbers of species can be brought to coexist in the model, limited only by the available computational resource.Typically, half the initial number S ini of species survive to co-exist in the model's quasi-steady state.Here we set S ini = 40, giving communities with about 20 co-existing species.To improve the statistical power of our analyses, data from three independent model communities were combined, yielding a total of 77 simulated S-R data sets.In view of the complexity of dynamics in these communities, it is legitimate to consider the 77 data sets, for the simple analyses performed here, as statistically independent.

Quantifying SSB and recruitment
With N i (m) denoting the density of individuals along the body mass (m) axis, we computed the SSB of each population as where Ω i (m) is the maturity ogive, defined in Eq. [1].All integrals are evaluated numerically as sums over log m bins after a corresponding change of the integration variable.
To obtain a measure of recruitment consistent with this definition of SSB, we computed where g i (m) is the somatic growth rate (Mass/Time) of species i at size m, and Ω i (m) = dΩ i (m)/dm is a hump-shaped function with mode near maturation body mass m * i .The value of R i is time-dependent through the time dependencies of both growth g i (m) and population structure N i (m).The definition c 2013 Crown copyright of R i in Eq. [4] is motivated by the observation that the expression g i (m)mN i (m) represents the mass per unit time of individuals reaching size m.Through the integration over Ω i (m), this value is averaged over the transition to maturity.The reason for averaging over the quantity g i (m)mN i (m) rather than g i (m)N i (m), the number per unit time reaching size m, is that its dependence on m is weaker (Rossberg 2012).In the idealization that maturation at m * i is a sharp transition and Ω i (m) a step function, R i /m * i equals the rate g i (m * i )N i (m * i ) at which individuals reach maturity, so that R i is proportional to recruitment by the conventional definition.Equation [4] simply takes into account that, in fact, maturation is a gradual process.

Simulated S-R relationships
As explained above, the simulated S-R relationship of each was obtained by increasing its added mortality µ a over time t, slow enough that the community always remained in a quasi steady state, taking the quasi steady state of the community without additional mortality as the starting point.Specifically, µ a was chosen as µ b ×0.02 yr −1 t, where µ b = 0.84 (M j /4) −1/4 yr −1 is the natural background mortality of species j (Hartvig et al. 2011).Eventually, a value µ a = µ lim was reached where the focal species became extirpated and simulations stopped.In our model the age at maturation depends on the somatic growth rate and therefore varies through time.Because of this complication, the delay between spawning and recruitment was not taken into account when constructing S-R graphs.
To make model fitting independent of assumptions regarding the distribution of residuals, 50% quantile regression (Koenker and Park 1994) in the implementation in R by Koenker (2012) was used.The algorithm determines the model parameters α, β, and γ such that, over any range in SSB, approximately 50% of data points lies above and below the fitted curve.To avoid numerical artifacts, only data points with SSB larger than 0.1% of the largest observed value were used for the model fits.

Heteroscedacity
The dependence of the variance of residuals of log R on SSB was investigated by Iles and Beverton (2000) and Minto et al. (2008) to probe for regulating mechanisms.For an easy comparison, we determined the heteroscedacity parameters η 1 of Minto et al. (2008) for our simulated data sets, following the method of Minto et al. (2008) in every detail; in particular, maximum likelihood fits assuming log-normally distributed R were used, rather than quantile regressions.The value of η 1 indicates whether variance in log R increases or decreases with SSB, and how much.Data points with low SSB were excluded as above.Minto et al. (2008) considered S-R curves corresponding to Eq. [5] with γ fixed at γ = −1, γ → 0, and γ = +1, but obtained very similar results.Therefore, only the case γ → 0 was investigated here.

Steepness
Steepness h is frequently used as a dimensionless index to characterize the shape of S-R curves.It is defined as the proportion of recruitment retained when SSB is reduced from the virgin biomass SSB 0 to 20% of this value.To obtain values of steepness comparable with published empirical results, we fitted the Beverton-Holt model to the simulated S-R data, assuming log-normal errors, and evaluated steepness as h = 0.2(1 + 0.2βSSB 0 )/(1 + βSSB 0 ).

Simulated S-R scatter plots
Simulations generated a wide variety of S-R relationships, of which the simple quantitative characterization discussed below do not capture all important features.Being unable to classify this variety of patterns automatically, we visually distinguished 10 different types according to characteristics explained in the caption of Fig. 1.The types were labeled from (a) to (j) in order of increasing visual evidence for compensatory recruitment.Figure 2 shows the prevalence of these types depending on asymptotic body mass M .
A large proportion (44%) of model S-R relationships had visual similarity with the Ricker model (types f, g, h, i.e. "Ricker", "blurred Ricker" or "Ricker-plus"), with a tendency for variability of recruitment at fixed SSB to increasing with M (f → g → h, Fig. 1).The types found for small (M < 1 g) model species, c, e, f, and g, exhibit only weak variability.This is likely an artifact of the simplified description of the spectrum of non-fish resources and the absence of competition for these in the model.The result is noteworthy nevertheless, as it indicates that the model's inherent oscillations of these species' predator populations only weakly affect recruitment.

Best-fitting Deriso-Schnute models
Of particular interest is the value of the parameter γ for fits of the Deriso-Schnute model to S-R curves, because this parameter determines the overall shape of the curves.Parameters α and β only scale them horizontally and vertically.As shown in Figure 3, we find a strong preferences for best-fitting values of γ near zero, confirming the visual impression of a tendency in our model for Ricker-type S-R relationships.In addition, however, the estimated distribution of γ has a pronounced shoulder around γ = −1, corresponding to the occasional occurrence of S-R relations better described by the Beverton-Holt model.Points represent simulation data, solid lines best-fitting curves, vertical dashed lines median virgin SSB.Both SSB and recruitment are normalized to the maximum found in simulations.Based on the visual impression, we labeled the types as follows, in approximate order of increasing evidence for non-linear compensation: (a) "cloud" when no clear structure was recognizable, (b) "depensatory" when recruitment increased faster than linear with SSB, (c) "linear" when recruitment rises nearly linearly with SSB, (d) "opening" when recruitment variability is strong and increases with SSB, (e) "maximizing" when recruitment has a local maximum near virgin SSB, (f) "Ricker" and (g) "blurred Ricker" for curves visually resembling a Ricker curve without and with additional variability, respectively (h) "Ricker-plus" for a variant of "Ricker" where variability abruptly becomes large for large SSB, (i) "curling" when SSB first increases, then decreases with added mortality, (j) "complex" when there is evidence for several regime shifts as mortality increases.

Heteroscedacity
While both Iles and Beverton (2000) and Minto et al. (2008) find indications for an increased variation of logarithmic recruitment with decreasing stock size (corresponding to η 1 < 0), this phenomenon is not reproduced by the present model.The distribution for the heteroscedastic coefficient η 1 found here is shown in Fig. 4. Next to a high preference for values of η 1 near 0, there is a strong contribution from large positive values.The latter correspond to an pronounced increased in the variability of log R for increasing SSB, as found, among others, for the "opening" S-R type (Fig. 1d).

Steepness
A histogram of the values of steepness h obtained from simulated S-R data is shown in Fig. 5a (histograms cope better with the constraint of possible h to the interval [0.2, 1] than kernel density estimation).Noteworthy is the significantly enhanced probability for values of h in the range [0.7, 0.8].As shown in Fig. 5b, S-R curves for stocks with smaller SSB, corresponding to smaller population densities, tend to exhibit lower steepness.
Steepness was independent of asymptotic body mass M in the model.The slope of the ordinary linear regression line of h against log 10 M was 0.016±0.018(s.e.), i.e., despite the small error indistinguishable from zero.

Coexistence
A first, basic finding of this study is that large model communities can indeed be built in which fully size-structured populations of fish interact and regulate their abundances exclusively through feeding interactions.The feasibility of community models of this type was hitherto unclear.The model originally proposed by Hartvig et al. (2011) was able to support at most nine co-existing populations.Variants of the model, employed by Blanchard et al. (2009), Andersen and Pedersen (2010) and Houle et al. (2012) use a hard-coded densitydependent fecundity or even model reproduction as a constant stream of larvae independent of spawner abundance (implying a strong negative dependence of fecundity on SSB).By building this additional density dependence into the FCSRM, the relative importance of population regulation and competition through trophic interactions is reduced, which explains the high species richness these variants can support.Other previous models of specious size-structured communities in which feeding was the only source of density dependence (Rossberg et al. 2008;Fung et al. 2013) did not resolve intraspecific size or age structure, and so were unable to represent the variable relationship between SSB and recruitment.OSMOSE (Shin andCury 2001, 2004) is one of a few other food-web models resolving intraspecific size structure where reproduction is modeled as a linear function of SSB; however, even OSMOSE contains additional sources of density dependence, e.g., a hard limit on the carrying capacity of non-piscivorous fish.

Stock-recruitment patterns
Perhaps not unexpectedly in view of the persistence of the unperturbed model communities, the simulated S-R data we extract from our model do indeed exhibit clear signs of densitydependent recruitment.Interestingly, density-independent recruitment variability is evident in the simulation data as well.This variability tends to be smaller than seen in empirical data, and this should not be surprising.The majority of densityindependent mechanisms known to control and modulate recruitment, such as climate-or oceanographic variability, are not represented in the model.Our results support the thesis that at least some of the observed recruitment variability at fixed SSB is attributable to variable abundances of the prey and predators of the focal species.
More interesting than the mere fact of density-dependent recruitment is the question how the structure of the simulated S-R data compares with structure identified in observation data.However, because of the relative brevity and enhanced variability of most empirical time series, it is difficult to decide which types of the variety of simulated S-R relations obtained in simulations (Fig. 1) have direct correspondences in nature.Yet, at least two overall patterns are in good agreement with observations.The first is the preference of the fitted Deriso-Schnute model for Ricker-and Beverton-Holt like S-R curves, with other conceivable shapes occurring less frequently (Fig. 3).Ricker-like curves (Fig. 1f,g,h) are found for species of all sizes (Fig. 2).Relations best described by Beverton-Holt-type models tend to be associated with the "opening" type (Fig. 1d) of S-R relations and arise predominantly for species with large maturation body size (Fig. 2).Empirically, the differences between Beverton-Holt and Ricker model are known to be difficult to resolve (Myers et al. 1994;Zhou 2007), and uncontrolled measurement errors in SSB can bias analyses towards flatter, Beverton-Holt type relationships (Walters and Ludwig 1981;Walters and Martell 2004).Whether the preference for Ricker that we find here is paralleled by observations or a model artifact is therefore difficult to say.Setting the distinction between these two models aside, the consistent preference for Ricker-or Beverton-Holt models in fisheries science itself can be taken as evidence that the other options offered by the c 2013 Crown copyright Deriso-Schnute model tend to yield less satisfactory descriptions of data, a finding that is consistent with our simulations.
The second pattern which we find to be in agreement with observations is the broad spread of steepness values, with a preference for steepness around 0.8.Shertzer and Conn (2012) demonstrated this preference in a meta-analysis covering 75 stocks of demersal fish, in good agreement with similar analyses by Myers et al. (2002).Contrary to the simulation results, Shertzer and Conn (2012) find no stocks with steepness in the range 0.2-0.3.The reason could be that the data set analyzed by Shertzer and Conn (2012) was inevitably biased towards wellstudied, rather abundant stocks, while steepness values near 0.2 are predominantly found for rare species in our simulations (Fig. 5b).Among the 28 most abundant socks, only one has steepness < 0.3.
This conformity with observations found for two quantitative characterizations is remarkable as no parameter in the model was adjusted for the purpose of the present study.Most parameter values were derived by Hartvig et al. (2011) from empirical data, except for the shape of the predator-prey sizeratio window, which was chosen following Houle et al. ( 2012)based on observations by Scharf et al. (2000), Pinnegar et al. (2003), andFloeter andTemming (2005)-and the parametrization of dietary preferences, Eq. [2], which follows from a quantitative analysis of fish-stomach data by Rossberg et al. (2011).
A third overall pattern, the broad variety of S-R relations found, conforms with observations on a qualitative level.There is broad agreement among fisheries scientists that a one-sizefits-all approach to S-R data is inappropriate, and our simulations support this view.Noteworthy are the two instances of the "depensatory" type (Fig. 1b) we find.The underlying mechanism must be very different from those usually invoked to explain this pattern, such as mate finding or shoal formation (Liermann and Hilborn 2001), because these mechanisms are not represented in the model.Remarkable are also patterns such as those seen in Figs.1h,j, which reflect qualitative transitions in stock dynamics with increasing added mortality.
Simulation results differ from the observation data analyzed by Iles and Beverton (2000) and Minto et al. (2008) in the predominant type of heteroscedacity.Partially, this might be explained by noting that empirical measurement errors in log R will be larger for small values of R, and so for small values of SSB.Neither Iles and Beverton (2000) nor Minto et al. (2008) had considered this explanation for their findings.Burrow et al. (2013) argue that data is insufficient to resolve heteroscedacity.Strong positive heteroscedacity (η 1 > 0) as seen in some simulated S-R relations could be masked by additional recruitment variability driven by effects not represent in the model.

Mechanisms
For both, the Beverton-Holt and the Ricker form of S-R relationships, a number of different mechanisms have been proposed.Some of these, such as contest or scramble competition, depend explicitly on spatial structure (Brännström and Sumpter 2005).However, when originally introduced, both forms were motivated by generic demographic models with mortality depending either on the size of the cohort of future recruits itself (Beverton and Holt 1957), on SSB, or number of eggs or hatchlings (Ricker 1954).The ecological mechanisms leading to these density dependencies remained unspecified.They may be purely trophic, as in the FCSRM, or involve other effects.Ricker (1954), for example, mentions density-dependent predation, parasitism, and density-dependent growth as potential regulating mechanisms.Ricker (1954) discusses cannibalism as the simplest mechanism, though probably relevant only for "a minority of populations".In Appendix A we analyze the regulating mechanism for a selected stock of "blurred Ricker" type (Fig. 1g), and find indeed no indication for regulation through density-dependent mortality (e.g., cannibalism).Only growth and reproduction exhibit clear density dependencies, demonstrating bottom-up regulation through food limitation.The strongest density-dependence is found for late juveniles and adults, compatible with the conclusions Lorenzen and Enberg (2002) and Lorenzen (2008) draw from survey data.For other types of S-R relations that we find, the regulating mechanisms are possible quite different, constrained only to derive from trophic interactions.
Our model does not resolve space, and so implicitly assumes a "well mixed" community.Because early life stages contribute only a small proportion to the total biomass of individuals in each size class (most individuals in each body-size class are adults, Rossberg 2012), the effect of early life stages on their environment is generally weak, and so the resulting direct densitydependent regulation of cohort size.This might explain the absence of clear Beverton-Holt-like patterns among the modeled S-R relations (Fig. 1).In reality there might be reasons that cohorts remain sufficiently localized to induce abundance regulation at early stages, despite the strong evolutionary pressure to disperse that this implies.
Another, independent argument for the dominant role of trophic interactions in determining abundances is the frequently observed power-law distributions of community biomass over body size, as represented by community size spectra (Sheldon et al. 1972).The FCSRM reproduces this phenomenon (Hartvig et al. 2011).The theories brought forward to explain it differ in detailed mechanisms they invoke (Rossberg 2012), but there is broad agreement that the observed coupling of abundances across size classes is achieved through trophic interactions, which implies a dominance of trophic density-dependencies over other mechanisms.To illustrate the relevance of these considerations for S-R models, consider the derivation of a Beverton-Holt S-R relationship from foraging-arena theory by Walters and Korman (1999).By this derivation, the characteristic scale of population abundances is essentially proportional to the mixing rate at which resources enter and leave the foraging arena, a parameter particular to this theory.Although several examples are known where the theory is likely to apply (Ahrens et al. 2012), its reconciliation with the phenomenon of power-law size spectra would require additional constraints, either on the mixing rate or on the number of co-existing species within a size class.

Implications
Our study demonstrates that the density-dependence of observed S-R relationships is possibly a consequence of densitydependent trophic interactions alone.Concepts such as foraging arenas (Walters and Juanes 1993;Walters and Martell 2004) or concentration effects (Iles and Beverton 2000), though biologically plausible, are not required to explain the observed density dependencies.A corollary of this possibility, if conc 2013 Crown copyright firmed, is that carrying capacity is not an intrinsic property of stocks, but the consequence of a complex interplay with their prey and predators.The notion that steepness is related to life-history traits of species (Myers et al. 1999) would need to be questioned.In the simulated communities steepness scatters broadly, even though life-history parameters differ among stocks only through differences in asymptotic size M , and this, on its own, turned out not to affect steepness.
If our hypothesis holds true, large changes in community structure are bound to lead to large changes in S-R relations, with corresponding implications for the maximum yield sustainable by stocks.Because in our model each stock is assigning its own resource size spectrum, the model might rather underestimate these interdependencies than overestimating them.Usage of population models with hard-coded S-R relationship in a management context should therefore remain mindful that these relations are likely to depend intrinsically on the structure of the entire community.Whether models with time-varying S-R parameters (Collie et al. 2012)

A. Mechanisms regulating abundance in a simulated stock with a Ricker-like S-R relationship
In this appendix we provide an analysis of the mechanism regulating abundance of a simulated stock with "blurred Ricker" S-R relationship (Fig. 2g).The analysis allows us to exclude cannibalism from the conceivable regulating mechanisms in this case.We selected for this example a stock for which the Ricker-type pattern was clear but overlayed by small additional variability (Fig. A1).

Methods
The growth rate of a population j is to a good approximation proportional to the (negative) deviation of the integral mortality, a quantity here defined as from its equilibrium value (see Eqs. ( 61), (63) of Rossberg 2012).In this expression, µ j (m) denotes the natural mortality rate of species j at size m, the constant µ a additional mortality applied here to vary stock sizes (µ a = 0 for only one species j at a time), and g j (m) and g r,j (m) are rates of growth and fecundity as explained in the main text.The expression in the integral over body sizes can intuitively be understood as the mortality at different body-sizes classes weighting by the duration of time spent in these size classes.
When µ a is slowly ramped up to reduce the size of population j in simulations, whilst maintaining a quasi steady state, there will be changes in size-dependent morality, growth rates and/or the fecundity of this population, compensating for the change in µ a so as to keep I j near its equilibrium value.Inspection of the integrand in Eq. [A1] as µ a varies therefore allows identifying the life stages and mechanisms playing the main role for this compensation.
To identify the relative roles of the regulation of growth and mortality at different life stages, we evaluated the natural integral mortality and the corresponding values obtained when replacing in Eq. [A2] either the natural mortality µ j (m) or the combined somatic and reproductive growth rate g j (m) + g r,j (m) by their averages µ j (m) and g j (m) + g r,j (m) over the entire ramp in µ a from Spawning Stock Biomass 0 µ a = 0 to µ lim , so isolating the effect of the other component, i.e. g j (m) + g r,j (m) or µ j (m)), respectively.
For numerical calculations, and also for an easy ecological interpretation of the functional form of the integrand, it is useful to express body mass m in Eq. [A2] in terms of a logarithmic size variable u as m(u) = exp(u), and to the change the variable of integration accordingly.This leads to the relation which is mathematically equivalent to Eq. [A2].
As µ a was varied, we recorded the value of the integral in Eq. [A3], its value when replacing either denominator or numerator by their temporal averages, and also the functional forms of the corresponding integrands, given by [m µ j (m)]/[g j (m)+ g r,j (m)] and [m µ j (m)]/[g j (m) + g r,j (m)].These quantities are closely related to physiological mortality rate (Beyer 1989) [m µ j (m)]/g j (m), the importance of which has been recognized by theorists (Cushing 1975;Werner and Gilliam 1984;Beyer 1989;Andersen and Beyer 2006;Rossberg 2012) and empiricists (Houde 1987;Houde and Zastrow 1993;Bohlin et al. 1994;Houde 1997;Nash et al. 2007;Houde 2008;Nash and Geffen 2012) alike.

Results and Discussion
Simulation outputs for integral mortality and the underlying size-dependent physiological mortality rates are shown in Figs.A2 and A3.
Two types of dynamics are overlayed in the time series of SSB recorded as added mortality µ a on this species is slowly increased (Fig. A2a).Variations in SSB on short time scales reflect variations in natural mortality (Fig. A2b) which affects the population growth rates mostly in the late juvenile and adult phase (Fig. A3a, dashed lines).Over longer time scales, SSB declines as a result of increasing added mortality, until at µ a = µ lim the species becomes extirpated.
Mortality (Fig. A2b) exhibits no notable density-dependent response.Changes in mortality seem to be largely driven by population cycles of consumers of the focal species that operate independent of it.By contrast, the effect of the densitydependent responses of growth and reproduction on natural integral mortality (Fig. A2b) is very clear.These responses, too, have their strongest impacts during the late juvenile and adult stages (Fig. A3b).The long-term trend in natural integral mortality (Fig. A2d) is largely explained by the changes due to density dependent growth and fecundity (Fig. A2c).The downward bending of the stock-recruitment curve in Fig. 1g is therefore largely due to a decline in adult growth and fecundity with increasing population size, obviously driven by a depletion of the food resources by the focal species.The bending down is here not attributable to changes in the mortality of larvae or juveniles, only the scatter of data seen for fixed SSB is.Cannibalism as a regulating mechanism can therefore be excluded in this example.
Figure A2e verifies the analytically motivated premise underlying the analysis above: despite substantial changes in the size-dependent life-history parameters, integral mortality, Eq. [A1], remain largely unchanged when including the effect of added mortality µ a .Changes in natural integral mortality, Eq. [A2], numerator of this expression by its temporal average.The ten lines correspond to ten moments in time evenly spaced as µa is slowly ramped up from 0 to µlim.In panel (a) lines show no particular temporal pattern (see Fig. A2b).In panel (b) lower-lying lines correspond to higher addition mortality µa, demonstrating a rise in growth and fecundity rates with decreasing population density.Dashed line show the same data for a 50-fold extended vertical scale (right side).The reason for choosing this particular representation of life-history parameters is that the area enclosed by two curves in a graph above is approximately proportional to hypothetical change in the corresponding natural populations growth rate if either growth and fecundity (panel a) or mortality (panel b) were kept fixed.To ease the biological interpretation of results, the thick curve above panel (a) indicates the maturation ogive.The focal species is that corresponding to Fig. A1.therefore reflect density-dependent responses to the artificial depletion by the added mortality µ a .
We caution that when conducting the kind of analysis demonstrated above for different stocks in the three simulated model communities studied in the main text, the graphs corresponding to Figs.A2 and A3 can attain quite different forms.Preliminary studies show that population regulation through densitydependent growth and fecundity in late juvenile and adult stages is frequent, but not the only pattern occurring in our simulations.Understanding the variety of patterns found and their implications for the mechanics of S-R relationship remain tasks for future studies.

Fig. 1 .
Fig. 1.Types of SR-relationships found in model simulations.Points represent simulation data, solid lines best-fitting curves, vertical dashed lines median virgin SSB.Both SSB and recruitment are normalized to the maximum found in simulations.Based on the visual impression, we labeled the types as follows, in approximate order of increasing evidence for non-linear compensation: (a) "cloud" when no clear structure was recognizable, (b) "depensatory" when recruitment increased faster than linear with SSB, (c) "linear" when recruitment rises nearly linearly with SSB, (d) "opening" when recruitment variability is strong and increases with SSB, (e) "maximizing" when recruitment has a local maximum near virgin SSB, (f) "Ricker" and (g) "blurred Ricker" for curves visually resembling a Ricker curve without and with additional variability, respectively (h) "Ricker-plus" for a variant of "Ricker" where variability abruptly becomes large for large SSB, (i) "curling" when SSB first increases, then decreases with added mortality, (j) "complex" when there is evidence for several regime shifts as mortality increases.

Fig. 2 .
Fig.2.Dependence of observed types of SR-curves on species size.Labels a-j refer to the types illustrated and defined in Fig.1.Some overlapping points (crosses) were slightly separated horizontally to improve visibility .

Fig. 4 .Fig. 5 .
Fig. 4. Probability distribution of the heteroscedastic coefficient η1 of Minto et al. (2008) estimated from simulation data.Positive η1 mean that the variance in log R found in S-R data increases with SSB.Density estimated by the adaptive kernel method followingPortnoy and Koenker (1989).

Fig
Fig. A1.Simulated S-R data for the stock analyzed in this appendix.Symbols used are the same as in Fig. 1 of the main text.

Fig
Fig. A2.Population regulation in response to external pressure.All curves are trajectories through time as added mortality µa is ramped up from 0 to µlim.Panel (a) shows the corresponding decline in SSB, (b) natural integral mortality, Eq. [A2], when fixing growth and fecundity, (c) natural integral mortality when fixing mortality, and (d) natural integral mortality when evaluating all rates at their momentary values.Panel (e) shows integral mortality including added mortality, Eq. [A1], demonstrating that this quantity remains largely unchanged, as theoretically expected for populations in equilibrium.The focal species is that corresponding to Fig. A1.
are the appropriate response to this challenge remains to be determined.