Niklas L P Lundström

Nicolas Loeuille

Xinzhu Meng

Mats Bodin

Åke Brännström

Meeting yield and conservation objectives by balancing harvesting of juveniles and adults

Keywords: fisheries management, maximum sustainable yield, pretty good yield, Pareto frontier, resilience, size-structure

Sustainable yields that are at least 80% of the maximum sustainable yield are sometimes referred to as pretty good yield (PGY). The range of PGY harvesting strategies is generally broad and thus leaves room to account for additional objectives besides high yield. Here, we analyze stage-dependent harvesting strategies that realize PGY with conservation as a second objective. We show that (1) PGY harvesting strategies can give large conservation benefits and (2) equal harvesting rates of juveniles and adults is often a good strategy. These conclusions are based on trade-off curves between yield and four measures of conservation that form in two established population models, one age-structured and one stage-structured model, when considering different harvesting rates of juveniles and adults. These conclusions hold for a broad range of parameter settings, though our investigation of robustness also reveals that (3) predictions of the age-structured model are more sensitive to variations in parameter values than those of the stage-structured model. Finally, we find that (4) measures of stability that are often quite difficult to assess in the field (e.g. basic reproduction ratio and resilience) are systematically negatively correlated with impacts on biomass and impact on size structure, so that these later quantities can provide integrative signals to detect possible collapses.

Introduction

Almost one third of the world's fished marine stocks are currently overexploited (FAO 2016). Some fish stocks have even collapsed, with examples including the Californian sardine (Sardinops sagax, Clupeidae) fishery in the 1950s [START_REF] Radovich | The collapse of the California sardine fishery -what have we learned?[END_REF], the Atlanto-Scandian herring (Clupea harengus, Clupeidae) fishery in the late 1960s [START_REF] Krovnin | Atlanto-Scandian herring: a case study[END_REF], the Peruvian anchovy (Engraulis ringens, Engraulidae) fishery in the 1970s [START_REF] Clark | The lessons of the Peruvian anchoveta fishery[END_REF], and the Northern cod (Gadus morhua, Gadidae) fishery off the east coast of Canada in the 1990s [START_REF] Hannesson | Fisheries Mismanagement: the case of the North Atlantic cod[END_REF][START_REF] Olsen | Maturation trends indicative of rapid evolution preceded the collapse of northern cod[END_REF]). The large proportion of overexploited marine fish stocks underscore the importance of implementing sustainable harvesting practices and for further improving modern fisheries-management methods.

Maximum sustainable yield (MSY) has long been a central concept in population ecology [START_REF] Smith | The gospel of maximum sustainable yield in fisheries management: birth, crucifixion and reincarnation[END_REF][START_REF] Hilborn | Defining success in fisheries and conflicts in objectives[END_REF][START_REF] Mesnil | The hesitant emergence of maximum sustainable yield (MSY) in fisheries policies in Europe[END_REF]. While maximization of yield from harvested populations is economically desirable, there is a rich scientific literature that criticizes the MSY concept and highlights its shortcomings, including the difficulty of correctly estimating MSY, the inappropriateness of long-term yield maximization as the single management objective, and the practical difficulty of accurately implementing the required level of harvesting effort [START_REF] Smith | The gospel of maximum sustainable yield in fisheries management: birth, crucifixion and reincarnation[END_REF]. MSY has further been criticized for its inability to prevent the collapse of important fisheries [START_REF] Beverton | On the dynamics of exploited fish populations[END_REF][START_REF] Larkin | An epitaph for the concept of maximum sustained yield[END_REF][START_REF] Mangel | Regime, phase and paradigm shifts: making community ecology the basic science for fisheries[END_REF][START_REF] Hilborn | Pretty good yield and exploited fishes[END_REF]. As an example, Alaska's Bering Sea Pollock fishery declined in 2009, and despite being known as a sustainable fishery which implements scientific recommendations, the management has been criticized for considering mainly MSY [START_REF] Morell | Can science keep Alaska's Bering sea pollock fishery healthy?[END_REF].

MacCall and Hilborn have introduced the concept pretty good yield (PGY; Hilborn 2010) for sustainable yields that are at least 80% of the MSY. In contrast to MSY harvesting-management objectives, PGY can be realized by a range of harvesting strategies and therefore leaves room to account for other desirable objectives in addition to the maximization of yield. The added value that PGY offers will depend on the extent to which the implemented harvesting strategies can successfully account for other desirable objectives beyond yield.

The aim of this paper is to investigate to which extent PGY harvesting strategies can simultaneously account for high yield and large conservation benefits. To increase the chances that our conclusions are valid over a broad range of circumstances, we base our study on two established population models. The first is an age-structured model (henceforth age model ) that is commonly used for modeling fish populations and evaluating fishing strategies [START_REF] Francis | Use of risk analysis to assess fishery management strategies: a case study using orange roughy (Hoplostethus atlanticus) on the Chatham Rise, New Zealand[END_REF][START_REF] Punt | Assessments of the stocks of Cape hake Merluccius spp. off South Africa[END_REF][START_REF] Punt | Stock assessment and risk analysis for the south Atlantic population of albacore Thunnus alalunga using an age-structured production model[END_REF][START_REF] Punt | Fisheries stock assessment and decision analysis: the Bayesian approach[END_REF][START_REF] Hilborn | Pretty good yield and exploited fishes[END_REF]. The second is a stagestructured consumer-resource population model (henceforth stage model ) that has been introduced by [START_REF] De Roos | Simplifying a physiologically structured population model to a stage-structured biomass model[END_REF]. Both models are capable of describing a range of aquatic and terrestrial animal populations. The age model belongs to a class of models that have a long history in fisheries science and that incorporates age-dependent fecundity, age-dependent survival, and density-dependent recruitment. The stage model is derived from a fully size-structured counterpart with food-dependent growth, fecundity and maturation, and accounts for feedbacks from resource depletion. In particular, it accounts for ontogenetic asymmetry, i.e., differential abilities of juveniles and adults to utilize available resources [START_REF] De Roos | Population and Community Ecology of Ontogenetic Development[END_REF]. With the age-model and stage-model being two fairly distinct representatives of contemporary population models, results on which they agree are likely to be fairly robust and results on which they differ are likely to be ones where careful description of the population ecology is important and may therefore differ from species to species. Using separate independently developed models to investigate robustness of findings should be a reasonable strategy [START_REF] Levins | The strategy of model building in population biology[END_REF].

We extend both models by introducing selective harvesting of juveniles and adults, giving wide ranges of possible harvesting strategies with different consequences for yield and conservation. While it is straightforward to quantify the yield of a harvesting strategy, it is less obvious how the conservations benefits should be measured. Here, we consider four different measures of conservation benefits: two measures that capture the direct impacts on the harvested population (the impact on population biomass and the impact on the population size structure) and two measures that capture the indirect risks of collapse due to changes in population dynamics (resilience and the basic reproduction ratio). We determine trade-offs between yield and conservation benefits by finding the so-called Pareto-efficient front; the set of strategies that cannot simultaneously be improved upon in both yield and conservation benefit. These trade-off curves allow us to assess how large conservation benefits can be gained while preserving PGY. Finally, we determine the relationship between the direct impact measures and the indirect risk measures, with the idea that the former are likely to be more easily observable in the field while the latter better reflect the risks of collapse. Taken together, our results show that there are large potential gains of using specific PGY harvesting strategies over traditional MSY strategies. Moreover, among PGY strategies, the ones that include equal harvesting of adults and juveniles often allow the best compromises between conservation and yield.

Methods

In this section we first present the two population models, one age-structured and one stagestructured. We extend both models by introducing selective harvesting of juveniles and adults, giving wide ranges of possible harvesting strategies with different consequences for the realized yield and for conservation. We next present our methods of stability analysis involving the impact measures and risk measures that we will use to evaluate different harvesting strategies. Finally, we recall the concept of maximum sustainable yield (MSY) and the economic concept of Pareto-efficiency which we will use to determine trade off curves between yield and conservation.

The age model

We adopt an age-structured population model that in different guises has been widely used when modeling fish populations and evaluating fishing strategies [START_REF] Francis | Use of risk analysis to assess fishery management strategies: a case study using orange roughy (Hoplostethus atlanticus) on the Chatham Rise, New Zealand[END_REF][START_REF] Punt | Assessments of the stocks of Cape hake Merluccius spp. off South Africa[END_REF][START_REF] Punt | Stock assessment and risk analysis for the south Atlantic population of albacore Thunnus alalunga using an age-structured production model[END_REF][START_REF] Punt | Fisheries stock assessment and decision analysis: the Bayesian approach[END_REF][START_REF] Hilborn | Pretty good yield and exploited fishes[END_REF]. The model incorporates densitydependent recruitment in the form of a Beverton-Holt stock-recruitment relationship with a tunable degree of random recruitment variability. Natural mortality is assumed to be independent of age and time, and age-specific harvesting is assumed constant over time. The central elements of the model, which are mainly derived from [START_REF] Hilborn | Pretty good yield and exploited fishes[END_REF], are described below.

We denote by N a,t the number of individuals of age a in year t and assume that individuals mature at age a mature after which they reproduce at an age-dependent rate proportional to their body size. Individuals younger than a mature are considered juveniles, while individuals older than or with age equal to a mature are considered adults. For simplicity, we assume that the population is made up entirely of female individuals, but as we show in the Appendix, our results are unchanged with a standard Fisherian sex-ratio of 50% females. The age-dependent fraction of mature females (m a ) and their corresponding egg production (f a ) are given by

m a = 0 if a ≤ a mature 1 if a > a mature f a = c s a , (1)
where c is a positive constant that scales the fecundity rate and s a is the mass of an individual at age a. We adopt a von [START_REF] Von Bertalanffy | Quantitative laws in metabolism and growth[END_REF] growth curve to describe individual length as a function of age, and assume that individual mass s a is proportional to the cube of individual length, i.e. s a = s max 1 -e -K(a-a 0 ) 3 ,

where s max is the asymptotic maximum body mass, K is a growth rate parameter and a 0 is a hypothetical negative age at which the individual has zero length. In Appendix Fig A 3 we illustrate von Bertalanffy growth curves for some parameter values. The total egg production in a year t = 0, 1, 2, . . . is found by summing over the offspring produced by mature females of different ages,

E t = amax a=0 m a f a N a,t , (3)
where a max is the maximum age of individuals and N a,t is the number of individuals, per unit of volume, of age a in year t.

The number of individuals in each age class changes from year to year according to

N 0,t = R t , and N a,t = N a-1,t-1 S(1 -γ a-1 ) for 1 ≤ a, (4)
where R t is the recruitment of newborn individuals in year t, t = 1, 2, 3, . . . , as described further below, and S(1 -γ a-1 ) is the probability that an individual survives from one year to the next. This survival probability is decomposed in survival from natural mortality S and from harvesting mortality (1 -γ a-1 ). Note that, as probabilities, these variables always take values from 0 to 1. We incorporate stage-selective harvesting by allowing separate constant fractions harvested of juveniles (F J ) and adults (F A ) and setting the vulnerability of individuals to

γ a = F J if a < a mature F A if a ≥ a mature . (5)
We assume Beverton-Holt recruitment [START_REF] Beverton | On the dynamics of exploited fish populations[END_REF],

