Abstract

Santos, M. B., González-Quirós, R., Riveiro, I., Cabanas, J. M., Porteiro, C., and Pierce, G. J. 2012. Cycles, trends, and residual variation in the Iberian sardine (Sardina pilchardus) recruitment series and their relationship with the environment. – ICES Journal of Marine Science, 69: 739–750.

Recruitment variability is an important component of the dynamics of Iberian sardine (Sardine pilchardus). Since 2006, poor recruitment has led to a decrease in stock biomass, the latest in a series of such crises for sardine fisheries. Understanding the mechanisms behind recruitment fluctuations has been the objective of many previous studies, and various relationships between recruitment and environmental variables have been proposed. However, such studies face several analytical challenges, including short time-series and autocorrelated data. A new analysis of empirical relationships with environmental series is presented, using statistical methods designed to cope with these issues, including dynamic factor analysis, generalized additive models, and mixed models. Relationships are identified between recruitment and global (number of sunspots), regional (NAOAutumn), and local [winter wind strength, sea surface temperature (SST), and upwelling] environmental variables. Separating these series into trend and noise components permitted further investigation of the nature of the relationships. Whereas the other three environmental variables were related to the trend in recruitment, SST was related to residual variation around the trend, providing stronger evidence for a causal link, possible mechanisms for which are discussed. After the removal of trend and cyclic components, residual variation in recruitment is also weakly related to the previous year's spawning-stock biomass.

Introduction

The Iberian sardine (Sardina pilchardus) exemplifies the life strategy of small pelagic fish around the world: short life, fast growth, high fecundity, and long spawning season (Carrera and Porteiro, 2003; Stratoudakis et al., 2007, and references therein). As is true for most marine fish (Cushing, 1996), recruitment of small pelagic species is highly variable and does not relate clearly to the abundance of the parental stock (measured as spawning-stock biomass, SSB), because of (among other factors) high and variable rates of mortality during the early life stages, which are thought to be strongly affected by environmental processes. In species with short lifespans, the small number of extant year classes in the population leads to a high dependence of SSB on recruitment. This, together with recruitment variability, presents great challenges to fishery managers because traditional approaches to management developed for fisheries on long-lived demersal species may be inappropriate. As part of the implementation of an ecosystem approach to fisheries management, European fishery scientists are required to take account of environmental forcing in stock assessment and forecasting (CE SEC, 2008), with the aim of achieving sustainable fisheries. Although sustainability has several dimensions (e.g. environmental, economic, and social), in terms of stock abundance, this means fisheries that would not collapse or suffer wide fluctuations in catches. For short-lived species, knowledge of recruitment dynamics and understanding of the mechanisms through which recruitment is controlled by the environment become even more important in ensuring sustainable and long-term exploitation (Rice, 2011).

Recognizing the need to incorporate environmental information into fisheries management is one thing; finding ways of doing so effectively is another. Marine ecosystem dynamics are notoriously complex, responding to physical processes at different spatial and temporal scales (Mann and Lazier, 1996). Several hypotheses have been proposed describing mechanisms by which the environment could be affecting recruitment or population variability in fish. These range from ideas based on the properties of the water column at short temporal scales (e.g. the stable ocean hypothesis; Lasker, 1975), mesoscale features (e.g. the optimal environmental window, OEW, hypothesis; Cury and Roy, 1989), and phenological processes (e.g. the match/mismatch hypothesis; Cushing, 1990) to theories based on ocean-scale, long-term climatic modes of variability (e.g. Chavez et al., 2003).

The Iberian sardine is considered to form a single stock distributed from the Strait of Gibraltar to the border between France and Spain (ICES, 2010; Figure 1). Although some exchange with other sardine populations outside the stock boundaries is likely, there is currently no evidence that this substantially influences stock dynamics (ICES, 2010). The Iberian stock has been assessed annually since 1978 by ICES, using data on Portuguese and Spanish fishery landings and the information obtained by fishery-independent acoustic and daily egg production method surveys.

Figure 1.

Map of the distribution area of the Iberian sardine stock. Limits of the ICES Subdivisions covered by the stock (VIIIc and IXa) are shown together with the main spawning locations (shaded areas). The dots indicate the position of the local environmental variables (SST, AT, wind strength, and the northerly and westerly components of wind direction) used in this study.

Figure 1.

Map of the distribution area of the Iberian sardine stock. Limits of the ICES Subdivisions covered by the stock (VIIIc and IXa) are shown together with the main spawning locations (shaded areas). The dots indicate the position of the local environmental variables (SST, AT, wind strength, and the northerly and westerly components of wind direction) used in this study.

The sardine fishery has been (and remains) economically important in Portugal, although it has decreased in importance in Galicia (NW Spain) over the years. In both countries, most of the catch is taken by purse-seiners, and landings show interannual and seasonal patterns reflecting the availability of the fish to the fleet. Periods of successive low recruitments in the past have led to “crises” in the fishery, most recently at the end of the 1990s (felt particularly in Galician waters), when the SSB reached a historical low (Wyatt and Porteiro, 2002; Carrera and Porteiro, 2003). Several studies have tried to identify possible mechanisms that could explain the variability of recruitment success of Iberian sardine and have explored the role of a series of environmental variables at global (e.g. sunspots as a proxy for solar irradiance; Guisande et al., 2004), regional (e.g. Atlantic Multidecadal Oscillation, AMO, and North Atlantic Oscillation, NAO; Guisande et al., 2001; Borges et al., 2003), and local scales [e.g. sea surface temperature (SST), upwelling, wind strength; Roy et al., 1995; Guisande et al., 2001; Santos et al., 2001; Borges et al., 2003].

Many of these studies have proposed that the environmental effects on recruitment success involve the modulation of the upwelling episodes taking place in the area. The Atlantic Iberian coast is characterized by a marked seasonality in oceanographic processes. There is a summer peak in the intensity of upwelling events responsible for the area's elevated primary production and high ecosystem productivity in general (Figueiras et al., 2002), and autumn sees the development of the Iberian Poleward Current (IPC), which flows north along the shelf break, is characterized by high temperature and salinity, and persists until spring (Pingree and Le Cann, 1990). That current generates frontal structures between coastal and offshore waters that affect plankton distribution and composition (Fernández et al., 1993) and may limit offshore advection of sardine larvae (González-Quirós et al., 2004; Santos et al., 2004). Sardine eggs, larvae, and juveniles could be vulnerable to high upwelling intensity that can cause offshore advection and increase mortality (Chesney and Alonso-Noval, 1989), but survival would also depend on the availability of suitable prey that will concentrate in waters enriched by mixing.

Sardine produce a large number of eggs (typically 50 000–60 000 eggs per female, or 300–400 eggs g−1; Muus and Nielsen, 1999; Zwolinski et al., 2001) over an extended spawning period. Spawning takes place from September to May, with two main peaks of activity, in winter (November–January; more important along the Atlantic coast; et al., 1990) and spring (March–May; more important in the Cantabrian Sea; Solá et al., 1992), and recruitment in a particular year integrates all this spatio-temporal variability in reproductive output (Cabanas et al., 2007). This life-history strategy offers obvious benefits for sardine survival by ensuring that at least part of the new generation will encounter favourable conditions but at the same time makes it more difficult to choose the most suitable variables for the study of environmental forcing on recruitment. For example, over which periods and geographic scales should one aggregate the environmental variables to arrive at a suitable index that captures the conditions in which good recruitment takes place?

