Mitigating fisheries-induced evolution in lacustrine brook charr (Salvelinus fontinalis) in southern Quebec, Canada

Size-selective mortality caused by fishing can impose strong selection on harvested fish populations, causing evolution in important life-history traits. Understanding and predicting harvest-induced evolutionary change can help maintain sustainable fisheries. We investigate the evolutionary sustainability of alternative management regimes for lacustrine brook charr (Salvelinus fontinalis) fisheries in southern Canada and aim to optimize these regimes with respect to the competing objectives of maximizing mean annual yield and minimizing evolutionary change in maturation schedules. Using a stochastic simulation model of brook charr populations consuming a dynamic resource, we investigate how harvesting affects brook charr maturation schedules. We show that when approximately 5% to 15% of the brook charr biomass is harvested, yields are high, and harvest-induced evolutionary changes remain small. Intensive harvesting (at approximately >15% of brook charr biomass) results in high average yields and little evolutionary change only when harvesting is restricted to brook charr larger than the size at 50% maturation probability at the age of 2 years. Otherwise, intensive harvesting lowers average yield and causes evolutionary change in the maturation schedule of brook charr. Our results indicate that intermediate harvesting efforts offer an acceptable compromise between avoiding harvest-induced evolutionary change and securing high average yields.


Introduction
Fishing profoundly affects the dynamics of fish populations and the ecological communities in which they are found. Although the demographic impacts of fishing on fish population dynamics are relatively well studied (e.g., Getz and Haight 1989;Jennings 2004), we are only beginning to appreciate the evolutionary consequences of intensive fishing, which arise when fishing mortality imposes strong selective pressures on the harvested fish populations (e.g., Law 2000;Olsen et al. 2004;Jørgensen et al. 2007;Kuparinen and Merilä 2007;Allendorf et al. 2008;Fenberg and Roy 2008;Heino and Dieckmann 2008;Heino and Dieckmann in press., and Hutchings and Fraser 2008). In particular, the size-selective removal of fish is likely to result in evolutionary changes in important life-history traits, such as the size at maturation, when such traits are heritable (e.g., Reznick and Endler 1982).
However, assessing the impact of different harvest regimes on the evolution of life-history traits poses a challenge, because owing to phenotypic plasticity the same genotypes often express different phenotypes depending on the environment an individual encounters. Indeed, observed changes in life-history traits can arise from a purely phenotypically plastic response to harvesting, rather than from genetic evolution (e.g., Nelson and Soule 1987;Rijnsdorp 1993). Thus, observed changes in maturation schedules resulting from harvesting could be within the range of phenotypes produced by the genes controlling those schedules. Maturation reaction norms (Heino et al. 2002) help disentangle the plastic and genetic impacts on maturation schedules (reviewed in Heino 2007 andHeino and. Much recent research has therefore focused on how harvesting affects the evolution of these maturation reaction norms (e.g., Ernande et al. 2004;Olsen et al. 2004;Dieckmann and Heino 2007;Dunlop et al. 2007;Dunlop et al. 2009a;Dunlop et al. 2009b;Enberg et al. 2009;Sharpe and Hendry 2009).
Lacustrine populations of the brook charr, Salvelinus fontinalis Mitchill, from the Canadian Shield represent a promising model system for studying the effects of selective harvesting on the evolution of life-history traits in fish populations. Previous modeling work on this species has shown that fishing can cause evolution in the migratory behavior in exploited populations (Thériault et al. 2008). Moreover, in a study comparing 17 populations, Magnan et al. (2005) found that charr from fished lakes mature significantly earlier than those from unfished lakes. As brook charr exhibit a short life cycle in these lakes, with lifespans of 3-7 years, even relatively recent harvesting could impose strong selective pressures on reproductive traits such as the size at maturation, leading to significant selection responses within a few generations. Genetic change caused by harvesting may thus explain the smaller sizes at maturation observed in fished lakes. However, previous studies on unfished brook charr populations found associations between early maturation and rapid growth rates (e.g., Hutchings 1993). It is also known that intraspecific competition from adults can depress juvenile growth rates, and therefore possibly delay maturation in fish populations (e.g., van Kooten et al. 2007). Hence, growth-mediated maturation plasticity might suffice to explain observed differences in maturation schedules between harvested and unharvested lakes.
To assess how fishing potentially impacts brook charr life history at the genetic level, we follow the example of earlier research (e.g., Olsen et al. 2004 andDunlop et al. 2007) by examining how fishing may cause evolution in the probabilistic maturation reaction norm (PMRN) for age and size at maturation. A PMRN for age and size at maturation describes the probability that an immature individual undergoes maturation within a given time interval, depending on its age and size (e.g., Heino et al. 2002).
Indeed, genetic change in the maturation reaction norm caused by fishing will often be more difficult to reverse by the simple cessation of fishing than growthmediated phenotypic plasticity in maturation (e.g., Grey 1989 andDunlop et al. 2009b). While fishery managers routinely seek to implement sustainable management regimes (e.g., Getz and Haight 1989;Fenichel et al. 2008), most such attempts currently do not explicitly address fisheries-induced evolution (but see Grey 1989 andHeino 1998).
Here we use the eco-genetic modeling approach Dunlop et al. 2009b) to describe the effects of harvesting on the genetic determinants of the PMRN in a harvested population of brook charr. Eco-genetic models seek to integrate key ecological processes, such as resource consumption and somatic growth, with an explicit treatment of changes in the distribution of genotypic traits.
Here we examine how such models could help develop advantageous management regimes for the brook charr fishery in Southern Quebec. We incorporate size-specific fishing mortality explicitly, and compare harvest regimes according to their mean annual yield, as well as to the amount of evolutionary change that they cause. Within this comparative framework, we seek to address how fisheries managers can regulate fishing effort to reduce future evolutionary change in fish populations, while also maintaining acceptably high annual yields. Although our model is tailored to brook charr in Southern Quebec, we believe that our approach and predictions will also be applicable to other fisheries.

Methods
Our model describes a size-structured population of brook charr consuming a dynamic biological resource. To study fisheries-induced maturation evolution in such a complex ecological setting, we use the eco-genetic modeling framework described in Dunlop et al. 2007 andDunlop et al. 2009b (see also, e.g., Thériault et al. 2008, and Dunlop et al. 2009a, and Enberg et al. 2009). Fish population dynamics are implemented in an individual-based model (e.g., DeAngelis and Mooij 2005), in which each brook charr i is characterized by its PMRN traits, sex S, and somatic mass W.
Somatic mass is divided into irreversible mass and reversible mass. An individual's irreversible, or structural, mass X is determined by components such as organ and skeletal tissue that cannot be starved away (de Roos and Persson 2001). In contrast, an individual's reversible mass is determined by energy reserves such as fat, other lipids, and gonadal tissue in mature individuals, that can be marshaled for basic metabolic functions during starvation and hence can be starved away. We partition reversible mass further into the mass of storage tissue, such as fat and other lipids, Y, and that of gonadal tissue, G (e.g., Broekhuizen et al. 1994;Persson et al. 1998;de Roos and Persson 2001). Thus, In males, gonadal mass G is interpreted as the amount of reversible mass expended on reproduction -for example, through the loss of somatic mass incurred by fighting other males over access to a redd, or spawning nest.
The model's dynamics are iterated on an annual time step. In each time step, the model cycles through all individuals to determine their fates, with processes occurring in the following sequence: 1 Reproduction 2 Natural mortality, resource consumption, and growth 3 When harvesting occurs, fishing mortality Model parameters were determined through analyzing data previously investigated in Magnan et al. (2005), calibrating the model to reproduce the observed data (Appendices A and B), and incorporating values available from the literature (Table 1).