R t+1 = E t α + βE t exp u t - σ 2 u 2 , ( 6
)
where u t are independent and normally distributed random variables with mean 0 and standard deviation σ u . The factor -σ 2 u /2 ensures that the expected number of recruits remains the same with varying σ 2 u because the lognormal random variable represented by the exponential has always mean 1. The parameters α and β are not used directly as they lack a direct ecological interpretation. Instead, they are determined from the expected total egg production at equilibrium in absence of harvesting (E 0 ) and from the steepness parameter (h) setting the sensitivity of the recruitment to the total egg production. The steepness is defined as the ratio of recruitment when egg production equals 20% of E 0 to recruitment at E 0 [START_REF] Mace | A generalised bioeconomic simulation model for fish population dynamics[END_REF][START_REF] Hilborn | Pretty good yield and exploited fishes[END_REF] and may take values between 0.2 and 1. If h is close to 1, recruitment is almost independent of the egg production, and if h is close to 0.2, recruitment is almost proportional to the egg production. In Appendix Fig. A 1 we illustrate how recruitment depends on E 0 and h and state their exact mathematical relationship with α and β.

For the age model, we use

R 0 = s max = c = 1, a max = 100, K = 0.23, a 0 = -2, a mature = 8, h = 0.7, S = 0.8, σ u = 0, (7)
as our default parameters values, with substantial motivations given in the Appendix. Here, R 0 is measured in number of individuals per unit of volume, s max and c -1 have an arbitrary mass unit, while a max , a 0 , a mature and K -1 are measured in years. However, after a rescaling of the equations, R 0 , s max and c (as well as s a , E t , R t and N a,t ) can be considered as non-dimensional and we can take R 0 = s max = c = 1 without loss of generality. See the Appendix for details.

While we believe that the parametrization in ( 7) is a reasonable choice, we have considered substantial variations and present how our results from the age model depend on parameter values in the result section. A systematic investigation of the robustness of our results with respect to variations of the parameter values in ( 7) are given in the Appendix.

The stage model

We adopt an archetypal consumer-resource model that has been introduced by de Roos et al. ( 2008) as a reliable approximation of a fully size-structured population model. The model is stage-structured and incorporates key aspects of individual life history such as fooddependent growth, maturation, and fecundity. In contrast to the age-model, the populationlevel feedback that results from resource depletion induces competition between life-history stages whenever juveniles and adults have differential abilities to utilize available resources.

Competition under such ontogenetic asymmetry can strongly influence the ecological dynamics (de Roos and Persson 2013) and may thus effect how harvesting affects population size structure. The central elements of this model are described below, with the detailed model formulation given in de [START_REF] De Roos | Simplifying a physiologically structured population model to a stage-structured biomass model[END_REF]. Individuals are composed into two stages, juveniles and adults, depending only on their size. Both juveniles and adults forage on a shared resource R = R(t). The juvenile biomass is denoted by J = J(t) while adult biomass is denoted by A = A(t). Juveniles are born with size s born and grow until they reach the size s max at which point they cease to grow, mature, and become adults. Juveniles use all available energy for growth and maturation, while adults do not grow and instead invest all their energy in reproduction. The juvenile growth rate and adult reproduction rate depend on resource abundance. In accordance with metabolic theory of ecology, foraging ability and metabolic requirements increase with individual body size [START_REF] Brown | Toward a metabolic theory of ecology[END_REF]). Juveniles and adults do not produce biomass when the energy intake is insufficient to cover maintenance requirements.

The rate at which the biomass of juveniles, adults, and available resources changes are given by three differential equations:

dJ dt = (w J (R) -v(w J (R)) -M -F J ) J + w A (R)A, dA dt = v(w J (R))J -(M + F A ) A, (8)
dR dt = r(R max -R) -I max R H + R (J + qA) .
Here, w J (R) and w A (R) are the net biomass production, per unit body mass, of juveniles and adults, respectively. The natural mortality is denoted by M while F J and F A are the respective stage-dependent harvesting rates of juveniles and adults. Continuing, v(w J (R)) is the resourcedependent rate at which juveniles mature and become adults, r is the resource turnover rate, and R max is the maximum resource density.

The net biomass production rates for juveniles and adults are assumed to equal the balance between ingestion and mass-specific metabolic rate T according to

w J (R) = max 0, σI max R H + R -T and w A (R) = max 0, σqI max R H + R -T .
Here, σ represents the efficiency of resource ingestion, and the maximum juvenile and adult ingestion rates per unit biomass equal I max and qI max , respectively, H is the half-saturation constant of consumers, and the factor q describes the difference in ingestion rates between juveniles and adults. The juvenile maturation rate depends on the net biomass ingestion and thus also on the resource density. It is derived from the fully size-structured counterpart by assuming that the population size structure is at equilibrium and by determining the rate at which juvenile individuals reach the maturation size s max . In the stage model, the juvenile maturation rate is given by

v(w J (R)) = w J (R) -M -F J 1 -(s born /s max ) 1-(M +F J )/w J (R)
The function v(w J (R)) results from a mathematical derivation and lacks a clear biologically interpretable form. When w J (R) = M + F J the function is undefined, and at this value it is defined by v(M + F J ) = -(M + F J )/ log(s born /s max ). In Appendix Fig.

A 7 we illustrate the maturation function v(w J (R)) as well as the net biomass production functions w J (R) and w A (R). Noting that the size at birth s born and the size at maturation s max appears only as the fraction s born /s max we can reduce the numbers of parameters by letting z = s born /s max .

For the stage model, we use

H = T = r = 1, R max = 2
, σ = 0.5, I max = 10, M = 0.1, z = 0.01, q = 0.85, (9) as our default parameters values. Here, H and R max are measured in biomass per unit of volume while T, r, I max and M are expressed per unit of time. However, after a rescaling of the equations all parameters, as well as the biomass densities J, A and R, can be considered as non-dimensional and we can take H = T = 1 without loss of generality. See the Appendix for details. The parameter values in (9) are inspired by [START_REF] De Roos | Simplifying a physiologically structured population model to a stage-structured biomass model[END_REF] and may be considered as archetypal. We show in the Appendix that our results are largely robust with respect to variations of these default values.

Model dynamics

Both the age model and the stage model are nonlinear dynamical systems and therefore complicated dynamics can not be ruled out a priori. However, extensive numerical investigations of basin of attractions indicate that solution trajectories end up, after sufficient time, at a globally stable equilibrium in both models. This equilibrium is therefore the only attractor which is either an interior (positive) equilibrium or an extinction equilibrium depending on the harvesting rates. Indeed, our nonlocal approach to resilience (presented below) tests the dynamics in both models for a large number of perturbations (inferred through initial conditions) and would detect any coexisting attractor with a high probability. We refer the reader to [START_REF] Meng | Dynamics and management of stage-structured fish stocks[END_REF] and [START_REF] De Roos | Simplifying a physiologically structured population model to a stage-structured biomass model[END_REF] for more on dynamic properties, mathematical analysis as well as numerical investigations of the stage model.

Stability analysis: measures of conservation

Stability of ecological systems is important for both conservation and harvesting purposes. In unstable systems, population dynamics may transiently go to low biomass values where the populations become vulnerable to demographic stochasticity or other factors. Hence, lack of stability promotes extinction. Stability is desirable also for harvest managers as it ensures stable yield. There are many definitions of stability, see e.g. [START_REF] Mccann | The diversity-stability debate[END_REF], and we propose here to study the consequences of harvesting on stability through four different measures of conservation. The first two are impact measures; impact of harvesting on the population biomass and impact of harvesting on the population size structure. The second two measures have a natural link to the risk of extinction and will be referred to as risk measures. These are the resilience and the basic reproduction ratio (the recovery potential in case of the stage model). We define the resilience as the reciprocal of the time needed for the population to recover from a perturbation, and we consider here both small and large perturbations. The basic reproduction ratio/recovery potential describes the population's rate of increase from very low abundances, and thus could be construed as the likelihood of population rebound, following a crash (e.g. due to a large disturbance).

Measures of impact on biomass and size structure

Let J * = J * (F J , F A ) and A * = A * (F J , F A ) denote the juvenile and adult biomass at equilibrium, respectively, of the harvested population. In case of the stage model J * and A * are given directly by the state variables at equilibrium. For the age model, J * and A * are obtained through the formulas

J * = amature-1 a=0 N * a s a and A * = amax a=amature N * a s a ,
where N * a denotes the number of individuals of age a at equilibrium. Furthermore, let J * u and A * u denote the juvenile and adult biomass at equilibrium in the absence of harvesting, i.e. J * u = J * (0, 0) and A * u = A * (0, 0). Moreover, let B * = J * + A * and B * u = J * u + A * u . We measure impact on biomass of harvesting through the expression

Impact on biomass = 1 - B * B * u .
Similarly, we consider impact on size-structure through the expression Impact on size-structure =

J * J * + A * J * u J * u + A * u -1
-1, which equals the relative change in the fraction of juvenile biomass following harvesting. If the impact on size-structure is positive (negative), then harvesting has increased (decreased) the fraction of juveniles in the population.

Resilience, basic reproduction ratio and recovery potential as risk measures

Resilience as a risk measure is increasingly used in ecology [START_REF] Pimm | Number of Trophic Levels in Ecological Communities[END_REF][START_REF] Loreau | Phenotypic diversity and stability of ecosystem processes[END_REF][START_REF] Petchey | Species richness, environmental fluctuations, and temporal change in total community biomass[END_REF][START_REF] Montoya | Ecological networks and their fragility[END_REF][START_REF] Loeuille | Influence of evolution on the stability of ecological communities[END_REF][START_REF] Valdovinos | Consequences of adaptive behaviour for the structure and dynamics of food webs[END_REF]. Resilience is now also increasingly discussed in a fishery management context [START_REF] Hsieh | Fishing elevates variability in the abundance of exploited species[END_REF][START_REF] Law | On balanced exploitation of marine ecosystems: results from dynamic size spectra[END_REF][START_REF] Fung | Why the size structure of marine communities can require decades to recover from fishing[END_REF]. The higher the resilience, the smaller the risk of extinction due to random drift. We consider resilience of the population by measuring the reciprocal of the time needed for the population to recover the positive equilibrium given a random perturbation. We do this by considering a large number of initial conditions. From each initial condition, we measure the time until the population (and also the resource in case of the stage model) returns to a small neighborhood of the equilibrium. The average value of this return time over the number of trials are then used to quantify the resilience: Resilience = 1 Average value of the return times .

Our resilience measure estimates the population's expected rate of return, given a random perturbation. In contrast to many other studies on resilience that assess resilience based on eigenvalues of the Jacobian matrix, our approach is not limited to the immediate neighborhood of the equilibrium, but can also tackle large disturbances, a point we will return to in the discussion section. The precise procedure by which we determine the resilience is described in the Appendix where we also present an alternative resilience measure, estimating the population's probability to return within a time limit, and reproduce some of our results using different magnitudes of disturbances.

We also consider the basic reproduction ratio as a risk measure, which represents the average number of offspring produced over the lifetime of an individual in the absence of densitydependent competition, i.e., when the population abundance is very low. For the age model, we derive the following expression for the basic reproduction ratio as functions of the harvesting rates F J and F A :

Basic reproduction ratio = (1 -F J ) amature × amax a=amature+1 s a S a (1 -F A ) a-amature 1 -h-0.2 0.8h × amax a=amature+1 s a S a . ( 10
)
A basic reproduction ratio larger than one ensures that the biomass of an initially small population increases on average, while a basic reproduction ratio less than one implies that the population will eventually become extinct. The derivation of expression (10) can be found in the Appendix.