The nature of the data available presents several statistical challenges, not always successfully addressed in previous studies, e.g. the short time-series, autocorrelated data, collinearity between putative explanatory variables, and the existence of non-linear relationships. We argue that this also justifies a re-evaluation of previously proposed environmental relationships in the Iberian sardine stock.

The aim of the present study was therefore to re-examine the pool of global, regional, and local environmental variables previously proposed as influencing sardine recruitment, using statistical methods (e.g. dynamic factor analysis, DFA; generalized linear models, GLMs; generalized additive models, GAMs; mixed modelling) designed to cope with the challenges presented by such data. Using these techniques, we make a new selection of relevant environmental variables and test their relationships with sardine recruitment. To deal with the issue that time-series showing simple trends may inevitably show empirical relationships with one another, each other, whether causal or coincidental, we decompose the recruitment and environmental time-series into their trend and residual components, with the aim of identifying which components of the series are related to each other. Owing to the uncertainties in the timing of the potential impact of environmental variables on the early life stages of sardine, we consider both contemporaneously measured environmental variables and those with a 1-year time-lag in relation to our recruitment series.

Material and methods

The Iberian sardine stock has been assessed by ICES since 1978. Recruitment (R, i.e. abundance of fish aged 0) and SSB for the whole stock are estimated annually using an age-structured model (AMCI, http://www.ices.dk/datacentre/software.asp) and were available up to 2009 (ICES, 2010), providing a time-series of 32 years.

Several studies have explored possible relationships between Iberian sardine recruitment and environmental variables at global, regional, and local scales (see above), so we selected environmental variables identified in those earlier studies. We have, however, avoided the use of smoothed series, e.g. based on moving averages, because the separation of trend and noise components is part of our subsequent analysis. Table 1 lists the environmental variables used, their periodicity, source, and details of any calculations applied.

Table 1.

Selected environmental variables, listing the periodicity with which data were available and their source and reference to how they were calculated, for local indices.

Abbreviation Name Periodicity Source Observations 
Global 
 avspots Average number of sunspots Monthly and annual http://solarscience.msfc.nasa.gov/greenwch/spot_num.txt  
Regional 
 NAO North Atlantic Oscillation Annual http://www.cru.uea.ac.uk/cru/data/nao.htm Atmospheric sea level pressure anomaly: difference in the normalized sea level pressure between Iceland and the Azores 
 NAOWinter Winter NAO Quarterly http://www.cru.uea.ac.uk/cru/data/nao.htm Averaged over December–February 
 NAOSpring Spring NAO Quarterly http://www.cru.uea.ac.uk/cru/data/nao.htm Averaged over March–May 
 NAOSummer Summer NAO Quarterly http://www.cru.uea.ac.uk/cru/data/nao.htm Averaged over June–August 
 NAOAutumn Autumn NAO Quarterly http://www.cru.uea.ac.uk/cru/data/nao.htm Averaged over September–November 
 AMO Atlantic Multidecadal Oscillation Annual http://www.esrl.noaa.gov/psd/data/timeseries/AMO/ Index of long-term SST in the North Atlantic Ocean 
 EA East Atlantic pattern Annual http://www.cpc.noaa.gov/data/teledoc/ea.shtml North–south dipole of anomaly centres spanning the North Atlantic from east to west 
 Gulf Gulf index Annual http://web.pml.ac.uk/gulfstream/data.htm Index of the variability of the position of the Gulf Stream 
Local 
 Iw_43_11 Upwelling index April–September Spanish Metereological Institute Average Ekman transport (m3 s−1 km−1) from April to September calculated at 43°N 11°W (Lavín et al., 1991
 IPC Poleward current index October–December Spanish Metereological Institute Average Ekman transport from October to December of previous year (Lavín et al., 1991
 SST Sea surface temperature Annual + 6 monthly http://dss.ucar.edu/datasets/ds540.1/ Annual and October–March and April–September averages (°C) at four locations: 40°N 10°W, 42°N 10°W, 44°N 8°W, and 44°N 6°W 
 AT Air temperature Annual + 6 monthly http://dss.ucar.edu/datasets/ds540.1/ Annual and October–March and April–September averages (°C) at four locations: 40°N 10°W, 42°N 10°W, 44°N 8°W, and 44°N 6°W 
 U Wind westerly component Annual + 6 monthly http://dss.ucar.edu/datasets/ds540.1/ Annual and October–March and April–September averages (m s−1) at four locations: 40°N 10°W, 42°N 10°W, 44°N 8°W, and 44°N 6°W 
 V Wind northerly component Annual + 6 monthly http://dss.ucar.edu/datasets/ds540.1/ Annual and October–March and April–September averages (m s−1) at four locations: 40°N 10°W, 42°N 10°W, 44°N 8°W, and 44°N 6°W 
 W Wind strength (scalar windspeed) Annual + 6 monthly http://dss.ucar.edu/datasets/ds540.1/ Annual and October–March and April–September averages (m s−1) at four locations: 40°N 10°W, 42°N 10°W, 44°N 8°W, and 44°N 6°W 
Abbreviation Name Periodicity Source Observations 
Global 
 avspots Average number of sunspots Monthly and annual http://solarscience.msfc.nasa.gov/greenwch/spot_num.txt  
Regional 
 NAO North Atlantic Oscillation Annual http://www.cru.uea.ac.uk/cru/data/nao.htm Atmospheric sea level pressure anomaly: difference in the normalized sea level pressure between Iceland and the Azores 
 NAOWinter Winter NAO Quarterly http://www.cru.uea.ac.uk/cru/data/nao.htm Averaged over December–February 
 NAOSpring Spring NAO Quarterly http://www.cru.uea.ac.uk/cru/data/nao.htm Averaged over March–May 
 NAOSummer Summer NAO Quarterly http://www.cru.uea.ac.uk/cru/data/nao.htm Averaged over June–August 
 NAOAutumn Autumn NAO Quarterly http://www.cru.uea.ac.uk/cru/data/nao.htm Averaged over September–November 
 AMO Atlantic Multidecadal Oscillation Annual http://www.esrl.noaa.gov/psd/data/timeseries/AMO/ Index of long-term SST in the North Atlantic Ocean 
 EA East Atlantic pattern Annual http://www.cpc.noaa.gov/data/teledoc/ea.shtml North–south dipole of anomaly centres spanning the North Atlantic from east to west 
 Gulf Gulf index Annual http://web.pml.ac.uk/gulfstream/data.htm Index of the variability of the position of the Gulf Stream 
Local 
 Iw_43_11 Upwelling index April–September Spanish Metereological Institute Average Ekman transport (m3 s−1 km−1) from April to September calculated at 43°N 11°W (Lavín et al., 1991
 IPC Poleward current index October–December Spanish Metereological Institute Average Ekman transport from October to December of previous year (Lavín et al., 1991
 SST Sea surface temperature Annual + 6 monthly http://dss.ucar.edu/datasets/ds540.1/ Annual and October–March and April–September averages (°C) at four locations: 40°N 10°W, 42°N 10°W, 44°N 8°W, and 44°N 6°W 
 AT Air temperature Annual + 6 monthly http://dss.ucar.edu/datasets/ds540.1/ Annual and October–March and April–September averages (°C) at four locations: 40°N 10°W, 42°N 10°W, 44°N 8°W, and 44°N 6°W 
 U Wind westerly component Annual + 6 monthly http://dss.ucar.edu/datasets/ds540.1/ Annual and October–March and April–September averages (m s−1) at four locations: 40°N 10°W, 42°N 10°W, 44°N 8°W, and 44°N 6°W 
 V Wind northerly component Annual + 6 monthly http://dss.ucar.edu/datasets/ds540.1/ Annual and October–March and April–September averages (m s−1) at four locations: 40°N 10°W, 42°N 10°W, 44°N 8°W, and 44°N 6°W 
 W Wind strength (scalar windspeed) Annual + 6 monthly http://dss.ucar.edu/datasets/ds540.1/ Annual and October–March and April–September averages (m s−1) at four locations: 40°N 10°W, 42°N 10°W, 44°N 8°W, and 44°N 6°W 

Data exploration

All dataseries were first explored for outliers, normality, homogeneity of variance, collinearity, etc., following the protocol proposed by Zuur et al. (2010). Logarithmic transformation was applied to the recruitment R series to obtain a better approximation to a normal distribution.

Because of the large number of potential explanatory variables, we first carried out DFA, a dimension-reduction technique specifically designed for time-series (Zuur et al., 2003), to try to identify common trends in the series of explanatory variables. First, for some variables, we had several alternative series available. Second, DFA was used to determine the presence of common trends over time among the different environmental variables. Where two different variables show similar trends, one would normally be dropped from subsequent models. As an example of the first process, seasonal and annual dataseries for average number of sunspots were available. We fitted models with 1, 2, and 3 common trends. The lowest Akaike information criterion (AIC) was obtained with one common trend, confirming the similarity of the 12 monthly series. We therefore considered it justified to choose the annual value as representative of the variable. This approach allowed for direct comparison with the results of previous studies.

DFA was also used to extract a common trend from the pool of finally selected environmental variables series (n = 26) with time-lags of 0 and 1 year relative to the R series. GAMs were then used to explore if these trends were related to the R series (log-transformed). Gaussian GAMs with an identity link were used, with the maximum number of degrees of freedom restricted to 3 (k = 4), to avoid overfitting.

Modelling approach

Time-series can be thought of as constituted by a trend, a cycle, an autocorrelation (AC) component, and residual variation (the variation that remains in the series after trend, cycle, and the AC have been removed). In the present work, this residual variation could include, for example, environmental/fishery/stock effects and unexplained noise. Given the relatively short series of R and SSB available, although evidence for cyclical patterns in all series using spectral analysis was sought, we did not attempt to remove cyclical patterns from the series. For AC, we took the view that in the response variables, AC could be “explained” by AC in the explanatory variables. However, if AC persisted in model residuals, mixed models were used in an attempt to control for it, using a specified variance structure (normally AR1). Hence, the focus of the analysis was the trend and noise components, following the argument that common trends could be coincidental, but that where response and explanatory variables show similar variation around trends, this is more likely to imply a causal link.

Trends, cycles, and autocorrelation

The existence of linear and simple polynomial trends in R and the explanatory variables was explored using GLMs, with time (year) as the explanatory variable. Models were fitted in which the highest order term was a linear, quadratic, or cubic function of time. The best model was chosen based on the lowest AIC. Spectral analysis was carried out to determine the presence and length of cycles in the detrended series of response and explanatory variables. Visual comparison was used to identify any cycles of equal length.

Modelling R as a function of environmental variables

GAMs were used to explore the relationships between R (response variable), SSB, and the suite of environmental explanatory variables (NAO, AMO, sunspots, etc.) for the same and previous years (i.e. time-lags 0 and 1 year in relation to the R series). Separate models were run for the values of environmental variables for the same year and for the previous year to take account of possible environmental influences not only on the survival of eggs and larvae, but also on the survival of juveniles. The rationale behind the use of SSB for the previous year as a predictor of R is that the effect of the abundance of the parental stock on recruitment success is likely to be mediated by the intensity of spawning, which occurs mainly in winter of the previous year. However, for completeness, we also tested for an effect of SSB in the current year, but this was never significant and is not discussed further.

Gaussian GAMs with an identity link were used always, with the maximum number of degrees of freedom restricted to 3 (k = 4), to avoid overfitting. Forward selection was used and, at each step, non-significant variables were dropped and one additional variable added. Candidate models were compared according to AIC values, and the model with the lowest AIC was selected as the basis for the next round of models. If final models contained non-significant terms, the consequence of removing these was tested using an F-test; they were retained if they improved the model fit significantly. Once final models were obtained, residuals were checked for patterns and for AC. If AC was discovered in the model residuals, mixed models (Zuur et al., 2009) with an autoregressive (AR) variance structure of time-lag 1 and the same fixed components were used instead. For relatively short time-series, more complex error structures are difficult or impossible to fit, so AC at lags >1 may require an alternative approach (see below).

As running separate models for contemporaneous (lag 0) and previous year's (lag 1) environmental variables separately is a somewhat artificial division, we also derived composite models by amalgamating the two final models and using backward selection to prune the model back to the new “best” model.

Modelling relationships between R and the environment for decomposed time-series

GLMs of the variable in question vs. time were used to decompose R and environmental explanatory variables (of the latter, only those remaining in final GAMs were used) into series of trend and residual (noise) values. Where no trend was found, the original environmental series was used. GAMs were then used to relate decomposed R and environmental series, including one environmental variable per model. As each variable could be used in its raw form, or as trend or noise components, this gives nine possible models for each environmental variable (Table 2). Of these, eight combinations are new, because GAMs of raw R vs. raw environmental series would have been fitted previously. However, in practice, several of these combinations are redundant. For environmental variables that showed no trend component and were therefore not decomposed, there would be just three models (two new ones). In addition, for environmental variables with trend components, the trend vs. trend model is trivial (it can always explain 100% of deviance in the recruitment trend), and if the trend in the environmental series is linear, both possible trend vs. residual models will a priori show no relationships.

Table 2.

Possible combinations of the raw and decomposed series.

Recruitment series Environmental series Comment 
Raw Raw Fitted previously 
Raw Trend 
Raw Noise 
Trend Raw  
Trend Trend 1, 2 
Trend Noise 1, 3 
Noise Raw  
Noise Trend 1, 3 
Recruitment series Environmental series Comment 
Raw Raw Fitted previously 
Raw Trend 
Raw Noise 
Trend Raw  
Trend Trend 1, 2 
Trend Noise 1, 3 
Noise Raw  
Noise Trend 1, 3 

1, null if no trend in environmental series; 2, trivial (100% of deviance will be explained); 3, no relationship a priori if trends in both series are linear.

Alternative decomposition of the R series using AR moving average

Where AC persisted in model residuals and mixed modelling failed to remove it, we explored the alternative of first fitting an AR moving average (ARMA) model to the R series, extracting trend and noise components, then running further GAMs using these components. All analyses were carried out using R (2.9.1) and the statistical programmes Brodgar (Highland Statistics Ltd) and Minitab (Minitab Inc.).

Results

Environmental effects on recruitment

Results from the GAMs for R (log-transformed) vs. the common trend in the environmental variables extracted with the DFA showed significant environmental effects on the R series. The model with the single common trend on contemporaneous (lag 0) environmental variables explained 26.4% of the deviance in the R series (p < 0.01). The model with the trend for the previous year's (lag 1) environmental data explained 37.2% of the deviance (p < 0.01). The inclusion of SSB as an additional explanatory variable did not improve either model.

Trends in the dataseries

Results from the GLM models for R (log-transformed) vs. year show a strong, linear, decreasing trend over time (p < 0.001; Figure 2). For higher order fits (quadratic, cubic), none of the terms were significant and AIC values were higher. SSB did not show a significant linear or quadratic trend, but all three terms in the cubic relationship were significant (p < 0.01) and the fitted trend captures the partial recovery of the stock after the all-time low estimated value at the end of the 1990s (Figure 2).

Figure 2.

Fitted trends superimposed in the (a) recruitment (log-transformed number of individuals) and (b) SSB (t) series for the Iberian sardine stock.

Figure 2.

Fitted trends superimposed in the (a) recruitment (log-transformed number of individuals) and (b) SSB (t) series for the Iberian sardine stock.

Among the environmental explanatory variables, no significant trends were found for upwelling or the NAOWinter or NAOSpring series. Most of the remaining environmental variables showed significant linear, quadratic, or cubic trends. The best fitting trends are summarized in Table 3.

Table 3.

Characteristics of the stock and environmental variables used in the analysis, showing the presence and type of linear and simple polynomial trends (arrows indicating the main direction of the trend), the presence and periodicity of cycles (years), and the presence and time-lag of partial autocorrelation (PAC).

Variable Cycle Trend PAC 
Ra Linear ↓ 
SSB 10 Cubic ↘ 
Average number of sunspots 11 Linear ↓ 1, 2, 3 
NAO 4, 8 Linear ↓ No AC 
NAOWinter 2–3, 8 None No AC 
NAOSpring None 
NAOSummer 4, 6 Linear ↓ No AC 
NAOAutumn 4-5 Quadratic ↘ No AC 
AMO Cubic ↗ 1, 9 
EA Linear ↑ 1, 2 
Upwelling Unclear None 1, 2, 3 
IPC Linear ↓ No AC 
Wind strength at 40°N10°Wb Linear ↑ 1, 2, 3 
SST at 40°N 10°Wb Unclear None 
AT at 40°N 10°Wb Unclear Quadratic ∩ 
Variable Cycle Trend PAC 
Ra Linear ↓ 
SSB 10 Cubic ↘ 
Average number of sunspots 11 Linear ↓ 1, 2, 3 
NAO 4, 8 Linear ↓ No AC 
NAOWinter 2–3, 8 None No AC 
NAOSpring None 
NAOSummer 4, 6 Linear ↓ No AC 
NAOAutumn 4-5 Quadratic ↘ No AC 
AMO Cubic ↗ 1, 9 
EA Linear ↑ 1, 2 
Upwelling Unclear None 1, 2, 3 
IPC Linear ↓ No AC 
Wind strength at 40°N10°Wb Linear ↑ 1, 2, 3 
SST at 40°N 10°Wb Unclear None 
AT at 40°N 10°Wb Unclear Quadratic ∩ 

aLog-transformed.

bWinter (October–March) averages.

Cycles in the dataseries

The R series shows an apparent main cycle with a periodicity of 4 years and secondary cycles with 2- and 10-year periodicities. For SSB, only a 10-year cycle is apparent from the spectral analysis (Figure 3). As the length of the time-series is just 32 years, identification of decadal cycles should be treated with caution.

Figure 3.

Cycles in detrended (a) recruitment (log-transformed number of individuals) and (b) SSB (t) series for the Iberian sardine stock.

Figure 3.

Cycles in detrended (a) recruitment (log-transformed number of individuals) and (b) SSB (t) series for the Iberian sardine stock.

Of the environmental explanatory variables investigated (Table 3), four of the NAO variables (annual, spring, summer, and autumn) and EA showed 4-year cycles, as seen for recruitment. Annual NAO also showed an 8-year cycle, and NAOSummer displayed an additional 6-year cycle. NAOWinter showed cycles with periodicities of 2–3 and 8 years. The average number of sunspots showed an 11-year cycle and AMO a 9-year cycle. No clear cycles were apparent in the upwelling index or SST and air temperature (AT) series.

Autocorrelation in the dataseries

Significant AC was present in both response series and the majority of the explanatory series. For R and SSB, partial AC (PAC) was significant only for a 1-year lag. For the environmental series, PAC was significant for both 1-year and some longer time-lags (Table 3). No significant AC was present in the series of NAO, three of the seasonal NAO components (excluding NAOSpring), or the index of the poleward current. The PAC for NAOSpring was significant for a 6-year time-lag.

GAMs for recruitment series

As explained above, based on the assumption that R strength should be a consequence of parent stock abundance and environmental conditions, GAMs were fitted to the R series using the SSB of the previous year (lag 1) and the environmental series for the same and previous years (lags 0 and 1) as explanatory variables. The previous year's SSB had no significant effect on recruitment (log-transformed), confirming the lack of a stock–recruitment relationship.

When considering the environmental variables at lag 0, the final model for R (log-transformed) included two local variables (winter wind strength and winter SST at 40°N 10°W), both of which showed negative linear effects, and the average number of sunspots, which showed a positive linear effect (Figure 4a). In models containing a single explanatory variable, these variables explained 20.2, 35.3, and 25.2% of deviance, respectively. Together in the same model, they explained 52.4% of deviance, with SST being the most important (p < 0.01). Significant AC was present in the residuals at lag 4. When using mixed models (GAMMs) in an attempt to remove AC, it was not possible to fit an ARMA (4, 0) variance structure to such a short series, and (as expected) a simple AR1, i.e. ARMA (1, 0), variance structure reduced, but did not remove, the AC in the model residuals.

Figure 4.

Final GAMs of recruitment (log-transformed number of individuals) in relation to environmental variables (lag 0): (a) smoothers for significant effects of wind strength (W40_10, in m s−1) and SST (SST40_10, °C) measured in winter at location 40°N 10°W and average annual number of sunspots (avspots), and (b) autocorrelograms (AC and PAC) showing that significant AC (time-lag 4) is present in the residuals.

Figure 4.

Final GAMs of recruitment (log-transformed number of individuals) in relation to environmental variables (lag 0): (a) smoothers for significant effects of wind strength (W40_10, in m s−1) and SST (SST40_10, °C) measured in winter at location 40°N 10°W and average annual number of sunspots (avspots), and (b) autocorrelograms (AC and PAC) showing that significant AC (time-lag 4) is present in the residuals.

The final model for R (log-transformed) vs. environmental variables for the previous year (lag 1) included positive linear effects of NAOAutumn and the average number of sunspots (positive linear effect; Figure 5a). Inclusion of the upwelling index improved the overall fit (as demonstrated by results of an F-test to compare models with and without this term), although its individual effect was not significant. In models containing a single explanatory variable, these variables explained 25.4, 21.5, and 26.7% of deviance, respectively; together in the same model, they explained 53.2% of deviance, with sunspots being the most important (p < 0.01). Model residuals showed no AC, and a mixed model with an AR1 variance structure had a higher AIC than the model with no variance structure.

Figure 5.

Final GAMs of recruitment (log-transformed number of individuals) in relation to environmental variables (lag 1): (a) smoothers for significant effects of average annual number of sunspots (avspots), NAOAutumn (averaged over September–November), and the upwelling index (lw_43_11, in m3s−1 km−1) measured at location 43°N 11°W, and (b) autocorrelograms (AC and PAC) showing no significant AC in the residuals.

Figure 5.

Final GAMs of recruitment (log-transformed number of individuals) in relation to environmental variables (lag 1): (a) smoothers for significant effects of average annual number of sunspots (avspots), NAOAutumn (averaged over September–November), and the upwelling index (lw_43_11, in m3s−1 km−1) measured at location 43°N 11°W, and (b) autocorrelograms (AC and PAC) showing no significant AC in the residuals.

Finally, a third approach was to use, in the same model, all the environmental variables at lags 0 and 1 year from the previous two final models and to generate a new best model by backward selection. The best model included positive linear effects of NAOAutumn (lag 1), the average number of sunspots, and winter SST at 40°N 10°W (lag 0). The upwelling index (lag 0) also had an effect, with the smoother indicating lower recruitment success at intermediate levels.

GAMs for decomposed time-series

As described above, GLMs vs. time were used to decompose R and explanatory variables (of the latter, only those remaining in final GAMs were used) into series of trend and residual values. Where no trend was found, the original environmental series was used. GAMs were then fitted for R vs. individual environmental variables for all meaningful combinations of raw, trend, and residual (noise) series (see Table 2 for an explanation).

Results of the GAMs run with the decomposed variables for R vs. environmental variables lagged at 0 and 1 year indicate that in the first case, the components that have significant effects on R were the linear trends of winter wind strength at 40°N 10°W, and the average number of sunspots (deviance explained 39.2%, p < 0.001 in both cases). The variation around the trend (residuals or noise) in R is not significantly related to either trend or residual components of winter wind strength at 40°N 10°W, or the average number of sunspots, but it is significantly related to winter SST at 40°N 10°W (deviance explained 16.3%, p < 0.05).

Similar results were obtained for the effects of the different components of the environmental variable time-series lagged by 1 year. Again, it is the linear trend of the average number of sunspots that is significantly related to R (deviance explained 39.2%, p < 0.001). The same applies to the effect of NAOAutumn on R, with the trend (this time a quadratic one) being the only component significantly related to R (deviance explained 36.5%, p < 0.001). The variation around the trend in R was not significantly related to either trend or residual components of the average number of sunspots or NAOAutumn.

Alternative decomposition using ARMA

The final model of R vs. environmental variables at lag 0 had AC at lag 4 in the residuals, which could not be removed using a mixed model. This suggests that the model had not captured the 4-year cycle previously observed in R. We therefore attempted an alternative decomposition of the R series, fitting an ARMA (4, 0) model to the series, and storing trends and residuals. Plots of the residuals and fitted values are shown in Figure 6. This process successfully removed AC. The “trend” component here captures both the previously identified trend and the 4-year cycle. We then tested these components for relationships with (i) the previously selected “best” explanatory variables, and (ii) the other explanatory variables available.

Figure 6.

ARMA decomposition of the recruitment series (log-transformed) into (a) trend, and (b) noise (residuals) components by fitting an ARMA (4,0) model. The autocorrelograms for (c) AC and (d) PAC show that the residuals of this model are free of AC.

Figure 6.

ARMA decomposition of the recruitment series (log-transformed) into (a) trend, and (b) noise (residuals) components by fitting an ARMA (4,0) model. The autocorrelograms for (c) AC and (d) PAC show that the residuals of this model are free of AC.

Running GAMs on the fitted values using the previously significant environmental variables (winter wind strength and SST, and the average number of sunspots), only wind strength showed a significant effect (deviance 34.2%, p < 0.001). None of these three variables were related to the residuals series. However, examining the other explanatory variables available, SSB of the previous year (lag 1) was significantly related to residual R (deviance 14.3%, p < 0.05).

Discussion

Understanding recruitment dynamics and the mechanisms responsible for its year-to-year variation has always been a target for fisheries management. For small pelagic fish such as sardine, devising effective management measures becomes more difficult if a clear relationship between stock and recruitment cannot be established, and reference limits such as Blim (defined as the SSB below which recruitment is impaired or the stock dynamics unknown) have to be approximated instead based on historical stock data, e.g. setting Blim = Bloss, the biomass corresponding to the lowest observed SSB (ICES, 1997). However, the implicit assumption of such an approach is that there is an underlying stock–recruitment relationship that is being obscured by other factors.

Since ICES began assessing the Iberian sardine stock in 1978, both SSB and R have fluctuated widely, and the past decade has seen both a relatively high peak (in 2006) and record lows. In 2000, SSB fell below the value estimated in 1978 for the first time, and in 2011, it reached an all-time low (ICES, 2011). It has been argued that Iberian sardine is particularly well adapted to thrive in a wider environmental window than other small pelagic fish species (e.g. anchovy; Borja et al., 2008), with a broad temperature tolerance recorded (both in terms of habitat and spawning distribution; Bernal et al., 2007), and a large number of eggs produced over an extended period (Stratoudakis et al., 2007). If environmental effects on early life stages (eggs, larvae, and juveniles) are the main drivers behind recruitment success, low population sizes could still produce strong recruitments. In 2000, the lowest sardine SSB on record (since 1978) produced just such a recruitment (ICES, 2010). However, this does not imply that the stock is in a healthy state. In the past decade, not only has recruitment continued to follow a downward trend, it has been more variable than recorded previously. The effect of pulses of good recruitment (another good one followed in 2004) on adult stock abundance appears to be of shorter duration now than was the case in the 1980s and (to a certain extent) in the 1990s.

Improving understanding of the mechanisms underlying sardine recruitment success (or lack thereof) has been the objective of several previous studies. Proposed environmental drivers include various global- to local-scale variables, integrated over the periods identified as the most critical to ensure the survival of eggs and larvae, by reducing the transport of eggs and larvae offshore. Indirect effects, e.g. on growth and condition through variations in food supply or water temperature, have been given less attention. The results from these studies show that environmental effects, although present, are often weak and in some cases contradictory. For example, upwelling intensity affects recruitment both positively and negatively (Dickson et al., 1988; Santos et al., 2007). Unfortunately, it is in the nature of short time-series that correlations tend to vary according to the length of the series available, and great care is needed to exclude coincidental relationships (Myers, 1998; Solow, 2002). Inconsistent results could also be attributable in part to the different analytical techniques used, e.g. whether non-linear relationships or time-lagged effects have been considered or whether AC has been controlled for. In the present study, there is no avoiding the fact that one is dealing with short time-series (32 years), where uncertainty in the recruitment estimates (arising from analytical assessments based on declared catches and hence subject to error; Myers, 1998) is potentially an issue.

Here, we have been less concerned with choosing the most suitable measure of particular environmental characteristics (e.g. our upwelling index could have been further subdivided into separate monthly values, or we could have tried to obtain local variables at different geographical scales or at time-lags of several months instead of 1 year), and we have concentrated instead on finding ways, within the empirical statistical modelling framework, to reveal more about the nature of the relationships, so providing some insight into mechanisms by which environmental variables could affect recruitment variability. Our final GAMs explained ∼50% of deviance in recruitment, with individual environmental variables explaining between 20 and 35% of deviance. Although these values are relatively high, the time-series are short and these are goodness-of-fit measures, not measures of predictive power. We have essentially erected hypotheses rather than formally tested them.

Decomposing the time-series into their constituent components (trends and noise) is helpful in that it reveals whether relationships arise from similar underlying trends or because the environmental variables explain residual variation around the trend. The former type of relationship carries more risk of being spurious. For example, in our analysis, the number of sunspots appeared to be one of the most important predictors of recruitment variability at both time-lag 0 and time-lag 1. However, R and sunspots share a linear trend, and there was no relationship between residuals around the respective trends (unsurprising given that sunspots show an 11-year cycle and recruitment includes a cyclic component of 4-year periodicity). The same is true of effects of wind strength (which show a 3-year cycle) and NAOAutumn (which has a cycle of 4–5 years). Nevertheless, in all three cases plausible mechanisms exist to link sardine recruitment to the environmental variable, as described previously by Borges et al. (2003) and Guisande et al. (2004), who also found significant relationships between wind indices, winter NAO (January–March), and sunspot cycle length and the averaged number of sunspots and sardine landings and recruitment.

The number of sunspots could be indirectly related to biological recruitment through a series of steps: more sunspots indicate higher solar irradiance (Wilson and Hudson, 1991), increasing the temperature of the stratosphere, affecting global wind patterns, and in turn modifying wind patterns at regional scales. The NAO controls the strength and direction of westerly winds and storm tracks across the North Atlantic (Wallace and Gutzler, 1981). Southerly and southwesterly winds dominate off the NW Iberian coast from October to March (Relvas et al., 2007), helping to develop the IPC, which is related in turn to the generation of convergence zones over the shelf (Frouin et al., 1990) that can act as retention areas for eggs, larvae, and their food (González-Quirós et al., 2004; Santos et al., 2004). Increased winter wind intensity could prolong and strengthen this current or could intensify the number and intensity of upwelling events depending on the wind component (e.g. northerly winds favour upwelling and are the dominant winds during summer in the area). Strong upwelling transports eggs and larvae offshore, but is simultaneously responsible for the elevated primary productivity in the area (Figueiras et al., 2002), potentially leading both to increased mortality and better growth among the survivors.

Compared with sunspots, NAO, and wind strength, SST is arguably a better candidate for a causal role in determining recruitment success, because it is the only variable in the final model that explained residual variation around the trend in recruitment (and individually, it explained a higher percentage of deviance). Winter SST, i.e. from October of the previous year to March of the current year, had a significant negative, linear effect on recruitment. However, in a model based on the previous year's environmental conditions (time-lag 1 year), there was no effect of (the previous year's) winter SST. These results are consistent with some of the hypothesized mechanisms for the effects of SST, but appear to contradict others. First, temperature has a direct effect on many physiological processes that may affect larval survival or adult reproduction. Higher winter temperatures may lead to reduced spawning success: winter SST refers to a 6-month period that includes one of the two important spawning peaks in the study area (Stratoudakis et al., 2007), i.e. autumn–winter. Spawning activity in sardine has previously been linked to SST, with an optimum at 14–15°C and avoidance below 12 or above 16°C (Stratoudakis et al., 2007). High summer temperatures persisting into the autumn–winter spawning period may therefore have a negative effect on total egg production of the stock. However, this hypothesis does not address the dynamics of early life stages, which are in theory responsible for much or even most of the variability of recruitment (Bailey and Houde, 1989). Therefore, higher SST could lead to increased larval growth rate, better survival, and consequently stronger recruitment (the opposite of the trend observed).

Such hypotheses about direct mechanisms in any case are probably overly simplistic. SST is an emergent property of a system's oceanography, and as such reflects the dynamic balance between a wide spectrum of different forcing factors (Cole and McGlade, 1998). On the West Iberian shelf from October to March, SST will be determined largely by the interplay between the upwelling of cold water, solar warming vs. atmospheric cooling of the surface layers, and the intensity of the IPC. During the autumn–winter period, characterized by low upwelling intensity and frequency over the West Iberian shelf (Fraga, 1981), upwelling events would be on the positive side of the upwelling–recruitment quadratic function defined by the OEW hypothesis (Cury and Roy, 1989), so favouring sardine recruitment. The IPC is characterized by warmer temperatures than the surrounding waters, and some authors have argued that it may favour larval retention over the shelf (González-Quirós et al., 2004; Santos et al., 2004). Therefore, temperature, as a proxy of the presence and intensity of the IPC, should be positively correlated with recruitment, again the opposite of the pattern inferred from the model.

The role of the IPC in determining food availability for early life stages of sardine is difficult to establish because, although the subtropical origin of this water mass could be associated with the prevalence of microbial assemblages (Fernández et al., 1991), its role in the development of spring blooms and frontal structures enhances heterogeneity in plankton composition, biomass, and distribution (Bode et al., 2002). Analysis of the role of these processes in determining SST, and the use of data at a higher temporal resolution (instead of a 6-month average), could help to disentangle the possible causal effect of winter SST on recruitment.

Further decomposition of the recruitment series using an ARMA (4, 0) model allowed extraction of a trend component that also captured the 4-year cycle, and with it all significant AC. It is interesting that the SSB of the previous year was significantly (but weakly) related to this residual series, indicating that indeed there may be a stock–recruitment relationship of sorts for Iberian sardine, though one that does not explain the main trend and cyclic components of the series.

As with previous studies, one can hypothesize plausible mechanisms to explain the observed empirical links, and there is no doubt that further research on the mechanisms by which the environment affects recruitment remains crucial to improving the understanding of the processes. However, such studies, if successful, may not provide a panacea: the environment in which fish live is intrinsically multidimensional, and it may simply be wrong to assume that a single process is the most important (Carrera and Porteiro, 2003).

In short-lived species, annual recruitment can represent a substantial component of the stock and, as such, it is likely that environmental signals are detectable not only in stock abundance trends, but also in sardine fishery landings, for which much longer time-series are available. In principle, landings will be a function of fishing effort and stock abundance. Therefore, we can also expect to be able to detect environmental signals in landings data (see Borges et al., 2003), especially if it is possible to control for variation in fishing effort. A recent investigation of environmental relationships for Mediterranean sardine and anchovy landings (Katara et al., 2011) indicates that such an approach may be worthwhile. The availability of longer series offers the prospect of both reducing the likelihood of finding coincidental relationships and allowing division of the series into model-construction and model-testing datasets.

Acknowledgements

We thank Portuguese and Spanish colleagues working on sardine, the crew and scientists in the acoustic surveys, those who collected landings data, and Alain Zuur (Highland Statistics) for statistical advice. MBS was supported by the Xunta de Galicia, Programa de Recursos Humanos (2010-104) and by the Spanish Ministry of Education, Programa Nacional de Movilidad de Recursos Humanos de Investigación (PR-2010-0518). These analyses were carried out as part of the LOTOFPEL project (Plan Nacional de I + D + I, CTM 2010-16053), and preliminary analysis was supported by the ANIMATE project MEXC-CT-2006-042337. Local environmental variables from this study are from the Research Data Archive (RDA) maintained by the Computational and Information Systems Laboratory (CISL) at the National Center for Atmospheric Research (NCAR); NCAR is sponsored by the National Science Foundation (NSF). Other sources of environmental data are acknowledged in Table 1. We thank the editors and two anonymous referees for their helpful comments on the manuscript.

References

Bailey
K.
Houde
E.
Predation on eggs and larvae of marine fishes and the recruitment problem
Advances in Marine Biology
 , 
1989
, vol. 
25
 (pg. 
1
-
83
)
Bernal
M.
Stratoudakis
Y.
Coombs
S. H.
Angelico
M. M.
Lago de Lanzos
A.
Porteiro
C.
Sagarminaga
Y.
, et al.  . 
Sardine spawning off the European Atlantic coast: characterization of and spatio-temporal variability in spawning habitat
Progress in Oceanography
 , 
2007
, vol. 
74
 (pg. 
210
-
227
)
Bode
A.
Varela
M.
Casas
B.
González
N.
Intrusions of eastern North Atlantic central waters and phytoplankton in the north and northwestern Iberian shelf during spring
Journal of Marine Systems
 , 
2002
, vol. 
36
 (pg. 
197
-
218
)
Borges
M. F.
Santos
A. M. P.
Crato
N.
Mendes
H.
Mota
B.
Sardine regime shifts off Portugal: a time-series analysis of catches and wind conditions
Scientia Marina
 , 
2003
, vol. 
67
 (pg. 
235
-
244
)
Borja
A.
Fontán
A.
Sáenz
J.
Valencia
V.
Climate, oceanography, and recruitment: the case of the Bay of Biscay anchovy (Engraulis encrasicolus)
Fisheries Oceanography
 , 
2008
, vol. 
17
 (pg. 
477
-
493
)
Cabanas
J. M.
Porteiro
C.
Carrera
P.
The effect of environmental changes in the NE Atlantic sardine fishery
 , 
2007
ICES Document CM 2007/ACFM: 25 (Annex 6)
pg. 
10 pp
 
Carrera
P.
Porteiro
C.
Stock dynamic of the Iberian sardine (Sardina pilchardus, W.) and its implication on the fishery off Galicia (NW Spain)
Scientia Marina
 , 
2003
, vol. 
67(Suppl. 1)
 (pg. 
245
-
258
)
CE SEC
Communication from the Commission to the Council and the European Parliament: the role of the CFP in implementing an ecosystem approach to marine management
 , 
2008
COM 2008/ 187
pg. 
10 pp
 
Chavez
F. P.
Ryan
J.
Lluch-Cota
S. E.
Ñiquen
M. C.
Climate: from anchovies to sardines and back: multidecadal change in the Pacific Ocean
Science
 , 
2003
, vol. 
299
 (pg. 
217
-
219
)
Chesney
E. J.
Alonso-Noval
M.
Coastal upwelling and the early life history of sardines (Sardina pilchardus) along the Galician coast of Spain
Rapports et Procès-Verbaux des Réunions du Conseil International pour l'Exploration de la Mer
 , 
1989
, vol. 
191
 (pg. 
63
-
69
)
Cole
J.
McGlade
J.
Clupeoid population variability, the environment and satellite imagery in coastal upwelling systems
Reviews in Fish Biology and Fisheries
 , 
1998
, vol. 
8
 (pg. 
445
-
471
)
Cury
P.
Roy
C.
Optimal environmental window and pelagic fish recruitment success in upwelling areas
Canadian Journal of Fisheries and Aquatic Sciences
 , 
1989
, vol. 
46
 (pg. 
670
-
680
)
Cushing
D. H.
Plankton production and year-class strength in fish populations: an update of the match/mismatch hypothesis
Advances in Marine Biology
 , 
1990
, vol. 
26
 (pg. 
249
-
293
)
Cushing
D. H.
Towards a Science of Recruitment in Fish Populations
 , 
1996
Oldendorf/Luhe, Germany
Ecology Institute
pg. 
175 pp
 
Dickson
R. R.
Kelly
P. M.
Colebrook
J. M.
Wooster
W. S.
Cushing
D. H.
North winds and production in the eastern North Atlantic
Journal of Plankton Research
 , 
1988
, vol. 
10
 (pg. 
151
-
169
)
Fernández
E.
Bode
A.
Botas
A.
Anadón
R.
Microplankton assemblages associated with saline fronts during a spring bloom in the central Cantabrian Sea: differences in trophic structure between water bodies
Journal of Plankton Research
 , 
1991
, vol. 
13
 (pg. 
1239
-
1256
)
Fernández
E.
Cabal
J.
Acuña
J. L.
Bode
A.
Botas
A.
García-Soto
C.
Plankton distribution across a slope current-induced front in the southern Bay of Biscay
Journal of Plankton Research
 , 
1993
, vol. 
15
 (pg. 
619
-
641
)
Figueiras
F. G.
Labarta
U.
Fernandez Reiviz
M. J.
Coastal upwelling, primary production and mussel growth in the Rias Baixas of Galicia
Hydrobiologia
 , 
2002
, vol. 
484
 (pg. 
121
-
131
)
Fraga
F.
Richards
F. A.
Upwelling off the Galician coast
Coastal Upwelling
 , 
1981
Washington, DC
American Geophysical Union
(pg. 
176
-
182
529 pp
Frouin
R.
Fiúza
A.
Ambar
I.
Boyd
T.
Observations of a poleward surface current off the coasts of Portugal and Spain during winter
Journal of Geophysical Research
 , 
1990
, vol. 
95
 (pg. 
679
-
691
)
González-Quirós
R.
Pascual
A.
Gomis
D.
Anadón
R.
Influence of mesoscale physical forcing on trophic pathways and fish larvae retention in the central Cantabrian Sea
Fisheries Oceanography
 , 
2004
, vol. 
13
 (pg. 
351
-
364
)
Guisande
C.
Cabanas
J. M.
Vergara
A. R.
Riveiro
I.
Effect of climate on recruitment success of Atlantic Iberian sardine Sardina pilchardus
Marine Ecology Progress Series
 , 
2001
, vol. 
223
 (pg. 
243
-
250
)
Guisande
C.
Vergara
A. R.
Riveiro
I.
Cabanas
J. M.
Climate change and abundance of the Atlantic Iberian sardine (Sardina pilchardus)
Fisheries Oceanography
 , 
2004
, vol. 
13
 (pg. 
91
-
101
)
ICES
Report of the Study Group on the Precautionary Approach to Fisheries Management
 , 
1997
ICES Document CM 1997/Assess: 07
pg. 
37 pp
 
ICES
Report of the Working Group on Anchovy and Sardine
 , 
2010
ICES Document CM 2010/ACOM: 16
pg. 
295 pp
 
ICES
Report of the Working Group on Anchovy and Sardine
 , 
2011
ICES Document CM 2011/ACOM: 16
pg. 
462 pp
 
Katara
I.
Pierce
G. J.
Illian
J.
Scott
B. E.
Environmental drivers of the anchovy/sardine complex in the eastern Mediterranean
Hydrobiologia
 , 
2011
, vol. 
670
 (pg. 
49
-
65
)
Lasker
R.
Field criteria for survival of anchovy larvae: the relation between inshore chlorophyll maximum layers and successful first feeding
Fishery Bulletin US
 , 
1975
, vol. 
73
 (pg. 
453
-
462
)
Lavín
A.
Díaz del Río
G.
Casas
G.
Cabanas
J. M.
Afloramiento en el noroeste de la Península Ibérica. Índices de afloramiento para el punto 43°N, 11°W, Período 1966–1989
 , 
1991
Informes Técnicos del Instituto Español de Oceanografía, 91
pg. 
40 pp
 
Mann
K. H.
Lazier
J. R. N.
Dynamics of Marine Ecosystems: Biological–Physical Interactions in the Oceans
 , 
1996
2nd edn
Oxford
Blackwell Science
pg. 
xxii + 394 pp
 
Muus
B. J.
Nielsen
J. G.
Sea Fish
 , 
1999
Hedehusene, Denmark
Scandinavian Fishing Year Book
pg. 
340 pp
 
Myers
R. A.
When do environment–recruitment correlations work?
Reviews in Fish Biology and Fisheries
 , 
1998
, vol. 
8
 (pg. 
285
-
305
)
Pingree
R. D.
Le Cann
B.
Structure, strength and seasonality of the slope currents in the Bay of Biscay region
Journal of the Marine Biological Association of the UK
 , 
1990
, vol. 
70
 (pg. 
857
-
885
)
P.
Cabral e Silva
R.
Cunha
E.
Farinha
A.
Meneses
I.
Moita
T.
Sardine spawning off Portugal
Boletín del Instituto Nacional de Investigación de Pescas, Lisboa
 , 
1990
, vol. 
15
 (pg. 
31
-
34
)
Relvas
P.
Barton
E. D.
Dubert
J.
Oliveira
P. B.
Peliz
A.
da Silva
J. C. B.
Santos
A. M. P.
Physical oceanography of the western Iberia ecosystem: latest views and challenges
Progress in Oceanography
 , 
2007
, vol. 
74
 (pg. 
149
-
173
)
Rice
J.
Managing fisheries well: delivering the promises of an ecosystem approach
Fish and Fisheries
 , 
2011
, vol. 
12
 (pg. 
209
-
231
)
Roy
C.
Porteiro
C.
Cabanas
J.
The optimal environmental window hypothesis in the ICES Area: the example of the Iberian sardine
ICES Cooperative Research Report
 , 
1995
, vol. 
206
 (pg. 
57
-
65
)
Santos
A. M. P.
Borges
M. F.
Groom
S.
Sardine and horse mackerel recruitment and upwelling off Portugal
ICES Journal of Marine Science
 , 
2001
, vol. 
58
 (pg. 
589
-
596
)
Santos
A. M. P.
Chícharo
A.
Dos Santos
A.
Miota
T.
Oliveira
P. B.
Peliz
A.
P.
Physical–biological interactions in the life history of small pelagic fish in the western Iberia upwelling ecosystem
Progress in Oceanography
 , 
2007
, vol. 
74
 (pg. 
192
-
209
)
Santos
A.
Peliz
A.
P.
Dubert
J.
Oliveira
P.
Angélico
M.
Impact of a winter upwelling event on the distribution and transport of sardine eggs and larvae off western Iberia: a retention mechanism
Continental Shelf Research
 , 
2004
, vol. 
24
 (pg. 
149
-
165
)
Solá
A.
Franco
C.
Lago de Lanzós
A.
Motos
L.
Temporal evolution of Sardina pilchardus (Walb.) spawning on the N–NW coast of the Iberian Peninsula
Boletín del Instituto Español de Oceanografía
 , 
1992
, vol. 
8
 (pg. 
97
-
114
)
Solow
A. R.
Fisheries recruitment and the North Atlantic Oscillation
Fisheries Research
 , 
2002
, vol. 
52
 (pg. 
295
-
297
)
Stratoudakis
Y.
Coombs
S. H.
Lago de Lanzós
A.
Halliday
N.
Costas
G.
Caneco
B.
Franco
C.
, et al.  . 
Sardine (Sardina pilchardus) spawning seasonality in European waters of the Northeast Atlantic
Marine Biology
 , 
2007
, vol. 
152
 (pg. 
201
-
212
)
Wallace
J. M.
Gutzler
D. S.
Teleconnections in the geopotential height field during the northern hemisphere winter
Monthly Weather Review
 , 
1981
, vol. 
109
 (pg. 
784
-
812
)
Wilson
R. C.
Hudson
H. S.
The suńs luminosity over a complete solar cycle
Nature
 , 
1991
, vol. 
351
 (pg. 
42
-
44
)
Wyatt
T.
Porteiro
C.
Sherman
K.
Skjoldal
H. R.
Iberian sardine fisheries: trends and crises
Large Marine Ecosystems of the North Atlantic: Changing States and Sustainability, Large Marine Ecosystems Series, 10
 , 
2002
Oxford
Elsevier Science
(pg. 
321
-
338
464 pp
Zuur
A. F.
Fryer
R. J.
Jolliffe
I. T.
Dekker
R.
Beukema
J. J.
Estimating common trends in multivariate time-series using dynamic factor analysis
Environment
 , 
2003
, vol. 
14
 (pg. 
665
-
685
)
Zuur
A. F.
Ieno
E. N.
Elphick
S.
A protocol for data exploration to avoid common statistical problems
Methods in Ecology and Evolution
 , 
2010
, vol. 
1
 (pg. 
3
-
14
)
Zuur
A. F.
Ieno
E. N.
Walker
N. J.
Saveliev
A. A.
Smith
G. M.
Mixed Effects Models and Extensions in Ecology with R
 , 
2009
Berlin
Springer
pg. 
574 pp
 
Zwolinski
J.
Stratoudakis
Y.
Soares
E.
Intra-annual variation in the batch fecundity of sardine off Portugal
Journal of Fish Biology
 , 
2001
, vol. 
58
 (pg. 
1633
-
1645
)

Author notes

Handling editor: Bill Turrell