Seasonal plankton dynamics along a cross-shelf gradient

Nils Chr Stenseth, Marcos Llope, Ricardo Anadón, Lorenzo Ciannelli, Kung-Sik Chan, Dag Ø Hjermann, Espen Bagøien, Geir Ottersen


Much interest has recently been devoted to reconstructing the dynamic structure of ecological systems on the basis of time-series data. Using 10 years of monthly data on phyto- and zooplankton abundance from the Bay of Biscay (coastal to shelf-break sites), we demonstrate that the interaction between these two plankton components is approximately linear, whereas the effects of environmental factors (nutrients, temperature, upwelling and photoperiod) on these two plankton population growth rates are nonlinear. With the inclusion of the environmental factors, the main observed seasonal and inter-annual dynamic patterns within the studied plankton assemblage also indicate the prevalence of bottom-up regulatory control.


1. Introduction

The structure of marine plankton food webs and therefore the linkages between trophic levels are related to the recurrence of mesoscale processes, especially when these cause the resupply of new nutrients from the subsurface reservoir (Legendre & Rassoulzadegan 1996; Falkowski et al. 1998). As a temperate sea, the deep winter mixing is the main annual input of nutrients to the upper layers in the southern Bay of Biscay. However, on the coast, pulses of upwelled nutrients are an important contribution during the summer stratification (Llope et al. in press). This extra input, along with the nutrient loading from rivers, establishes a cross-shelf gradient in the elemental stoichiometric ratios, indicating the development of different communities under different environmental conditions (Llope et al. submitted). The classical diatom–zooplankton–fish (DZF) food web would be favoured towards the coast during most of the year, but more restricted to the spring bloom offshore (Kiørboe 1993). Conversely, the persistent summer depletion of nutrients enables the recycling microbial food web to have a competitive advantage over the DZF in the open ocean (Azam 1998; Legendre & Rivkin 2002; Mahaffey et al. 2004).

Under this complex scenario of physical and biogeochemical variability, we address the issue of plankton growth seasonal dynamics as well as the type of coupling between trophic levels in three stations along a coast–ocean gradient (figure 1). We used real phyto- and mesozooplankton biomass estimates to assess the effect of environmental factors, such as the amount of various nutrients, temperature, photoperiod and upwelling. Our approach differs from the previous nutrient–phytoplankton–zooplankton models, in that we do not pre-define the formulation of the biological processes (e.g. Gibson et al. 2005 and references therein), but let the data ‘tell us’ how the change in plankton biomass is determined by various environmental factors, and then use the information available in the literature to explain the obtained relations on the basis of possible mechanisms.

Figure 1

The sampling area, Cantabrian Sea (South Bay of Biscay) and the plankton dynamics. Three permanent stations off the Asturian coast. st. 1 (43°36′ N, 06°08′ W, 65 m), st. 2 (43°42′ N, 06°09′ W, 135 m) and st. 3 (43°46′ N, 06°10′ W, 870 m). Phytoplankton is shown in green (right axis, circles) and zooplankton in blue (left axis, triangles).

Our analysis builds upon earlier methodological work within the tradition of a statistical modelling of the time-series data, and on this basis deduces the possible underlying ecological model (Royama 1992; Bjørnstad et al. 1996; Kendall et al. 1998; Stenseth et al. 1999).

2. Material and methods

(a) The time-series

We sampled phytoplankton (measured as chlorophyll a) and mesozooplankton (more than 200 μm) biomass at three stations along an inshore–offshore gradient. Station 1 (st. 1) is a shallow coastal station influenced by freshwater discharges, tidal currents and frequent wind-driven upwelling during summer. Station 2 (st. 2) is located on the continental shelf, and as such is also affected by upwelling events, although less intensively than st. 1. Station 3 (st. 3) is near to the slope and is the most oceanic, as it is only marginally affected by coastal processes, except for the indirect effect of upwellings, through inshore–offshore advection (Llope et al. in press). The sampling methodology is in the electronic supplementary material.

(b) Statistics and model structure: a generalized additive model approach