In case of the stage model we use the recovery potential introduced in [START_REF] Meng | Dynamics and management of stage-structured fish stocks[END_REF],

Recovery potential = w A (R max ) M + F A × v(w J (R max )) v(w J (R max )) -w J (R max ) + M + F J .
The recovery potential is the generational net biomass production (per unit body mass) in a pristine environment (free from density-dependent mortality) and is therefore closely related to the basic reproduction ratio. Similar as for the basic reproduction ratio, a recovery potential larger than one ensures that the biomass of an initially small population increases on average, while a recovery potential less than one implies that the population will eventually become extinct. The basic reproduction ratio, as well as the recovery potential, are directly linked to the probability of surviving a period of low population abundance during which random drift caused by demographic stochasticity can lead to extinction. We further discuss this fact, as well as giving overall motivations of our choices of conservation measures, in the discussion section.

Maximum sustainable yield and trade-off through Pareto efficiency

Recalling that J * and A * denote the juvenile and adult biomass at equilibrium, for any given harvesting rates F J ≥ 0, F A ≥ 0, the yield objective function is given by

Yield = F J J * + F A A * . (11)
Moreover, the maximum sustainable yield (MSY) is obtained by taking the maximum of the yield objective function across all harvesting strategies (F J , F A ). In addition to the yield function we are, in case of both the age model and the stage model, armed with four measures of conservation as functions of the harvesting rates (F J , F A ). Using these objective functions we can calculate both the yield and the conservation for given harvesting strategies, see Fig. 1 in the Results section.

To determine the trade-off between the two objectives yield and conservation, we plot the yield as a function of each conservation measure in the results section and apply the economic concept of Pareto efficiency to evaluate different harvesting strategies. A harvesting strategy is Pareto efficient if it cannot be improved upon without trading off one of the considered objectives against the other, see e.g. Karpagam (1999, page 11). The Pareto front is the set of all Pareto efficient harvesting strategies. Hence, managers can restrict the choice of harvesting strategy to this set, rather than considering the full range of possible harvesting strategies. The closer a strategy is to the Pareto front, the more efficient it is.

Results

Figure 1 shows how the four measures of conservation and the yield changes with harvesting intensity for equal harvesting rates of juveniles and adults (henceforth equal harvesting), i.e. F J = F A . As harvesting pressure increases, the yield first increases after which it decreases as the population becomes "overexploited". The impact on biomass and the impact on sizestructure increase with harvesting pressure, while the basic reproduction ratio and the recovery potential decrease. The resilience decreases with harvest pressure in case of the age model, but first increases to a maximum and then decreases in case of the stage model. Note that due to the different nature of the age model and the stage model, the values of the harvesting rates in the two models may not be immediately compared.

We are now ready to present the trade-offs between yield and the four measures of conservation. Figure 2 represents results from the age model, while Fig. 3 gives the corresponding results for the stage model.

Pretty good yield allows large conservation benefits

Focusing on the age model, Figs. 2 (b)-(d) show that the basic reproduction ratio is relatively low and the impact on size structure is also relatively large at MSY, while the resilience is relatively high at MSY. Focusing on the stage model, Figs. 3 (c) and (d) show slightly different results; harvesting for MSY (which is obtained by harvesting only adults) gives a resilience and a recovery potential that is only a tiny fraction of the unexploited state and is close to the boundary of extinction. Hence, harvesting for MSY may substantially increase the risk of stock collapse. Fig. 3 (b) shows also that the impact on size structure is at a maximum at MSY.

Common for both models and all four measures is, see Figs.

Stage model

Figure 1: The four measures of conservation and yield as functions of the harvesting rates considering equal harvesting (F J = F A ) in case of the age model and the stage model. Yield (black, solid), impact on biomass (blue, dash-dot), impact on size-structure (green, dotted), resilience (red, solid) and basic reproduction ratio (recovery potential) (yellow, dashed). To visualize all objectives in the same plot graphs show (left) 14.5×Yield, 50×Resilience, 0.1×Basic reproduction ratio and (right) 3.3×Yield, 5×Resilience, 0.13 × log(Recovery potential).

with nearly half the impacts on biomass and half the impact on size structure, and also with nearly twice the basic reproduction ratio (age model) and a much higher recovery potential (stage model). Resilience can also be improved in case of both models, though the difference in resilience is most impressive for the stage model, see Figs. 3 (c). Hence, both the age model and stage model give the result that PGY allows for large conservation benefits. Varying parameter values show that this conclusion is robust in both models (see Appendix). However, stepping back in yield into the range of PGY does not automatically ensure conservation in terms of any of the measures we consider. To exclude the non-optimal harvesting strategies and to find the best ones within the range of PGY, we apply the economic concept of Pareto efficiency, as introduced in the previous section. Following the Pareto front (the set of all Pareto efficient harvesting strategies shown as the green curves in Figs. 2 and3) reveals these preferable harvesting strategies. In the following, we will discuss simple harvesting strategies which are relatively close to the Pareto fronts in all cases.

Equal harvesting rates on juveniles and adults is often a good strategy Figs. 2 and3 show that equal harvesting (F J = F A ), performs well with respect to both models and all four measures of conservation. In particular, in the range of PGY, the black curves come rather close to the Pareto front in all subfigures (especially in a neighborhood of the black dots). Therefore, we can harvest juveniles and adults at equal rates, which should be strategies that are rather easy to implement, without losing too much yield or conservation. The black dots in Figs. 2 and3 show one such strategy. Indeed, harvesting only adults is costly on some aspects, particulary in terms of resilience (stage model) and the impact on size structure as well as basic reproduction/recovery potential (both models). Harvesting only juveniles is costly in terms of resilience (both models) and in terms of impact on biomass (age model).

Varying parameter values show that this conclusion is very robust in the stage model, where it seems to remain in the wide ranges .5, 100]. We refer the reader to the Appendix for substantial investigations (both numerical and analytical) of robustness with respect to variations of parameter values. Additional trade-off curves, as those presented in Fig. 3, are given for six different parametrizations in Figs. A 8, A 9 and A 10. However, equal harvesting is often but not always suggested by the age model. Here, the efficient strategies seem to depend on the fraction of juveniles at the unharvested equilibrium, J * u /(A * u + J * u ), as well as on the survival from natural mortality, S. We proceed by investigating this dependence by comparing pure adult harvesting (henceforth adult harvesting), equal harvesting and pure juvenile harvesting (henceforth juvenile harvesting) for a wide range of parameter values in the age model. Figure 4 gives an approximation of regions in which the age model suggest adult harvesting, equal harvesting and juvenile harvesting. Juvenile or adult harvesting is suggested only if such strategies are the most Pareto-efficient once, within the range of PGY, with respect to all four conservation measures. The borders in Figure 4 are approximations which are produced by examining a large number of variants of Figure 2 for parameter values in the intervals

