Implications of fisheries-induced evolution for stock rebuilding and recovery

Worldwide depletion of fish stocks has led fisheries managers to become increasingly concerned about rebuilding and recovery planning. To succeed, factors affecting recovery dynamics need to be understood, including the role of fisheries-induced evolution. Here we investigate a stock's response to fishing followed by a harvest moratorium by analyzing an individual-based evolutionary model parameterized for Atlantic cod Gadus morhua from its northern range, representative of long-lived, late-maturing species. The model allows evolution of life-history processes including maturation, reproduction, and growth. It also incorporates environmental variability, phenotypic plasticity, and density-dependent feedbacks. Fisheries-induced evolution affects recovery in several ways. The first decades of recovery were dominated by demographic and density-dependent processes. Biomass rebuilding was only lightly influenced by fisheries-induced evolution, whereas other stock characteristics such as maturation age, spawning stock biomass, and recruitment were substantially affected, recovering to new demographic equilibria below their preharvest levels. This is because genetic traits took thousands of years to evolve back to preharvest levels, indicating that natural selection driving recovery of these traits is weaker than fisheries-induced selection was. Our results strengthen the case for proactive management of fisheries-induced evolution, as the restoration of genetic traits altered by fishing is slow and may even be impractical.


Introduction
One quarter of the world's fish stocks are overexploited, depleted, or recovering, according to the United Nations' Food and Agriculture Organization (FAO 2009). Although biological extinctions are very rare (Dulvy et al. 2003;Swain and Chouinard 2008), several of these declines have lead to collapses of fishing activities to a state of 'commercial extinction,' where targeted fisheries are no longer commercially viable (Myers et al. 1996). Infamous examples of commercial extinctions of major fisheries targets include stocks of sardine Sardinops sagax off California and Japan in the late 1940s, the Peruvian and Chilean stocks of anchovy Engraulis ringens in 1972 (Csirke 1977;Murphy 1977), the Norwegian spring-spawning stock of herring Clupea harengus in 1968 (Toresen and Østvedt 2000), and the Newfoundland-Labrador stock of cod Gadus morhua through the 1980s (Hutchings 1996). Classic theory of fishing suggests rapid population recovery if fishing is ceased, but in practice recovery rates have been much slower than expected, and in some cases the expected recovery has not taken place at all (Hutchings 2000a).
Stock collapses have enormous social and economic costs (e.g., Haedrich and Hamilton 2000), and painful experiences have led politicians and fisheries management institutions to gradually shift focus from 'how much to catch' to 'how to make sure there is something to catch'.
For example, the United Nations 2002 World Summit on Sustainable Development declared that all fish stocks should be restored to levels that produce maximum sustainable yield by the year 2015. In the United States, the Magnuson-Stevens Fishery Conservation and Management Act of 1996 mandates that overfished stocks should be rebuilt within <10 years, unless some circumstance such as species biology dictates a longer time frame (Safina et al. 2005;Rosenberg et al. 2006).
To rebuild a stock successfully, it is crucial to have a good understanding of the factors that influence recovery dynamics and the associated timescales. Recovery can be interpreted as a reversal of processes involved in fishing a population down to low abundance. Hence, if we first understand what happens to fish stocks when we harvest them, it then becomes easier to address the subsequent question of how various harvest-induced changes are reverted if fishing is reduced or ceased. The dynamics a population exhibits in response to fishing involves changes at several levels of system organization. Below we review five of these. The classic theory of fishing considers the first two levels. The most obvious effect of fishing is a reduction in population abundance and biomass. This effect is represented in all fishery models. However, not all biomass is the same: at a second level, fishing changes the demographic composition of a stock toward a dominance of younger and smaller fish. A truncated age and size structure may have consequences for population dynamics, and has been shown to reduce reproductive potential (Marteinsdottir and Thorarinsson 1998;Murawski et al. 2001), increase variability in recruitment or population abundance (Longhurst 2002;Hsieh et al. 2006), and make a stock more vulnerable to environmental fluctuations (Ottersen et al. 2006;Hsieh et al. 2008).
At a third level, fishing may change the biotic environment, triggering changes in phenotypically plastic traits of individuals. Fishing may reduce intraspecific competition, thereby promoting phenotypic plasticity in the form of increased growth, which often leads to earlier maturation (Trippel 1995;Lorenzen and Enberg 2002;Kell and Bromley 2004). Fishing may also change interspecific interactions, an appreciation that has sparked a move toward the ecosystem approach to fisheries management (Pikitch et al. 2004;Francis et al. 2007).
At a fourth level, an exploited population may be evolutionarily adapting to the new mortality regime (Rutter 1902;Law and Grey 1989). A growing body of research suggests that such fisheries-induced evolution is taking place in a number of fish species and stocks worldwide (reviewed in e.g. Jørgensen et al. 2007; Kuparinen and Merilä 2007;Allendorf et al. 2008;Fenberg and Roy 2008;Heino and Dieckmann 2008;Hutchings and Fraser 2008;Heino and Dieckmann in press). The possibility that these changes are genetic means that fishing may be changing harvested species more fundamentally than previously thought. There are theoretical reasons to believe that such evolutionary changes may be slow and, within practical timescales, even impossible to reverse (Law and Grey 1989;). Yet, a recent laboratory study on Atlantic silverside Menidia menidia showed that the reversal of at least one key life-history trait (body size) was possible and in some cases relatively fast, although the recovery rate would depend on the form of selection that led to the changes in the first place (Conover et al. 2009).
At a fifth level, because of several of the processes above, harvesting can lead to changes in ecosystem structure, potentially causing ecological regime shifts and alternative stable states (Jackson et al. 2001;Scheffer et al. 2005). Frank et al. (2005) have suggested that removal of top predators from the ecosystem may lead to cascading effects that affect the whole ecosystem. Such trophic cascades may have long-lasting consequences, and also influence recovery processes of harvested populations (Frank et al. 2005).
The extents to which overfishing affects stock biomass, population structure, phenotypic plasticity, adaptive evolution, or ecosystem structure have important implications for recovery. In particular, the failure of some fish stocks to recover after their collapse (Hutchings 2000a;Hutchings and Reynolds 2004) has raised concerns as to whether fisheries-induced evolution is contributing to this lack of recovery (Hutchings 2004(Hutchings , 2005Walsh et al. 2006). For example, before the infamous collapse of the northern cod off Newfoundland and Labrador, marked changes occurred in key life-history traits such as maturation schedule (Olsen et al. 2004). Around the same time, other cod stocks in the region showed changes in growth that are interpreted as being indicative of evolution (Swain et al. 2007. The design and implementation of a rebuilding plan depend on the nature of the changes in population-level and individual-level characteristics that occurred through fishing, and on the degree to which these changes are anticipated to be reversible. Our aim in this study is to investigate the rebuilding of stocks following a period of exploitation, focusing on the role of fisheries-induced evolution and using a model that includes four of the five biological response levels described above (excluding ecosystem-level effects). The life history of our model population resembles that of a slow-growing and latematuring fish species such as Atlantic cod Gadus morhua in the northern parts of its range. By comparing the evolutionary model to a model version in which evolution is 'turned off,' both during the exploitation phase and throughout the subsequent recovery, we address how fisheries-induced evolution affects the biological dynamics and examine implications for recovery. Focal questions are: How will adaptive evolution change the expected rate and extent of recovery? And, when fisheries management is oblivious of fisheries-induced evolution, what errors will be made?

Model description
We use an individual-based model that combines the quantitative genetics of evolving life-history traits with individual-level ecological processes of growth, survival, and reproduction and with population-level ecological processes such as density dependence and environmental variability. Our modeling methodology extends earlier individual-based evolutionary models (Holland 1992;Huse et al. 1999) and falls within the framework of eco-genetic modeling Höök and Wang 2009;Okamoto et al. 2009).
Briefly described, each individual carries genetic traits affecting its life history through growth rate, maturation schedule, and reproductive investment. An individual's genetic traits are expressed imperfectly, to allow for the chance environmental variation of phenotypic traits around genetic traits. The individual's genetic traits together with its environment thus determine its life history of growth and reproduction, and its risk of dying from natural causes or fishing. Population dynamics emerge when many such interacting individuals are coupled through a shared environment. The environmental influences are manifested as density dependence and stochastic fluctuations acting on recruitment and growth conditions, and fitness emerges as the success of individuals to grow, survive, mate, and reproduce-thus producing offspring that transmit their genetic traits to future generations. We model only female life histories. However, the 'females' in our model reproduce sexually and mate with each other, and can therefore be considered as hermaphrodites. The model thus includes sexual reproduction and evolves toward evolutionarily stable strategies for the female sex; it is not meant to predict how evolution might differentially affect males.

Environmental variability
A model in which genotypes are inherited and selection acts on phenotypes requires careful consideration of the noise processes affecting the link from genotypes to phenotypes. Without environmental variation confounding this link, the coupling of genotypes to selection pressures would be too strong, the quantitative traits would have unrealistically high heritabilities, and the speed of evolution would be exaggerated. Adding environmental variation obscures the link from genotypes to phenotypes, and thus to selection pressures, and helps bring heritabilities down to levels often observed in nature (0.2-0.3 for life-history traits ;Gjedrem 1983;Carlson and Seamons 2008). In living organisms, the 'noise' in the correlation between genotypes and selection pressures is a conglomerate of different processes, most of which are not fully understood. We have chosen to infuse such environmental variability at different stages through the following noise components: (i) Inheritance of genetic traits. For each genetic trait, an offspring's trait value typically deviates from the mid-parental value, reflecting the effects of mutation, segregation, and recombination. (ii) Phenotypic expression of genetic traits. At birth and for each genetic trait, the individual's expressed phenotypic trait deviates from its inherited genetic trait, reflecting micro-environmental variation as well as chance effects of epistatis and dominance. (iii) Population-level inter-annual variation in growth conditions. This reflects temporal fluctuations in the abiotic and biotic environment, for example, in temperature or resource availability. (iv) Individual-level inter-annual variation in growth conditions. This reflects variability between individuals because of chance events, for example, in resource acquisition or environmental exposure. (v) Population-level variation around recruitment function. This reflects a stochastic element of population dynamics and affects growth conditions through density dependence. (vi) Demographic stochasticity in mortality. An individual's probability of dying from natural causes or harvesting is determined by its growth strategy and size, but whether it actually dies is a random event. (vii) Demographic stochasticity in mating. Similarly, an individual's probability of being a parent depends on its gonad size, but random parents are drawn according to this probability. We have used normally distributed random variables for (i), (ii), and (iv). For (iii) and (v), lognormal distributions were used. A lognormal distribution was chosen for (iii) because ecological data is often characterized by distributions with a long tail (Hilborn and Mangel 1997) and for (v) following several studies that consider this process lognormal in nature (Hennemuth et al. 1980;Caputi 1988;Fogarty 1993a,b;Myers et al. 1995).

Genetic traits and their expression
Each individual possesses four inherited quantitative genetic traits: a growth coefficient g that affect its resource acquisition, two traits that specify its maturation schedule through the slope s and intercept y of a linear probabilistic maturation reaction norm (PMRN; Heino et al. 2002a;Dieckmann and Heino 2007) and a trait r that quantifies its reproductive investment in terms of its gonado-somatic index (gonad mass/somatic mass) and thus governs resource allocation to reproduction from maturation onward. The population's distribution of these genetic traits determines their additive genetic variances. During expression, phenotypic traits (denoted by G, S, Y, and R) are subject to different noise processes as described above. Here and in the following, e with different subscripts denotes a random number drawn from a normal distribution with mean 1 and given standard deviation ( Table 1). The phenotypic trait S(i) for individual i's genetic PMRN slope s(i) is thus where e s (i) describes the expression noise and is drawn once per lifetime (the argument i in e s (i) indicates that different values are drawn for each individual). Analogously, we have YðiÞ ¼ e y ðiÞ Á yðiÞ and RðiÞ ¼ e r ðiÞ Á rðiÞ (see Table 2 for units).
The phenotypic growth coefficient G is influenced by four processes. First is the expression noise e g (i). Second is an individual-level inter-annual noise, e i (i, t). Third is a population-level inter-annual noise e k E ðtÞ , where k E ðtÞ is a normally distributed random deviate (Table 1) so that e kEðtÞ is lognormally distributed. This noise component represents environmental fluctuations that influence all individuals in a similar way, for example, through resource availability or temperature, which affect the growth of many fish species, including Atlantic cod (Hansson et al. 1996). Fourth is the population-level density dependence D(t) specified under 'density-dependent growth' below (Eqn 14). Thus, Gði; tÞ ¼ e g ðiÞ Á e i ði; tÞ Á e kEðtÞ Á DðtÞ Á gðiÞ: The variance in the expression noise for each trait is set such that the total expressed variance r 2 E is related to the additive genetic variance r 2 A as r 2 E ¼ r 2 A h À2 À 1 ð Þ , with r 2 A determined by an assumed initial genetic coefficient of variation CV G and by the initial mean trait values (Table 1; for further details see ). Our assumed value of 6% for CV G is on the conservative end of estimates from empirical work (Houle 1992;see Dunlop et al. 2007 for the sensitivity of evolutionary rate to assumed initial genetic coefficient of variation; see also ). See Table 1 for means and standard deviations of the different noise processes. Each listed random variable is drawn from a normal distribution with the shown mean and standard deviation. -, dimensionless parameters. *Equation applies by analogy.

Life-history processes
An individual's phenotype consists of its expressed genetic traits S(i), Y(i), R(i), and G(i, t), its length L(i, t), its age A(i), and its maturity status. The time step in our model is 1 year. The annual modeling cycle can be divided into four processes: (1) sexual maturation, (2) growth of soma and gonads, (3) natural and fishing mortality, and (4) reproduction and inheritance.

Sexual maturation
We model the maturation process using PMRNs (Fig. 1, Heino et al. 2002a;reviewed in Dieckmann and Heino 2007). Whether an immature individual is likely to mature in a given year probabilistically depends on its age, size, and its PMRN. We assume a linear PMRN with constant width. The PMRN is thus determined by the phenotypic intercept Y(i) and slope S(i). The probability to mature is described by the logistic regression: where L(i, t) is the length of individual i in year t, and L p50 (i, t) is the length at 50% maturation probability, calculated as L p50 ði; tÞ ¼ YðiÞ þ SðiÞ Á Aði; tÞ. The parameter d is determined by the PMRN width w as with logit p ¼ lnðp=ð1 À pÞÞ, where p l and p u specify the lower and upper probability bounds, respectively, chosen for defining the PMRN width (in our model, we choose quartiles, p l = 0.25 and p u = 0.75; Fig. 1A). Maturation is modeled as a stochastic process of Bernoulli trials, and takes place if a number randomly drawn from a uniform distribution between 0 and 1 is smaller than p(i, t).

Growth of soma and gonads
For growth, we use a model that generalizes the model by Lester et al. (2004), flexibly treating the allometric scaling exponents as parameters (D. S. Boukal and U. Dieckmann, unpublished). Length L was assumed to scale with weight W as W ¼ a 1 L b1 , and gonadic and

Fisheries-induced evolution and stock recovery
Enberg et al.
somatic tissue were assumed to be energetically equivalent. Resource acquisition scales with weight as d dt W ¼ Gði; tÞW 1Àb 2 year À1 . The growth coefficient G(i, t) thus specifies the amount of resources an individual has available, which it can invest into the growth of soma or gonads. Individual length growth is thus determined by the phenotypic growth coefficient G(i, t) and reproductive investment R(i) according to Before maturation, R(i) = 0, as juveniles use all acquired resources for somatic growth. Equation (5) implies a maximum possible gonado-somatic index of G i; t ð ÞL i; t ð Þ Àb 1 b 2 a Àb 2 1 , at which R(i) hence is capped. Individual fecundity Q is given by gonad weight divided by the weight W e of a single egg, The amount of available resources is variable, so that under unfavorable resource conditions, mature individuals may not have enough resources to grow as well as to reproduce. In such cases, individuals will prioritize reproduction over growth.
Natural and fishing mortality Natural mortality consists of three components: (i) size-independent mortality m 0 because of, for example, diseases and parasites; (ii) size-dependent predation mortality m p during activities other than foraging (e.g., resting, migrating, and hiding); and (iii) size-dependent predation mortality m f related to foraging (C. Jørgensen and Ø. Fiksen, unpublished data). We base the size dependence of mortality on observations in marine systems showing that mortality scales with length as L Àd 2 with an allometric exponent of about d 2 = 0.75 (Peterson and Wroblewski 1984;Brown et al. 2004). Because length can change substantially over 1 year, we use the average length in a year to determine the instantaneous rate of predation mortality, While smaller individuals are generally more vulnerable to predators, all individuals can accept a higher foraging mortality to achieve a higher resource intake and, consequently, a higher growth rate (Walters and Juanes 1993;Biro et al. 2006;Biro and Post 2008;C. Jørgensen and Ø. Fiksen, unpublished data). In our model, the higher resource intake enabled by higher risk-exposure thus implies higher foraging mortality, The foraging mortality scales with the overall sizedependent predation mortality m p , and is thus higher for smaller fish, which is in line with observations of juvenile fish spending much of their time hiding from predation and trying to minimize the risk associated with foraging (Walters and Juanes 1993 and references therein). This risk associated with foraging can depend, for example, on the total time spent foraging, which, when increased, results in higher encounter rates with predators.
In addition to the natural mortality components, individuals are potentially subject to fishing mortality at an instantaneous rate F. Fishing mortality is size-dependent, and we use a sigmoid selectivity curve as follows, where g determines the steepness of the selectivity curve, 1 2 ðLði; tÞ þ Lði; t þ 1ÞÞ is the mean length of individual i in year t, and L 50 is the length at which an individual has a 50% probability of being captured relative to the asymptotic maximum capture probability at large sizes (the maximum slope of selectivity as a function of mean length occurs at L 50 and equals 1 4 g). The instantaneous fishing mortality rate depends on the selectivity curve and on the harvest rate f max at sizes at which fish are fully vulnerable to the fishery, The total instantaneous mortality rate is , and individual i's resultant annual probability of dying is P(i, t) = 1)e )Z(i,t)AEyear , which the model again realized through Bernoulli trials. The parameter values chosen (Table 3) produce natural mortality rates that are comparable with estimates from field studies and that are used in assessment work of Atlantic cod (Sinclair 2001;ICES 2003), giving a total instantaneous natural mortality rate of m 0 + m f (i, t) + m p (i, t) % 0.25AEyear )1 for an individual of the average age and size at maturation.

Reproduction and inheritance
The number n 0 (i, t) of offspring produced by parent i in year t is proportional to that parent's share of total population fecundity, Here, the sum extends over the entire mature population and N 0 (t) is the total number of recruits in year t as determined by Beverton-Holt recruitment as explained under 'density-dependent recruitment' below (Eqn 13). The parent produces each offspring with a randomly selected partner. The partner is chosen with a probability proportional to its gonad size. The use of multiple partners is motivated by the many marine fish that are batch spawners. For example, Atlantic cod can produce 20+ batches within a month (Kjesbu et al. 1996) and there is thus a high probability that the offspring of one female are sired by several partners.
In our model, trait inheritance follows quantitative genetics theory (Roughgarden 1979;Falconer and Mackay 1996;Lynch and Walsh 1998), as life-history traits are usually highly polygenic quantitative characters determined by many loci (Roff 1992;Conner and Hartl 2004). An offspring inherits the genetic traits from its parents, and we randomly draw the offspring's trait values from a normal distribution with a mean given by the trait's midparental value. Thus, the offspring o of parents i and j will have its genetic traits, here specified for the growth trait, determined by where h g is randomly drawn from a normal distribution with zero mean that reflects the effects of mutation, segregation, and recombination of the underlying loci; its standard deviation is specified for each trait separately (Table 1) and corresponds to a coefficient of variation of 7.1% in the population prior to fishing, implying a constant mutation-segregation-recombination kernel (Roughgarden 1979). The emergent heritability for age at maturation is around 0.2 at equilibrium before fishing, a conservative value within the range typically observed for life-history traits in general (Gjedrem 1983;Law 2000;Carlson and Seamons 2008) and for the proportion of mature 2-year-old Atlantic cod in particular (Kolstad et al. 2006).  Gjedrem (1983), Mousseau and Roff (1987), Houle (1992), and Carlson and Seamons (2008); (2) values chosen such that the life-history characteristics resemble those of Atlantic cod in its northern range (e.g., Rose and Driscoll 2002;Heino et al. 2002b;McIntyre and Hutchings 2003;Marshall et al. 2004;Olsen et al. 2004Olsen et al. , 2005. (3)  Fisheries-induced evolution and stock recovery Enberg et al.

Density regulation
Density-dependent recruitment The number of newborns in a given year is density-dependent and determined by a Beverton-Holt recruitment function (Beverton and Holt 1957) that depends on the total fecundity P j Qðj; tÞ of the population (Fig. 1C), Here, e k R ðtÞ describes population-level inter-annual environmental variability in recruitment modeled as a lognormal process, where k R ðtÞ is a normally distributed random deviate (Table 1, see Fig. 1F for a time series of recruitment). The parameter a measures the survival of recruits when total fecundity is low, while b specifies the strength of density dependence; the asymptotic number of recruits when total fecundity is high is given by a/b (Quinn and Deriso 1999).

Density-dependent growth
Conspecific density may affect the growth of fish through resource competition (Lorenzen and Enberg 2002). Realized growth is thus influenced also by population biomass (Fig. 1B). We express the effect of density dependence through the factor D(t) that influences growth multiplicatively (Eqn 2), where d 1 specifies the strength of density dependence, B(t) is the total biomass of all individuals aged 1 year or older in year t, and we choose B so that the time average of D(t) equaled 1 at the stochastic equilibrium before harvesting.

Model parameterization and model runs
The model was parameterized to describe a population resembling Atlantic cod Gadus morhua in the northern part of its range. We expect this parameterization to be roughly representative also of other slow-growing and late-maturing fish species. Where available, parameters for cod were taken from the literature (Table 3). However, some parameters are unknown and cannot readily be estimated from available data. Thus, the unknown parameters were chosen following a pattern-oriented modeling strategy, which ensures that the emergent model properties qualitatively and quantitatively resemble the observed natural patterns (Grimm et al. 2005). To achieve this, the model was initially run with likely parameter values and its output compared with data available in the literature and in stock-assessment reports. Parameters responsible for discrepancies were adjusted. This was repeated until the modeled patterns-such as growth curves, age and size distributions, natural mortality levels, and fecundity measures-resembled the natural patterns observed for Atlantic cod. All model parameters are listed in Tables 1 and 3. Figure 1D shows sample time series of the population-level inter-annual noise e kEðtÞ affecting resource acquisition. The corresponding sample growth trajectories and annual fecundities are shown in Fig. 1E, while Fig. 1F shows population biomass B(t) and the number N 0 (t) of recruits. The average PMRN of the preharvest population is shown in Fig. 1A.
Before harvesting was started, the modeled population was allowed to reach a stochastic evolutionary and ecological equilibrium, so that all genetic traits and the correlations among them had converged close to an evolutionarily stable strategy. This was achieved by running the model for several hundred thousand years before saving a population from which the simulations were started. During each of 50 replicate runs with different random seeds, we let the model again equilibrate for 100 years before harvesting was started. In individualbased models, evolution is less influenced by genetic drift and individual-level stochasticity when a population's size is large. Because the populations we used were large (around 220 000 individuals), genetic drift and other forms of historical evolutionary contingency are less relevant for our results than they would be for smaller populations. Ideally one might have wanted to start simulation from several populations, but because of the long computing time needed for reaching the stochastic evolutionary equilibrium, we created only one such population. In all figures except for Fig. 2 the harvest period was 100 years. There was no fishing during the subsequent recovery period.

Nonevolutionary model
To assess the consequences of evolutionary change for recovery and rebuilding, we created a nonevolutionary version of our model in which the genetic traits were not allowed to evolve. For the initial comparison in Fig. 2, we did this simply by giving each individual born into the nonevolving population the genetic traits of a random individual that was alive at the time when fishing was initiated. We first constructed a library of about 200 000 individuals assembled in the last year of the long stabilization period without fishing; then, during fishing and the subsequent recovery, offspring were assigned genetic traits from this library. In this way, the distribution of genetic traits at birth was prevented from evolving under fishing. Nevertheless, the distribution of genetic traits later in life could change under fishing, through the differential vulnerability of genotypes to fishing. The distribution of phenotypic traits later in life changed under fishing for two further reasons: through the effects of fishing on density and thus on the phenotypic expression of the density-dependent growth coefficient, and through the knock-on effects of density-dependent growth and recruitment altered by fishing on a stock's length structure and thus on the differential vulnerability of phenotypes to fishing.
A related challenge arises because the ecological conditions of the evolutionary model differ from those in the nonevolutionary model as soon as the populations are fished. This prevents the comparison above from satisfactorily isolating the effects of evolutionary changes on the recovery process, because age, size, and maturity distributions, and hence density-dependent feedbacks on growth and recruitment, differ among populations starting their recovery from different initial biomasses. To Figure 2 Decrease and subsequent recovery of population biomass (of individuals aged 1 year and older) in dependence on harvest duration (increasing from left to right) and instantaneous harvest rate (increasing from top to bottom). Black curves: evolutionary model; gray curves: nonevolutionary model; gray shading: harvesting period. For the highest harvest rate and the longest harvest duration (lower right corner), the nonevolutionary population contains <100 individuals at its lowest population size, so that increasing the harvest rate further would lead to its extinction.

Fisheries-induced evolution and stock recovery
Enberg et al.
better isolate the effects of evolution on recovery (in Figs 4 and 5, as opposed to Fig. 2), we therefore used four steps to scale the stock biomass of the nonevolving population at the beginning of the recovery period to the corresponding level of the evolving population, while ensuring the former population's demographic and genetic composition matches its adjusted biomass. First, we determined the evolving population's biomass B(T) at the time T at which fishing was stopped. Second, we fished the nonevolving population with the same harvest rate and identified the year s just before it reached the evolving population's target biomass B(T). Third, we determined the nonevolving population's biomasses B N (s) and B N (s + 1) and separately stored all individuals it contained in those 2 years. Fourth, we randomly drew individuals from these two stored populations-complete with genetic traits and other individual traits such as length, age, and maturity status-until the population assembled in this way reached the target population biomass B(T). During this assembly, the probability that an individual was drawn from the population at time s was while the remainder was drawn from the population at time s + 1. In this way, we separately constructed a new nonevolving population for 50 replicate model runs.

Results
To quantify the impact of evolution on recovery, we compared responses in the biomasses of evolving and nonevolving populations to fishing and a subsequent moratorium (Fig. 2), investigated the corresponding evolution of genetic traits (Fig. 3), and analyzed the differential dynamics of evolving and nonevolving populations that were of equal biomass at the recovery's start (Fig. 4), with a special focus on identifying differences in recovery times (Fig. 5).
We found that increased harvest, in terms of intensity or duration, magnified the evolutionary response of the harvested stock. This can be seen by comparing time (E) age at maturation; (F) length at maturation; (G) length increment at age 2 years; (H) lengths at ages 3 years (below) and 10 years (above). The preharvest average population biomass, spawning stock biomass, and eggs per spawner are scaled to 1.

Fisheries-induced evolution and stock recovery
series of total biomass for the evolving population with its hypothetical nonevolving counterpart (Fig. 2). The evolutionary response involved adaptations in the genetic traits that allowed the evolving population to better withstand fishing during the harvested phase; these adaptations persisted after harvesting had ceased. Since the evolving population adapted to the harvest, its biomass began to rebound after the initial drop, whereas the nonevolving population declined monotonically as long as it was harvested.
One consequence of fisheries-induced evolution is that populations that adapted to the new mortality regime could withstand considerably higher fishing pressures than if evolution was not occurring. While the hypothetical nonevolving populations went extinct for harvest rates exceeding 0.5 year )1 , the evolving populations acquired a capacity for withstanding considerably higher harvest rates. However, even though fisheries-induced evolution appeared to make populations more resistant to extinction, the flip side of adaptation was observed when the harvest pressure was removed. Although starting recovery from a higher biomass, the evolving populations could not fully recover to the preharvest level within the recovery period of 350 years shown in Fig. 2. This lag effect was the more pronounced the more intense and prolonged a population's exposure to fishing had been. The capacity of evolving populations to withstand higher harvest pressure resulted from evolution of their genetic traits. The most prominent evolutionary change took place in the maturation schedule, and in the PMRN intercept in particular (Fig. 3A). The PMRN also evolved to become more steeply inclined during harvesting (Fig. 3B), the gonado-somatic index evolved to higher values (meaning that individuals invested progressively more of their resources into reproduction; Fig. 3C), and the growth coefficient evolved to lower values (Fig. 3D). The higher the harvest pressure, the larger the magnitude of the evolutionary response. In all genetic traits, recovery was slower than the preceding harvest-induced changes. Since the heritabilities of genetic traits did not change significantly over time, this implies that the natural selection pressures reverting the genetic traits during the moratorium were weaker than the preceding fisheries-induced selection pressures.
Harvesting changed a range of stock characteristics in the evolving and nonevolving populations, including those commonly used for quantifying biomass, recruitment, maturation, and growth (Fig. 4). Although the genetic growth coefficient evolved to lower values (Fig. 3D), the phenotypic growth rate increased because of relaxed density regulation (Fig. 4G). This finding underscores the importance of considering phenotypic plasticity and density dependence when modeling fisheries-induced evolution, or when interpreting observed empirical changes in the light of fisheries-induced evolution, as phenotypic plasticity and density dependence can both mask and exaggerate fisheries-induced changes in a stock's genetic composition.
As evolution changes a population's genetic traits, it also affects population dynamics and biomass variation (Fig. 4). A naïve comparison of evolving and nonevolving populations, such as in Fig. 2, thus cannot separate the effects of evolution from the effects of phenotypic plasticity and density dependence, because the recovery processes of the evolving and nonevolving populations start from different population biomasses after harvesting. We eliminated this confounding factor by rescaling the nonevolving population's biomass to the same level observed for the evolving population at the end of the harvesting period (see 'Nonevolutionary model' under 'Methods' above). Initiating the recovery of evolving and nonevolving populations from the same biomass revealed that fisheries-induced evolution had little influence on the recovery of population biomass during the moratorium's first 15 years (Fig. 4A). However, after this initial period, the evolving population took hundreds of years to fully recover to its preharvest biomass level. In contrast, the hypothetical nonevolving populations were fully recovered within about 50 years.
The recovery of some stock characteristics was faster for populations that had undergone fisheries-induced evolution. For example, spawning stock biomass increased more rapidly during the moratorium, because individuals were maturing earlier as a result of fisheries-induced evolution. Likewise, when harvesting ceased, spawning stock biomass in the evolving population exceeded its preharvest level, because a larger part of the population was mature and thus contributed to the spawning stock (Fig. 4B). At the same time, the number of eggs per spawner and the number of recruits per spawner decreased considerably due to fisheries-induced evolution, because mature individuals were on average smaller (Fig. 4C,D). Furthermore, the recovery of these metrics was extremely slow, and there was little short-term recovery. In contrast, the nonevolving populations showed a demographic increase in the number of eggs and recruits per spawner at the beginning of the recovery period, because the number of old and large individuals increased when harvesting was ceased. With time, density dependence began to kick in and reproduction thus fell to preharvest levels. The initial demographic increase and subsequent density-dependent decrease in reproductive traits were also present, albeit less pronounced, in the evolving populations.
The average age at maturation is one of the stock characteristics that responded most strongly to fishing. The reason for this is threefold: first, evolution of the PMRN increases the probability of maturation for smaller individuals at younger ages; second, released density dependence increases growth and thus allows fish to reach sizes at which they are more likely to mature early in life; and third, population-level averages of age at maturation are based on live fish and therefore are biased toward early maturation, because harvest removes fish that otherwise would have matured late in life. In Fig. 4, the evolving populations exhibited all three effects, whereas the nonevolving populations underwent only the latter two effects (Fig. 4E). Even though the average age at maturation decreased in the nonevolving populations, the decline was smaller and the recovery to the preharvest level was considerably quicker (Fig. 4E).
The average length at maturation decreased dramatically because of fisheries-induced evolution and recovered very slowly thereafter, whereas in the nonevolving populations it not only decreased less during fishing, but also rebounded to its preharvest level within just a few years Fisheries-induced evolution and stock recovery Enberg et al.
( Fig. 4F). For a PMRN with significantly negative slope, one would expect that faster growth leads to earlier maturation at larger size when everything else is unchanged. In the nonevolving populations, however, (i) average maturation length at the end of the fishing period lies below its preharvest value even though individuals grow faster during the fishing period and (ii) average maturation length increases during the early phase of recovery even though growth rates go down concomitantly. Both points seemingly contradict expectations based on the PMRN slope. The contradiction is only apparent, as late-and large-maturing individuals were more exposed to fishing, resulting in a bias toward fish that are smaller at maturation. When fishing pressure was released, this bias vanished, and the average length at maturation quickly increased. Thereafter, the average length at maturation slowly decreased again together with the decreasing average growth rate and thus in line with expectations based on the PMRN slope. The average length increment of 2-year-old individuals was only minimally influenced by fisheries-induced evolution (Fig. 4G). This early age was chosen to avoid the confounding effects resulting from maturation evolution and the associated changes in resource allocation to gonads. Although the growth coefficient evolved to lower values during harvesting (Fig. 3D), the length increment of 2-year-old individuals increased slightly and decreased relatively quickly to the preharvest level as population density increased. This weak effect of growth evolution on body length can also be observed at other ages: length at age 3 years decreased little because of evolution, and although length at age 10 years showed a dramatic response, this was mainly because of maturation evolution (Fig. 4H). In the nonevolving populations, the average length at age 10 years increased during the harvesting period because of the release of density dependence acting on growth. In contrast, fisheries-induced evolution led to a smaller length at age 10 years, because individuals evolved to mature earlier in life, so that their growth slowed down as part of their resources were used for reproduction (Fig. 1E). The recovery of the length at age 10 years was very slow in the evolved populations, reflecting the slow evolutionary recovery of the maturation schedule.
Increased harvest pressure led to longer recovery times. We quantified this as the time it took a population to reach 50%, 70%, and 90% of its preharvest levels of biomass, number of recruits, and average age at maturation (Fig. 5). Biomass recovery to 50% and 70% of the preharvest level was not strongly influenced by fisheriesinduced evolution, whereas recovery to 90% of the preharvest level took roughly 10 times longer when the population had adapted to fishing (Fig. 5A,B). Recruit-ment, on the other hand, recovered faster when the population had evolved: the fish matured earlier and made higher reproductive investments, and the resultant number of recruits was also larger. As a consequence, the recovery time to preharvest levels of recruitment was shorter for evolved populations than for nonevolved populations (Fig. 5C,D). Recovery of the average age at maturation was strongly affected by fisheries-induced evolutionary changes (Fig. 5E,F). Recovery to 70% of the preharvest maturation age took more than 300 years when harvest rates had exceeded 0.5 year )1 , and recovery to 90% of the preharvest maturation age could take more than a thousand years (Fig. 5E). In the hypothetical nonevolving populations, age at maturation recovered to 90% of its preharvest level within <40 years (Fig. 5F), and it never fell below 80% of its preharvest level (not shown).

