New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Cycles and trends in cod populations

Communicated by Robert T. Paine, University of Washington, Seattle, WA (received for review March 9, 1998)
Abstract
Yeartoyear fluctuations in fish stocks are usually attributed to variability in recruitment, competition, predation, and changes in catchability. Trends in abundance, in contrast, are usually ascribed to human exploitation and largescale environmental changes. In this study, we demonstrate, through statistical modeling of survey data (1921–1994) of cod from the Norwegian Skagerrak coast, that both short and longterm variability may arise from the same set of agestructured interactions. Asymmetric competition and cannibalism between cohorts generate alternating years of high and low abundance. Intercohort interactions also resonate the recruitment variability so that longterm trends are induced. The coupling of agestructure and variable recruitment should, therefore, be considered when explaining both the short and longterm fluctuations displayed by the coastal cod populations. Resonant effects may occur in many marine populations that exhibit this combination of traits.
 agestructured dynamics
 Gadus morhua
 Markov chain Monte Carlo estimation
 nonlinear stochastic dynamics
 stock recruitment
Most commercial fish stocks vary greatly in abundance and exhibit a combination of shortterm (yeartoyear) fluctuations and longer term trends. Different factors are usually thought to give rise to these two types of changes in population size. Trends are usually related to external forcing on the populations, especially from human exploitation (1, 2) and largescale environmental variation, such as changes in atmospheric and oceanographic circulation (3) or temperature (4). The yeartoyear variation, in contrast, is usually explained by two different classes of hypotheses, both of which are internal to the populations in the sense of affecting vital rates, either in a densityindependent or densitydependent manner. First, variable (stochastic) recruitment caused by changes in winds, currents, temperature, or resources can generate substantial variation in juvenile survivorship (3–9). For instance, the relative timing of the spring algal bloom and egg hatching is a critical determinant of food availability to the larvae and postlarvae. Because the spring bloom is itself determined by environmental conditions, this “match–mismatch” hypothesis offers an explicit mechanism for the effect of climate on year class strength (6, 10, 11). Second, densitydependent interactions between and within cohorts (in the form of cannibalism and competition) can generate population fluctuations (12–14). Such interactions can lead to cycles or quasiperiodic oscillations (15–17).
Fromentin et al. (18, 19) described patterns of abundance in a longterm survey (1919–1994) of cod (Gadus morhua L.) from the Norwegian Skagerrak coast. They demonstrated both 2 to 3year cycles and longterm trends in these time series. The cycles were hypothesized to arise from interactions within and among cohorts (18, 20). The longterm trends appeared as a regional feature but were not clearly linked to any largescale changes in either zooplankton abundance or climatic variables (19). In the present study we built a dynamic model for the Skagerrak cod survey data to explore the cycles and trends. The model includes stochastic reproduction and agestructured interactions and is designed to test whether this set of regulatory and disruptive processes can account for the variability in abundance. We asked whether a stochastic agestructured model alone can account for the dual pattern of cycles and trends in cod abundance. The parameters of the model were estimated from the timeseries data with a Bayesian Markov chain Monte Carlo method (21, 22). Markov chain Monte Carlo allows us to estimate the parameters in the nonlinear agestructured model, including the variance and the mean of the per capita reproduction rate.
MATERIALS AND METHODS
The Data.
Recruitment has long been recognized as a key parameter in year class strength of fish (5, 6). To investigate the recruitment dynamics of local cod stocks, a research survey was initiated along the Norwegian Skagerrak coast in 1919. Standardized beach seines were used to target the two juvenile cohorts of cod: the 0group and the 1group (see ref. 18 for details). Since then, the program has been running according to the same sampling protocol at more than 250 locations along the Southern Skagerrak coast of Norway. In the Southern part of the study area (between Kristiansand and Porsgrunn), records have been made every year for the period 1921–1994 (except for a gap during World War II) at 38 locations (see ref. 18 for details). In the Northern part of the study area (between Porsgrunn and the NorwegianSwedish border), 31 stations provide complete records for the period 1936–1994 (except for a gap during the World War II). The total counts of the 0group and the 1group are shown in Fig. 1A for the Southern area and Fig. 1B for the Northern area. The temporal variability in abundance is extensive at all locations. In our analyses, we will focus on the two sets of aggregate time series. The first set is the summed count of the 0group and the summed count of the 1group across the 38 Southern locations. The second set is the summed series across the 31 Northern locations. The median for the 0group and the 1group counts in the Southern set is 458 (range, 5–1,902; maximum to minimum ratio, 380) and 39 (range, 1–164; maximum to minimum ratio, 164), respectively. The corresponding medians for the Northern set are 227 (range, 13–2,325; maximum to minimum ratio, 179) and 35 (range, 1–301; maximum to minimum ratio, 301), respectively. The four preWorld War II observations for the Northern set represent such a short continuous period that they were not used in the estimation.
Descriptive analyses of the 38 time series from the Southern area have identified significant shortterm cycles and longterm trends in most populations (18). The two aggregate time series considered here (Fig. 1) reflect this dual pattern. Significant 2 to 3year cycles and significant longterm trends (here taken to mean variability in the lowest frequencies observable, corresponding to 25+ years) are seen in both the Southern (for the contiguous postWorld War II period) and the Northern area. The evidence for cycles is stronger in the 1group [permutation test of the periodogram (ref. 23, Chap. 11) by using 10,000 permutations: Southern, P = 0.03; Northern, P < 0.01] than in the 0group data (Southern P = 0.08; Northern P < 0.01). In contrast, the trends are more apparent in the 0group (Southern, P < 0.01; Northern, P = 0.02) than in the 1group (Southern, P = 0.16; Northern, P = 0.36). The combined high and lowfrequency variability lead to bimodal (Ushaped) spectra.
AgeStructured Model: The Cod Life Cycle.
Adults and juveniles are parapatric and do not interact strongly (24–26). The two juvenile age classes (the 0group, X and the 1group, Y), in contrast, have similar habitat preferences (18, 24). An essential component of the dynamics is believed to be competitive and cannibalistic interactions within and between individuals of these two classes. The reduction in survival caused by juvenile interactions is thought to be density dependent. The functional form for the density dependence in cod survival appears to be approximately log linear (20, 27). The number of individuals surviving from the 0group to the 1group (survival during the first year) may, therefore, be approximated as 1 where β (<0) represents densitydependent reduction in survival caused by intracohort competition and γ (<0) represents density dependence caused by intercohort competition and cannibalism (Fig. 2). The former, thus, embraces the interaction among 0group individuals, and the latter embraces the interaction between 0group individuals and 1group individuals. The constant c (<0) represents the level of annual densityindependent survival.
Individuals mature at 2–3 years, and the annual survival rate, λ, of adults appears to be independent of density and relatively constant across years (25). The average annual adult survival was estimated to be 0.34 ±0.02 (range, 0.44–0.32) in a mark–recapture study involving 40,000 marked individuals during 1986–1989 (28). Fishing mortality is known to have been higher during the late 1980s (28); an adult survival of approximately 0.4 is therefore likely to be representative for these populations. A constant adult survival of 0.4 implies that the group 2, Z, and group 3, W, together make up more than 80% of the spawning stock. In agreement with this, Gjøsæter et al. (25) found these two age groups to make up approximately 75% of the spawning stock of local populations. Therefore, the two dominant adult classes are given by 2 where λ is set to 0.4, and the reproductive stock in year t may be approximated as Z_{t} + W_{t}.
Spawning is usually in March in the Norwegian Skagerrak cod, and the eggs hatch into larvae after a few weeks. The larvae metamorphose into young fish (the 0group individuals) in May/June. After a pelagic phase, the 0group individuals settle and feed on the bottom when they are about 3–5 cm long (around August). Thus, the numbers of 0group individuals counted in the fall arise from reproduction by the adult stock in the spring of the same year. Approximating the stock by the 2 and 3year olds, the number of the 0group recruits, X, will be 3 where e^{α}t is the reproduction rate. We define the reproductive process as the transition from eggs through fish larvae to pelagic and finally demersal 0group individuals (thus including the settlement on the bottom). This extended reproductive process is highly variable between years. To allow for variability between years in reproduction, we assume that α_{t} is drawn from a Gaussian distribution with mean ᾱ and variance σ_{α}^{2}. This leads to a per capita reproductive rate, e^{α}t, that is lognormally distributed across years (27, 29, 30). To understand the dynamic effects and importance of reproductive variability, we also consider the contrasting situation: α_{t} is constant across years (denoted α_{const}), and estimate the parameters under this assumption.
Under the heuristic assumption of constant recruitment, equations 1–3 describe a deterministic nonlinear dynamic system that produces cyclic or dampened cyclic dynamics, depending on the absolute and relative strength of inter and intracohort interactions. Limit cycles will arise if the intercohort interaction is substantially larger than intracohort interaction. For adult survival of 0.4, limit cycles arise if the intercohort interaction γ < −0.6 + 0.4β (unpublished numerical results). With deterministic reproduction, all variability in the data is because of measurement errors inherent in the counting of individuals (approximated by Poisson variability) and because of the intrinsic dynamic (deterministic) variability (31). Statistically speaking, the census counts arise from a nonlinear and nonGaussian process with fixed coefficients (32): α_{const}, β, γ, and c.
When the reproductive rate is affected by changing environmental conditions, there is an additional level of variability: the processinherent stochasticity (see, for example, refs. 30 and 33). Under the stochastic recruitment model, the census counts arise from a stochastic nonlinear and nonGaussian process with one random coefficient (32), α_{t} (which we assume to be approximately normally distributed with mean ᾱ and variance σ_{α}^{2}).
For a known λ, the parameters α_{const}, ᾱ, σ_{α}^{2}, β, γ, and c will be fully specified by the census data for 0group and 1group (Fig. 1). The parameters can be estimated by using the Markov chain Monte Carlo method implemented in bugs, version 0.6 (21, 22, 34). We base our estimates (and 95% credible intervals (ci)—the Bayesian alternative to the confidence interval) on 100,000 iterations of the Gibbs sampler after a burnin of 10,000 iterations.
RESULTS
The parameter estimates demonstrate that there is both inter and intracohort density dependence in the survival of the juvenile cod (Table 1). The parameters β and γ are negative and their 95% ci does not include 0. The parameters for the Southern and Northern areas are very similar. The baseline juvenile survival, e^{c}, does not differ significantly between the two areas and is estimated to be about 0.6 (Southern, 0.56; Northern, 0.68; Table 1). The estimated agestructured interactions greatly reduce the survival. In years of median abundance of juveniles (for comparison, we use the grand median, which is 37 for 1group and 333 for the 0group), the survival of juveniles is estimated to be approximately 0.1 (Southern, 0.10, ci = {0.09, 0.11}; Northern: 0.15, ci = {0.14, 0.15}). Fig. 2 shows the model of the 1group abundance for the Southern area. When the four most outlying observations are omitted from the analyses, the betweencohort interaction is of comparable strength to the withincohort interaction (β ≈ γ ≈ −0.2).
The per capita reproductive parameter in the constant recruitment model, e^{α}const, predicts that an adult produces approximately 10 eggs that successfully hatch, metamorphose, and enter the 0group (Southern, 15; Northern, 8). The constant recruitment model predicts dampened oscillations in abundance, because the intercohort interactions are too weak to induce stable limit cycles (see “AgeStructured Model”). The spectral density of the model (including the transient period) is Lshaped (with a peak at 2 years). Arguably, this model may be seen as capturing the shortterm fluctuations in cod abundance (the 2 to 3year periodicity). It does, however, fail to account for any lowfrequency behavior. Some essential part of the cod dynamics appears to be missing from the deterministic model (35).
Allowing the reproduction to be variable does not affect estimates of the interaction strengths. It does, however, demonstrate significant interannual variability in the per capita reproductive rate (Table 1). The mean reproductive rate does not differ significantly from the constant reproduction model. The reproductive variance, σ_{α}^{2}, is approximately unity and is significantly larger than 0 in both areas (Table 1). The coefficient of variation in α_{t} is estimated to be 0.41 (with ci = {0.34, 0.51}) for the Southern area and 0.63 (with ci = {0.48, 0.84}) for the Northern area. Typical years may be considered those that fall within 1 SD of the mean. The estimated variability is such that the per capita production of eggs that successfully hatch, metamorphose, and enter the 0group may vary by a factor of 5 between two such typical years (ᾱ ± σ_{α}: Southern, {4.6, 38.0}; Northern, {3.1, 15.0}).
The stochastic model has a spectral density that differs nontrivially from its deterministic counterpart (Fig. 3A). As in the deterministic model, there is a strong cyclic component corresponding to a 2 to 3year oscillation. However, there is an additional lowfrequency component to the dynamics. The bimodal spectrum of realizations of the stochastic model (Fig. 3A) bears resemblance to that of the data. Repeated simulations (across 50 iterations corresponding to the numbers of contiguous years in the data) and spectral estimation from the stochastic model reveal that the model can account for the bimodal spectrum of the coastal cod populations. Thus, both the cyclicity and the lowfrequency variability observed in the cod dynamics can be predicted from the model. The observed spectrum in the data is enveloped by the 95% distribution of spectra from the stochastic model (Fig. 3). When the model is run for 5,000 iterations with the Northern parameter estimates (Table 1), the dominant period is always the lowfrequency behavior with an apparent period of more than 100 years (median across 1,000 simulations, 175 years; interquartile range, {101, 370}). The subdominant period is the 2year cycle (median, 2.1 years; interquartile range, {2.0, 2.2}). The model for the Southern area has a similar behavior (median dominant period, 153 years; interquartile range, {88, 313}; median subdominant period, 2.2 years; interquartile range, {2.1, 2.3}).
DISCUSSION
We studied a simple model for cod dynamics that included densitydependent interactions and interannual variation in reproductive rates. The parameters in the model were estimated from survey data from the Norwegian Skagerrak coast between 1921 and 1994. We demonstrated here a significant inter and intracohort density dependence in juvenile mortality, and estimated the magnitude of the recruitment variation. By using the parameter estimates, we showed that the model will give rise to shortterm cycles in abundance superimposed on longterm trends when there is reproductive variability.
The estimated strength of the density dependence in juvenile survival of cod (Table 1) is comparable to estimates previously reported for other cod populations (27). The estimated per capita production of juveniles (approximately 8–15 0group juveniles per adult) also compares well with previous reports (36, 37). One of the main features of the demography of the Norwegian Skagerrak cod is a low age of maturation (2–3 years) and a short lifespan (<5 years) (18, 25). A cycle in abundance caused by interactions in agestructured populations is typically related to time to maturation (ref. 38, but see ref. 17). Any cycles are therefore likely to have longer period in Atlantic cod populations that have longer lifespans and slower maturation rates. The qualitative aspects of the dynamics, i.e., the possibility of a bimodal spectrum with both cycles and trends, are, however relatively robust to time to maturation and lifespan (unpublished numerical analyses). Loglinear density dependence in the mean reproductive rate or adult survival rate does not affect the type of dynamics observed (26). The estimation is insensitive to the assumption of an adult survival rate of 0.4. Assuming higher or lower survival in the estimation will result in altered (lower or higher, respectively) estimates of the mean per capita reproductive rates, but no other estimates are affected. The predicted dynamics are unchanged. The critical feature appears to be the combination of variability in reproduction and asymmetric intercohort interactions.
Simple deterministic (35) or stochastic (30, 33) models of fluctuating populations typically predict that spectra will be dominated by high frequencies (to be “blue shifted”) (35, 39, 40). The discrepancy between the prediction of these simple models and the “multicolored spectra” (i.e., also containing lowfrequency signals) of real populations has received recent attention (35, 39–41). Spatial dynamics (40) or longterm stochastic forcing (refs. 42 and 43, but see ref. 41) has been invoked to account for the discrepancy. This study shows that stochastic (white noise) recruitment can be echoed through the age structure to produce a twopeaked spectrum (cycle and trend) from a single set of densitydependent mechanisms. Disruptive forces (such as stochastic reproduction) have previously been shown to be capable of sustaining transient periodicity in population models (44–46). Here we show that disruptive forces can interact with regulatory mechanisms to produce dynamics that differ qualitatively from both the purely stochastic and the purely deterministic analogs.
The effect that (uncorrelated) stochasticity in reproductive rates can generate lowfrequency variability in dynamics may appear counterintuitive because the deterministic agestructured model exhibits highfrequency cycles only. Stenseth et al. (26) investigated a simplified version of the model (equations 1–3) with only one adult age class and showed that this may be written in delay coordinates to obtain a model of the form: 4 where, x_{t} = log(X_{t}), η is a constant, and the other symbols are as defined above. The random component of this model, α_{t} − γα_{t−1}, represents an explicit formulation of how random reproduction, {α_{t}}, is dynamically echoed in time (recall that γ < 0). Ecologically speaking, the effect is caused by the intercohort interaction (γ) in the agestructured dynamics operating as a “resonator” for the variability in recruitment.** Statistically speaking, however, the random component represents a noncentered moving average process. Moving average processes are generally capable of inducing lowfrequency behavior (49). Therefore, this effect is the dynamic cause of the presence of lowfrequency signals in the model trajectory. Because the intercohort interaction is estimated to be fairly weak and the reproductive variance fairly high, the resonant effect is relatively weak in the model realizations (Fig. 3). To illustrate the effect more fully, we simulate a system with strong interactions (β = −0.3, γ = −0.9) and low reproductive variance (ᾱ = 2, σ_{α}^{2} = 0.01). This set of parameters yields a low coefficient of variation in the per capita reproduction rate (the coefficient of variation of e^{α}t is 0.1). Fig. 4 shows 100 years of simulation (after a transient period of 50 years) of two stochastic realizations and one deterministic realization (α_{const} = 2). The lowfrequency oscillations arising from the resonance is subdominant in power to the 2year cycle for this set of parameters but nevertheless highly conspicuous.
The agestructured model may account for both the cycles and trends in abundance of the Skagerrak cod populations. However, external factors such as pollution, harvesting, and climatic changes should not be ruled out. Such factors have already been shown to account for some of the variability in fish populations (1–3, 19, 50). Community interactions are also likely to be of importance (13). Our main aim here is to point to an interesting alternative source of lowfrequency behavior that deserves future study. Future studies must be undertaken to understand which of these causal mechanisms are the most important in governing the abundance of the cod populations along the Norwegian Skagerrak cod.
Asymmetric interactions between age or stage classes are common in many animal life cycles (51) and are a major source of nonlinearity and of dynamic complexity in deterministic population models (16, 52, 53). We show that when agestructured dynamics occur in a varying environment, both cycles and longterm trends may arise. Nonlinear interactions combined with stochastic reproduction can induce resonant effects that represent a neglected level of dynamic richness. Separating anthropogenic and environmentally induced trends from this intrinsic lowfrequency variability is a challenge that will require more attention.
Acknowledgments
The late Rangvald Løversen and Aadne Sollie carried out the survey every year for more than 70 years. Tore Johannessen provided help with the data. Rolf A. Ims, Fiorenza Micheli, and three anonymous reviewers commented on the manuscript. We received financial support from the Norwegian Science Foundation (O.N.B. and N.C.S.), the University of Oslo (N.C.S.), and the Marine Sciences and Technologies program of the European Union (J.M.F.).
Footnotes