Reproduction
The probability that an immature individual matures during a year is determined by its PMRN in conjunction with its age A and total body length l. An individual's body length l(X) ¼ pX u is allometrically determined by its irreversible mass X. Estimation of the parameters p and u is described in Appendix A.1. An individual's PMRN describes the length-and age-specific probabilities of maturation between one season and the next. We assume that these PMRNs have logistic shape (as illustrated in Fig. 1A-C) and that their widths are constant across ages (as illustrated in Fig. 1D), where d is the PMRN width, measuring the difference between lengths leading to 25% and 75% maturation probability. We assume that the length l p50,A at which the maturation probability equals 50% at age A can be described by the linear function ( Fig. 1D) with a PMRN intercept m and a PMRN slope r that are specific to each individual Dunlop et al. 2009b). Brook charr breed once per year, and the model assumes random mating between males and females, conditioned on gonadal mass. The number F of fertilized eggs in the population is proportional to the total gonadal mass of the female population, where n f,t is the total number of mature females in year t, G j is the gonadal mass of female j, and W 0 is the average mass of an egg (Appendix A.1). For each egg, the mother and father are drawn at random from mature males and females in the population. The probability that the egg comes from mature female i is a nonlinear function of the female's relative gonadal mass G i in the population, The nonlinearity in reproductive value induced by t f > 1 reflects the positive correlation between gonadal mass and body size, as well as the superior ability of larger females to ensure offspring survivorship by, for example, identifying superior redd sites or remaining on spawning grounds longer to find suitable mates (e.g., Blanchfield and Ridgway 1997). Moreover, larger females also produce larger eggs, which in turn improve offspring survival (Mann and Mills 1985;Sehgal and Toor 1991;Maruyama et al. 2003).
Similarly, the probability that mature male i fertilizes a given egg is a function of its structural and gonadal mass, X i + G i , relative to the mass of other mature males in the population. To reflect the common observation that body length is positively correlated with reproductive value in males (e.g., Power 1980), a male's structural mass is incorporated into calculating its probability of fertilization, where n m,t is the total number of mature males and t m determines how strongly a male's reproductive value increases with its relative irreversible mass.

Natural mortality
In any given year, the number U of newborns that survive to the growing season is related to the population's total egg production (Power and Power 1995).
Here we model the recruitment U according to a Beverton-Holt recruitment function (e.g., Beverton and Holt 1957), After the surplus newborns have died, individual survivorship from one year to the next is assumed to be related to total body mass W and condition Y/W at the beginning of the year, Mitigating fisheries-induced evolution in lacustrine brook charr Okamoto et al.
Prðindividual i survivesjX; This functional form describes an exponential increase in survival with improved condition, and a sigmoidal increase in survival with larger body mass, and thus captures the type-III survival relationships consistently found for salmonids (e.g., Power 1980). The estimation of parameters b 1 and b 2 is described in Appendix A.1, while the value of b s was taken from the literature (Table 1).

Resource consumption
To describe resource competition among brook charr, we focused on zooplankton, one of the brook charr's major prey in some situations, such as when competitors (e.g., the white sucker Catostomus commersoni) are present, or when benthic organisms are rare (e.g., Magnan 1988;Tremblay andMagnan 1991, andMagnan et al. 2005). We also focused on zooplankton because they were the prey item for which sufficient data were available to parameterize our model. The resource-consumption rate E g,i,s of individual i at time s during a year is a function of its body weight W s and of the resource (zooplankton) number R s , where h(R s ) describes the proportion by which an individual's resource consumption is diminished when the resource density R s falls short of the maximum daily consumption rate, based on a type-II functional response (e.g., Koski and Johnson 2002), where J 1 is the maximum number of zooplankton that can be eaten in a day, J 2 is the half saturation constant, and L is the volume of the habitat shared by zooplankton and fish. We assume that the resource changes according to semi-chemostat dynamics (e.g., Claessen et al. 2000) and predation by charr, where k is the product of the average weight of an individual zooplankton (e.g., Gamble et al. 2006) and a conversion factor that takes into account assimilation efficiency (e.g., Claessen et al. 2000), E g,i,0 is the consumption rate of zooplankton by individual i at the beginning of the year, and n t is the population size of brook charr at the beginning of year t. We assume that the resource's population dynamics are much faster than the charr's population dynamics and that the impact of the charr population on the resource is roughly constant throughout the year. Thus, in each year t the resource density quickly attains an equilibrium, where d describes the brook charr's metabolic rate. Thus, we assume that when the mass derived from consumption exceeds metabolic costs, the surplus mass is invested in somatic growth. Since R s is assumed to be constant throughout the year (once it has equilibrated to R t ) and since c < 1, the growth increment DW(R t ) from year t to year t + 1 can be obtained as where f is a normally distributed random variable with mean 0 and standard deviation r f , and n ¼ exp½Àð1 À cÞdT ð 15Þ and gðR t Þ 1Àc ¼ hðR t Þad À1 ð1 À nÞ ð 16Þ (Hiyama and Kitahara 1993), where T is the length of the year. Individual brook charr are assumed to be born with the maximum ratio q J between reversible mass and irreversible mass. As long as resources are not limiting, juvenile brook charr will maintain this ratio of reversible to irreversible mass. The fraction F A (X t ,Y t ) of DW(R t ) that is allocated to irreversible mass is an empirically derived function of irreversible and reversible mass at the beginning of the growing season (Appendix C).
In mature individuals, a fixed fraction of (1)F A (X t ,Y t ))DW(R t ) is set aside for reproduction. In females, this takes the form of allocation to gonadal mass G, while in males this includes, for example, the loss in mass associated with searching and securing redds. For simplicity, we assume that gonadal mass in juveniles is negligible. Thus, at the end of the growing season, reversible, irreversible, and gonadal mass are given by where I f equals 0 if the individual is immature, and equals 1 if the individual is mature, and [x] + ¼ max(x,0). The growth equations above are based on the assumption that all gonadal weight from the previous year is spent on reproduction. The estimation of parameters of the growth model is described in Appendix A.1. A Mann-Whitney test indicated a good fit between the model-predicted size distribution and the observed size-frequency data (Appendix A.2), suggesting that the model's structural assumptions and parameters are appropriate for describing the modeled populations of brook charr.

Genetics
We adopt the eco-genetic modeling approach, and hence its reliance on the principles of quantitative genetics, to characterize the evolution of the PMRN as a result of harvesting in our simulated brook charr populations Dunlop et al. 2009b). We model both males and females, with sex being determined randomly assuming an even sex ratio at birth and in the initial population.