The time-series data were analysed using generalized additive models (GAMs), as implemented in the mgcv library of R (Hastie & Tibshirani 1990; Wood 2000). GAMs are non-parametric regression techniques and therefore do not require an a priori specification of the underlying functional forms between dependent and independent variables. The GAM regression technique consists of fitting smooth additive functions for each covariate included in the model structure. The smooth functions are linear combinations of a finite number of basis cubic spline functions, with the smoothness of the function estimated by minimizing the generalized cross-validation criterion (Wood 2000) that balances the goodness-of-fit and the smoothness of the functions. To avoid overfitting, we constrained the number of basis functions to be three, at most. We regressed per-month variations of phyto- or zooplankton biomass (i.e. roughly a proxy for net population growth rates) against population biomass (i.e. density dependence), abundance of the other trophic level and a selection of environmental covariates recorded at the time of plankton sampling.

The general model structure, used in the phytoplankton analysis, wasEmbedded Image(2.1)where Embedded Image is the phytoplankton biomass increase/decrease calculated as the difference in the logarithmic abundance of two consecutive phytoplankton abundances (e.g. Pt and Pt+1); fp, hp and gp,j, non-parametric, smoothing functions, specifying the effect of the phytoplankton abundance (i.e. density dependence), the zooplankton density and the environmental forcing of the covariate Ej on the phytoplankton increase rate Embedded Image, respectively; bp, an intercept; and Embedded Image, the noise term.

A similar model structure was used to study the zooplankton (Z) dynamics. The corresponding functions are fz, hz and gz,j.

(c) Covariates

We used the amount of various nutrients, temperature, photoperiod and an upwelling index as covariates. Temperature and upwelling are proxies for water-column status and offshore advection. The interpretation of nutrients and photoperiod varies depending on the trophic level.

Thus, for phytoplankton, these indicate energy availability. We introduced nitrate, phosphate and silicate as nutrients, and let each model select among them. Only phosphate and silicate entered the best models. However, the replacement of these by nitrate led to minor differences, suggesting that the preference for one or another was not essential. But this was not the case at the coastal station, where the best model made a firm selection for phosphate.

Unlike phytoplankton, zooplankton biomass may increase through processes other than the multiplication of individuals, such as the growth (in size) of individuals or the seasonal input of the meroplankton (larvae of benthic organisms or fishes). Thus, for zooplankton, photoperiod would act as a cyclic proxy accounting for some of these processes.

The concentration of nutrients decreases progressively, and not stepwise, from winter to spring (Llope et al. submitted, see also the silicate dynamics in fig. 6 of the electronic supplementary material). This is a consequence of active nutrient uptake during the whole winter, rather than only at the spring bloom period. These winter blooms last for only a few days, being too short-lived to be accurately detected by our sampling. In this sense, nutrient dynamics capture, in a smooth fashion, the more pulsing dynamics of primary production. Therefore, high nutrient concentrations may also be seen as resources for mesozooplankton. In fact, nutrients are better proxies for phytoplankton production than measurements of chlorophyll (on the scale that matters to zooplankton), and at the same time they also account for the dynamics of other non-chlorophyll resources that may be more affected by nutrients than are phytoplankton.

(d) Density dependence

The model structure includes a term for density dependence that proved to be important in all the models. In order to avoid the fallacy of overestimating the degree of within-population density dependence, we also analysed models with both phyto- and mesozooplankton abundances as the dependent variables (Solow 2001). The model structure is the same as reported in §3 (shown in figure 2); however, as expected, the effects of lagged biomass were weaker on all stations and non-significant for both phyto- and zooplankton at the coastal and shelf stations. Thus, despite most phyto- and zooplankton having fast generation times, the accumulation or loss of biomass shows an apparent density-dependent regulation, which is significant in the case of the oceanic station that is most nutrient-limited.

Figure 2

The best (i.e. selected) phyto- and zooplankton growth models for each of the three stations. The y-axis indicates the partial additive effect that the term on the x-axis has on the response variable (phyto- or zooplankton growth). The numbers in parentheses on the y-axis indicate the estimated degrees of freedom, which are also used in table 1. The final value of the response variable is given by the sum of all partial effects plus a constant. The density-dependent effect is shown first (the f-functions; see equation (2.1)), then the response of the covariates entering the models (i.e. the gj functions) and the effect of the other trophic level (when detected, i.e. the h-functions). Pt is the abundance of the phytoplankton at time t, and Zt is the corresponding abundance of the zooplankton. The terms used for the covariates are the same as those used in table 1. The broken lines give the 95% (point-wise) confidence envelopes of the predicted curves.