M ∈ [0, 0.5], z ∈ [0.0001, 0.2], σI max ∈ [3, 100], q ∈ [0.6, 2], R max ∈ [0
a 0 ∈ [-3, -0.2], K ∈ [0.1, 1], a mature ∈ [3, 15], h ∈ [0.3, 0.9], σ u ∈ [0, 0.5].
Indeed, we varied each parameter at a time, keeping the others at the default values given in ( 7), and tested at least 10 values in each interval. Further parameter combinations have also been tested in order to refine the borders in Figure 4. Points P 0 -P 8 in Fig. 4 correspond to different parametrizations of the age model. The default parametrization in (7) gives J * u /(A * u + J * u ) ≈ 0.6, S = 0.8 and is marked with P 0 . In Appendix Figs A 2, A 4, A 5 and A 6 we present trade-off curves, similar to those in Fig. 2, for different parametrizations corresponding to the remaining eight points P 1 -P 8 . The Appendix also contains motivations and explanations for the dependence shown in Fig. 4.

It turns out that if it is possible to obtain PGY for a wide range of harvesting strategies (including adult, juvenile and equal harvesting), then our conservation measures are in favor of equal harvesting. When adult or juvenile harvesting performs better than equal harvesting, it is usually because equal harvesting can not give a yield in the range of PGY.

The age model is more sensitive to variations in parameter values than the stage model

Focusing on the age model we first note that for the parameter values used in Figs. 1 and2 we have a survival from natural mortality of S = 0.8 [START_REF] Mills | Abundance, annual survival, and recruitment of unexploited and exploited lake charr, Salvelinus namaycush, populations at the Experimental Lakes Area, northwestern Ontario[END_REF] and the fraction of juveniles in an unharvested population, J * u /(A * u +J * u ) ≈ 0.6. We conclude that in this case equal harvesting is a good strategy. Varying the parameter values, it turns out that an increase in the fraction of juveniles implies an increase in the yield obtained when harvesting only juveniles, i.e. the blue curves will be lifted in Fig. 2. Similarly, a decrease in the fraction of juveniles implies an increase in the yield obtained when harvesting only adults, i.e. the red curves will be lifted in Fig. 2. This dependence, which is expected and natural, can be observed in both models, but it is much stronger in the age model. While Fig. 4 gives an approximation of the borders between adult harvesting, equal harvesting and juvenile harvesting, a similar investigation on the stage model gives a much larger region suggesting equal harvesting. In particular, in the stage model, the most Pareto-efficient strategies, within the range of PGY, seems to be dominated by equal harvesting as long as 0.1 < J * u /(A * u + J * u ) < 0.9. 4: The harvesting strategy suggested by the results of the age model depends on S and the fraction of juveniles in the unharvested population. In the dark-grey region, equal harvesting is suggested by the age model, while adult harvesting is better for low fraction of juveniles and juvenile harvesting is to recommend when the fraction of juveniles is high. The light-grey regions refer to borderline cases. Juvenile or adult harvesting is suggested only if such strategies are the most Pareto-efficient once, within the range of PGY, with respect to all four conservation measures. Point P 0 corresponds to the default parametrization, while points P 1 -P 8 correspond to parametrizations considered in the Appendix.

the age model and the stage model agree on similar results. For populations outside of this region the age model suggests adult harvesting, or, for some rare parameter settings, juvenile harvesting.

Impact on size structure and impact on biomass serve as warning signals

As neither the resilience nor the basic reproduction ratio (recovery potential) can be directly measured in the field, it is important to identify reliable proxies for conservation management that can be measured in field surveys. Figs 5 and6 show that a harvesting strategy with a high impact on population size structure, or a high impact on biomass, implies a low basic reproduction ratio (recovery potential) and a low resilience and hence a high risk of collapse. Indeed, we find that resilience and basic reproduction ratio (recovery potential) are systematically negatively correlated with impacts on biomass and size structure, so that these later quantities, which should be relatively easy to measure in field surveys, can provide integrative signals to detect possible collapses.

Discussion

We have investigated how well stage-dependent harvesting strategies that qualify for pretty good yield (PGY) can account for conservation as a second objective. To increase the chances that our results apply to a broad range of populations, we have studied two established population models and reported conclusions that are common to both. We have also investigated a wide range of parameter values for both models. To incorporate conservation as a second objective for our optimization procedure, we have used four different measures of conservation applied to both the age model and the stage model. First, this extended analysis allows us to conclude strong robustness of the results when all measures agree for both models; e.g., that there are large potential gains of using specific PGY harvesting strategies that often, but not always, correspond to equal harvesting rates of juveniles and adults. Second, we are able to discuss and compare both the two models as well as the four measures of conservation with each other.

Implications for management of harvested populations

Our study supports the implementation of PGY. Furthermore, our results support implementation through equal harvesting of juveniles and adults, in conjunction with regular surveys that aim to detect changes in population biomass and size structure. Managers aiming to implement optimal regulations may want to parameterize the age and stage model (or other suitable population models) for the specific species in question. A similar analysis as the one presented here can then be carried out and will give the specific harvesting strategy that maximizes conservation benefits, e.g. as described by the four conservation measures considered here, for a given target-value of sustainable yield.

Managers relying on other approaches may still be interested in assessing changes in the size structure of a population, as well as changes in biomass, as these are strongly linked to our risk measures and may thus serve as warning signals for an impending collapse. In fisheries management, changes in size structure and biomass can be measured through trial fishing, reinforcing our conclusion that size-structure and biomass are appropriate proxies for the risk of collapse and possible extinction.

Stage model

Figure 6: The relation between impact measures and risk measures for the stage model. We observe that a large impact on biomass implies a low basic reproduction ratio (recovery potential) and also a low resilience. The same is true for impact on size structure. Curves, green squares and parameters are as in Fig. 3.

One should note that equal harvesting should be relatively easy to implement. Indeed, in a single-species setting an equal harvesting strategy is implemented by setting the same harvesting rate over all sizes of individuals. This joint rate should then be tuned against the smallest value giving the desired yield. In a multi-species setting, the harvesting rate should still be the same over all ages/sizes within each species, but the preferable rate may differ between species. Naturally, the rate should be higher for species with a higher productivity, and lower for species having lower productivity. Related to this is the concept of balanced harvesting, which has attracted considerable attention recently, and aims to distribute "a moderate mortality from fishing across the widest possible range of species, stocks, and sizes in an ecosystem, in proportion to their natural productivity, so that the relative size and species composition is maintained" [START_REF] Garcia | Reconsidering the consequences of selective fisheries[END_REF]. While implementing balanced harvesting is difficult since such strategy may be selective within each species as well (since productivity may depend on age/size), see e.g. [START_REF] Reid | Implementing balanced harvesting: practical challenges and other implications[END_REF], our results show that size-structure can be preserved fairly well by implementing the simple strategy of equal harvesting.

We have optimized for yield and conservation, not for economic yield. Therefore, depending on the market (price of small fish versus price of large fish), managers may obtain different preferable harvesting strategies if aiming for economic yield.

Why harvest juveniles? Differences and similarities between the age model and the stage model

While harvesting individuals before they mature is a debated topic, we have seen that both the age model and the stage model give arguments for equal harvesting rates of juveniles and adults. Indeed, relying on the stage model this argument is robust with respect to variations in parameters values. The age model is more sensitive and the suggested harvesting strategy varies between mainly equal harvesting and adult harvesting as a function of parameter values.

To understand these results we first recall (see Results) that if it is possible to obtain PGY for a wide range of harvesting strategies, then our conservation measures are in favor of equal harvesting. Therefore, we can focus the following discussion on when and why the models allow for such wide range of harvesting strategies.

By extensive numerical experiments we illustrate this dependence for the age model in Figure 4. Varying the parameters values in the age model, it turns out that an increase in the fraction J * u /(A * u +J * u ) of juveniles implies an increase in the yield obtained when harvesting only juveniles, and that a decrease in the fraction of juveniles implies an increase in the yield obtained when harvesting only adults. This is natural; when considering harvesting of a population consisting of mainly juveniles it is not possible to obtain a good yield by harvesting only adults. From Fig. 4 we also see that as the survival from natural mortality S increases, the recommendation goes towards including more juveniles in the harvesting strategy. A reason for this is that for small S additional mortality through harvesting on young individuals implies that too few individuals survive and become adults and the population declines.

A corresponding parameter dependence, as described in Fig. 4, is much weaker in case of the stage model. Indeed, as mentioned in the results section, the stage model suggest equal harvesting for wide ranges of parameter values, see the Appendix for more on this. One reason for this difference in sensitivity of the recommended harvesting strategies between the two models, with respect to parameter values, is as follows. The stage model explicitly models the resource R(t) through the third equation in (8), and reproduction, growth and maturity are assumed to be increasing functions of the resource. Therefore, removing adult or juvenile biomass through harvesting results in more resource available for the remaining population, which in turn increases biomass production through all three mechanisms reproduction, growth and maturity of juveniles. This feedback implies that the dynamics of the stage model allows for wide ranges of efficient harvesting strategies.

On the other hand, the age model incorporates the Beverton-Holt spawner recruit curve in (6) for reproduction, and, independent of the recruitment, individuals are assumed to grow following the Bertalanffy growth curve in (2). Growth and recruitment are thus assumed to be independent in the age model, while they are dependent through the resource in the stage model. This means that if some juveniles are removed by harvesting it will not be in favor of the recruitment of newborns in case of the age model, as it would be in case of the stage model. Thus, it is more costly to harvest juveniles in the age model than in the stage model, and therefore the age model more often suggests to leave small individuals, let them grow, and catch them as adults.

In conclusion, the more extensive population-level feedbacks in the stage model makes the population productive for a wider range of harvesting strategies than the age model does, and the age model is more restrictive to juvenile harvesting than the stage model. This explains why equal harvesting performs well through wider ranges of parameter values in the stage model, than in the age model.

Importance of preserving population size structure

Our advice is based on our finding that large impacts on size-structure generally implies a high risk of collapse as captured by our risk measures, see Fig. 5 and 6. To reduce the impact of harvesting on population size structure, it seems advisable to harvest juveniles as well as adults, see Fig. 2 and3. Thus, equal harvesting is more likely to preserve the size structure than single-stage harvesting. (A similar conclusion was reached by [START_REF] Jacobsen | The consequences of balanced harvesting of fish communities[END_REF] We have shown that large impacts on size structure typically indicate unfavorable readings of our risk measures. Our work thus reinforces the conclusions from a large and growing number of studies (considering both ecological and evolutionary aspects) that argue for the importance of preserving the size structure of harvested populations. These studies, which we discuss below, reinforce the importance of including impact on size structure explicitly as an important conservation measure when discussing harvesting strategies. In fact, not accounting for the impact on size structure explicitly in our analysis means that we should find the recommended harvesting strategies from Figs. 2 and3 From an ecological point of view, our analysis that size structure largely impacts the structure and functioning of the system is in agreement with previous works. [START_REF] Anderson | Why fishing magnifies fluctuations in fish abundance[END_REF] show that populations with a larger fraction of juveniles have less stable population dynamics because of changes in demographic parameters, and, therefore, suffers a larger risk of extinction. (In this context, see also [START_REF] Wikström | The role of harvesting in age-structured populations: Disentangling dynamic and age truncation effects[END_REF] Moreover, changes in the size-structure of the population affect the balance of intra-and interspecific competition [START_REF] Loreau | Biodiversity and ecosystem stability: a synthesis of underlying mechanisms[END_REF].

From an evolutionary point of view, affecting the size-structure of a population can potentially induce changes in biological traits such as size-at-age and age-at-maturation. One reason is that harvesting only large individuals creates a large mortality selective pressure so that only adults that reproduce early, and at small size, pass their genes [START_REF] Dunlop | From genes to populations: how fisheriesinduced evolution alters stock productivity[END_REF]. A conse-quence may be evolution toward small individuals reproducing early, which is generally not desirable from an ecological (low reproduction) nor from an economical (too small to be valuable) point of view [START_REF] Grift | Fisheries-induced trends in reaction norms for maturation in North Sea plaice[END_REF][START_REF] Olsen | Maturation trends indicative of rapid evolution preceded the collapse of northern cod[END_REF]. Evolutionary changes may have large impact on economic profit and future management [START_REF] Conover | Sustaining fisheries yields over evolutionary time scales[END_REF][START_REF] Jorgensen | Managing evolving fish stocks[END_REF][START_REF] Belgrano | How fisheries affect evolution[END_REF], and may be difficult to reverse. Studying the collapse of the northern cod (Gadus morhua, Gadidae), it has been shown that, before government imposed a moratorium, the life history shifted towards maturation at earlier ages and at smaller sizes [START_REF] Olsen | Maturation trends indicative of rapid evolution preceded the collapse of northern cod[END_REF]), suggesting fisheries-induced evolution of maturation patterns. Moreover, a recent study provides experimental evidence for rapid evolution induced by changes in the population size-structure of a fished population [START_REF] Van Wijk | Experimental harvesting of fish populations drives genetically based shifts in body size and maturation[END_REF]). Significant genetic variation for production-related traits is also present in fished populations [START_REF] Law | Fishing, selection, and phenotypic evolution[END_REF], and [START_REF] Cameron | Eco-evolutionary dynamics in response to selection on life-history[END_REF] experimentally demonstrate evolutionary changes, in response to harvesting juveniles or adults. In [START_REF] Kuparinen | Detecting and managing fisheries-induced evolution[END_REF] the authors argue that we should stop targeting only large individuals to avoid evolutionary impact on fisheries. See also [START_REF] Garcia | Reconsidering the consequences of selective fisheries[END_REF] and [START_REF] Law | On balanced exploitation of marine ecosystems: results from dynamic size spectra[END_REF] for more arguments for harvesting preserving population size structure.

In conclusion, we recommend that managers consider the impact on size structure and that they avoid large deviations from the size structure of a pristine, unharvested population.

Relations between the four measures of conservation

We have considered conservation as a second objective, beyond yield, in our optimization procedure. To quantify conservation we have chosen two impact measures; impact on biomass and impact on size structure, and two "risk" measures; resilience and basic reproduction ratio (recovery potential). It is not obviously true, even though it is expected, that the impact measures relate simply to the risk measures. Therefore, we present Figs 5 and 6 which show that a large impact on biomass, or size structure, implies a low basic reproduction ratio (recovery potential) and also a low resilience. From this fact we concluded that it is important to preserve size structure and biomass in order to preserve stability of the population, and that impact on biomass and impact on size structure work as warning signals for a collapse.

From Figs. 5 and 6 it is clear that the relations between the four measures of conservation are not simple. This is clearest from Fig. 6 (a) and (c) representing the stage model; resilience can vastly differ from the other measures by increasing with harvest pressure for some harvesting strategies. This phenomena, which bears resemblances to the paradox of enrichment [START_REF] Rosenzweig | Paradox of enrichment: destabilization of exploitation ecosystems in ecological time[END_REF][START_REF] Rip | Cross-ecosystem differences in stability and the principle of energy flux[END_REF], deserves attention since it is very strong in the stage model under equal harvesting, but not present at all under adult harvesting. This thus provides a substantial argument in favour of equal harvesting. Using an alternative resilience measure, we give further illustrations and explanations of this behaviour of the stage model in the Appendix, see Figs. A 11 and A 12.

Comparing resilience simulations in Figs. 2 (c) and 3 (c) we conclude that harvesting of only adults is among the best strategies in the age model, while such strategy is among the worst using the stage model (considering resilience only). The resilience in the stage model instead suggest equal harvesting rates. This is because in the age model, the population will return fast to the equilibrium when harvesting only adults, after a given perturbation. In the stage model however, the population returns very slowly when harvesting only adults, compared to the case of equal harvesting. Hence, transient behavior, and therefore the resilience, behaves different in the models.

Concepts of stability are numerous in ecology [START_REF] Mccann | The diversity-stability debate[END_REF], and how the different stability measures relate to one another is considered a timely and important question [START_REF] Donohue | On the dimensionality of ecological stability[END_REF]. Indeed, measures of stability in ecological systems is today an active research area, see e.g. [START_REF] Neubert | Alternatives to resilience for measuring the responses of ecological systems to perturbations[END_REF] for alternatives to resilience, [START_REF] Nimmo | Vive la résistance: reviving resistance for 21st century conservation[END_REF] for discussions of resistance and resilience, and [START_REF] Isbell | Biodiversity increases the resistance of ecosystem productivity to climate extremes[END_REF] as well as [START_REF] Dunne | Network structure and biodiversity loss in food webs: robustness increases with connectance[END_REF] for stability and its relations to biodiversity. Importantly, our model suggests that our different conservation measures covary, and may be usefully assessed through changes in biomass and size structure.

Motivations of our choice of conservation measures

As our results may depend on our chosen conservation measures, we consider here additional motivations and discussions concerning this topic. First, the measures impact on biomass and impact on size structure are important to consider simply since they can be measured in reality.

Second, these measures are natural, simple and easy to interpret and a large impact on biomass would certainly imply impact on the surrounding ecosystem. Moreover, in the subsection Importance of preserving population size structure, we further motivated the impact on size structure as a central measure, based on the fact that population size structure is important to preserve from both an ecological and evolutionary point of view.

To motivate the basic reproduction ratio and the recovery potential as a risk measure, we note that, as already mentioned in the methods section, these measures are directly linked to the probability of surviving a period of low population abundance during which random drift caused by demographic stochasticity can lead to extinction. To see this, consider a small population in a pristine environment in which all individuals are, for simplicity, assumed to be identical. In this case, the basic reproduction ratio (or the recovery potential) is simply the ratio between birth-rate b and death-rate d, that is Θ = b/d. As proved by Grimmett and Stirzaker (1992, page 272), the probability of avoiding extinction through random drift is given by 1 -1/Θ if Θ > 1, and zero if Θ ≤ 1. Hence, there is a direct link between the basic reproduction ratio (recovery potential) and the probability of surviving a period of low population abundance; a high basic reproduction ratio (recovery potential) ensures a high probability of surviving. To justify the investigation of the effects of large disturbances that bring the population to small numbers, so that density dependence can be ignored, placing individuals in a pristine environment, we mention mass mortality events [START_REF] Fey | Recent shifts in the occurrence, cause, and magnitude of animal mass mortality events[END_REF], drastic climate variability such as heat waves, storms, and floods [START_REF] Reusch | Ecosystem recovery after climatic extremes enhanced by genotypic diversity[END_REF], and heavily exploited ecosystems [START_REF] Jones | Rapid recovery of damaged ecosystems[END_REF].

To motivate our choice of resilience as a risk measure we first note that, when dealing with nonlinear models, many works considers only local stability and local resilience measures (based on eigenvalues of the Jacobian matrix). However, such approach gives only information arbitrarily close to the equilibrium, saying little about the basin of attraction (the set of initial conditions attracted by the equilibrium). If the equilibrium is locally stable but the basin of attraction is small, then even a small perturbation can force the dynamics to jump to another attractor, having possibly dangerous behaviour. A large and convex basin of attraction, with the equilibrium in the middle, ensures that the population recovers a perturbation with a high probability. Therefore, it is natural to use both the size and the shape of the basin of attraction as stability/risk measures (Lundström and Aidanpää 2007[START_REF] Menck | How basin stability complements the linear-stability paradigm[END_REF][START_REF] Lundström | How to find simple nonlocal stability and resilience measures[END_REF]. However, such "nonlocal" stability measure does not deliver any information in our case because for both models the equilibrium is the unique globally stable attractor (The basin of attraction is the whole positive space, see Methods). We proceed toward a nonlocal resilience measure by invoking the next natural candidate, the return time to equilibrium given a perturbation, and define resilience as the reciprocal of the expected time needed for the population to retain the equilibrium (see Methods). In contrary to many works on resilience, our nonlocal approach can invoke effects from small as well as from large perturbations which is in line with classical definitions of ecological resilience [START_REF] Walker | Stability of semi-arid savanna grazing systems[END_REF][START_REF] Holling | Resilience and stability of ecological systems[END_REF]. Our resilience measure estimates the population's expected rate of return, given a random perturbation. In the Appendix we present an alternative nonlocal resilience measure, the basintime resilience, which is based on the size of a subset of the basin of attraction from which trajectories return fast. Basin-time resilience estimates the probability that the population recovers the equilibrium within a time limit. Results from this measure strengthen our previous conclusions and are illustrated in Appendix Figs. A 11 and A 12.

In general, our approach to resilience is applicable to advanced nonlinear models (with complicated dynamics involving multiple attractors) as well as to simple linear models with one unique stable equilibrium. We have chosen to impose perturbations by random sampling from a uniform distribution, but any set of perturbations may be considered, e.g. normally distributed from equilibrium or a deterministic choice. One may also consider perturbations only in the juvenile-, adult-or the resource dimension. Our resilience measures link the widely used local approach (analyzed through eigenvalues) with the nonlocal one that is usually considered relevant by ecologists (accounting for basins of attraction and large disturbances) as we may consider various ranges of disturbances. We refer the reader to [START_REF] Lundström | How to find simple nonlocal stability and resilience measures[END_REF] for further discussions and constructions of nonlocal stability and resilience measures as well as their relations to local measures. For discussion on the use of local resilience in ecology and the fact that it can be difficult to assess from an empirical point of view, see [START_REF] Haegeman | Resilience, invariability, and ecological stability across levels of organization[END_REF].

Topics for future research

In both the age model and the stage model we considered individuals in only two stages, as juveniles or as adults. We then considered harvesting strategies that allow for different mortality rates in these two stages. Using slightly generalized versions of the age model and the stage model, many more possible harvesting strategies can be explored. A natural first step is to consider harvesting on a size interval, and this can later be extended to include several size intervals as well as more realistic descriptions of harvesting mortality as a function of size. Classical works of [START_REF] Beverton | On the dynamics of exploited fish populations[END_REF] and [START_REF] Holt | The evaluation of fisheries resources by the dynamic analysis of stocks, and notes on the time factors involved[END_REF] consider separate harvesting rates on each year/size class and show that given a fixed harvesting effort, the yield is maximized if fish are caught at the size or age where cohort biomass is maximum. Extending our modeling to allow for different harvesting rates on each year/size class would allow for evaluating such result with respect to our suggested measures of conservation.

We point out that even though the stage model is purely deterministic, one strength of our work is to tackle, through the nonlocal resilience measure and the recovery potential, how the population recover from small as well as large stochastic perturbations, and how the population may survive demographic stochasticity at low density levels. To account for more stochastic effects, a possible direction would be to expand the modeling towards demographic and environmental stochasticity, since the interplay between stochasticity of demographic parameters and deterministic nonlinearity is important [START_REF] Sugihara | Are exploited fish populations stable?[END_REF]. A study in this direction by [START_REF] Engen | Spatial distribution and optimal harvesting of an agestructured population in a fluctuating environment[END_REF] has found that adding environmental stochasticity may change predictions of which harvesting strategy (adult, juvenile, or mixed harvesting) that gives the highest yield.

They show that even when a deterministic model gave the highest yield from adult or juvenile harvesting, adding environmental stochasticity caused mixed harvesting to give the highest yield in many cases. This indicates that our result "equal harvesting is often a good strategy" might be further strengthened when environmental stochasticity is given further account.

Another promising extension of the work presented here is to move beyond single-species management towards ecosystem-based management. We believe in a trend from single-objective towards multi-objective approaches (i.e. optimizing for yield and conservation, not only yield), strengthened by the present paper. This trend may evolve towards multi-objective approaches using multi-species population models, that is, towards ecosystem management. Studies in this direction already exists, see e.g. [START_REF] White | The value of coordinated management of interacting ecosystem services[END_REF], [START_REF] Tromeur | Balancing yield with resilience and conservation objectives in harvested predator-prey communities[END_REF] and [START_REF] Jacobsen | Efficiency of fisheries is increasing at the ecosystem level[END_REF]. Our four measures of conservation can be extended to more general multispecies settings, and the present method using Pareto frontiers to find sustainable harvesting strategies can then be applied also in such general settings. For example, harvesting on a set of k species in a food web with respective harvesting rates (F 1 , F 2 , ..., F k ), we can for any desired yield determine the harvesting strategy that offers the highest conservation benefits. Hence, the methods presented here open a door for reconciling economic and conservation issues in ecosystem management and can be extended to more complex scenarios including for example management of multiple fisheries and maintaining species diversity.

A concept which has attracted considerable attention recently is balanced harvesting (see e.g. [START_REF] Garcia | Reconsidering the consequences of selective fisheries[END_REF], see also the beginning of Discussion. Balanced harvesting strategies should preserve ecosystems' relative size and species composition, and thus harvesting rates may need to be adjusted in proportion to the productivity of individuals. As productivity differs among species, and also within a single species, a balanced harvesting strategy is probably selective and nontrivial to find and implement. [START_REF] Law | Squaring the circle: reconciling fishing and conservation of aquatic ecosystems[END_REF] argue that switching from size-at-entry regulations to balanced harvesting can increase both yield and conservation. Noting that equal harvesting is a more balanced strategy than adult harvesting, their result is in line with the present paper. As our approach can potentially be extended to a general multi-species settings it can also be applied to evaluate general balanced harvesting strategies in the framework of advanced population models.

Appendix-supporting research for the manuscript Meeting yield and conservation objectives by harvesting of both juveniles and adults

The appendix gives additional motivations, explanations and details on the main manuscript, as well as additional numerical and analytical results which strengthen the main findings of the main text. We begin by motivating parameter values as well as investigating robustness with respect to variations of parameter values in the age model, and proceed with a similar section for the stage model. Next, we describe our resilience measure in detail and give some additional resilience investigations of the stage model using an alternative resilience measure. We end by deriving the analytical expression for the basic reproduction ratio for the age model.

Motivations and variations of parameter values in the age model

Concerning the parametrization of the age model, we have used

R 0 = s max = c = 1, a max = 100, K = 0.23, a 0 = -2, (
A1) a mature = 8, h = 0.7, S = 0.8, σ u = 0, as default values in the main text, and we have considered substantial variations from these particular values. When varying parameter values in the age model, we have seen that the preferable harvesting strategies depend on the parametrization in a way that can be well-explained by Fig. 4 in the main text. That is, the suggested harvesting strategy (when comparing adult harvesting, equal harvesting and juvenile harvesting) depends on the survival from natural mortality, S, as well as on the fraction of juveniles in the unharvested case, J * u /(J * u + A * u ). The results in Fig. 4 seem to be robust whenever parameters take on values in the wide intervals

K ∈ [0.1, 1], a 0 ∈ [-3, -0.2], a mature ∈ [3, 15], h ∈ [0.3, 0.9], σ u ∈ [0, 0.5]. (A2)
Indeed, we have varied each parameter at a time, keeping the others at the values given in (A1) and tested at least 10 values in each interval. Further parameter combinations in the above intervals have also been tested.

In the following we will give motivations for the parameter values in (A1) as well as for the considered intervals in (A2). We will also give some explanations of how and why our results depend, or not depend, on the parameter values. Finally, we present additional trade-off curves (similar to those in Fig. 2 in the main text) but for parametrizations corresponding to points P 1 -P 8 in Fig. 4 in the main text. These curves are given in Figs. A 2, A 4, A 5 and A 6. Careful reading of these figures should convince the reader that our main findings are robust with respect to variations of parameter values. We do not present additional versions of main text Fig. 5 (showing relation between measures). However, the correlation can be seen from the trade-off figures by first focusing on the MSY point and then follow curves, for different measures, in the direction of either increasing or decreasing harvesting pressure.

Reduction of parameters

We begin by motivating that without loss of generality we can choose R 0 = s max = c = 1, as well as a max = 100, and we will therefore not consider variations of these parameter values. In particular, by careful investigation of the equations in the main text we see that if we divide f a = cs a in eq. ( 1) by c, eq. ( 2) by s max and eq. ( 4) by R 0 , then R 0 , s max , c, s a , E t , R t and N a,t can be considered as non-dimensional. We obtain non-dimensional equations, similar to the old once but with R 0 = s max = c = 1. This shows why these parameters can be set to 1 without loss of generality. Indeed, the results for arbitrary R 0 and s max can be obtained by multiplying the yield obtained for R 0 = s max = 1, given by main text eq. ( 11), by arbitrary values of these parameters. The parameter R 0 scales the number of individuals in the population, per unit volume, while s max scales the mass of each individual. The constant c sets the number of eggs per individual, and the assumed Beverton-Holt recruitment relation (A3), giving the relation between egg production and offspring, saturates in relation to the number of eggs at the unharvested equilibrium. As our results are independent of the number of eggs in the lake, the value of the constant c is unimportant. This explains also that we can take m a = 1/2 in place of m a = 1 and consider half the population as adults; that would only correspond to considering half the number of eggs per individual, i.e. c = 1/2.

Finally, we will take the maximum age of an individual to be so large that it only cut away a negligible amount of biomass, i.e., very old fish. For our needs, a max = 100 is enough.

Steepness h in the Beverton-Holt recruitment, points P 1 and P 2

We recall the assumed Beverton-Holt recruitment

R t+1 = E t α + βE t exp u t - σ 2 u 2 , ( A3
)
in which u t are independent and normally distributed random variables with mean 0 and standard deviation σ u . The parameters α and β are given by

α = E 0 R 0 1 - h -0.2 0.8h , β = h -0.2 0.8hR 0 ,
where E 0 and R 0 are, respectively, the average egg production and recruitment at equilibrium in the absence of harvesting mortality. The parameter h is the steepness which sets the sensitivity of recruitment with respect to egg production and may take values between 0.2 and 1. The steepness is defined as the ratio of recruitment when egg production equals 20% of E 0 to recruitment at E 0 [START_REF] Mace | A generalised bioeconomic simulation model for fish population dynamics[END_REF][START_REF] Hilborn | Pretty good yield and exploited fishes[END_REF]).

The relation between E 0 and R 0 yields

E 0 = R 0 amax a=0 m a f a S a . (A4)
To derive this relation, put R t = R 0 and γ a-1 = 0 in main text eq. ( 4). It then follows that numbers of individuals at equilibrium, in the absence of harvesting, are

N 0,t = R 0 , N 1,t = R 0 S, N 2,t = R 0 S 2 , . . . , N amax,t = R 0 S amax .
By summing up these individuals total egg production (recall main text eq. ( 3)) we obtain relation (A4).

We have already chosen R 0 = 1 and the value of E 0 follows from m a , f a , s a via (A4). Moreover, by numerical investigations we have convinced ourselves that there is little effect of varying σ u . Therefore, in the following we focus on the steepness parameter h. [START_REF] Myers | Maximum reproductive rate of fish at low population sizes[END_REF] reviewed the steepness of 244 stocks of fish and found that steepness values varies mainly in the range from 0.3 and 0.9. An intermediate value of 0.7 and slightly greater steepness values are common [START_REF] Myers | Maximum reproductive rate of fish at low population sizes[END_REF], Table 1). We have considered variations of steepness in the range h ∈ [0.3, 0.9], and to clarify the dependence of steepness we present additional trade-off curves in Fig. A 2 for h = 0.5 (point P 1 in Fig. 4) and h = 0.9 (point P 2 in Fig. 4). It turns out that our main results are not sensitive to the steepness. We can see in Fig. A 2 that (i) our conclusion that PGY allows large conservation benefits becomes more clear when steepness increases, and that (ii) equal harvesting performs well rather independent of steepness.

To understand why (i) and (ii) should hold in general, and how our conservation measures depend on the steepness, we proceed by noting the following.

If h = 1 then α = 0, β = 1/R 0 E 0 Egg production, E t R 0 Recruitment, R t h = 1 h = 0.9 h = 0.7 h = 0.5 h = 0.2

Beverton-Holt recruitment curves

Figure A 1: Recruitment as functions of egg production following (A3) for different values of the steepness parameter h. When h = 0.2 then recruitment is proportional to egg production, and when h = 1 then recruitment is constant, independent of egg production. Remaining parameters are as in (A1).

and

R t = R 0 (constant recruitment). If h = 0.2 then α = E 0 /R 0 , β = 0 and R t = E t • R 0 /E 0 (proportional recruitment), see Fig. A 1.
Thus, populations with high steepness values can be assigned heavy harvesting pressures, implying small abundances, and still give high yield because the reproduction stays high also for small egg production. To explain this, we note that at unharvested equilibrium we have R t = R 0 and E t = E 0 , i.e. where all curves intersect in Fig. A 1. Imposing harvesting on the population will reduce the population biomass, and hence reduce the total egg production, and we therefore move to the left of the point of intersection in Fig. A 1. If steepness is high (close to 1), then recruitment is still close to R 0 , but if steepness is low (close to 0.2), then recruitment decrease linearly. From this reasoning we realize that high steepness and high harvesting pressure give a situation with (relatively) high reproduction, high harvesting mortality and high yield. Therefore, young individuals should be dominating, and hence the impact on size structure and the impact on biomass will be high. Moreover, the basic reproduction ratio as well as the resilience will be high. These phenomena can be seen in Fig

. A 2.
From the above reasoning we also realize that (i) should hold in general. In fact, when steepness increases PGY can be obtained for wider ranges of harvesting strategies, giving more possibilities to chose harvesting strategies (giving high yield) from. See [START_REF] Hilborn | Pretty good yield and exploited fishes[END_REF] for more on this.

To understand why (ii) should hold in general, we recall that at the unharvested equilibrium we have R t = R 0 and E t = E 0 independent of the steepness h. Hence, the fraction of juveniles at the unharvested equilibrium, J * u /(J * u + A * u ), as well as the location of the parametrization in Fig. 4 in the main text, are also independent of the steepness. As the suggested harvesting strategy mainly depends on J * u /(J * u + A * u ) and S, it follows that it should be rather independent of the steepness h. 4 and h = 0.9 (thick curves) giving point P 2 . Juvenile harvesting (blue, dashed), equal harvesting (black, solid), adult harvesting (red, dash-dot). Remaining parameters are as in (A1) and yield is normalized as in Fig. 2 in the main text. Age at maturity a mature , point P 5 and P 6

The quotient of length of first maturity to maximum length of fish may vary in the wide range from 0.2 to 1, but most species seem to be in the range 0.3 to 0.9 (Fishbase Fig. 43, http://www.fishbase.org/manual). Recalling the von Bertalanffy growth curves in Fig. A 3, we realize that it is satisfactory that our results hold for age at maturity in the interval a mature ∈ [3, 15].

We present additional trade-off curves in Fig.

A 5 for a mature = 5 (point P 5 in Fig. 4) and a mature = 11 (point P 6 in Fig. 4). If a mature = 5 then J * u /(J * u + A * u ) ≈ 35% and adult harvesting is the only strategy on the Pareto front, in the range of PGY, for all conservation measures. If a mature = 11 then J * u /(J * u + A * u ) ≈ 78% and we are at the border where equal harvesting and juvenile harvesting perform rather equally. Thus, a mature has a large effect on the fraction of juveniles in the population, and, therefore, a large impact on the preferable harvesting strategies. Survival from natural mortality S, points P 7 and P 8

Motivated by e.g. [START_REF] Mills | Abundance, annual survival, and recruitment of unexploited and exploited lake charr, Salvelinus namaycush, populations at the Experimental Lakes Area, northwestern Ontario[END_REF] we used the default value S = 0.8 for the survival from natural mortality. We have seen that results from the age model are rather dependent of S, and this dependence has been investigated and presented in Fig. 4 in the main text. shows additional trade-off curves for S = 0.65 (point P 7 in Fig. 4) and S = 0.9 (point P 8 in Fig. 4). If S = 0.65 then J * u /(J * u + A * u ) ≈ 88% and we are at the border where equal harvesting and juvenile harvesting perform rather equally. If S = 0.9 then J * u /(J * u + A * u ) ≈ 32% and adult harvesting is the only strategy on the Pareto front, in the range of PGY, for all conservation measures.

Motivations and variations of parameter values in the stage model

We adopt the following values of the model parameters, H = T = r = 1, R max = 2, σ = 0.5, I max = 10, M = 0.1, z = 0.01, q = 0.85. (A6) These values may be considered as archetypal, and motivations can be found in de [START_REF] De Roos | Simplifying a physiologically structured population model to a stage-structured biomass model[END_REF]. However, these values were derived in the [START_REF] De Roos | Simplifying a physiologically structured population model to a stage-structured biomass model[END_REF] using not only data on fish species. To ensure the validity of results in the setting of fisheries, we give below a rather lengthy section concerning motivations and sensitivity of our results with respect to parameter values.

First, we recall that we have considered substantial variations of these values and concluded, by numerically calculating trade-off curves, that our results from the stage model are robust for the wide ranges of parameter values given by

M ∈ [0, 0.5], z ∈ [0.0001, 0.2], σI max ∈ [3, 100], q ∈ [0.6, 2], R max ∈ [0.5, 100]. (A7)
Indeed, we varied each parameter at a time, keeping the others at the values given in (A6) and tested at least 10 values in each interval. Several further parameter combinations have also been tested, e.g. letting each parameter take on values at the boundary of the intervals in (A7). The upper bound of 100 on the intervals for R max and σI max does not seem necessary; the structures of the population, relevant for our study, seem to remain for much larger values as well. The lower limits of these intervals ensure that the population has a positive stable equilibrium.

In the following we will give, in addition to the motivations in de Roos et al. ( 2008), discussions and motivations for the parameter values in (A6) as well as for the considered intervals in (A7). We will also give some explanations of how and why our results depend, or not depend, on the parameter values. Finally, we present additional trade-off curves (similar to those in Fig. 3 in the main text) but for other parametrizations. These curves are given in Figs.

A 8, A 9 and A 10. Careful reading of these figures should convince the reader that our main findings are robust with respect to variations of parameter values. We do not present additional versions of main text Fig. 6 (showing relation between measures). However, the correlation can be seen from the trade-off figures by first focusing on the MSY point and then follow curves, for different measures, in the direction of either increasing or decreasing harvesting pressure.

Reduction of parameters

We begin by motivating that we can take H = T = 1 without loss of generality. The halfsaturation constant H represents a resource density and is measured in biomass per unit volume. Changing its value can be considered as changes in the volume in which we express the densities R, J and A. Without loss of generality we can choose H = 1. This fixes the environmental volume in which the population is assumed to live and thus scales the biomass densities R, J and A. Mathematically, we rescale the stage model by introducing the new dimensionless variables J ′ = J/H, A ′ = A/H, R ′ = R/H and R ′ max = R ′ max /H, resulting in that H goes away from the set of equations. We also rescale the time variable t with the mass-specific metabolic rate parameter T , measured per unit of time, by introducing the dimensionless time t ′ = tT and expressing the dynamics of the resource and the consumers as functions of this Biomass production functions

w J (R) w A (R) v(w J (R)) v(w J (R))
Figure A 7: Biomass production function for juveniles, w J (R), and for adults, w A (R). The maturation function v(w J (R)) for M + F J = 0.1 (blue, dashed) and for M + F J = 1 (green, dotted). Remaining parameters are as in (A6).

From (A12) we note that the root R * will be independent of both r and R max , and from (A11) we see that both the juvenile and adult biomass, J * and A * , increase linearly with both r and R max , and so is the yield. Moreover, from the expressions in (A11) we see that

J * J * + A * = M + F A M + F A + v(w J (R * ))
and ( A13)

Impact on size-structure = M + F A M + F A + v(w J (R * )) M M + v(w J (R * u )) -1 -1.
Therefore, the fraction of juveniles in the population, as well as the population size structure, are independent of r and R max (because R * is). We can also conclude, from the expressions for the recovery potential and the impact on biomass given in the main text, that both these measures are independent of r. We may also note that the recovery potential becomes independent of R max as R max becomes very large as the ingestion for a single individual saturates for very large resource densities. These arguments clarifies why our results are more or less independent of r and R max (whenever these parameters take on values giving a non-extinct population at equilibrium). Therefore, we feel comfortable with letting r = 1 and R max = 2 and put focus on the remaining parameters.

Ingestion parameters σ and I max

Numerical investigations strongly indicate that our results are very robust with respect to variations of the ingestion parameters σ and I max . To understand why this is the case, we expand on the above analytical reasoning and prove that at equilibrium, the fraction of juvenile biomass in the population and the impact on size structure, see (A13), as well as w A (R * ), w J (R * ) and v(w J (R * )) are in fact independent of r, R max , σ and I max . Indeed, they depend only on the parameters M + F J , M + F A , q, and z. To prove this, consider equation (A12) from which we solve the resource equilibrium biomass. In this equation, R * , σ and I max appear only implicitly in the functions w A (R * ) and w J (R * ), see expressions (A9). Thus, by setting

α * = σI max R * 1 + R *
we can solve (A12) for α * , independent of the parameters σ and I max . Therefore, varying σ and I max will not change w A (R * ), w J (R * ) and v(w J (R * )), and thus not change the fraction of juvenile biomass in the population nor the impact on size structure. This shows why our results are very robust with respect to σ and I max , and from here we fix σ = 0.5 and I max = 10 and focus on the remaining parameters.

Ratio of size at birth to maximum size z De Roos et al. ( 2008) use the value z = 0.01 for the ratio of size at birth to maximum size, z = s born /s max . To further motivate this value in a fisheries context, we note that Punt et al. (1995, page 290) studying the albacore (Thunnus alalunga, Scombridae) use K = 0.23 and a 0 = -1 in the von Bertalanffy growth growth curves, recall Fig. A 3. This gives a value of z ≈ 0.009 which is very close to 0.01. (Changing a 0 to -3 gives z ≈ 0.12, and a 0 = -0.2 gives z ≈ 0.00009.)

By numerically calculating trade-off curves as those in Fig. 3 in the main text, we have seen that our results are robust with respect to variations of z in the wide interval z ∈ [0.0001, 0.2]. Indeed, results from the stage model suggest that equal harvesting performs well rather independent of z as long as it takes on reasonable values. We present additional trade-off curves in Fig. A 8 for the values z = 0.001 and z = 0.1. Here, we clearly see that the value of z has little effect on how juvenile-, adult-and equal harvesting performs.

Difference in ingestion between juveniles and adults q

It seems hard to motivate a specific value for the parameter q describing the difference in ingestion rates between juveniles and adults: "An estimate for this parameter can hardly be derived from experimental data as it is only a phenomenological representation of stage-specific differences in resource availability and resource use between juveniles and adults." [START_REF] De Roos | Simplifying a physiologically structured population model to a stage-structured biomass model[END_REF]. By substantial numerical investigations we have seen that our results hold at least as long as q ∈ [0.6, 2]. We present additional trade-off curves in Fig. A 9 for q = 0.6 and q = 1.5, showing that our main findings hold for these parametrizations as well. The conservation measures stays rather stable also when q is goes out of the bounds 0.6 and 2, but then the difference in yield between the harvesting strategies increase. Indeed, smaller q implies that adult harvesting is better, and larger q implies that juvenile harvesting is better.

Natural mortality rate M

Following De [START_REF] De Roos | Simplifying a physiologically structured population model to a stage-structured biomass model[END_REF] we have adopted the default value of M = 0.1 for the natural mortality rate. Fig. A 10 shows additional trade-off curves for the values M = 0.05 and M = 0.4. The main structure of the curves remains and thus our main results are robust also to changes in M. Let us note that M is simply added to the harvesting rates F J and F A in the stage model, recall (A8). Therefore, an increase in M should force the trade-off curves for juvenile harvesting and adult harvesting toward the trade-off curves for equal harvesting (if M is high, then there is then always high mortality on both juveniles and adults). This can be

Yield

Variation of parameter q We can also see that the impact on size structure is much smaller for higher M, which may be understood by the fact that now the unharvested state already has a high death rate. We can also observe the expected decrease in yield when death rate increases from M = 0.05 to M = 0.4.

Detailed description of our resilience measure

We have considered resilience of the population by measuring the reciprocal of the time needed for the population to recover the positive equilibrium given a random perturbation. To describe the procedure in detail, let (J * , A * , R * ) denote the equilibrium of biomass in the stage model, and let N * a , a = 0, 1, 2, . . . , a max denote the equilibrium of number of individuals in the age model. For a given constant κ > 0 that scales the maximum displacement of the population from the equilibrium, we start a trajectory from a random point uniformly distributed in the cube (0, κJ

* ) × (0, κA * ) × (0, κR * ) (stage model), (A14) (0, κN * 0 ) × (0, κN * 1 ) × (0, κN * 2 ) × • • • × (0, κN * amax ) (age model).
We then find the return time as the time needed for this trajectory to be close enough to the equilibrium in the sense that

J * -J(t) J * 2 + A * -A(t) A * 2 + R * -R(t) R * 2 1/2 ≤ ǫ (stage model), (A15) amax a=0 N * a -N a,t N * a 2 1/2 ≤ ǫ (age model),
for some small ǫ > 0. After repeating this procedure n times we find the resilience, as a function of the harvesting rates F J and F A , by taking the reciprocal of the average of the corresponding return times, i.e.

Resilience = 1 Average of the n return times .

Our resilience measure estimates the population's expected rate of return, given a random perturbation. We have used the parameter values κ = 2, ǫ = 0.01 and n = 100 as default values, but main results are not very sensitive to the magnitude of the imposed perturbations (κ) nor the size of the small neighborhood (ǫ) as long as κ >> ǫ take on reasonable values. We present some results on κ = 1 in Fig. A 12 below, which means that we consider only perturbations that remove biomass. It is worth noting that when implementing our approach to resilience, one tests the model for a large number of initial conditions. Therefore, one may find out if the system is bistable, i.e. if other attractors (e.g. stable states or periodic solutions) are present, which is valuable information. Indeed, if this is the case, a tested initial condition may produce a trajectory that never return to the equilibrium. This will be revealed by our simulation procedure. If applying our approach to bistable systems, one should define another stopping criteria for the other attractors. The resilience can then be defined as the average of the reciprocal of the return times, setting the reciprocal to zero whenever an initial condition reached another attractor. In that way initial conditions (perturbations) that escape to another states contributes with resilience zero, c.f. measure R in [START_REF] Lundström | How to find simple nonlocal stability and resilience measures[END_REF].

Basin-time resilience and further investigations of the stage model

From our results on the stage model we have seen that resilience, in contrary to the other conservation measures, may vastly increase with harvesting pressure for some harvesting strategies. This is true for equal harvesting and juvenile harvesting, but not at all for adult harvesting, recall Fig. 3 in the main text as well as Figs. A 8, A 9 and A 10. Indeed, adult harvesting gives the worst resilience while equal harvesting gives the best resilience (which can be up to 20 times higher in the range of PGY). As this phenomenon provides a substantial argument in favour of equal harvesting compared to adult harvesting, we investigate it further here in order to increase the credibility of our results.

Each dot in Figs. A 11 (a)-(c) show an imposed perturbation as an initial condition from which we integrated the trajectory until it reached the small neighborhood of the equilibrium given by (A15) with ǫ = 0.1. Tested initial conditions are distributed uniformly in (0, κJ * ) × (0, κA * ) with R = R * and κ = 2. Grey dots return to the equilibrium within a time limit which we set to t ≤ 5, while red dots need more time to recover the equilibrium. There is no harvesting pressure in Fig. A 11 (a), adult harvesting F J = 0, F A = 35 in (b) and equal harvesting F J = F A = 0.8 in (c). The harvesting strategies in (b) and (c) give similar yields (close to the maximum possible yield). The fraction of grey dots is largest in case of equal harvesting, while it is smallest in case of adult harvesting, which means that the population returns fast from more perturbations (initial conditions) in case of equal harvesting than for adult harvesting. This reflects the higher resilience in case of equal harvesting. Indeed, the fraction of grey dots estimates the probability that the population recover the equilibrium in a time limit, given a random perturbation, and constitutes a natural candidate for measuring resilience.

In general, let basin-time be a subset of the basin of attraction from which trajectories return to (a small neighborhood of) the attractor within a time limit, t ≤ τ . A large and convex basin-time set should then reflect high resilience. Therefore, both the size and the shape of the basin-time constitute natural candidates for measuring resilience [START_REF] Lundström | How to find simple nonlocal stability and resilience measures[END_REF]). We will consider the size of the basin-time set as a resilience measure. To do so we impose perturbations as described in (A14) and integrate each trajectory until it reaches the neighborhood defined by (A15). We further let N τ return denote the number of initial conditions from which trajectories return in time τ , and N tot denote the total number of tested initial conditions. We measure resilience as the size of the basin-time set through the estimate b) and (c) may indicate that the higher resilience in case of equal harvesting is due to fast return of trajectories starting at high biomass. Indeed, the basin-time, estimated by the fraction of grey dots, is largest in Fig. A 11 (c) and includes all points corresponding to larger initial biomass than at equilibrium. Therefore, we complement the given resilience estimations by a simulation in which we consider only perturbations that remove biomass. In particular, we take κ = 1 and reproduce the trade-off curves between resilience (both measures) and yield, see Fig. A 12. We can observe that equal harvesting is now challenged by juvenile harvesting, but both resilience measures still agree on the result that equal harvesting is much better than adult harvesting.