Trait inheritance
An individual's expressed size at maturation is a function of genetic and environmental effects. Each individual charr i possesses two quantitative traits that determine its PMRN, the PMRN slope r i and the PMRN intercept m i . These traits undergo diploid inheritance and expression through a sequence of two steps.
First, each parent produces a haploid gamete that is envisaged to contain, at random, half the parental alleles. As these alleles are not explicitly modeled in quantitative genetics, deviations in gametic genetic values from parental genetic values as a result of the combined effects of recombination, segregation, and mutation are described by random deviates. The genetic values of a gamete produced by parent i thus equal r i + . r and m i + . m , where . r is randomly drawn from the normal distribution N(0,(zr i ) 2 ) with probability l and equals r i with probability 1 ) l, where l describes the probability that the gamete's genetic value is recognizably different from the parent's genetic value; . m is drawn analogously. Although z, the coefficient of variation of the gametic genetic value from the parental genetic value is assumed to be constant throughout time, the magnitude of the recombinationsegregation-mutation effect varies as the parental genetic values evolve.
Second, the offspring's genetic values are obtained as the midparental values resulting from the union of two gametes, given by the arithmetic means of the gamete's genetic values. When these midparental values are phenotypically expressed, environmental variance is added in the form of a normally distributed random deviate with a mean of 0 and a standard deviation equal to a fraction k of the midparental value. The coefficient k of environmental variation was chosen to yield heritabilities of 0.15, consistent with common observations for numerous lifehistory traits (e.g., Mousseau and Roff 1987;Roff 1997). Values for z, l, and k are given in Table 2.

Initial genetic structure
We chose the initial population-level means r 0 and m 0 (Table 1) of the two genetic PMRN traits in accordance with the average PMRN of the modeled brook charr populations in unharvested lakes, using the estimation method described by Barot et al. (2004) and matching the estimation results to a linear PMRN with constant width (Appendix B). Figure 1 shows the PMRN thus obtained. The effect of age on maturation probability is relatively small compared to the effect of body length, indicating that the latter is the dominant indicator of maturation probability for brook charr from unharvested lakes.
To describe genetic variation in the PMRN traits,which allows selection to occur with and without harvesting, we initialized populations with combinations of the two genetic traits with variances (Cr) 2 and (Cm) 2 around the population's initial mean PMRN slope and intercept, r 0 and m 0 , respectively. The standard deviations of these two normal random distributions were set to a fifth of the respective mean genetic trait values (C ¼ 0.2; Table 2) to ensure that the initial genetic coefficients of variation (e.g., Houle 1992 andBürger 2000) were of an order of magnitude that is comparable to values reported for lifehistory traits in Houle (1992). These genetic coefficients of variation determine a population's ability to respond to selection. The values given in Houle (1992) were doubled to reduce the effect of the initial distribution of breeding values on subsequent evolutionary trajectories. Moreover, by increasing the genetic coefficient of variation, we also sought to ensure that by the time fishing began, there was still sufficient genetic variation in the PMRN traits on which selection could act. Most importantly, when fishing began, the genetic coefficient of variation in the model (mean coefficient of variation of 14% with a standard deviation of 2% for both PMRN traits) approached values reported in Houle (1992) for other life-history traits.

Genetic assumptions
The model choices described above are consistent with empirical work on the genetics of life-history traits in other taxa. For instance, life-history traits are believed to generally have low additive genetic variances (e.g., Mousseau and Roff 1987;Price andSchluter 1991, andBürger 2000). Thus, we assume the proportion of genetic variance in a particular trait's genetic value attributable to the effects of recombination and segregation to be small. Moreover, existing research suggests that most mutations have negligible effects (e.g., Lynch and Walsh 1998), indicating that the contribution of mutation to genetic variation is restricted. Most mutational variance appears to be attributable to a few mutations of large effect in, for example, Drosophila (e.g., Bürger 2000 and references therein). Given this limited potential for genetic differences between parental genetic values and gametic genetic values for lifehistory traits to arise through recombination, mutation, and segregation, we introduced the parameter l to control for how frequently one can expect the combined effect of these processes to produce gametic genetic values of a trait that are recognizably distinct from the corresponding parental genetic value.
Our model also assumes no interaction effects and free recombination between the loci controlling the PMRN intercept and slope, even though some degree of pleiotropy or linkage may well exist between these quantitative traits. In addition, genetic correlations and genetic or developmental constraints could limit the degree of evolutionary change in the modeled brook charr populations. However, in the absence of studies revealing the degrees of genetic interactions and genetic correlations between these two traits, we feel that the case of no interactions and free recombination provides a basis on which to develop further work examining the effects of relaxing these simplifying assumptions by investigating more complicated genetic settings. Indeed, our model, phrased in terms of individual quantitative traits, is formally equivalent to a single-locus, infinite-allele model for each trait. On that basis, one could readily extend our analyses to include complications such as pleiotropy and linkage between the PMRN traits.
For simplicity, we also assumed that life-history traits other than the two PMRN traits -such as offspring size, investment in gonads, and rates of somatic growth -are not subject to significant evolutionary change over the modeled time frame. This could apply because of low selection pressures or low genetic variation, with the former assumption being supported by results in Dunlop et al. (2009b), Dunlop et al. (2009a)), and Enberg et al. (2009).

Harvesting
The goal of this study is to identify harvest regimes that best mitigate harvest-induced evolutionary change. We compare a range of potential harvest regimes, and characterize each by three parameters, a, M s , and b H . The first two parameters describe the size-selectivity of harvesting. Specifically, the selectivity to which a fish of mass W is subjected is assumed to be given by I lðWÞ ða; að1ÀMsÞ Ms Þ, where I x (a,b) describes the cumulative distribution function in x of a beta distribution with shape parameters (a, b) (Fig. 2) and M s describes the length, as a fraction of the maximum observed length for brook charr (700 mm), at which selectivity equals 50%. For comparison, the length at which maturation probability equals 50% for 2-year-olds is approximately 30% of the maximum length in the population. The parameter a controls the steepness of the selectivity curve around M s (Fig. 2), and thus, in effect, how size-selective the harvest regime is. For example, if the 50th percentile of the size-selectivity function I lðWÞ ða; að1 À M s Þ Ms Þ is interpreted as the minimum catch size, a can be interpreted as the stringency with which the minimum catch size is enforced. The third parameter, b H , governs the density dependence of harvesting. In particular, the total allowable catch in year t is assumed to be given by (e.g., Hilborn and Walters 1992), where is the total harvestable biomass in year t, with the sum extending over all n t brook charr in year t and a H is the total allowable catch for H t ¼ 0. The small non-zero constant a H describes rare poaching and by-catch mortality (e.g., Hilborn and Walters 1992). We assumed that the fraction of biomass harvested when b H ¼ 0 was minimal, and therefore a H was set equal to 50 g, the approximate mass of a single charr at full condition with a maturation probability of 50%. For virtually all harvest regimes considered here, this meant that the density-dependent component of harvesting, b H , approximately equaled the harvest probability, i.e., the fraction of the total harvestable biomass H t designated as total allowable catch Y t in year t. The parameters used to characterize the harvest regimes are summarized in Table 3. The fishing season is assumed to run concurrently with the growing season, but the probability that an individual charr is harvested depends on its size at the beginning of the fishing season. Thus, the annual probability l H with which a brook charr at weight W i is harvested is given by Figure 2 Illustration of different selectivity curves describing how the exposure to fishing varies with the length of fish. The horizontal axis shows the length of fish as a fraction of the maximum observed length for brook charr. The vertical axis shows the probability that an individual charr of a given length will be harvested. The three shown curves correspond to different choices of a, which controls the steepness of the selectivity curve. In all three cases, M s ¼ 0.5, so that selectivity equals 50% for fish possessing half the maximum length.
We highlight that this probabilistic treatment of harvesting, when applied to a finite brook charr population, implies sampling variation in annual yield Y t .