↵‡ email: bjornsta{at}nceas.ucsb.edu.

↵¶ To whom reprint requests should be addressed. email: n.c.stenseth{at}bio.uio.no.

↵** This mechanism is not a form of stochastic resonance in which the signal to noise ratio increases with increasing noise levels (47, 48). The signal to noise ratio of the present model decreases monotonically with reproductive variance.
ABBREVIATION
 ci,
 Bayesian credible intervals
 Received March 9, 1998.
 Accepted January 21, 1999.
 Copyright © 1999, The National Academy of Sciences
References
 ↵
 ↵
 ↵
 Dickson R R,
 Brander K
 ↵
 ↵
 Hjort J
 ↵
 Cushing D H
 ↵
 Sundby S,
 Bjorke H,
 Soldal A V,
 Olsen S
 ↵
 Ellersten B,
 Fossum P,
 Solemdal P,
 Sundby S
 ↵
 Sundby S,
 Fossum P
 ↵
 Blaxter J H S
 May R C
 ↵
 ↵
 Fortier L,
 Villeneuve A
 ↵
 ↵
 Bailey K M,
 Houde E D
 ↵
 Cushing J M
 ↵
 Cushing J M
 ↵
 ↵
 Fromentin JM,
 Stenseth N C,
 Gjøsæter J,
 Bjørnstad O,
 Falk W,
 Johannessen T
 ↵
 Fromentin JM,
 Stenseth N C,
 Gjøsæter J,
 Johannesen T,
 Planque B
 ↵
 Bjørnstad O N,
 Fromentin JM,
 Stenseth N C,
 Gjøsæter J
 ↵
 Best N G,
 Spiegelhalter D J,
 Thomas A,
 Brayne C E G
 ↵
