Fisheries Research 109 (2011) 131–139 Contents lists available at ScienceDirect Fisheries Research journal homepage: www.elsevier.com/locate/fishres Six decades of pike and perch population dynamics in Windermere Øystein Langangen a,∗ , Eric Edeline b , Jan Ohlberger a , Ian J. Winfield c , Janice M. Fletcher c , J. Ben James c , Nils Chr. Stenseth a , L. Asbjørn Vøllestad a a Centre for Ecological and Evolutionary Synthesis (CEES), University of Oslo, Department of Biology, P.O. Box 1066 Blindern, N-0316 Oslo, Norway UPMC-Paris 6, UMR (CNRS) 7618, Laboratoire Biogéochimie et Ecologie des Milieux Continentaux, 46 rue d’Ulm, F-75230 Paris Cedex 05, France c Centre for Ecology & Hydrology, Lancaster Environment Centre, Library Avenue, Bailrigg, Lancaster, Lancashire LA1 4AP, UK b a r t i c l e i n f o Article history: Received 14 September 2010 Received in revised form 26 January 2011 Accepted 26 January 2011 Keywords: Age-structured population models Fish stock assessment Hyperdepletion Intraguild predation Population dynamics a b s t r a c t We use six decades of catch-at-age data for perch (Perca fluviatilis) and pike (Esox lucius) in Windermere (UK) to estimate age- and sex-specific population sizes, natural mortalities and catchabilities in both species. Population sizes are estimated by fitting age structured population models to the catch-at-age data using standard maximum likelihood methods. We validate our methods using data simulations, and use our estimates of vital rates (natural mortality, recruitment and catahability) to address important aspects of fisheries biology. Our model indicates that strong fishery selection against male perch apparently triggered a population collapse, highlighting that sex-selective fisheries can be harmful even at a reasonable exploitation rate (here ≤30%). Recruitment (R) increased with the abundance of spawners (S) in both species, but it also responded to both abiotic and other biotic factors. In particular, increased predator (pike) abundance induced a change from compensation to depensation in the prey (perch) SR relationship, thus favouring the occurrence of an Allee effect. Our study provides reference points for the effective exploitation of pike and perch populations, and underscores the need for ecosystem-based harvesting management. © 2011 Elsevier B.V. All rights reserved. 1. Introduction Estimation of population demographic parameters is at the very core of ecological research and management (Walters and Martell, 2004). Fisheries based data and fish stock assessment models provide unique opportunities to estimate demographic parameters that often remain unknown in unexploited species. Here, we use fisheries data and fish stock assessment models to estimate ageand sex-structured population size, natural and fishery mortality, catchability, and relationships between spawner abundance (S) and subsequent recruitment (R) in pike (Esox lucius) and perch (Perca fluviatilis), two species widely distributed around the northern hemisphere (Nelson, 2006). Perch and pike are often keystone species in ecosystems, and both species support widespread commercial and recreational fisheries. Hence, studies of perch and pike population dynamics are of great value for both researchers and managers. We base our approach on fitting age-structured population models to 60-year long catch-at-age data series from a scientific sampling on pike (the top-predator) and perch (its main prey) in Windermere, a glacial valley lake comprised of two (north and ∗ Corresponding author. Tel.: +47 22 85 44 00; fax: +47 22 85 40 01. E-mail address: oysteol@bio.uio.no (Ø. Langangen). 0165-7836/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.fishres.2011.01.029 south) basins in the English Lake District (UK). Detailed description of the lake, its fish and its fisheries can be found in Le Cren (1987), Le Cren (2001) and Winfield et al. (2008). Paxton et al. (2004, 2009) give full sampling details for the pike and perch, respectively. Both pike (e.g. Kipling and Frost, 1970; Kipling, 1984) and perch (e.g. Le Cren et al., 1977; Kipling, 1984; Le Cren, 1987) have been subjected to a large number of studies over time, some of which have dealt with estimating their population numbers. The models used in this paper are an improved (different models) and extended (more data) version of the unpublished work of des Clers et al. (1994). Moreover, des Clers et al. (1994) did not validate their model suggesting that it may be overparameterized; thus, we validated our model using simulations. One of our goals is to make the model parsimonious with respect to the data, by (i) adding 10 years of data to the analysis and (ii) testing the performance of different fisheries-based models using synthetic data. Our results are generally consistent with previously reported population trends, but often return very different estimates. In order to robustly support our results we extensively evaluate our models with simulated data, and technical aspects thus represent a substantial part of this paper. We further our analyses by exploring perch and pike stock–recruitment (SR) relationships using newly estimated population numbers. Previous SR analyses either used crude catch per unit of effort (CPUE, sum of fish at age/number of traps used 6 weeks on the spawning ground) (e.g. 2.1. Study site Windermere is situated (54◦ 22′ N, 2◦ 56′ W; altitude 39 m) in the English Lake District, U.K. It comprises a mesotrophic north basin (surface area 8.1 km2 , maximum depth 64 m) and a eutrophic south basin (surface area 6.7 km2 , maximum depth 44 m), although cultural eutrophication has impacted both parts of the lake (Winfield et al., 2008). The fish communities of the two basins are qualitatively similar and contain only seven species of current or former numerical importance, i.e. Arctic charr (Salvelinus alpinus), Atlantic salmon (Salmo salar), brown trout (Salmo trutta), European eel (Anguilla anguilla), perch, pike and, in recent years, roach (Rutilus rutilus), although at least a further nine species of minor importance are also present (Pickering, 2001). 1000 6000 5000 100 4000 3000 10 2000 1000 0 1 7000 1000 6000 5000 100 4000 3000 10 2000 1000 2.2. Fish sampling 0 1960 1970 1980 1990 1 2000 Year Fig. 1. Total perch catch (dashed lines) in the north (upper panel) and the south basin (lower panel) of Windermere are shown. The total (commercial and scientific) effort (dotted line with diamonds) and scientific effort (solid line with diamonds) are also shown. Note that the scientific monitoring effort accounts for the whole effort after 1948 in the north basin and after 1964 in the south basin. Catch [total number of fish] Scientific monitoring for perch and pike in Windermere started in 1943 for perch and in 1944 for pike and was initiated with a commercial-like fishery for perch (hereafter commercial fishery) and a pike removal program to enhance perch numbers. The scientific sampling continues today with effectively no change in gear type and fishing methods. Perch are sampled on some of their spawning grounds usually from mid April to early June, using the same traps as employed in the commercial fishery. Perch traps are unselective for individuals of size 90–300 mm (Le Cren, 1987, meaning the traps retain equally well individuals in this size range) and capture mainly mature males because such fish remain on the spawning grounds for long periods while females visit only to spawn (Le Cren et al., 1977; Craig, 1987). In the scientific monitoring, a percentage (sometimes 100%) of the catch was measured for length and weight and the left opercular bone is removed to determine the age of the fish in years (Le Cren, 1947). Sex is determined by internal examination. Catch from the commercial fisheries was reported as total weight and effort. A summary of perch effort and catch is shown in Fig. 1. Pike are sampled in winter (usually from mid October to January) in shallow areas of the lake using 64 mm gill nets. The gill nets are size selective for pike larger than 550 mm (pike smaller than 550 mm are rarely caught), with a slight reduction in catchability above 650 mm (Frost and Kipling, 1967; Carlson et al., 2007). The selectivity of the gillnet means that we mainly catch pike of age 2 and older. The gear has remained the same over the entire study period. Before 1990 nets were set at a location and kept there until yield decreased (usually in one or two weeks), but since 1990 the netting time has been kept fixed at 240 net-days per year within the period from mid October to late December and site depletion no longer drives site-specific timings. As for the perch, pike are measured for length and weight, and the left opercular bone is used to determine the age of the fish in years (Frost and Kipling, 1959). Sex was again determined from internal examination. A summary of pike catch and effort is shown in Fig. 2. 1950 1000 10000 800 8000 600 6000 400 4000 200 2000 0 Effort [net days] 2. Materials and methods 7000 Catch [Total number of fish] Mills and Hurley, 1990; Paxton et al., 2004, 2009) or the previous estimates from des Clers et al. (1994) (Edeline et al., 2008). Compared to this previous work, our SR analysis confirms that temperature is a key parameter for the success of both pike and perch recruitment, and brings new insights about predator–prey interactions. We finally examine the causes for perch dynamics under commercial exploitation, and propose sex-selective fishing as a possible cause for a transient perch collapse in the late 1940s. Effort [#traps 6 weeks on spawning ground] Effort [#traps 6 weeks on spawning ground] Ø. Langangen et al. / Fisheries Research 109 (2011) 131–139 Catch [Total number of fish] 132 0 1950 1960 1970 1980 1990 2000 Year Fig. 2. Pike gill net effort (dotted line) and total catch (solid line) for the combined north and south basins of Windermere. 2.3. Single species models of perch and pike Single species models are widely used in fish stock assessments (e.g. Myers et al., 1997; Millar and Meyer, 2000; Aanes et al., 2007). The main reason is simplicity (e.g. independence of cohorts, fewer parameters, etc.) but also it has been questioned if multi-species models are an improvement over single-species Ø. Langangen et al. / Fisheries Research 109 (2011) 131–139 models in fisheries management (Hollowed et al., 2000). As an exploratory analysis we tested a simple lag recruitment survival and growth model (LRSG, see Hilborn and Mangel, 1997) with a Ricker stock–recruitment model for pike. This test showed that such models probably underestimate the population size (all fish in the model are caught in some years) and miss the yearly fluctuations in abundance (predicts only smooth trends). Another drawback of such models is that the form of the stock–recruitment is a priori defined, possibly resulting in errors in recruitment that can propagate through the model and give unacceptable errors in the numbers at age. We therefore use a more flexible single species age-structured model including fishing mortality as well as natural mortality (see e.g. Hilborn and Walters, 1992, chapters 10–11). The Windermere perch and pike data are particularly well suited for fitting age-structured population models because they do not have the limitations commonly found in fisheries datasets (Quinn and Deriso, 1999): (i) data are scientifically collected, (ii) they cover an exceptionally long time period with almost no change in fishing gear and methods, and (iii) the fishing effort has varied greatly and is accurately documented over the study period. 2.3.1. The pike model Pike disperse freely between the two basins of Windermere and we therefore modeled them as one population (Haugen et al., 2006, 2007). We allow for seven year-classes (age 3–9) of pike to be included in the analyses, and we assume that all pike older than age 9 die in our model. Our pike model is a standard age-structured model 133 2.3.2. The perch model Perch populations are considered to be independent between the north and the south basins of Windermere according to capturemark-recapture and genetic data (Kipling and Le Cren, 1984; Bodaly et al., 1989). We therefore fit two models both having the same core structure as our pike model (Eqs. (1) and (2)). Only 5 year classes (age 2–6) appeared frequently enough to be included in the analysis. Hence, we assume that in all perch die after age 6 in the model. We have assumed equal catchability in the commercial fishery and scientific sampling due to the fact that equal traps and similar fishing grounds (spawning grounds) were used both in the commercial fishery and in the scientific sampling. In the model, we allowed for (i) sex-specific, maturation-dependent catchability of 2 vs. 3+ perch (traps are set on spawning grounds only), (ii) sexspecific, stepwise constant natural mortality before, during (1976) and after a severe disease (Bucke et al., 1979), (iii) equal recruitment for both males and females, (iv) age-dependent fishing mortality due to strong maturity selection effects in catchability, (v) ageindependence of natural mortality. A model with yearly temporal variability in natural mortality was less robust, possibly due to over parameterization. Each perch model estimated 10 catchability and natural mortality parameters, plus recruitment for each year. 2.4. Parameter estimation We assume lognormal sampling errors in the catch data. The objective function then becomes LL∼  2 obs sim (ln(Ct,a,sex + 1) − ln(Ct,a,sex + 1)) , (3) t,a,sex Nt+1,a+1 = (Nt,a − Ct,a ) · e−Mt , (1) where Ct,a is the catch-at-age a in time t, and we have assumed no age variation in the natural mortality (Mt ). In this model the recruitment, fishery removals and natural mortality act sequentially on the population vulnerable to the fisheries. The population first receives the recruits (Rt = Nt,a=min where min is the age at recruitment to the fisheries, age 3 for pike), then the fishery removal is applied before the natural mortality is applied to the population. The sequential model is acceptable for the system we study in this paper because the fishing seasons are rather short (approximately 8 weeks for pike and approximately 6 weeks for perch). The fishery catch equation is ˇ Ct,a = Nt,a (1 − e−qa Et ), (2) where qa is the catchability and Et is the effort (where fishing mortality F ≡ qE). ˇ is a parameter describing the non-linear relationship between the population and the catch (see e.g. Hilborn and Walters, 1992; Wilberg et al., 2010), which was not present in the model of des Clers et al. (1994). Indeed, pike were sampled more intensively in the zones of highest CPUE in a commercial fishery-like manner (see Le Cren, 2001). Furthermore, the pike has a non-homogeneous density distribution in Windermere (Paxton and Winfield, 2000). Both effects are likely to contribute to a nonlinear relationship between catch and number of pike in gill net fisheries (see Hilborn and Walters, 1992, chapter 5). We allowed for (i) sex- and age-specific (age 3 vs. age 4+, testing models with more age classes did not perform well under model validation, see Section 2.5) catchability because both factors strongly influence the probability to be recruited to the fishery, and (ii) sex-specific and constant natural mortality. The pike model contained 6 catchability and natural mortality parameters to be estimated, ˇ, plus recruitment for each year (equal for males and females at age 3). where the (+1) is added to avoid infinities in the occasional events of zero catch of an age group. The catch used in Eq. (3) is the observed and simulated scientific catch, and the subscripts t, a and sex represent time, age and sex, respectively. In the objective function (Eq. (3)) we use the modeled scientific catch to compare with the observed scientific catch. In the model we obtain the scientific catch by multiplying the catch with the fraction of scientific effort to total effort C s = C tot Es Es tot · N(1 − e−qE ) = s tot E + Ec E (4) where the superscripts s, c and tot stand for scientific, commercial and total, respectively. Model parameters (catchability qa , natural mortality Mt and recruitment Rt ) were estimated by numerical optimization of Eq. (3) using a genetic algorithm (PIKAIA, Charbonneau, 1995). The genetic algorithm is constructed to find the global extrema by mimicking the process of natural selection. The method is inspired by concepts such as heredity, mutations and fitness. In short, the algorithm uses the objective function (Eq. (3)) as a measure of the fitness of a specific set of parameter values. For each iteration (or generations), a set of parameters (typically 100 different parameter values) is constructed based on the best-fit parameters from the last generation (typically 10% of the best fit parameters are allowed to contribute to the next generation). Note that the first generation of parameters is randomly selected from the search interval. The new set of parameters is given random fluctuations by allowing crossover (random crossing of the digits in two surviving parameter values) and random “mutations” (in the digits of a parameter). The crossover and mutation give the algorithm the necessary flexibility to move through the fitness landscape and eventually find (with a certain probability) the global extrema. The algorithm is used to minimize the objective function for a given model to obtain point estimates for the parameters. The catchability and the natural mortality are estimated simultaneously by minimizing the total (over 134 Ø. Langangen et al. / Fisheries Research 109 (2011) 131–139 all cohorts, as in Eq. (3)) objective function. For each proposed value for the natural mortality and the catchability we calculate a cohort specific objective function value (similar to Eq. (3), but only summing over the relevant cohort). The objective function over each cohort is minimized (Brent’s algorithm for 1 dimensional minimization, see e.g. Press et al., 1992) to fix the yearly recruitment. To summarize, the method makes a “guess” on the values of the natural mortality and catchability, which are then used to minimize the objective function over each cohort with respect to the recruitment. The sum of all cohort objective functions is used to evaluate the “guess” values of the natural mortality and catchability. Our models could not simultaneously estimate natural mortality, catchability and ˇ due to confounding effects, and we therefore manually tested model fit with different fixed ˇ values (from 0.8 to 2.0 by steps of 0.1). The result favoured a ˇ larger than one in the pike model (i.e. hyperdepletion), with 1.6 giving the best fit, similar to previously reported values for pike by Pierce and Tomcko (2003). In the two perch models we found a value ˇ close to 1, hence we use ˇ = 1 in the analysis of the perch data (see Section 2.5 for more detail on the testing of ˇ). We tested the convergence properties of each model over 10 runs, each run with different random initial starting values for the model parameters. The 10 runs converged to the same parameter values (with an error in the third digit due to the randomness of the genetic algorithm) and we use the best fit result as the point estimate. The correlation between the observed and the modeled catches is shown in supplementary material (Figs. Siii and Siv). 2.5. Model validation We extensively tested our models with simulated data (following the methods in Chapter 11, Hilborn and Mangel, 1997) of the same size as the real perch data, but generated from known, widely varied parameter values (see supplementary material, Table S1). The estimation method was applied to the simulated data to test whether our models could successfully estimate the known input parameters and were robust to parameter variation. These tests have shown that our models performed very well (a few percent deviation in parameter values). Furthermore, we tested the robustness to observational error by adding random noise to the data set (Monte Carlo tests redone 100 times). The returned parameter estimates formed bell-shaped curves centred at (or very close to) the input values of the parameters. Hence, we conclude that our method is very robust to observational noise. We also specifically tested the ˇ estimation method. A synthetic data set (similar to the pike data) was constructed from random recruitment and fixed natural mortality and catchability (close to, but not similar to the point estimates from the model), while ˇ was set to 1.5. Attempts to estimate ˇ as a free parameter resulted in underestimation (and overestimation of the population) due to strong confounding effects with catachability. Furthermore, we tested the method of fixing ˇ, with favourable results. However, this method does not allow us to include the error in ˇ in the final population estimates. We tested the error in the population estimates by using different ˇ values on the real data set (ˇ values of 1.5, 1.6 and 1.7, close to the ones used in the pike model). An error in the ˇ mainly introduces a total bias (over the full time series) of about 40% for ˇ = 1.5 (over estimation) and 10% for ˇ = 1.7 (under estimation). The time series of catch-at-age in our model and observed catches correlates well (correlation coefficient above 0.9, see supplementary materials (Figs. Siii and Siv)), but due to lack of independent measure of catch error it is hard to give an objective measure of the model fit. We have also compared the reported commercial catch with the predictions of the model. In general the model performs well, but with some discrepancies (see Appendix A for more details). 2.6. Parameters confidence intervals Parameter confidence intervals were estimated using parametric bootstrapping (Efron and Tibshirani, 1998; Press et al., 1992), assuming only sampling errors and independent and identical distribution (IID) of the observations. We construct a simulated best-fit dataset from the point estimates, and new simulated data sets are drawn from this data set by adding a random observational error to each of the observed catch-at-age data. The observational standard deviation is obtained by fitting a Gauss curve to the log-transformed catch residuals from the point estimates (standard deviations of 0.61, 0.62, and 0.59 for pike, north basin perch and the south basin perch, respectively). We test how many samples are necessary to obtain reasonable ranges for the parameter confidence intervals by running 1000 samples for the north basin perch model. From this test we conclude that about 200 samples give reasonable confidence limits (less than 5% deviation from the full sample), and we use this number in the south basin as well as for the pike model. Removing the 5 lowest and 5 highest estimates from the sample, we obtain the 95% confidence limits. The residuals show some marginally significant autocorrelations (in about half of the time series, correlations typically around 0.3, see Fig. Sii). However, we conclude that the IID assumption is reasonable, based on simulated data using observed autocorrelation and showing only a slight increase in parameter confidence bands (around 20%). We have tested the effect of including process error on the estimated confidence limits. We performed the test by allowing the natural mortality and catchability to temporally fluctuate around a mean value in a simulated dataset (random Gaussian distribution with a variance of 20% of the point value, see Section 2.5 for more details on synthetic data). Furthermore, we have allowed for a factor 2 increase in the catchability after the disease in the simulated data. The method estimated values close to the mean values and the bias is very small (a few percent), hence indicating that the method was robust to these process errors. From the above tests we conclude that the method performs satisfactorily on the type of data presented in this paper. 2.7. Stock–recruitment analysis SR relationships describe the form of density-dependence in populations, and are thus critical elements of fisheries models and management plans (Myers, 2001). From our population estimates we explored the relationship between S (age 3–9 pike, age 2–6 perch) in year t (from 1944 in pike and from 1943 in perch) and R (number of age 3 pike, age 2 perch) in year t + 3 (pike) and year t + 2 (perch). We modeled perch and pike SR using Gamma (i.e., Gamma distributed R) with a log-link (i.e., link function between the mean of the distribution function and the linear predictor) GAMs in the mgcv library of R (Wood, 2006a,b). These models are equivalent to Gamma (with a log-link) GLMs, except that nonparametric smoothing functions (piecewise polynomials) can be placed into the linear predictor (here on the S effect). Hence, in contrast to classical approaches (e.g. Ricker or Beverton and Holt SR relationships; Myers, 2001), we did not assume any a priori functional form of the SR relationship. In our SR GAMs we tested the additive and nonadditive effects of a set of environmental factors including the whole combinations of available climatic (yearly mean water temperature, sum of degree days above 14 ◦ C, winter and summer North Atlantic Oscillation index – NAO index) and biotic variables (recruit and adult numbers of perch or pike, and mean body size of the stock) during first year of life, second year of life, third year of life (pike Ø. Langangen et al. / Fisheries Research 109 (2011) 131–139 135 North basin . Total number of perch age 2-6 Total number of Pike age 3-9 2.0 10 4 . 1.5 10 4 . 1.0 10 4 . 5.0 10 3 1.5•106 1.0•106 5.0•105 0 1950 1960 1970 1980 1990 0 2000 1950 Year 1970 1980 1990 2000 1990 2000 Year Table 1 Parameters estimated in the pike model for the two basins combined with 2.5 and 97.5 confidence limits deduced from the 200 bootstrapping samples with ˇ = 1.6. Male M and Female M are the male and female natural mortality and Male qage and Female qage are the sex specific catchabilities, respectively. 2.5% Pointest. 97.5% 0.221 0.318 0.435 0.810 0.268 1.35 0.257 0.371 0.0640 1.05 0.368 1.75 0.283 0.396 0.955 1.57 0.550 2.51 South basin Total number of perch age 2-6 Fig. 3. Total estimated pike population (age 3–9 years) in the combined north and south basins of Windermere. Values are given as point estimates (solid line) with 95% confidence limits (grey area). Male M Female M Male q3 (×106 ) Male q4+ (×106 ) Female q3 (×106 ) Female q4+ (×106 ) 1960 1.5•106 1.0•106 5.0•105 0 1950 1960 1970 1980 Year only), or from birth to recruitment. All these variables are likely to influence recruitment success (Myers, 2001). We found that total phosphorus concentration was highly correlated with pike number (Pearson correlation = 0.82), precluding using both predictors simultaneously. We chose to discard total phosphorus concentration and to keep pike numbers in our perch SR analyses based on the assumption that perch was more likely top-down than bottomup controlled. This assumption is based on theoretical predictions from (Lotka–Volterra) predator–prey models in which increased primary productivity translates into increased predator biomass but not prey biomass (Case, 2000). We quantified the relative contributions of environmental variables to recruitment variance using Brian Ripley’s Tree function for the R software (Venables and Ripley, 2002, binary recursive partitioning). At this stage we discarded NAO index (both winter and summer), mean water temperature and mean body size of spawners because these variables were weakly informative. We then incorporated the remaining, most informative climatic and biotic predictors in GAMs, and selected the best model structure as the one minimizing model’s AIC (i.e., the AIC of the lme component of the model in a GAMM; Wood, 2006a). 3. Results and discussion 3.1. Pike demography The total estimated population with 95% confidence bands is shown in Fig. 3 and the estimated parameters are presented in Table 1. Note the rather sharp decline in the total pike population at the very end of the time-series in Fig. 3. As this decline appears only at the edge of the time-series it has an appreciable level of uncertainty due to truncated cohorts. The estimated natural mortality parameters are very close to previously reported values Fig. 4. Total estimated perch populations in the north (upper panel) and the south (lower panel) basins of Windermere. Values are given as point estimates (solid lines) with 95% confidence limits (grey areas). for Windermere pike (Craig, 1996; Haugen et al., 2007) of about 0.3 [year−1 ]. Furthermore, the higher natural mortality in female pike reported in this paper was also reported by Haugen et al. (2007) from extensive capture-mark-recapture data from the same population. The fact that we estimate natural mortalities close to previously reported values derived from capture-mark-recapture analysis indicates that our method is robust in estimating natural mortality. 3.2. Perch demography The resulting point estimated populations with 95% confidence bands for the two basins are shown in Fig. 4. Parameter estimates for the parameters are presented in Table 2. The estimated natural mortalities are quite large, but they are in fact smaller than earlier reported values for Windermere perch of 0.41 [year−1 ] and larger (see e.g. Craig, 1987, and references within). The perch stock severely declined during the beginning of the time-series (Fig. 4), especially in the north basin on which we will focus (Fig. 5). The perch decline in the north basin was possibly due to the fishery because it corresponds to a period of intense fishing effort. The number of fish caught in the scientific and commercial fisheries was never more than 30% of the total population (Fig. 5), suggesting at first sight that the fishery could not be responsible for the perch decline. However, the perch decline was driven by a failing recruitment (Fig. 5), which corresponds to a strong asymmetry in the ratio of mature males to mature females and a much stronger selectivity (more than 10 times larger) of mature males in the fishery. Hence, we specu- Ø. Langangen et al. / Fisheries Research 109 (2011) 131–139 136 Table 2 Point estimates (P.est.) of the parameters in the two basin-specific (North and South) perch models with 2.5% and 97.5% confidence limits deduced from the 200 bootstrapping samples. Male Myear and Female Myear denote male and female natural mortality in the relevant time period, respectively, and Male qage and Female qage denotes the catchability (×103 ) for males and females for the relevant age group, respectively. Parameter North Male M1976− Male M1976 Male M1976+ Female M1976− Female M1976 Female M1976+ Male q2 (×103 ) Male q3+ (×103 ) Female q2 (×103 ) Female q3+ (×103 ) South 2.5% P.est. 97.5% 2.5% P.est. 97.5% 0.277 5.03 2.91 0.0991 2.72 2.26 0.590 3.16 0.0569 0.199 0.381 6.21 3.12 0.212 4.10 2.51 0.896 3.79 0.0944 0.308 0.508 11.9 3.54 0.332 11.5 3.01 1.28 4.56 0.130 0.452 0.168 7.62 2.79 0.234 4.54 1.77 0.505 2.00 0.0289 0.0967 0.299 11.7 2.87 0.304 10.3 1.86 0.677 2.47 0.0424 0.139 0.404 11.9 3.19 0.402 11.9 2.20 1.01 3.25 0.0608 0.219 late that the “recruitment collapse” could be due to large portions of unfertilized eggs due to overfishing of male perch. A similar pattern occurred in the south basin but less severely, possibly due to a lower catchability of mature males in the south basin (Table 1). This failing recruitment was probably also affected by other factors since recruitment did not recover immediately after the sex ratio had recovered (Fig. 5), but our results illustrate a possible but often neglected aspect of gear selectivity in fisheries research. Sex-selective fisheries are one of the possible selections that might profoundly affect ecosystem based fisheries management (Zhou et al., 2010). Nevertheless, the effect of skewed sex ratio on the reproductive output of fish has only recently been investigated in some detail (see e.g. Ballon et al., 2008) and such results are usually scarce. Our results therefore add to our understanding of the population dynamical effects of sex-selective fisheries. 6•105 Number of perch Perch AIC 52.0 9.7 3.3 1.3 0.4 0.0 113.1 93.6 53.3 8.3 5.6 0.0 *Equivalent degree of freedom of the smoother(s) = 1, indicating a linear relationship. R = ˇ0 + f1 (S) + ˇ1 T 2 + ˇ2 Pe + ˇ3 P + ˇ4 Pe × P + e (5) where ˇ0 is model intercept and the other ˇi are slopes of the linear predictors, f2 is a natural cubic spline with 3 knots (Wood, 2006b), Pe stands for numbers of age 2–6 perch experienced by pike from birth to recruitment, T2 is number of degree-days above 14 ◦ C during second year of life, and et are normally distributed, autoregressive errors 5 3•10 2•105 1•105 0 et = et−1 + at 1.0 4•10 5 0.8 5 0.6 2•10 5 0.4 Sex ratio Number of perch Pike Model formulation R = ˇ0 + f1 (S) + ε R = ˇ0 + f1 (T2) * + ε R = ˇ0 + f1 (S) + ˇ1 T2 + ˇ2 Pe + ε R = ˇ0 + f1 (S) + ˇ1 T2 + ˇ2 Pe + ˇ3 P + ε R = ˇ0 + f1 (S) + ˇ1 T2 + ε R = ˇ0 + f1 (S) + ˇ1 T2 + ˇ2 Pe + ˇ3 P + ˇ4 Pe × P + ε R = ˇ5 + f2 (S) * + ε R = ˇ5 + f2 (T1) * + ε R = ˇ5 + f2 (S, Pi) + ˇ6 T1 + ε R = ˇ5 + f2 (S) * + f3 (Pi) * + f4 (T1) * + ˇ6 P + ˇ7 Bas + ε R = ˇ5 + f2 (S, T1) + ˇ6 Pi + ˇ7 P + ˇ8 Bas + ε R = ˇ5 + f2 (S, Pi) + ˇ6 P + ˇ7 T1 + ˇ8 Bas + ε In pike, the best supported GAM was: 4•105 1•10 Species 3.3. Pike stock/recruitment relationship 5•105 3•10 Table 3 SR model formulations and associated AIC for a relevant subset of formulations tested during the model selection procedure. AIC provides AIC difference with the AIC of the best model. 5 0.2 0 0.0 1945 1950 1955 1960 Year Fig. 5. Total perch population in the north basin (upper panel, thick solid line) of Windermere, number of mature males (upper panel, dotted line) and females (upper panel, dashed line), recruitment (lower panel, dashed line), total catch in number of fish (lower panel, dotted lines with diamonds) and mature male to female ratio are shown (lower panel, thick solid line). (6) where  is the autoregressive parameter and at is a homoscedastic noise term centred at 0 and assumed to be independent of et−1 . The whole model incorporating Eqs. (5) and (6) for Gamma-distributed R was fit by penalized quasi likelihood and returned approximate maximum likelihood estimates (Wood, 2006a). A relevant subset of the tested models and associated AIC is provided in Table 3. The  parameter estimate was 0.497. In our analysis we deleted the last two data points of the time series (2001 and 2002, n = 54) that were strong outliers, possibly due to population estimation problems at this extreme of the time series. Model parameter estimates for pike SR are provided in Table 4. Pike SR relationship is almost linear but shows weak densitydependent compensation (Fig. 6). Increased water temperature stimulated pike recruitment (Table 4, T2 effect), as previously reported for pike (Le Cren, 1987) and marine species (Myers, 2001). The underlying mechanisms of this positive thermal effect often remain overlooked. The explanation might involve changes in primary productivity and size-dependent metabolic constraints, Ø. Langangen et al. / Fisheries Research 109 (2011) 131–139 137 Table 4 SR models, approximate maximum likelihood parameter estimates and their summary statistics for Windermere perch and pike. S: stock of spawners, P: Pathogen factor corresponding to the absence (period 1 before 1976) or presence (period 2 after 1975) of the perch pathogen, T1 or T2: water temperature (sum of degree days above 14 ◦ C) during first or second year of life, Pi: age 3–8 pike numbers, Pe: age 2–6 perch numbers. The main effects of Pe, and P in the pike model were estimated from a model with no Pe × P interaction. N Species Equation number Pike 54 (5) Perch 114 (7) Element Parameter Estimate SE of the estimate p value Intercept S T2 Pe P (2 vs. 1) Pe × P Intercept S × Pi T1 P (2 vs. 1) Bas (south vs. north) ˇ4 f2 ˇ5 ˇ6 ˇ7 ˇ8 ˇ0 f1 ˇ1 ˇ2 ˇ3 6.90E+00 NA 2.53E−03 −4.82E−07 2.75E−01 3.00E−06 9.21E+00 NA 7.66E−03 1.68E+00 5.55E−01 1.40E−01 NA 3.28E−04 8.84E−07 1.86E−01 1.60E−06 2.71E−01 NA 8.65E−04 4.27E−01 2.40E−01 <0.0001 3.67E−03 6.63E−10 5.88E−01 1.46E−01 6.76E−02 0.00E+00 5.36E−04 2.32E−14 1.55E−04 2.30E−02 which have been well explored in perch (see below). The effect of perch on pike recruitment was disease-dependent, possibly due to size-dependent ecological interactions (Edeline et al., 2008). Indeed, perch had a weakly negative effect on pike recruitment before the disease, but a positive effect after the disease (Fig. 7). This perch × disease interaction corresponds to perch age- and sizetruncation by the pathogen (Edeline et al., 2008), suggesting that large perch compete with small pike for juvenile perch and other resources, and that perch and pike interactions involve a mixture of predation and competition typical of intraguild predation (Holt and Polis, 1997; Mylius et al., 2001). 3.4. Perch stock/recruitment relationship In perch we followed the same model selection procedure as for pike. The best-supported model structure was: R = ˇ5 + f2 (S, Pi) + ˇ6 P + ˇ7 T 1 + ˇ8 Bas + ε, (7) recruitment, P is a pathogen effect (P1: pre-disease period before 1976, P2: post-disease period after 1975), T1 is number of degreedays above 14 ◦ C during first year of life, Bas is an effect of lake basin on mean R level (we found no support for a basin-specific form of the SR relationship). The tensor product approach was here possible due to the larger number of data points (n = 114). The  parameter estimate was 0.345. Our model detected a nonlinear, threshold predator–prey interaction between pike and perch (Fig. 8, f1 effect in Table 4). When pike density is low, perch SR clearly shows compensation but hardly shows overcompensation, suggesting Beverton and Holt dynamics. When pike density is high, perch recruitment sharply decreases and the SR relationship changes to a positive exponential-like shape indicative of depensation, also known as the Allee effect (Fig. 8). Therefore, our results suggest that pike might drive perch to depensation, and that the pike–perch predator–prey system might be prone to alternative stable states with catastrophic collapses when 3000 2500 2000 500 −0.5 1000 1500 Pike recruit numbers 0.5 0.0 Pike recruits 1.0 3500 4000 where f2 is a tensor product smooth term with 3 knots, Pi stands for mean age 3–9 pike number experienced by perch from birth to 0 50000 100000 150000 200000 250000 Perch numbers 2000 4000 6000 8000 10000 12000 14000 Pike stock Fig. 6. Nonlinear effect of pike stock (numbers of age 3–9 pike in years t) on pike recruitment (numbers of age 3 pike in years t + 3) with 95% confidence limits, in partial residual units. Dots represent partial residuals from the model, i.e., residuals that would be obtained by dropping the focal term from the model while leaving all other estimates fixed. Equivalent degrees of freedom of the smoother = 1.636. Fig. 7. Linear interaction between the perch disease and the perch (average numbers of age 2–6 perch from year t to year t + 3) on pike recruitment (numbers of age 3 pike in years t + 3). Dots represent the raw data (filled circles: before the disease, open squares: after the disease), while lines represent predicted effects from model (4) with 95% confidence limits (dashed lines: before the disease, solid lines: after the disease). Predictions were obtained by letting the focal variables varying from their minimum to their maximum value while fixing other variables to their mean value. Apparent lack of fit is due to correction for the effects of other covariates by the model. Ø. Langangen et al. / Fisheries Research 109 (2011) 131–139 138 AIC was increased from 303.7 to 307.7. Therefore, we chose to keep the most parsimonious model in our SR analysis. 4. Conclusions ruits h rec Perc Pe rch rs c sto be k ke um n Pi Fig. 8. Surface plot for the nonlinear interaction between perch stock (numbers of age 2–6 perch in year t) and pike numbers (mean numbers of age 3–9 pike in years t and t + 1) on perch recruitment (numbers of age 2 perch in year t + 2). The surface was predicted from model (6) by letting focal variables varying from their minimum to their maximum value while fixing other variables to their mean value. Equivalent degrees of freedom of the tensor product smoother = 6.079. pike density exceeds a threshold value. Predators can drive an Allee effect in their prey if (i) they have a type I (linear) or type II (saturating) functional response without a type III (sigmoid) aggregative response (predators aggregate on patches of high prey density) or vice versa, (ii) predation is the main driver of prey dynamics, and (iii) prey have little spatial or temporal refuge from predation (Gascoigne and Lipcius, 2004). Pike functional and aggregative responses are not known, but Windermere perch probably have indeed little spatial and temporal refuge from predation by pike. The form of the recruitment surface in the middle area fits only a very few data points and its shape should not be taken as reflecting any biological process. The disease had a positive effect on perch recruitment (P effect, Table 4), suggesting that increased mortality from the pathogen (i) selected for increased energy allocation to reproduction (Edeline et al., 2008; Ohlberger et al., 2011a) and/or (ii) decreased competition and cannibalism among perch. Increased temperatures had a positive effect on perch recruitment (T1 effect, Table 4), in line with previous literature (Craig and Kipling, 1983; Kipling, 1984; Le Cren, 1987; Paxton et al., 2004). This suggests that increased temperatures may (i) increase lake primary productivity resulting in increased adult fecundity and/or decreased juvenile mortality, and/or (ii) increase competitive ability of juvenile against adult perch. Indeed, high temperatures favour small perch because small individuals become more efficient at exploiting a common resource with increasing temperature, which makes them competitively superior over large ones (Ohlberger et al., 2011b and references therein). Finally, we found that perch recruitment was higher in the eutrophic south basin than in the colder, deeper, mesotrophic north basin (Bas effect, Table 4), possibly highlighting some bottomup control on perch recruitment that was not captured by our data and models. Taken together, our results suggest that perch is not only top-down but also bottom-up controlled, possibly due to a lag in the predator (pike) response to prey (perch) increase. Inclusion of a sex-ratio effect in our perch SR model (Eq. (7)) showed a positive effect of an increased proportion of males on R (p = 0.0487), supporting our hypothesis. However, this sex-ratio effect did not change the form of perch SR surface (Fig. 8), and model Age-structured models are commonly used in fish stock assessment (Megrey, 1989). Fitting such models to catch-at-age time-series by estimation of the parameters is sometimes considered problematic due to confounding effects between fishing mortality, natural mortality, catchability and recruitment (Quinn and Deriso, 1999). Common objections are usually based on poor data quality and not on intrinsic problems with the methodology (Quinn and Deriso, 1999). The robust fitting of age-structured models requires long term data, a contrast in population abundance, low recruitment variability, high contrast in harvest rate and high natural mortality (Yin and Sampson, 2004; Magnusson and Hilborn, 2007; Aanes et al., 2007). These requirements are all met by our data. Additionally, we have thoroughly tested the method using simulated data and resampling methods to further back up our results. Our population estimates expand our knowledge of perch and pike dynamics, provide reference points for management, and allowed us to show that perch and pike recruitment variability in Windermere is simultaneously driven by density-dependent (stock of spawners, pathogens, competitors, predators and prey) and density-independent factors (temperature) through both bottom-up (primary productivity, basin effect) and size-dependent top-down (predation, competition) processes. In particular, our models show that predators can induce depensation in their prey, i.e. can favour prey collapse and make the predator–prey system fundamentally unstable. This result bears important management consequences because if depensatory dynamics exist in a population, reduced exploitation may be insufficient to elicit stock recovery (Case, 2000). Finally, our perch model suggests that a fishery-induced drop of the male to female ratio is at least partially responsible for a stock decline in perch in the north basin in the late 1940s. This result highlights an aspect of fishing-induced selection that is largely neglected by researchers and managers who usually focus on size selection. In total, our results underscore the need for ecosystem-based fisheries management in which the biotic (particularly predator–prey relationships) and abiotic factors are accounted for along with the trait selectivity of the fishery. Acknowledgments We are grateful to the many individuals who have participated in the Windermere data collection over the years. We are also grateful to the Freshwater Biological Association for their joint stewardship of these invaluable data. This research received support from the Research Council of Norway. Support from the Natural Environment Research Council to the Centre for Ecology and Hydrology is also acknowledged. We would also like to thank Jaime O. Villar, L.C. Stige, T.O. Haugen, D. Hjermann, G. Ottersen and T. Schweder for valuable input during this research. Appendix A. Since the perch model does not use the reported commercial catch (in kg), we have calculated the commercial catch from our model (by taking the commercial catch from the fraction of the commercial effort, similar to Eq. (4)) and compared these with the reported catch (see Fig. Si). The model performs well for the north basin and for the early years in the south basin. However, the model has a tendency to weakly underestimate the commercial catch in some year in the south basin. Possible reasons for basin specific Ø. Langangen et al. / Fisheries Research 109 (2011) 131–139 differences could include: slight time variation in the catchability and basin interactions not taken into account in the model (e.g. migration of perch between basins, changing predation pressure due to pike migration). We have included three figures (Figs. Sii, Siii and Siv) showing the performance of the models. Appendix B. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.fishres.2011.01.029. References Aanes, S., Engen, S., Saether, B.E., Aanes, R., 2007. Estimation of the parameters of fish stock dynamics from catch-at-age data and indices of abundance: can natural and fishing mortality be separated? Can. J. Fish. Aquat. Sci. 64, 1130–1142. Ballon, M., Wosnitza-Mendo, C., Guevara-Carrasco, R., Bertrand, A., 2008. The impact of overfishing and El Nino on the condition factor and reproductive success of Peruvian hake, Merluccius gayi peruanus. Prog. Oceanogr. 79, 300–307. Bodaly, R.A., Ward, R.D., Mills, C.A., 1989. A genetic stock study of perch, Perca fluviatilis L., in Windermere. J. Fish Biol. 34, 965–967. Bucke, D., Cawley, G.D., Craig, J.F., Pickering, A.D., Willoughby, L.G., 1979. Further studies of an epizootic of perch Perca fluviatilis L., of uncertain aetiology. J. Fish Dis. 2, 297–311. Carlson, S.M., Edeline, E., Vollestad, L.A., Haugen, T.O., Winfield, I.J., Fletcher, J.M., James, J.B., Stenseth, N.C., 2007. Four decades of opposing natural and humaninduced artificial selection acting on Windermere pike (Esox lucius). Ecol. Lett. 10, 512–521. Case, T.J., 2000. An Illustrated Guide to Theoretical Ecology. Oxford University Press, Oxford. Charbonneau, P., 1995. Genetic algorithms in astronomy and astrophysics. Astrophys. J., Suppl. Ser. 101, 309–334. Craig, J.F., Kipling, C., 1983. Reproduction effort versus the environment; case histories of Windermere perch, Perca fluviatilis L., and pike, Esox lucius L. J. Fish Biol. 22, 713–727. Craig, J., 1987. The Biology of Perch and Related Fish. Croom Helm Ltd., Beckenham, Kent. Craig, J.F., 1996. Population dynamics, predation and role in the community. Fish Fish. Ser. 19, 201–217. des Clers, S., Fletcher, J.M., Winfield, I.J., Kirkwood, G.P., Cubby, P.R., Beddington, J.R., 1994. Long-term Changes in Windermere Perch and Pike Populations at the Basin Level. Institute of Freshwater Ecology, IFE Report WI/T11050d5/4. Edeline, E., Ben Ari, T., Vollestad, L.A., Winfield, I.J., Fletcher, J.M., James, J.B., Stenseth, N.C., 2008. Antagonistic selection from predators and pathogens alters foodweb structure. Proc. Natl. Acad. Sci. U. S. A. 105, 19792–19796. Efron, B., Tibshirani, R., 1998. An Introduction to the Bootstrap. Chapman & Hall/CRC, Boca Raton, p. 436. Frost, W.E., Kipling, C., 1959. The determination of the age and growth of pike (Esox lucius L.) from scales and opercular bones. ICES J. Mar. Sci. 24, 314–341. Frost, W.E., Kipling, C., 1967. A study of reproduction early life weight–length relationship and growth of pike Esox lucius L. in Windermere. J. Anim. Ecol. 36, 651–693. Gascoigne, J.C., Lipcius, R.N., 2004. Allee effects driven by predation. J. Appl. Ecol. 41, 801–810. Haugen, T.O., Winfield, I.J., Vollestad, L.A., Fletcher, J.M., James, J.B., Stenseth, N.C., 2006. The ideal free pike: 50 years of fitness-maximizing dispersal in Windermere. Proc. R. Soc. Lond. B: Biol. Sci. 273, 2917–2924. Haugen, T.O., Winfield, I.J., Vollestad, L.A., Fletcher, J.M., James, J.B., Stenseth, N.C., 2007. Density dependence and density independence in the demography and dispersal of pike over four decades. Ecol. Monogr. 77, 483–502. Hilborn, R., Walters, C.J., 1992. Quantitative Fisheries Stock Assessment Choice Dynamics and Uncertainty. Chapman and Hall, Inc., New York. Hilborn, R., Mangel, M., 1997. The ecological detective: confronting models with data. Monogr. Popul. Biol. 28. Hollowed, A.B., Bax, N., Beamish, R., Collie, J., Fogarty, M., Livingston, P., Pope, J., Rice, J.C., 2000. Are multispecies models an improvement on single-species models for measuring fishing impacts on marine ecosystems? ICES J. Mar. Sci. 57, 707–719. Holt, R.D., Polis, G.A., 1997. A theoretical framework for intraguild predation. Am. Nat. 149, 745–764. Kipling, C., Frost, W.E., 1970. A study of mortality, population numbers, year class strengths, production and food consumption of pike, Esox lucius L., in Windermere from 1944 to 1962. J. Anim. Ecol. 39, 115–157. Kipling, C., 1984. A study of perch (Perca fluviatilis L.) and pike (Esox lucius L.) in Windermere from 1941 to 1982. ICES J. Mar. Sci. 41, 259–267. 139 Kipling, C., Le Cren, E.D., 1984. Mark-recapture experiments on fish in Windermere, 1943–1982. J. Fish Biol. 24, 395–414. Le Cren, E.D., 1947. The determination of the age and growth of the perch (Perca fluviatilis) from the opercular bone. J. Anim. Ecol. 16, 188–204. Le Cren, E.D., Kipling, C., McCormack, J.C., 1977. Study of numbers, biomass and yearclass strengths of perch (Perca fluviatilis L.) in Windermere from 1941 to 1966. J. Anim. Ecol. 46, 281–307. Le Cren, E.D., 1987. Perch (Perca fluviatilis) and pike (Esox lucius) in Windermere from 1940 to 1985 – studies in population dynamics. Can. J. Fish. Aquat. Sci. 44, 216–228. Le Cren, D., 2001. The Windermere perch and pike project: an historical review. Freshwater Forum 15, 3–34. Magnusson, A., Hilborn, R., 2007. What makes fisheries data informative? Fish. Fish. 8, 337–358. Megrey, B.A., 1989. Review and comparison of age-structured stock assessment models from theoretical and applied points of view. In: American Fisheries Society Symposium , pp. 8–48. Millar, R.B., Meyer, R., 2000. Bayesian state-space modeling of age-structured data: fitting a model is just the beginning. Can. J. Fish. Aquat. Sci. 57, 43–50. Mills, C.A., Hurley, M.A., 1990. Long-term studies on the Windermere populations of perch (Perca fluviatilis), pike (Esox lucius) and arctic charr (Salvelinus alpinus). Freshwater Biol. 23, 119–136. Myers, R.A., Hutchings, J.A., Barrowman, N.J., 1997. Why do fish stocks collapse? The example of cod in Atlantic Canada. Ecol. Appl. 7, 91–106. Myers, R.A., 2001. Stock and recruitment: generalizations about maximum reproductive rate, density dependence, and variability using meta-analytic approaches. ICES J. Mar. Sci. 58, 937–951. Mylius, S.D., Klumpers, K., de Roos, A.M., Persson, L., 2001. Impact of intraguild predation and stage structure on simple communities along a productivity gradient. Am. Nat. 158, 259–276. Nelson, J.S., 2006. Fishes of the World, 4th ed. John Wiley & sons, Hoboken, NJ. Ohlberger, J., Langangen, Ø., Edeline, E., Moland Olsen, E., Winfield, I.J., Fletcher, J.M., Ben James, J., Stenseth, N.C., Vøllestad, L.A., 2011a. Pathogen-induced rapid evolution in a vertebrate life-history trait. Proc. R. Soc. Lond. B: Biol. Sci. 278, 35–41. Ohlberger, J., Edeline, E., Vøllestad, L.A., Stenseth, N.C., Claessen, D., 2011b. Temperature driven regime shifts in the dynamics of size-structured populations. Am. Nat. 177, 211–223. Paxton, C., Winfield, I.J., 2000. Some statistical aspects of the long-term gill net monitoring programme for pike Esox lucius in Windermere (English Lake District). Freshwater Forum 14, 35–50. Paxton, C.G.M., Winfield, I.J., Fletcher, J.M., George, D.G., Hewitt, D.P., 2004. Biotic and abiotic influences on the recruitment of male perch in Windermere, U.K. J. Fish Biol. 65, 1622–1642. Paxton, C.G.M., Winfield, I.J., Fletcher, J.M., George, D.G., Hewitt, D.P., 2009. Investigation of first year biotic and abiotic influences on the recruitment of pike Esox lucius over 48 years in Windermere, UK. J. Fish Biol. 74, 2279–2298. Pickering, A.D., 2001. Windermere: Restoring the Health of England’s Largest Lake. Freshwater Biological Association Special Publication No. 11. Freshwater Biological Association, Ambleside. Pierce, R.B., Tomcko, C.M., 2003. Variation in gill-net and angling catchability with changing density of northern pike in a small Minnesota lake. Trans. Am. Fish. Soc. 132, 771–779. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1992. Numerical Recipes in FORTRAN. The Art of Scientific Computing, 2nd ed. University Press, Cambridge. Quinn II, T.J., Deriso, R.B., 1999. Quantitative Fish Dynamics. Oxford University Press. Venables, W.N., Ripley, B.D., 2002. Modern Applied Statistics with S, 4th ed. Springer, New York. Walters, C.J., Martell, S.J.D., 2004. Fisheries Ecology and Management. Princeton University Press, Princeton and Oxford. Wilberg, M.J., Thorson, J.T., Linton, B.C., Berkson, J., 2010. Incorporating timevarying catchability into population stock assessment models. Rev. Fish. Sci. 18, 7–24. Winfield, I.J., James, J.B., Fletcher, J.M., 2008. Northern pike (Esox lucius) in a warming lake: changes in population size and individual condition in relation to prey abundance. Hydrobiologia 601, 29–40. Wood, S.N., 2006a. Texts in Statistical Science, Generalized Additive Models. Boca Raton, Chapman & Hall, p. 391. Wood, S.N., 2006b. Low rank scale invariant tensor product smooths for generalized additive mixed models. Biometrics 62, 1025–1036. Yin, Y.S., Sampson, D.B., 2004. Bias and precision of estimates from an age-structured stock assessment program in relation to stock and data characteristics. North Am. J. Fish Manage. 24, 865–879. Zhou, S.J., Smith, A.D.M., Punt, A.E., Richardson, A.J., Gibbs, M., Fulton, E.A., Pascoe, S., Bulman, C., Bayliss, P., Sainsbury, K., 2010. Ecosystem-based fisheries management requires a change to the selective fishing philosophy. Proc. Natl. Acad. Sci. U. S. A. 107, 9485–9489.