In addition to the above resilience estimations, we have also considered normally distributed perturbations (centered at the equilibrium) as well as the "standard" local resilience measure given by the real part of the eigenvalue of the Jacobian matrix at equilibrium with largest real part. Results from these investigations also agree very well with our findings.

For further ideas on how to use the basin of attraction, the basin-time as well as the return time to measure stability and resilience, we refer the reader to [START_REF] Lundström | How to find simple nonlocal stability and resilience measures[END_REF] and references therein.

Derivation of the basic reproduction ratio for the age model

The basic reproduction ratio represents the average number of offsprings produced over the lifetime of an individual in the absence of density-dependent competition. We consider an individual satisfying the dynamics of the age model and calculate an expression for this individuals expected offsprings during its lifetime. For this purpose, we let Ẽt denote the egg production and Rt the recruitment from this individual, at time t. That no density-dependent competition is present lead us to assume that the Beverton-Holt recruitment is completely unsaturated. Such situation is achieved by taking β = 0 (or R 0 = ∞) in (A3), giving Rt+1 = Ẽt

α exp u t - σ 2 u 2 in which α = 1 - h -0.2 0.8h amax a=0 m a f a S a , ( A16
)
where we used relation (A4) to find the expression for α. Now, since u t -σ 2 u /2 is normally distributed with mean -σ 2 u /2 and standard deviation σ u , we obtain, where E denotes expectation,