Discussion
Our results show that evolutionary life-history changes caused by fishing do affect recovery, as has previously been suggested by some authors (Hutchings 2005;Walsh et al. 2006;Stenseth and Dunlop 2009). Below we discuss these impacts, explain their origin, consider limitations of our analysis, and highlight management implications.

Evolutionary impacts on recovery
A fish stock's recovery after a period of exploitation is a multifaceted process involving numerous traits, during which the different traits recover at different rates and to different levels. Even though the rebuilding of population biomass is not strongly influenced by fisheries-induced evolution during the first 10-15 years of recovery, it can take dramatically longer to reach preharvest levels of biomass when evolution has taken place. On the other hand, our model predictions suggest that some traits, such as spawning stock biomass and recruitment, recover faster after evolution, and that the adaptations to high mortality rates make a stock less prone to extinction during intense fishing.
Atlantic cod is an illustrative example of a species that has experienced periods of declines and recoveries over its whole distributional range (reviewed by Lilly et al. 2008). Having been fished down to fractions of their pristine biomass, some stocks appear to experience great difficulty in recovering from their depleted state. Most infamous is the collapse and nonrecovery of the northern cod off Newfoundland (e.g., Haedrich and Hamilton 2000;Shelton et al. 2006), but also the southern Gulf of St Lawrence cod and North Sea cod are at historical lows, with good news being few and far between. For southern Gulf of St Lawrence cod, the situation is so grave that this Enberg et al.
Fisheries-induced evolution and stock recovery stock has been predicted to be extirpated within mere decades (Swain and Chouinard 2008). All of these cod stocks have shown life-history changes that have partly been attributed to fisheries-induced evolution (Olsen et al. 2004;Yoneda and Wright 2004;Olsen et al. 2005;Swain et al. 2007; see also Pérez-Rodríguez et al. 2009 for Flemish cap cod). Similar evolutionary changes induced by fishing have taken place in numerous other fish populations worldwide (reviewed in e.g., Jørgensen et al. 2007;Kuparinen and Merilä 2007;Fenberg and Roy 2008;Hutchings and Fraser 2008). Before we can use our model to examine what implications such fisheries-induced evolutionary changes have for stock recovery, we need to ascertain that the predictions of harvest-induced changes from this model are consistent with empirical observations and theoretical results from other models. That fishing leads to earlier maturation has been suggested or shown in several models, with the earliest examples being Borisov (1978) and Law and Grey (1989), while the model with a methodology most similar to ours is . Our study predicts a drop in age at maturation of 3 years after 100 years of harvesting at a rate of 0.3 year )1 and of 5 years for a harvest rate of 0.7 year )1 . In comparison, the change in age at maturation observed in Northeast Arctic cod is 2.5 years over 70 years, under harvest rates of around 0.5 year )1 (Heino et al. 2002b). Further studies of the evolution of maturation age and size were reviewed by Dieckmann and Heino (2007), but there appear to be no other examples of late-maturing stocks with sufficiently longtime series to enable similar comparisons. However, the rate of change in Northeast Arctic cod, as measured in darwins, is within the range of rates observed in numerous other species worldwide (reviewed in Jørgensen et al. 2007). The concomitant evolution of increased reproductive investment is also supported by model results (Jørgensen et al. 2006;) and empirical observations (Yoneda and Wright 2004;Rijnsdorp et al. 2005;Wright 2005). Reduced growth rates in response to harvesting have been shown most unambiguously in pink salmon, where maturation age is constant (Ricker 1981(Ricker , 1995, and experimentally in Atlantic silversides (Conover and Munch 2002), but a weaker effect more in line with our model's predictions has been reported in data on plaice and cod (Rijnsdorp et al. 2005;Swain et al. 2007) and in other models (Favro et al. 1979;Brown et al. 2008;Hilborn and Minte-Vera 2008;). Based on these comparisons with theory and observations, we conclude that the harvestinduced changes our model predicts are within the expected range and therefore pertinent to considerations of empirical recovery processes.
The effect of evolution on observed phenotypes and population dynamics depends on the timescale of interest. At short and medium timescales (years to decades), the primary role of evolutionary trait changes is that they alter population dynamics and thereby the 'rules of the game' for stock recovery. Some stock characteristics recover faster, some slower, and some incompletely, because of the evolutionary changes that took place while the stock was harvested. On longer timescales (decades to centuries), the evolution of genetic traits back to their preharvest levels may take considerably longer, as reverse selection pressures will often be weak.
The interconnected processes and diverse dynamics of different stock characteristics paint a more nuanced picture than the words 'recovery' or 'rebuilding' suggest, especially when considering the way these concepts are used in the nonscientific literature. For example, the Johannesburg Declaration states that to achieve sustainable fisheries, the requirement is to 'maintain or restore stocks to levels that can produce the maximum sustainable yield with the aim of achieving these goals for depleted stocks on an urgent basis and where possible not later than 2015' (UN 2002). The text is, perhaps deliberately, unclear on how the 'levels that can produce the maximum sustainable yield' are to be defined. For fisheries scientists and managers devising rebuilding plans and monitoring recoveries, an operational definition of this objective is needed. Below, we point to some of the processes at work during recovery and rebuilding that impinge on this question, and show how expectations for recovery depend on the stock characteristic at focus and on the processes that led to depletion or overfishing.

Three processes with different timescales
Three processes determine the recovery of different stock characteristics. These processes act on top of each other, and each one of them dominates on different timescales. Distinguishing these processes helps us to understand the biological dynamics that influence stock recoveries from a depleted state.
First, when fishing is ceased, fish that previously would have been harvested will survive and influence population dynamics. This is observed as an increase in biomass, but also as rapid increases in recruits per spawner, mean age and size, and mean age and size at maturation. These effects result from the restoration of an age and size structure in which cohorts are dying off at a slower rate. This restoration is observed in terms of very rapid initial recoveries of most of the stock characteristics shown in Fig. 4.
Second, as population density and/or total biomass increase, density-dependent effects begin to alter the Fisheries-induced evolution and stock recovery Enberg et al.
phenotypic composition of the population by affecting, in particular, individual rates of growth and reproduction. It can be instructive to think of this effect as having two components: ageing individuals that survived until fishing ceased gradually find themselves in a denser population that decelerates their growth and reproduction, while new cohorts spawned under higher densities will grow slower also early in life, thus exhibiting different adult characteristics than their parents. Such density dependence and phenotypic plasticity cause several of the stock characteristics in Fig. 4 to show a slow decline after a fast increase: for biomass this is hardly visible, but for reproductive traits, age and size at maturation, and length at age, this decrease can be pronounced. For the number of recruits per spawner, the influence of increasing density dependence is even stronger than the initial increase because of the restoration of age and size structure. When monitoring a stock's recovery, an initial rapid and promising restoration of population structure can thus be overturned as population dynamics become dominated by density regulation. Third, the slowest process is evolution of the genetic traits toward their preharvest levels. Some of the genetic trait changes that were induced by fishing evolve much slower in the opposite direction. This was suggested by Law and Grey (1989), who noted that selection pressures toward early maturation during fishing can be very strong, as most late-maturing individuals die before they can reproduce, whereas many early-maturing individuals can reproduce at least once. When fishing is ceased, both early-and late-maturing individuals can reproduce, but the late-maturing phenotypes will produce slightly more offspring, as they follow fine-tuned resource allocation strategies that maximize their lifetime reproductive success under conditions of natural mortality. Although evolution by natural selection thus gradually moves the population toward later maturation, this process can be very slow. A first quantitative corroboration of this asymmetry in the context of a detailed eco-genetic model was provided by . Our results show that similar asymmetries between the rates of fisheries-induced evolution and reverse evolution when fishing is relaxed can occur also for other traits. In our model, it took natural selection thousands of years without fishing to undo genetic changes caused by only a century of harvesting.
Our model does not, however, support predictions by de Roos et al. (2006), who found that fisheries-induced evolution toward smaller maturation size could be irreversible and even continue despite a cessation of fishing. We believe their results could be a consequence of phenology in seasonal environments, where fitness valleys associated with annual cycles may prevent genetic traits from recovering. A decade-long experimental study by Conover et al. (2009) supports the reversibility of evolutionary changes in life-history traits caused by harvesting. Their study showed, however, that the rate of evolutionary recovery depended on the harvest regime: largeharvested populations, which evolved toward smaller body size during a selective-harvesting period of five generations, showed significant recovery during a subsequent nonharvesting period of five generations. In contrast, small-harvested populations, which evolved toward larger body size during selection, showed no significant evolutionary reversal after harvesting was ceased (Conover et al. 2009).

Model limitations
Our model population represents a long-lived and latematuring species exposed to a specific harvest regime. How well it can capture the essence of recovery and rebuilding dynamics in populations harvested differently (e.g., less size-selectively) or in populations with different life histories (e.g., short-lived and early-maturing species) remains to be assessed by future work. Our model also included only a limited number of evolving traits. In contrast, Walsh et al. (2006) found in laboratory experiments with Atlantic silversides that when large individuals were selectively harvested, a suite of physiological, behavioral, developmental, and life-history traits were altered, many of which can affect a stock's recovery potential. For example, egg volume was reduced, leading to smaller larval size at hatching and to lower viability of the larvae. Also, food consumption rate and conversion efficiency decreased. In addition, models have studied traits that are difficult to manipulate in the lab, including migration patterns (Jørgensen et al. 2008;Thériault et al. 2008) and skipped spawning (Jørgensen et al. 2006). Including more traits into a model might lead to different evolutionary responses during the harvesting phase, with other consequences for recovery dynamics when harvesting is ceased.
Four mechanisms that were not included in our model but are likely to affect recovery and rebuilding are parental effects, Allee effects, sexual selection, and trophic interactions. Younger and smaller females are often inferior spawners compared with their older and bigger counterparts, with a lower hatching rate of eggs and reduced offspring survival because of smaller egg size (Trippel 1998;Berkeley et al. 2004). Parental care has been shown to affect fisheries-induced maturation evolution in an eco-genetic model of smallmouth bass ). Allee effects may arise from intraspecific or interspecific interactions and may cause delays in recovery (Shelton and Healey 1999). Allee effects might also emerge from fisheries-induced changes in food-web structure (Van Leeuwen et al. 2008). Sexual selection has been hypothesized to influence the course and rate of fisheries-induced evolution (Hutchings and Rowe 2008a,b;Urbach and Cotton 2008), and the inclusion of males and of more elaborate mating structures in the model might indeed affect the ecological and evolutionary recovery process. Removal of top predators may have led to trophic cascades in formerly cod-dominated ecosystems in the Northwest Atlantic (Frank et al. 2005), and changes in ecosystem structure have also been suggested to contribute to the nonrecovery of the cod populations there. Our model ignores interspecific interactions, and as a consequence, the rate of population recovery observed in our study might differ from what could be expected if a whole ecosystem or a subset of it were modeled (Gårdmark et al. 2003).

Management implications of slow evolutionary recovery
The fact that reverse evolution is slower than fisheriesinduced evolution-in our model it may take 20-30 times longer to bring a trait back to its preharvest level-has been referred to as causing a 'Darwinian debt' we impose on our descendants through current fishing practices (U. Dieckmann in an interview with the Financial Times 29 August 2004, p. 1). If these asymmetric rates in our model are representative, we might have to accept effects of past fisheries-induced evolution as unavoidable characteristics of a new reality. Below we discuss three implications.
If fisheries-induced evolution is deemed undesirable, proactive management should attempt to prevent future fisheries-induced evolution (for considerations of fisheries-induced evolution caused by different gear types see Hutchings 2009;Jørgensen et al. 2009). Implementing such proactive management requires a more informed and rigorous understanding of potential trait changes and their effects in specific stocks. It would also require a management process that can identify undesirable outcomes and enact and enforce regulations that prevent them from happening. In reality, recovery planning is sometimes needed, and in these cases acknowledging the role of evolution can facilitate the setting of realistic goals. These goals may need to be specified for different time horizons, taking into account demography, phenotypic plasticity, and evolutionary change.
A second implication stems from practical limitations of fisheries management, in a world in which widespread bycatch, illegal landings, and high-grading cannot be disregarded. For example, in northern cod, fishing mortality went up when a fishing moratorium was implemented, apparently because of low stock size, a sentinel fishery, bycatch, and range contraction (e.g., Shelton et al. 2006).
The slow recovery of genetic traits discussed above results from much weaker selection pressures for reverse evolution, which implies that even relatively low fishing mortalities are still likely to outweigh those weak pressures and hence hinder recovery. Fishing regulations that engender reverse evolution may thus not be a practical option. Models are needed to address this question for more realistic fisheries scenarios including bycatch and unreported landings.
The third implication is a need for understanding how fisheries-induced evolution changes the biological dynamics of stocks. The pronounced differences between evolved and nonevolved populations found in our study illustrate how many stock characteristics are bound to change as a result of fisheries-induced evolution. This may limit the utility of older observations for managing current stock dynamics. In consequence, management targets and reference points may be changing continuously. To counteract this problem, monitoring programs could be modified to include other or more stock characteristics than traditional survey protocols suggest. Successfully tackling the three challenges outlined above would require evolutionary ecology to become more strongly integrated with fisheries science and management (Hutchings 2000b).
To summarize, our study underlines that the recovery and rebuilding of fish stocks are influenced by both ecological and evolutionary processes. Although evolution has little direct effect during initial recovery, its indirect influences are important, because the traits that evolved in response to fishing affect demography and phenotypic plasticity. In the longer term, evolution itself also plays an important role, as full evolutionary recovery to original trait values can be very slow or even impractical. The slow rate of evolutionary recovery is reflected in the rebuilding of population biomass, leading to incomplete biomass recovery on intermediate timescales.
and Technology Foundation. This manuscript was finalized during a research stay in the Mangel Lab at the University of California Santa Cruz by KE and CJ.