(e) Simulation

The skeleton (i.e. the deterministic part) of the statistically deduced models was simulated using only the observed values of the environmental covariates as given, and letting the phyto- and zooplankton be proper dynamic variables predicted by the model for each time-step, except zooplankton at station 1. Each simulation was initiated using the first observation of the focus variable as a starting point, and then updated using the model-predicted value from the previous time-step and the observed values of environmental covariates. When an environmental covariate was lacking for one or more months, we used the linear interpolation based on the last month before and the first month after the ‘break’ (no such breaks were longer than three months). We also performed simulations for an ‘average year’ based on the monthly mean values of the environmental covariates (figure 3, right columns). The 95% confidence intervals for the average year were calculated from bootstrapping. For each station, we made 1000 bootstrap samples from the data, keeping the number of observations per month equal to that of the original dataset. For each bootstrap sample, we re-estimated the GAM, recalculated the monthly means and ran new simulations as explained earlier.

Figure 3

Seasonal structure of the phyto- and mesozooplankton growth simulated from the skeleton of the selected best models. In the left column, simulations were started using the first observation of the series (in 1993), while for the rest of the series, densities were predicted based on the last month's predicted population growth (as well as observed abiotic predictor variables). During months when some predictor variables were missing, linearly interpolated covariate values were used. In the right column, we have performed simulations for an ‘average year’, starting with the mean density for January, and predicting for the remaining months using the last month's predicted population growth, as well as using the monthly mean values of the predictor variables. The broken lines in the zooplankton panels for st. 1 represent the models with the additional seasonal component (found to be significant in the statistical analysis). Solid lines represent the best model without an explicit seasonal component (these were the overall best models for all series, except for the zooplankton of st. 1). Dots represent the real values. The confidence intervals, obtained by bootstrapping, are in all cases based on the models without seasonal components. The fraction of variance explained is also given.

3. Results and discussion

Figure 2 and table 1 summarize the resulting model structures for the three stations (details regarding model selection are provided in table 2 of the electronic supplementary material). The statistically selected best models explain a major fraction of the observed variance. The apparent density dependence may represent the seasonal evolution of the environment and the consequent succession of different species/communities with different growth capabilities (Huisman et al. 2004). The highest phytoplankton biomass increase appears when there is little phyto-biomass and the level of resources is still high. These typical winter–spring conditions are more favourable for the phytoplankton to grow, especially for the diatoms. On the other hand, high biomass (i.e. the spring bloom) usually led to nutrient depletion, and therefore to future biomass reduction. Intermediate variations could correspond to summer communities, upon which the seasonal stratification imposes low change rates. A similar explanation can be given for zooplankton: periods of high biomass cannot be prolonged when the phytoplankton diatom-based community changes towards a more heterotrophic structure.

View this table:
Table 1

The structure of the statistically selected best models; Pt is the abundance of the phytoplankton at time t and Zt is the corresponding abundance of the mesozooplankton (both analysed on log-scale). (The state variables (as opposed to the covariates) are bold-faced. Approximate significance of smoothed terms is given as the p-values. Edf is the estimated degree of freedom of the examined covariate. Edf equal to 1 implies a linear effect and values greater than 1 indicate a progressively stronger nonlinear effect. An additive cosine function of 1-year period was included in the model at st. 1 for zooplankton, in order to account for some slight seasonal variability. phos, phosphate; tmp, temperature; phot, photoperiod; nit, nitrate; sil, silicate; upw, upwelling index.)

Apart from the changes in growth rate, advection must also be seen as an important process controlling the plankton biomass dynamics, since at the outermost station, both phyto- and zooplankton models include upwelling as an explanatory variable. In our analyses, we do make the assumption that we can equate change in biomass to the change in population abundance.

(a) Phytoplankton dynamics

For phytoplankton dynamics, we found either a direct or an indirect (through changes in nutrients) effect of upwelling at all stations. At the oceanic station, the effect of upwelling—detected in our analysis—may be a consequence of phytoplankton advection (Marañón & Fernández 1995) rather than an in situ enhancement of primary production; positive values of upwelling imply advection from the coast while negative values (‘downwelling’) signify advection towards the coast. Only phosphate and silicate entered the best models; however, the replacement of these by nitrate (excepting st. 1) led to minor differences among models.