E exp u t - σ 2 u 2 = 1.
Therefore, the expected number of offsprings in year t becomes where α is given in (A16). It remains to calculate the total number of eggs produced, i.e. the sum in the right hand side of (A17). Indeed, following (3) and (4) in the main text we conclude that

amax t=0 Ẽt = m 0 f 0 + m 1 f 1 S(1 -γ 0 ) + m 2 f 2 S 2 (1 -γ 0 )(1 -γ 1 ) + . . . = m 0 f 0 + amax a=1 m a f a S a a-1 i=0 (1 -γ i ).
By the expression for m a in (1) in the main text and by γ a in (5) in the main text we find that the above expression yields Inserting main text eq. ( 3) and (A18) into (A17) and simplifying yield R0 (F J , F A ) = (1 -F J ) amature × amax a=amature+1 s a S a (1 -F A ) a-amature 1 -h-0.2 0.8h × amax a=amature+1 s a S a . This expression is identical to eq. (10) in the main text which gives the basic reproduction ratio for the age model.

  2 (a)-(d) and Figs. 3 (a)-(d), that by stepping back in yield by 20% into the range of PGY, we can find harvesting strategies

Figure 2 :

 2 Figure2: Trade-offs between yield and the four conservation measures in the age model. The gray regions show "all possible" combinations that can be realized when varying the harvesting rates on juveniles and adults. The solid green curves represent the Pareto front, while the dotted grey lines give the border for PGY, i.e. 80% of MSY. The yellow dots represent MSY and the green squares give the unfished state. We observe that within the range of PGY, equal harvesting performs well with respect to all measures. The black dots represent a suggested harvesting strategy, within the range of PGY, produced by F A = F J = 0.06. Parameter values are as in (7) and yield normalization is as in Figure1.