Evaluation of harvest regimes
Choosing suitable time horizons Each model run was initiated with 1500 individuals. To evaluate the effect and desirability of different harvest regimes, the model described above was run for 100 years without harvesting (e.g., Tenhumberg et al. 2004) to allow the population to reach a demographic steady state and allow the build-up of an endogenous genetic structure. In the 101st year of model runs, harvesting began, and the brook charr population became subject to the harvesting mortality l H . The model runs proceeded for another 50 years, or until the brook charr population went extinct. Thus, we investigated a management time frame of 50 years during which the same harvest regime was consistently applied. We focused on combinations (a,M s ,b H ) of parameters of the harvest regime that did not lead to deterministic extinction of the brook charr population during the first 50 years of harvesting. Focusing on a 50-year time horizon for harvesting is desirable for providing insights to decision makers: the model is parameterized to reflect the current life history of the brook charr populations and aims at providing decision makers with relevant information about the effects of different harvest regimes within a time frame that can be deemed relevant for current management.
Although the current genetic distributions of the PMRN coefficients in brook charr populations are unknown, using the distribution of genotypes after running the model for 100 years can be justified on a number of grounds. First, the coefficients of genetic variation for the PMRN coefficients resembled values of the genetic coefficients of variation observed for life-history traits in other taxa (e.g., Houle 1992 andHoule 1996). Second, the distribution of genotypes after 100 years was unimodal, approximately symmetric, and qualitatively similar to the distribution of observed life-history traits in other taxa, as well as to the normal distribution of traits assumed in many other studies (e.g., Houle 1996;Lynch and Walsh 1998;Bürger 2000). Finally, the mean trait values after the model was run without harvesting for 100 years, or even for 15 000 years, were not significantly different from values measured in the field. Thus, we feel the distribution of genotypes after 100 years provided a reasonable approximation of natural genetic variation in the brook charr populations.
Indeed, an earlier study that also examined the impact of different management regimes on harvest-induced evolution found that, for an explicit multi-locus model, simulating 100 years without harvesting allowed the distribution of allele frequencies to stabilize (Tenhumberg et al. 2004). Because our model does not involve complications that result from explicit multi-locus genetics (such as linkage disequilibrium between loci resulting from the model's initialization), we decided that running the model for 100 years prior to harvesting was sufficient to lessen the effects of transient population dynamics and of the initial genetic distribution of PMRN traits.
To assess the plausibility of the model's evolutionary equilibrium, we also ran 100 replicates of the model without harvesting for 15 000 years, i.e., for the approximate number of years since the brook charr colonized the study area (e.g., Angers and Bernatchez 1998). The average values of the predicted PMRN traits after 15 000 years without harvesting (283 mm and 7.15 mm year )1 for the PMRN intercept and slope, respectively) were well within two standard deviations of the PMRN coefficients estimated for unharvested populations in the field. We therefore conclude that our model evolved to reasonable values for unharvested brook charr.
Because it could possibly take longer than 15 000 years for the PMRN traits to reach evolutionary equilibrium, we sought to assess how quickly evolutionary change was occurring in the PMRN traits in the absence of harvesting. For the 100 replicates of the model that were run without harvesting for 15 000 years, we fit an exponential decay model (e.g., Ritz and Streibig 2005) to the difference between the PMRN traits in a given year and the PMRN traits after 15 000 years. The resultant fits suggest a deterministic rate of relative evolutionary change, in the absence of harvesting, of 0.08% per year for the PMRN intercept and of 0.03% per year for the PMRN slope.
Assessing harvest-induced evolution under residual trends and genetic drift Residual deterministic trends and stochastic genetic drift affecting the evolving traits under natural selection cause uncertainty in assessing the evolutionary effects of harvesting. In particular, genetic drift could predominate when intensive harvesting strongly reduces population abundance. Even in large marine stocks, the effective population size of exploited fish stocks can be several orders of magnitude smaller than their census population size (e.g., Hauser et al. 2002). The potential for genetic drift resulting from harvesting is exacerbated in freshwater stocks, where census population sizes are considerably smaller than in marine stocks even in the absence of harvesting. For instance, Fraser et al. (2004) used microsatellites to estimate the number of breeders in seven populations of brook charr, and found that results ranged between 57 and 200 individuals in six of the seven populations. Because the brook charr populations considered in our study inhabit relatively small lakes, with an average surface area of approximately 13 ha, these stocks have relatively low population sizes compared to many marine stocks. As harvesting can further diminish these population sizes, it is desirable to adopt a probabilistic framework for comparing evolutionary trends in harvested and unharvested populations. Indeed, quantifying the probability that a particular harvest regime causes evolutionary change addresses this inherent uncertainty. Comparing the magnitude of evolutionary changes under a particular harvest regime to the average magnitude of changes expected in the absence of harvesting provides one approach to achieving this.
To quantify how much evolutionary change we expect in the absence of harvesting, we ran 2500 replicates of our model for the full 150 years without harvesting. For each replicate model run without harvesting, we calculated the differences in the population's mean values of the PMRN traits between the end (year 150) and the beginning (year 0) of the simulation. We thus determined the empirical distribution of the amount of evolutionary change that would occur in the PMRN traits over 150 years in the absence of harvesting (Fig. 3).
For a particular harvest regime H, we evaluated in each model run the differences in the population's means of the evolving traits between the simulation's end (r H,150 and m H,150 for the PMRN intercept and slope, respectively) and its beginning (r H,0 and m H,0 for the PMRN intercept and slope, respectively). We then calculated their two-sided empirical P-values, that is, the probability that the magnitude of evolutionary change in the absence of harvesting would be at least as large as the predicted magnitude of evolutionary change resulting from a particular harvest regime. These P-values were calculated by comparing the changes in the population means of the PMRN traits for a particular harvest regime with the distribution of such changes without harvesting. The measured amount of evolutionary change a harvest regime caused in a given model run was characterized by the smaller of the two empirical P-values thus obtained for changes in the two PMRN traits. The goal here was not to categorically accept or reject the null hypothesis of no evolutionary change. Rather, we sought to use the smaller of the two empirical P-values to quantify the probability of a given harvest regime causing evolutionary change as extreme as or more extreme than would be expected in the absence of harvesting. Thus, if r W,t and m W,t describe the mean PMRN traits in populations Based on this construction, 1 ) E can be interpreted as the confidence with which an observer can conclude that harvest-induced evolution is occurring. The probabilities in the expression above were obtained by integrating the relative areas within the corresponding tails of the frequency distributions in Fig. 3. We used a large sample of 2500 model runs to adequately characterize the distribution of genotypes in the absence of harvesting. However, for each harvest regime, we were not ultimately interested in the distribution of genetic values, but rather in their means. Trial and error indicated that about 15 replicate model runs were sufficient to estimate these means consistently. For each harvest regime, we also evaluated the mean annual yield during the years in which harvesting took place. We ran 15 replicates of each harvest regime, and characterized the evolutionary change caused by a harvest regime as the mean of E across the 15 replicates. Similarly, the yield of a given harvest regime was estimated by taking the average of the mean annual yield across the 15 replicates. In this manner, our assessment of the evolutionary effect of a harvest regime accounts not only for residual evolutionary trends potentially occurring in the absence of harvesting, but also for uncertainty in predicted evolutionary outcomes in the absence and presence of harvesting.