Gilks, W. R., Richardson, S. & Spiegelhalter, D. J. (1996) (Chapman & Hall, London).
 ↵
 Manly B F J
 ↵
 ↵
 Gjøsæter J,
 Enersen K,
 Enersen S E
 ↵
Stenseth, N. C., Bjørnstad, O. N., Falck, W., Fromentin, J.M., Gjøsæter, J. & Gray, J. S. (1999) Proc. R. Soc. London Ser. B, submitted.
 ↵
 ↵
Julliard, R., Stenseth, N. C., Gjøsæter, J., Lekve, K., Fromentin, J.M. & Danielsen, D. S. (1999) Ecol. Appl., submitted.
 ↵
 Cushing J M,
 Dennis B,
 Desharnais R A,
 Costantino R F
 ↵
 ↵
 ↵
 Longford N T
 ↵
 ↵
 Spiegelhalter D J,
 Thomas A,
 Best N G,
 Wilks W R
 ↵
 ↵
Myers, R. A., Mertz, G. & Barrowman, N. J. (1996) Int. Counc. Explor. Sea Tech. Rep. D 11.
 ↵
 Myers R A,
 Barrowman N J,
 Hutchings J A,
 Rosenberg A A
 ↵
 ↵
 Sugihara G
 ↵
 White A,
 Begon M,
 Bowers R G
 ↵
 Kaitala V,
 Ylikarjula J,
 Ranta E,
 Lundberg P
 ↵
 ↵
 ↵
 Kaitala V,
 Ranta E,
 Lindstrom J

 Higgins K,
 Hastings A,
 Sarvela J N,
 Botsford L W
 ↵
 Stenseth N C,
 Bjørnstad O N,
 Saitoh T
 ↵
 ↵
 ↵
 Priestley M B
 ↵
 ↵
 Tuljapurkar S,
 Caswell H
 ↵
 ↵
 Costantino R F,
 Desharnais R A,
 Cushing J M,
 Dennis B