Figure 3 :

 3 Figure 3: Trade-offs between yield and the four conservation measures in the stage model. The gray region, curves, dots, and squares are as in Fig. 2. We observe that within the range of PGY, equal harvesting performs well with respect to all measures. The black dots represent a suggested harvesting strategy, within the range of PGY, produced by F A = F J = 0.8. Parameter values are as in (9) and yield is normalized as in Fig. 1.

Figure 5 :

 5 Figure5: The relation between impact measures and risk measures for the age model. We observe that a large impact on biomass implies a low basic reproduction ratio (recovery potential) and also a low resilience. The same is true for impact on size structure. Curves, green squares and parameters are as in Fig.2.

  (a), (c) and (d) only, not including the Pareto efficiency in Figs. 2 (b) and 3 (b). This would result in a shift towards recommending adult harvesting, especially in case of the age model.

Figure A 2 :

 2 Figure A 2: Age model for h = 0.5 (thin curves) giving point P 1 in Fig.4and h = 0.9 (thick curves) giving point P 2 . Juvenile harvesting (blue, dashed), equal harvesting (black, solid), adult harvesting (red, dash-dot). Remaining parameters are as in (A1) and yield is normalized as in Fig.2in the main text.

Figure A 3 :

 3 Figure A 3: Growth of an individual as functions of age following the von Bertalanffy growth curves (A5) for different values of the parameters a 0 and K. (a) variation of a 0 , (b) variation of K. Remaining parameters are as in (A1).

  Figure A 4: Age model for a 0 = -3 (thin curves) giving point P 3 in Fig. 4, and a 0 = -1 (thick curves) giving point P 4 . Remaining description is as in Fig. A 2.

  Figure A 5: Age model for a mature = 5 (thin curves) giving point P 5 in Fig. 4, and a mature = 11 (thick curves) giving point P 6 . Remaining description is as in Fig. A 2.