Evaluating the desirability of harvest regimes
We expect fisheries managers and other stakeholders to vary in their concern for the risk of fisheries-induced evolution in a particular fishery. We therefore introduce a parameter u to describe the relative importance that the avoidance of harvest-induced evolution has to management decision making. Low values of u describe a situation in which managers are comfortable implementing a harvest regime that risks significant evolutionary change, whereas high values of u magnify the relative importance of evolutionary change to management decision making, and thus describe situations in which managers or other stakeholders are averse to inducing evolutionary change.
On this basis, we characterize the desirability of a particular harvest regime by the product of mean annual yield with E u . Consequently, regimes that cause a lot of evolutionary change (corresponding to consistently low E) have their mean annual yields penalized considerably relative to regimes that cause little evolutionary change. Our goal is to identify those regimes that maintain a high yield while simultaneously limiting the amount of harvest-induced evolutionary change.

Results
Evolutionary changes induced by harvest regimes Figure 4 shows how evolutionary changes varied with the three parameters a, M s , and b H characterizing the harvest regimes. Qualitatively, the evolutionary changes predicted by the model were independent of the steepness a of the selectivity curve. Evolutionary changes were largest when harvesting was intense (b H > 0.15, where b H scales the density-dependent component of harvesting) and smaller brook charr were exposed to fishing (M s < 0.30, where M s is the fraction of the maximum observed length at which selectivity equals 50%). These results show that when harvesting is intense (b H > 0.15), harvesting brook charr below the size at 50% maturation probability for 2-year-old brook charr will almost certainly result in significant evolutionary change. Even when the size-selectivity of harvesting was extremely diffuse (a ¼ 0.5), intense harvesting extending to small sizes caused evolutionary change. As expected, light harvesting (b H < 0.05) caused the least amount of evolutionary change. Evolutionary change was also mitigated rather well by restricting the harvesting effort to brook charr larger than the size at 50% maturation probability for 2-year-old brook charr.
Annual yields achieved by harvest regimes Figure 5 shows how mean annual yields varied with the three harvest-regime parameters a, M s , and b H . Less intense harvesting (b H < 0.05) resulted in poor mean annual yields, while intensive harvesting, even when restricted to larger brook charr, typically improved mean annual yields. It is worth pointing out that intensive harvesting of smaller brook charr generally reduced the mean annual yield. Such regimes had the effect of producing a very large yield on average in the first year they were applied, but subsequent yields suffered, as the brook charr population recovered slowly. Examples of this effect are shown in Fig. 6. Figure 7 exemplifies how the intensive harvesting of smaller fish, which are less likely to be mature, can immediately have a pronounced effect on brook charr population dynamics even before a strong evolutionary effect is induced. We also ran 40 additional replicates of the model for the highest harvesting intensity (b H ¼ 0.2), to determine the size of the brook charr population during the course of harvesting. Figure 7 shows that the number of spawners during the last 10 years of harvesting (years 140-150) dropped to an average of only about 41 individuals, which illustrates that harvesting intensity could not be increased further without driving the population extinct. In contrast, harvesting brook charr that are quite likely mature (i.e., M s > 0.30, with the latter value corresponding to the length at 50% maturation probability for 2-year-old brook charr) ensured a continuous supply of recruits to maintain relatively high yields. Figure 8 shows how the desirabilities of harvest regimes varied with the three harvest-regime parameters a, M s , and b H . Here, the desirability of each harvest regime was measured by its mean annual yield weighted by E (implying u, measuring the relative importance that harvestinduced evolutionary change has to management decision making, ¼1). Results are comparable to those for annual yields (Fig. 5), except for two effects: first, there is less disparity between the different harvest regimes once evolutionary change is taken into account, and second, exposing smaller brook charr to high harvest intensities is even less desirable than an evaluation based on yield alone suggests. The high average yields obtained when intensive fishing is restricted to larger brook charr, in conjunction with their relatively weak evolutionary impacts, meant that these regimes were frequently identified as desirable, although when harvesting intensity was intermediate (0.05 < b H < 0.15), the size at 50% selectivity had little discernible effect on the desirability of a harvest regime. Finally, the steepness of the selectivity curve, measured by a, had a statistically significant but minimal effect on the desirability of the harvest regimes (Spearman's rank correlation coefficient, calculated based on 1000 bootstrapped replicates, equaled 0.03 with P < 0.01). Figure 9 shows how the desirabilities of harvest regimes varied with u, which measures the relative importance  This held true for managers ascribing differing levels of importance to avoiding harvest-induced evolution. Nevertheless, Fig. 9D shows that when avoiding harvestinduced evolutionary change is an important consideration (u ¼ 2) for managers, the range of intermediate harvesting intensities that produce desirable outcomes shrinks.

Discussion
Many theoretical and empirical studies have supported the notion that harvesting natural fish populations can induce evolutionary changes in key life-history traits, especially in traits governing maturation (e.g., Law 1991;Heino 1998;Law 2000;Ernande et al. 2004;Olsen et al. 2004;Jørgensen et al. 2007;Kuparinen and Merilä 2007;Fenberg and Roy 2008;Heino and Dieckmann 2008;Hutchings and Fraser 2008;Heino and Dieckmann in press.). Geared to the management of brook charr in Canadian lakes, this study assessed how harvest regimes can best be designed to respect the competing objectives of maximizing mean annual yield and minimizing harvest-induced maturation evolution.
The magnitudes of harvest-induced evolutionary changes in probabilistic maturation reaction norms (PMRNs) found  in our study were broadly consistent with qualitative predictions made in the existing literature on harvest-induced maturation evolution. In particular, as in Ernande et al. (2004); Dunlop et al. (2007), and Dunlop et al. (2009b), we found that harvesting resulted in evolutionary shifts of PMRN traits that caused earlier maturation at smaller size. Our results suggest that intense harvesting of brook charr extending below the size at 50% maturation of 2-year-old fish readily causes evolutionary changes in the PMRN. As we systematically averaged results across replicate model runs and extensively sampled the range of potential harvest regimes, we believe that these results are robust. In particular, were random genetic drift, rather than directional selection, a major driver of the modelpredicted evolutionary changes, differences between harvested and unharvested populations would not exhibit directional trends. The directional evolutionary trends we found in PMRN traits under a range of intense harvest regimes and across many replicate model runs therefore suggest that genetic drift is not a primary cause of the described harvest-induced evolution of PMRN traits.
Our results suggest that to reliably avoid harvestinduced maturation evolution in lacustrine brook charr, harvest regimes should be designed to limit the harvest to less than approximately 15% of the population's biomass, or restrict harvesting to individuals larger than 50% of the observed maximum length. Moreover, intensive harvesting extending to smaller brook charr proved undesirable from the perspectives of minimizing harvest-induced evolution and maximizing annual yield. One feature of our results that underscores the trade-off inherent in the choice of harvest regimes is that, regardless of the harvest's size-selectivity, less intense harvesting caused minimal evolutionary change (Fig. 4), but also resulted in low yields (Fig. 6). The trade-off between annual yield and evolutionary change is aggravated when managers (or other stakeholders) are very concerned about harvestinduced evolution (i.e., when u, measuring the relative Figure 9 Variation in the desirabilities of different harvest regimes. The three panels correspond to increasing values of u, which describes the importance managers attach to avoiding harvest-induced evolutionary change: (A) u ¼ 0.1, (B) u ¼ 0.5, (C) u ¼ 1.5, and (D) u ¼ 2. The steepness of selectivity curves is constant across panels, a ¼ 5.9. The vertical axis shows the length at 50% selectivity, measured as a fraction of the charr's maximum length, while b H , which scales the density-dependent component of harvesting, is shown along the horizontal axis. The desirabilities of different harvest regimes are color-coded, with darker blue corresponding to lower desirabilities, and thus to less successful harvest regimes, while red indicates high desirabilities. Values shown for each harvest regime are means of 15 replicate model runs.  E (equation 21), i.e., by the smaller of the two P-values measuring the probability of harvest-induced evolutionary change in the two probabilistic maturation reaction norm traits. The three panels correspond to selectivity curves with increasing steepness: (A) a ¼ 0.5, (B) a ¼ 5.9, and (C) a ¼ 25. The vertical axis shows the length at 50% selectivity, measured as a fraction of the charr's maximum length, while b H , which scales the density-dependent component of harvesting, is shown along the horizontal axis. The desirabilities of different harvest regimes are color-coded, with dark blue corresponding to lower desirabilities, and thus to less successful harvest regimes, while red indicates high desirabilities. Values shown for each harvest regime are means of 15 replicate model runs.
importance of avoiding harvest-induced evolution to management decision making, is greater than 1; Fig. 9C,D). Conversely, the trade-off is relaxed when managers are relatively unconcerned about harvest-induced evolution (u < 1, Fig. 9A,B). In the latter case, risking harvest-induced evolution by intensively harvesting smaller individuals becomes an acceptable outcome.
Our results also underscore that from the perspective of maximizing mean annual yield, intensively harvesting immature brook charr is, in any event, undesirable. Low yields need not be a consequence of the evolutionary effects of harvesting, as models that do not include evolution also show that increased juvenile mortality can result in much reduced population biomass (e.g., Chesson 1998) and, hence, in reduced yields. Our model accounts for the length-dependence of maturation probabilities, so that, as in reality, harvesting ever smaller individuals implies that more juveniles are harvested. Even in the absence of any evolutionary effects of harvesting, yields can be lowered through recruitment overfishing, i.e., through the demographic consequences of increased juvenile mortality. Indeed, it is widely understood that additional mortality imposed on small, immature fish can have strongly negative effects on the sustainability of fisheries, and therefore ultimately also on yields (e.g., Beverton and Holt 1957;Getz and Haight 1989;Beck et al. 2001;Roberts et al. 2005). Hence, there are good reasons, apart from the risk of harvest-induced evolution, to eschew harvest regimes that intensely exploit small, immature fish. It is therefore important that our approach to evaluating the desirability of alternative harvest regimes accounted for the demographic as well as the evolutionary effects of harvesting.
We chose to commence harvesting after running our model for 100 years without harvesting, and we initialized the model with the PMRN traits observed in the field. As, in the model as well as in nature, residual evolutionary transients may exist as a result of natural directional selection pressures that have not yet had enough time to run their course, and as resultant trends in genetic traits may thus be superimposed on the response of populations to harvesting, we devised a general probabilistic approach to quantifying the likelihood of harvest-induced evolutionary change. We achieved this by comparing the magnitudes of evolutionary changes in unharvested populations with those in harvested populations across multiple model runs, so as to estimate the probabilities, separately for each evolving trait, that the former exceed the latter. We then took the smallest of these probabilities to quantify the confidence (Equation 21) with which an observer may conclude that the evolutionary changes encountered during the harvest period are merely the consequence of genetic drift and/or residual evolutionary transients. In this manner, our probabilistic approach can also deal with environmental trends, such as with those implied by climate change, and their evolutionary consequences. We therefore expect this approach to be applicable under a broad range of management scenarios.
While we have assumed throughout our analysis that managers will seek to avoid harvest-induced maturation evolution, some fisheries managers may regard such evolution as desirable under certain circumstances. For example, harvest regimes could potentially select for increased growth rates (e.g., Conover and Munch 2002), or fisheries may target mature individuals to select for later maturation at larger size (e.g., Heino 1998). In some cases, such change may be considered beneficial if it enhances a stock's resilience to harvesting in the short term (e.g., Enberg et al. 2009). We expect our approach to be applicable to situations in which harvest-induced evolution is valued differently. Technically, this may be achieved as easily as by replacing E with 1 ) E in the definition of a harvesting regime's desirability.
The model used in this study did not consider gene flow from unharvested populations, which can counteract local adaptation in harvested lakes towards smaller maturation size (e.g., Kawecki and Ebert 2004). Gene flow between oligotrophic lakes in a boreal forest is probably quite limited on the short temporal scale considered here. Nevertheless, depending on the strength of selection, even extremely limited amounts of gene flow can homogenize populations (e.g., Hartl and Clark 1989). Transplanting brook charr from unharvested to harvested lakes could therefore be a viable option for mitigating harvestinduced evolution.
Some other simplifying assumptions could also affect our estimation of the short-term effects of harvesting regimes. For instance, if there were strong interactions between the evolutionary effects of harvesting and residual evolutionary transients, our comparative approach could misestimate the former. For example, if strong stabilizing selection operates on the PMRN traits independently of harvesting, and if this stabilizing selection has not yet run its course when harvesting commences, the effects of even mild changes in the harvesting regime on the PMRN traits after 50 years could be amplified by the residual stabilizing selection. By initializing the model with trait values observed in the field and by commencing harvesting after 100 years of natural selection, we might therefore be underestimating the predicted evolutionary effects of milder harvesting regimes. Furthermore, the relatively high genetic coefficient of variation used to initialize the model elevates the propensity of the PMRN traits to respond to harvest-induced selection. If that coefficient were lower, the evolutionary response to harvesting over the time frame considered here would not be as pronounced (e.g., Dunlop et al. 2007).
In spite of some simplifying assumptions adopted in this study, we expect that the framework developed here can help guide the design of harvest regimes in brook charr. The comparative evaluation of different harvest regimes is mandated by a precautionary approach to fisheries management (e.g., Fenichel et al. 2008). The defining feature of the precautionary approach is the development of risk-averse objectives (Food and Agricultural Organization of the United Nations 1995; Peterman 2004). Moreover, harvest regimes may need to be designed to comply with legal mandates or ethical imperatives to conserve standing genetic variation (e.g., Humphries et al. 1995 andBalmford et al. 2005). Additionally, fisheries-induced evolution could impede recovery efforts for a fishery or have cascading ecosystem effects (e.g., Enberg et al. 2009). Hence, when harvestinduced evolution poses unacceptable risks to yields, fish stocks, or ecosystems, or when there are pressing legal or ethical guidelines to minimize fisheries-induced evolution, stakeholders should carefully assess 'worst-case' scenarios that are particularly likely to induce evolutionary change. Indeed, a sound, transparent, and well-communicated understanding of such scenarios can contribute to avoiding them.
We expect the results of our modeling work to be generalizable to formulating management strategies for other harvested fish species with similar life histories. By integrating size-specific life-history processes with elements of harvest-induced life-history evolution and management strategy evaluation, we expect our approach and results to shed new light on the causes and consequences of harvest-induced evolution in fish and thereby aid the development of sustainable harvest regimes. on the evolution of alternative life-history tactics in brook charr. Evolutionary Applications 1:409-423. Tremblay, S., and P. Magnan. 1991 Brook charr were captured with experimental multifilament gillnets set randomly perpendicular to the shore. A minimum of 100 brook charr were collected over a span of 1-5 days. Data were available from two field seasons, and were pooled for the purposes of the present analysis. For a detailed description of the sampling methodology and the data collection methods employed, see Magnan et al. (2005). Specifically, the following parameters are based on this data set: (i) the coefficient p and exponent u of the length-mass relationship, (ii) the average mass W 0 of an egg, (iii) the strength b 1 of size-unspecific effects on mortality, the strength b 2 of size-specific effects on mortality, (iv) the metabolic rate d, the size-specific consumption coefficient a, the size-specific consumption exponent c, the standard deviation r f of the stochastic component of the growth equation, and (v) the fraction of energy intake allocated to reversible mass that is used for reproduction. The probabilistic maturation reaction norm (PMRN) parameters d,m, and r are also estimated from this data set, as described in Appendix B.
(i) Coefficient p and exponent u of the length-mass relationship l(X) ¼ pX u . The parameters p and u were obtained as coefficients of a least-squares linear regression of log-transformed irreversible mass on log-transformed lengths. Irreversible mass were estimated as follows. Because sampling was conducted at the end of the growth season (mid-September through early October), we assumed that the ratio between reversible mass and irreversible mass was at or near its maximum, where q is the maximum ratio. For mature females, mass before and after gonad removal were available. We assumed that once gonadal mass is subtracted from the total mass of mature females, their maximum ratio between reversible and irreversible mass is comparable to the maximum ratio for juveniles. Thus, the mass without gonads for mature females, as well as the total mass for immature individuals, were divided by 1 + q J . For mature males, we calculated their irreversible mass as W i / (1 + q A ), where q A is defined in Appendix C. For the regression analysis, the irreversible masses of immature and mature individuals of both sexes were pooled, yielding a total sample of 793 individuals. Our demographic model also uses the length-mass relationship to calculate the survival of eggs and newborn individuals. Because our data set did not contain individuals from this group, we included published information on the length-mass distribution of recently hatched alevins (described in Table 1 in Curry et al. 1993, p.133) in the regression analysis.
(ii) Average mass W 0 of an egg. Eggs were sampled from the gonads of mature females and were weighed in the laboratory.
(iii) Strength b 1 of size-unspecific effects and strength b 2 of size-specific effects on mortality. Our demographic model was initialized with values of b 1 and b 2 resulting in a logistic increase of survival with body mass. The model was then run for 100 years without genetic variability (l ¼ 0) and without harvesting. For each model run, the model's predicted mass distribution at the end of the run was compared to the mass distribution observed in the data. The parameters b 1 and b 2 were then simultaneously adjusted through trial and error until the average two-sampled Kolmogorov-Smirnov test statistic (e.g., Press et al. 2007) between the predicted mass distributions of a run (based on a uniform sampling of the simulated brook charr from the size range caught in the nets) and the observed mass distribution pooled with backcalculated mass was minimized. Figure A1 shows the resultant natural survival functions.
(iv) Metabolic rate d, size-specific consumption coefficient a, size-specific consumption exponent c, and standard deviation r f of the stochastic component of the growth equation. The parameters d and c were estimated following the method described in Hiyama and Kitahara (1993). Briefly, this method required fitting a von Bertalanffy growth equation to the length-at-age data using nonlinear least-squares regression. The predicted lengthat-age data were then transformed back to mass-at-age data. For females, the log-transformed gonadal mass was regressed on log-transformed length. The surplus energy (measured in terms of mass gain), DW, accumulated between years t and t + 1 was calculated according to Hiyama and Kitahara (1993) Because the white sucker Catostomus commersonii was absent in all but two of the lakes used for this analysis, we assumed that brook charr generally focused on zoobenthos as their primary prey (e.g., Magnan 1988 andTremblay andMagnan 1991). Consequently, predation pressure imposed on zooplankton by brook charr was treated as being generally light, so that h(R t ) % 1. The parameters n,c, and g could then be estimated by applying nonlinear least-squares regression to the equation, and solving for a and d in Equations (15) and (16). r f was inferred from the residual SE of this regression. The nonlinearity inherent in Equation (A1) ensures that as individuals reach their maximum observed size, DW approaches 0.
(v) Fraction of energy allocated to reversible mass used for reproduction. The estimation procedure for was similar to that used to estimate the mortality parameters b 1 and b 2 . The demographic model was initialized with a value of between 0 and 1. The model predicted the ratios between female gonadal mass and total mass. These ratios were compared to the ratios between female gonadal mass and total mass observed in the data. The value of was adjusted through trial and error until it appeared that the value that minimized the mean squared difference between these ratios was obtained. This value was then selected for use in the model. Figure A2 shows that our demographic model of brook charr population dynamics in the absence of fishing successfully recovered the size distribution observed in unharvested lakes in the field (Magnan et al. 2005). We confirmed this through two different tests. First, we used a Mann-Whitney test to determine how well the predicted size distribution (based on a uniform sampling from the size range caught in the nets) matched the empirically observed size distribution pooled with backcalculated mass, finding a mean P-value of P ¼ 0.20 and a median P-value of P ¼ 0.08 over 50 replicate model runs. Second, we also split the mass distribution into 20 size classes and used linear regression

Appendix A.2. Model validation
to calculate how well the model-predicted counts x i for size class i matched the empirically observed counts y i , finding an average R 2 of 0.96 over 50 replicate model runs. We could also confirm that a never significantly differed from 0 and that the average value of b(1.05) was within two standard deviations of one. This good fit seems largely driven by the fact that our demographic model performs well for small individuals (£ 50 g) and large individuals ( ‡ 1000 g), even though it tends to underestimate counts in the range between 50 g and 150 g, and slightly overestimates counts in the range between 250 g and 950 g (Fig. A2).

Appendix B. Estimating the probabilistic maturation reaction norm
We used the bootstrapping method described by Barot et al. (2004) to estimate the PMRN from field data, Figure A1 (A) Annual natural survival as a result of processes other than starvation risk as a function of irreversible mass and; (B) annual natural survival as a function of reversible and irreversible mass. As the ratio between reversible and irreversible mass declines, survival declines. In contrast, for a given ratio between reversible and irreversible mass, survival rapidly improves as individuals grow from avelins (with an irreversible mass of %0.004 g) to 100 g.

Figure A2
Comparison of the model-predicted size structure (continuous curve) with the empirical size structure observed in the field (open circles).
taking the following five steps for each bootstrap replicate. First, we estimated maturity ogives (defined by the probability o(A,l) of a charr's being mature given its age A and length l) and annual growth increments Dl(A) based on back-calculated lengths for brook charr from unharvested lakes (yielding a total sample of 711 individuals). Second, we estimated the probability of being mature at age A and length l as ( Barot et al. 2004) mðA; lÞ ¼ oðA; lÞ À oðA À 1; l À DlðAÞÞ 1 À oðA À 1; l À DlðAÞÞ : Third, we used linear regression to describe the obtained m(A,l) by a logistic length dependence at each age A, lnð mðA; lÞ 1 À mðA; lÞ Fourth, denoting the length at which m(A,l) ¼ 0.5 by l p50,A ¼ )d 0,A /d 1,A , we estimated the PMRN intercept m and the PMRN slope r by fitting the linear model Fifth, we estimated the PMRN width d, which measures the length difference between m(A,l) ¼ 0.25 and m(A,l) ¼ 0.75 for all ages A, as the mean of l p75;A À l p25;A lnð0:75=0:25Þ À lnð0:25=0:75Þ (Dunlop et al. 2009) for A ¼ 2,3, and 4, where l p25,A and l p75,A denote the lengths at which m(A,l) ¼ 0.25 and m(A,l) ¼ 0.75. We justify using a PMRN width that is independent of age A by noting that the mean bootstrapped values for l p25,A and l p75,A for each age were within one standard deviation of the mean bootstrapped values l p25,A and l p75,A , respectively, for the other two ages. Finally, we averaged the values of m 0 , r 0 , and d over 1000 bootstrap replicates of this procedure to obtain the values reported in Table 1.
To increase sample sizes, we pooled males and females. Although maturation patterns in fish populations often differ between males and females, it is unknown whether some or all loci coding for the PMRN are indeed sexlinked. Sexual dimorphism in maturation patterns is consistent with a PMRN shared by the sexes when combined with the common observation that males and females grow at different rates. Because currently almost nothing is known about the loci controlling the PMRN, we feel that the case of the loci underlying the PMRN being autosomal, resulting in a common PMRN for both sexes, will provide a basis on which to develop further work relaxing this simplifying assumption once further empirical data can clarify this question.
The assumption that PMRN loci are autosomal is supported by an earlier study quantifying the PMRN of another salmonid, Oncorhynchus keta, which found no sex-specific differences in the estimated PMRN (Morita and Fukuwaka 2007). O. keta is anadromous and spends much of its life at sea, whilst our populations of Salvelinus fontinalis are inland populations that remain resident around the year. Although our sample size, when split across the sexes, does not provide sufficient statistical power to draw clear conclusions about whether PMRN traits differ between males and females, to the extent that PMRN loci are conserved across salmonids, the results of Morita and Fukuwaka (2007) suggest that such loci are indeed autosomal.
Appendix C. Estimating the allocation of energy intake to irreversible and reversible mass The function F A (X t ,Y t ) describes the fraction of energy intake that is allocated to irreversible mass throughout year t, in dependence on an individual's irreversible mass X and reversible mass Y at the start of year t. Here we estimate this function as a quadratic polynomial in three steps described below. Notice that individuals are born with a maximum ratio between reversible mass and irreversible mass; as long as individuals do not starve, juveniles will retain this maximum ratio.
Since F A (X t ,Y t ) is based on a continuous-time process governing growth in reversible and irreversible mass, we begin by numerically integrating the following set of differential equations over a specified range of sampled initial values and input parameters (Table C2), where s denotes time during the year and ranges from 0 to T, C(X,Y) describes the instantaneous allocation of incoming energy to irreversible mass, E A is a constant describing the food consumed during the year less the metabolic costs (i.e., the surplus energy sensu Hiyama and Kitahara 1993 -for the functional form of C(X,Y), see, e.g., Broekhuizen et al. 1994 andPersson et al. 1998), and 0 £ m s (M) £ 1 is a function with a parameter M describing how evenly the incoming energy is temporally distributed throughout the year. The procedure we used to obtain samples of m s (M) is as follows: 1 We define 2 We draw 74 random deviates, each representing an interval of approximately 5 days out of the year, from a uniform distribution with a minimum of 0 and a maximum of 20.  (Table C1). Equation (C1) are integrated separately over each of the 74 intervals. Repeating the sampling of points from m(M) did not affect subsequent analyses.
In a second step, we assume that the instantaneous allocation of energy intake is described by where q is the maximum feasible ratio between reversible and irreversible mass. This functional form ensures that when Ys Xs ¼ q the proportion of incoming energy allocated to reversible and irreversible mass maintains their ratio, so that dYs=Xs ds ¼ 0. Considerable empirical work has resulted in general agreement that q for non-reproducing fish (q J ) equals approximately 1.6 across a range of fish taxa (Broekhuizen et al. (1994) and references therein). For reproducing fish, we describe q(q A ) as the sum of the maximum gonadosomatic index (%0.75) and q J . Thus, upon maturation, the maximum ratio between reversible and irreversible mass increases, because individuals, in addition to fat and lipid reserves, now allocate energy intake to gonads.
Following the integration of Equation (C1) from s ¼ 0 to s ¼ T, we fit to the results a quadratic statistical model of the form where Q i denote quadratic terms of the five predictors X(0),Y(0)/X(0),E A ,M and q,I j denote all possible twoway interaction terms between the predictors, and e is a normally distributed stochastic error term.  Figure A4 (A) Annual energy intake and; (B) fraction F A of this intake allocated to irreversible mass, when the resource is at its carrying capacity (K ¼ 1.8 · 10 9 zooplankton individuals). As the ratio between reversible mass and irreversible mass increases, individuals can allocate a higher proportion of their energy intake to the growth of irreversible mass, while individuals with a low ratio must allocate a lower fraction of their energy intake to the growth of irreversible mass to maintain the maximum ratio between reversible and irreversible mass. The boundaries in (B) arise from two distinct biological constraints. The upward sloping boundary on the left side of the panel demarcates biologically feasible ratios between reversible and irreversible mass. For example, the model does not allow an individual with an irreversible mass of 10 g to have a reversible mass of 3000 g. The downward sloping boundary on the right side of the panel results from the model assumption that for fish with a combined irreversible and reversible mass above a certain threshold (%4000 g), metabolic costs become so high that there is no surplus energy intake at the end of the year. Panel (B) applies to an individual that matures at a length of 215 mm (corresponding to a mass of %30 g), with q J ¼ 1.61 and q A ¼ 2.35. X(0),E A ,M and q. For each predictor, 1028 equidistant points were sampled over the ranges given in Table C2. The points were sampled in such a way that the maximum correlation between any two predictors was 0.0004, while also ensuring more than sufficient data points were available to accurately estimate the coefficients in Equation (C4) (e.g., Cioppa and Lucas 2007). As by definition F A ðX; YÞ ¼ XðTÞ À Xð0Þ E A , the most parsimonious model was inferred by the stepwise elimination of terms in Equation (C4), and the final linear regression model was given by Xð0Þ where e a is normally distributed random deviate with mean 0 and standard deviation r e . The estimated values of B 1 to B 12 and r e are reported in Table C3. This linear regression model had an adjusted R 2 of 0.99, and indicated that temporal heterogeneity in resource intake has little effect on the final pattern of resource allocation. Figure A4 depicts the estimated function F A (X,Y) assuming no resource scarcity. The heavier an individual, the more resources must be allocated towards maintaining reversible mass, at the expense of growth in irreversible mass.  Parameter Description Sampled range Unit Comments on range X(0) Initial structural mass 4.5-5000 g Chosen to be well inclusive of the range of X(0) described in Broekhuizen et al. (1994)

Y(0)
Initial reversible mass 0.019-9990 g Chosen to be sufficiently wider than the range of X(0) E A Expected resource intake 19.5-200 000 g Chosen to cover the range of intra-annual surplus energy across fish taxa (e.g., Hiyama and Kitahara 1993  Effect of initial structural mass )6.02 · 10 )5 g )1 B 2 Effect of initial ratio between reversible and structural mass 0.0714 -B 3 Effect of energy intake 1.70 · 10 )5 g )1 B 4 Effect of maximum ratio between reversible and structural mass 1.00 -B 5 Quadratic effect of initial structural mass 1.18 · 10 )9 g )2 B 6 Quadratic effect of initial ratio between reversible and structural mass )7.81 · 10 )3 -B 7 Quadratic effect of maximum ratio between reversible and structural mass )0.587 -B 8 Effect of initial reversible mass 3.48 · 10 )5 g )1 B 9 Interaction effect of initial structural mass and total energy intake 3.89 · 10 )5 g )2 B 10 Interaction effect of energy intake and initial ratio between reversible and structural mass 1.07 · 10 )5 g )1 B 11 Interaction effect of initial and maximum ratio between reversible and structural mass 0.178 -B 12 Interaction effect of energy intake and maximum ratio between reversible and structural mass 1.27 · 10 )5 g )1 r e Standard deviation of error term in Equation (15) 0.0300 g See Appendix C for details.