Temperature enters the phytoplankton models at both st. 1 and st. 2 in a somewhat different manner. There is a positive linear effect of temperature on phytoplankton at the continental shelf station (st. 2), whereas in the coastal station (st. 1), temperature shows a U-shaped relationship with a minimum at about 15 °C. The combined inclusion of two nutrients in the model may indeed explain the linear effect of temperature at the middle station (st. 2). The interpretation of the U-shape effect of temperature at the coast is more complex and may be related to two processes, both of which enhance the phytoplankton growth. Towards the end of winter and in early spring (cool temperatures), episodic periods of high phytoplankton growth may occur in these shallow waters where the upper mixed layer would remain always above Sverdrup's critical depth (Sverdrup 1953; Huisman et al. 1999). Moreover, salinity-driven stratification owing to river runoff typically occurs during the winter season, generating a nutrient-rich mixed layer. Together, these processes might explain the slight elevation with regard to decreasing temperature of the GAM function for values below the turning point. During summer (warm temperatures), upwelling will frequently occur along the coast. Upwellings are short-lived processes in this region. However, their effect on the phytoplankton community persists longer and can be detected as biomass increases, even when water conditions have reverted to typically summer values (González et al. 2003). As neither silicate nor nitrate (upwelling related) enters the model, this might explain the positive temperature–phytoplankton relation for the higher temperatures.

The only nutrient entering the phytoplankton model at the coast is phosphate. At this station, the N : P ratio differs from the other two stations (see fig. 5 of the electronic supplementary material). During winter, there is a very small residue of nitrate (0.027 μg kg−1) at st. 1 when phosphate is depleted, while it is phosphate which is in excess at st. 2 and st. 3 when there is no nitrate (0.021 and 0.055 μg kg−1, respectively). Using year-round values, this difference is less important. This suggests that the role of phosphate would be more important on the coast, probably owing to river discharges loading nitrate in excess (Tyrrell 1999).

Interestingly, our analysis may suggest that the phytoplankton dynamics at st. 2 (the most completely sampled station) are non-additive (see fig. 6 of the electronic supplementary material). By using silicate as a threshhold variable, phytoplankton growth shows a higher degree of density dependence at high silicate levels, i.e. from January to April when stratification is low and transient blooms can result in sharp growth changes. The reduced density dependence during summer could be owing to nutrient restriction and the different communities that develop under such conditions, which are less prone to sudden variations.

(b) Zooplankton dynamics

Our analysis indicates that photoperiod, temperature and upwelling are the key environmental variables affecting zooplankton dynamics (figure 2). Mesozooplankton biomass is linearly related to photoperiod at st. 1, and asymptotically at st. 2 and st. 3. Temperature enters the zooplankton growth model negatively at st. 2, reflecting water stratification and related conditions. Upwelling enters the model for st. 3, again most probably as a result of advection from coastal water. Phytoplankton biomass positively affects zooplankton growth at the coastal station (figure 2), indicating that direct transfer of biomass—via the classical food web—occurs year-round. However, this link between the two trophic levels was not detected at the outermost stations where the severe summer stratification leads to a more oligotrophic community (Serret et al. 1999). At this time, mesozooplankton may switch their feeding towards omnivory and be less dependent on phytoplankton. Therefore, it seems reasonable to interpret silicate and nitrate levels as proxies for food resources in general, i.e. phytoplankton and the other intermediate trophic levels (see fig. 5 of the electronic supplementary material). Thus, high nutrient levels enhance zooplankton growth the following month as these are incorporated into the organic matter, while low nutrients imply low new production to be used.

(c) Trophic level interaction