Figure A 6 :

 6 Figure A 6: Age model for S = 0.65 (thin curves) giving point P 7 in Fig. 4, and S = 0.9 (thick curves) giving point P 8 . Remaining description is as in Fig. A 2.

  Fig. A 6

Figure A 9 :Figure A 10 :

 910 Figure A 9: Stage model for q = 0.6 (thin curves) and q = 1.5 (thick curves), remaining description is as in Fig. A 8.

Basin

Figure A 11

 11 Figure A 11 (d) shows trade-off curves between basin-time resilience and yield. (parameters are κ = 2, n = 1000, ǫ = 0.1 and τ = 5.) It is clear that also the basin-time resilience measure strengthen the result of the stage model saying that equal harvesting performs best while adult harvesting gives the worst resilience; these trade-off curves are very similar to those in Fig.3(c) in the main text.Figs. A 11 (a), (b) and (c) may indicate that the higher resilience in case of equal harvesting is due to fast return of trajectories starting at high biomass. Indeed, the basin-time, estimated by the fraction of grey dots, is largest in Fig.A 11 (c) and includes all points corresponding to larger initial biomass than at equilibrium. Therefore, we complement the given resilience estimations by a simulation in which we consider only perturbations that remove biomass. In particular, we take κ = 1 and reproduce the trade-off curves between resilience (both measures) and yield, see Fig.A12. We can observe that equal harvesting is now challenged by juvenile harvesting, but both resilience measures still agree on the result that equal harvesting is much better than adult harvesting.In addition to the above resilience estimations, we have also considered normally distributed perturbations (centered at the equilibrium) as well as the "standard" local resilience measure given by the real part of the eigenvalue of the Jacobian matrix at equilibrium with largest real part. Results from these investigations also agree very well with our findings.For further ideas on how to use the basin of attraction, the basin-time as well as the return

Figure A 12 :

 12 Figure A 11: (a)-(c): Population recover fast from most perturbations in case of equal harvesting (c) and from fewest perturbations in case of adult harvesting (b).The grey dots represent initial conditions that return to the equilibrium within time 5, while the red dots need more time to recover the equilibrium. There is no harvesting pressure in (a). (d): Trade-off between basintime resilience and yield. Juvenile harvesting (blue, dashed), equal harvesting (black, solid), adult harvesting (red, dash-dot). Model parameters are as in (A6) and yield is normalized as in Fig.3in the main text.

  expected total number of offsprings from the single individual during its lifetime, i.e. the basic reproduction ratio R0 (F J , F A ), becomes R0 (F J , F A ) = E

  F J ) amature × amax a=amature+1 f a S a (1 -F A ) a-amature .We next use the relation f a = cs a given in (1) in the main text to deriveamax t=0 Ẽt = c (1 -F J ) amature × amax a=amature+1 s a S a (1 -F A ) a-amature .Using (1) in the main text once more we see that α can be writtenα = c 1 -h -

Parameters a 0 and K in the von Bertalanffy growth, points P 3 and P 4

We have assumed von [START_REF] Von Bertalanffy | Quantitative laws in metabolism and growth[END_REF] growth to describe individual length as a function of age, and that individual mass is proportional to the cube of individual length, i.e.

where s a is the mass of an individual at age a, s max is the asymptotic maximum body mass, K is a growth rate parameter and a 0 is the hypothetical negative age at which the individual has zero length. We have already motivated s max = 1. Punt et al. (1995, page 290) studying the albacore (Thunnus alalunga, Scombridae) motivate us to take a 0 = -1 and K = 0.23. However, the value of a 0 may be decreased to -2 as well since the age model is discrete and calculates the size of all individuals in the beginning of the year, but harvesting and egg production naturally occur continuously during the year. We have considered variations in the ranges a 0 ∈ [-3, -0.2] and K ∈ [0.1, 1]. Figure A 3 shows the corresponding growth curves for some values of a 0 and K.

We present additional trade-off curves in Fig. A 4 for a 0 = -3 (point P 3 in Fig. 4) and a 0 = -1 (point P 4 in Fig. 4). When a 0 = -1 then especially young individuals will have lower biomass and therefore the fraction of juveniles in the population will be small, J * u /(J * u + A * u ) ≈ 53%. As a result, adult harvesting performs better, with respect to all conservation measures, than equal harvesting and juvenile harvesting. When a 0 = -3 then especially young individuals will be of more biomass and therefore the fraction of juveniles in the population will be larger, J * u /(J * u + A * u ) ≈ 66%. In this case equal harvesting is on the Pareto front within the range of PGY for impact on biomass and resilience, and, therefore we recommend equal harvesting.

Considering variations in the growth parameter K gives a similar effect: A large K implies that individuals grow fast and this affects especially young individuals. If K = 1 then J * u /(J * u + A * u ) ≈ 81% and juvenile harvesting is the only strategy on the Pareto front, in the range of PGY, for all conservation measures. If K = 0.1 then J * u /(J * u + A * u ) ≈ 41% and adult harvesting is the only strategy on the Pareto front, in the range of PGY, for all conservation measures. new time variable. Moreover, we introduce the new dimensionless rate parameters r ′ = r/T , I ′ max = I max /T , M ′ = M/T , F ′ J = F J /T and F ′ A = F A /T , resulting in that T goes away from the set of equations. This shows that we can choose H = T = 1 without loss of generality, and that results for other values of the parameters H and T can be obtained by multiplications from the results for H = T = 1. We therefore omit further investigations in H and T and put focus on the remaining parameters. For simplicity, we will omit the "prime"-notation in the following even though we work with the dimensionless quantities. The rescaled stage model equations are identical to main text eqs. ( 8) but with H = T = 1. In particular,

where the net biomass production by a juvenile and an adult are, respectively,

and the maturation is given by

for x = M + F J and v(M + F J ) = -(M + F J )/ log(z).

For clarity we plot the net biomass production functions for juveniles and adults, w J (R) and w A (R), as well as the maturation function v(w J (R)), in Fig. A 7. All functions are zero (juveniles and adults do not produce biomass, and juveniles do not mature) until energy intake is sufficient to cover maintenance requirements, after which they increase with resource abundance.

Resource parameters r and R max Numerical investigations strongly indicate that our results are very robust with respect to variations of the resource parameters r and R max . We will now give a theoretical explanation of this. Let J * , A * and R * denote the juvenile, adult and resource biomass at equilibrium, respectively. Let also J * u , A * u and R * u denote these quantities in case of no harvesting pressure, i.e. when F J = F A = 0. Setting derivatives to zero in the stage model (A8) and solving the corresponding set of equations for equilibria we find that

where R * is given as the unique root of the equation