Both bottom-up and top-down regulatory mechanisms (Oksanen 1991; Sinclair et al. 2003; Worm & Myers 2003) have been assumed to be responsible for the evolution of current life histories, morphologies and behaviours of pelagic organisms (Verity & Smetacek 1996; Lehman et al. 2004). However, direct evidence of the predominant type of control is rarely found in field studies, based on quantitative long-term estimates of trophic-level abundances (i.e. biomass; Aebischer et al. 1990; Micheli 1999; Richardson & Schoeman 2004). Our analysis (figure 2) indicates that the coastal phytoplankton–mesozooplankton system is mainly bottom-up regulated. The classic grazing food chain seems to be the prevailing pathway of primary production transfer to higher trophic levels at the coast. The outer sampling sites exhibit a more complex structure. Copepods, which constitute the most abundant group of mesozooplankton, are known to feed not only on phytoplankton cells, but also on other non-chlorophyll resources, such as the microzooplankton or marine snow. Microzooplankton have recently received much attention as important grazers, as well as for the role they play as an intermediate trophic level between primary producers and omnivorous mesozooplankton (Quevedo & Anadón 2001; Calbet & Landry 2004). This omnivorous feeding in the mesozooplankton is necessary to meet their metabolic demands, and it has been found experimentally that herbivorous feeding suffices only in upwelling-affected areas (Isla et al. 2004). We found no evidence of mesozooplankton exerting a significantly negative effect on the phytoplankton growth, so there is no support for the top-down control.

(d) The simulation of the modelled dynamics

The statistically deduced models were used to simulate the dynamics of the phyto- and zooplankton (figure 3). A good agreement between observations and simulations was obtained, especially for zooplankton whose models showed higher R2 (figure 3, left columns). Moreover, when simulating the average year, the models are able to capture the seasonal patterns typically observed for both phyto- and zooplankton (figure 3, right columns). The phytoplankton simulation exhibits the two classical periods of high growth rate corresponding to the transition between the two water-column states (spring and autumn). At the most oceanic station, the summer period of low growth imposed by stratification (i.e. oligotrophic conditions) broadens out to almost include the whole second half of the year. At the coast, this period is restricted to June. The simulation shows a clear zooplankton peak in spring at st. 2 and st. 3 following that of phytoplankton. At st. 1, the simulation shows a more extensive zooplankton peak, probably owing to more efficient fuelling by phytoplankton during most of the spring–summer. Altogether, these results suggest that the environmental factors included in our models jointly generate the observed dynamics in the studied marine system in an adequate way, a conclusion supported by the lack of remaining structure in the residuals (see fig. 4 of the electronic supplementary material). It is feasible that data collected at a higher temporal resolution (days) would have revealed a more variable interaction between phyto- and zooplankton, an interaction where biotic forces may have occupied a greater role. However, at a seasonal level, well captured by our sampling resolution, the observed dynamics are adequately explained by abiotic factors. As such, our deduced model structure offers a viable hypothesis worthy of testing against future data.

4. Conclusions

The non-parametric statistical modelling approach (based on mid-term time-series sampling) adopted in this paper is clearly able to capture the complexity of the biological dynamics. It is worth noting that the models differ along a coast–ocean gradient. Our results demonstrated that the seasonal dynamics of plankton in the southern Bay of Biscay are primarily driven by abiotic factors. The derived models show a clear evidence of bottom-up regulation. Our novel approach permits us to further investigate the dynamic of the phyto- and zooplankton communities, and to better understand which abiotic factors affect their dynamics along the ecological gradient. Our models open the possibility of predicting how the dynamics of plankton might have been affected if climate should change—predictions which will be highly valuable when trying to manage these coastal ecosystems in a sustainable way.


This work was initiated during visits of M.L. to the CEES. M.L. was supported by a grant from the Ministerio de Ciencia y Tecnología (FPI fellowship). The Asturian time-series sampling program was funded by the project ‘Control a largo plazo de las condiciones químico-biológicas de la Plataforma Continental de Asturias’ (SV-02-IEO y SV-97-IEO-1; Universidad de Oviedo, Instituto Español de Ocenografía) and Excelencia Investigadora del Principado de Asturias ‘Ecología de ecosistemas acuáticos’ (PR-01-GE-3) 2001–2004. Meteorological data were obtained from the Spanish Instituto Nacional de Meteorología. Four anonymous reviewers are thanked for their valuable input on an earlier version of this paper.


  • Present address: Programme for Plankton Biology, Department of Biology, University of Oslo, PO Box 1050 Blindern, 0316 Oslo, Norway.

  • Present address: Centre for Ecological and Evolutionary Synthesis (CEES), Department of Biology, University of Oslo, PO Box 1050 Blindern, 0316 Oslo, Norway.

  • The electronic supplementary material is available at or via

    • Received June 17, 2006.
    • Accepted June 23, 2006.


View Abstract