The use of GEREM estimates to analyze correlations between large-scale environmental factors and glass eel recruitment at the European scale

Get Complete Project Material File(s) Now! »

Recruitment time series data available at the European level Glass eel migration process

Various recruitment time series data exist throughout Europe, they are derived from both fishery-dependent sources and fishery-independent surveys (Dekker, 2002b; Dekker, 2002c). Ideally, the strength of the year-class of youngest larvae should be measured annually, but since a monitoring in the Sargasso Sea is not or hardly possible, the nearest surrogate corresponds to glass eel annually recruiting to continental waters. The earliest stage measurable coincides with the recruitment from the sea into a river (estuarine recruitment) contrary to the later migration from the estuary into the upstream river (fluvial recruitment).
The natural migration process provides several opportunities for the development of catching methods and derivation of recruitment time series. Once glass eels arrive on continental shelves, they subsequently migrate towards the coast, using a mechanism known as selective tidal transport: during flood tide, glass eel swim in the water column and passively drift towards the coast, while during ebb tide; they hide near the bottom to resist being washed back to the ocean. This phase of selective tidal transport is followed by a resting phase in the estuary: glass eels or elvers cannot progress further upstream until the temperature reaches a certain level (12-15°C) resulting in large concentrations in early spring, especially at the upstream tidal limit. Once temperature has risen further, active migration into the river begins, where eels swim against the river flow (Harrison et al., 2014).
As a consequence, a wide range of fishing gears (Tela, glass eel dipnet, Fykenets) is used by commercial fishermen to filter glass eels during the flood tide. Moreover, the resting phase causes a great concentration of glass eels in time and space allowing active fishing by commercial fisheries at the upstream tidal limit. In some places, the presence of weirs, sluices or ship locks at the marine freshwater interface results in concentration of glass eels providing another opportunity for fishing. Lastly, glass eels or young yellow eels can be collected in traps when they swim actively upstream and encounter obstacles, later in spring (Dekker, 2002a).
In this study, the different recruitment time series are considered to be indicators of the recruitment produced by the spawning process and the subsequent ocean migration. As such, we seek to avoid the effects of local conditions on recruitment levels. However, the survival of leptocephalus larvae is affected by ocean and atmospheric conditions at their spawning grounds and during their oceanic migration on both large and regional scales (Bonhommeau et al., 2007; Bonhommeau et al., 2009; Bonhommeau et al., 2008). On the other hand, local factors may play a role in determining the migration of glass eels to estuarine habitats after metamorphosis (Arribas et al., 2012). In this way, local conditions may also contribute to the year to year variability in glass eel recruitment to estuaries. To address this issue, we should prioritize the most outward monitoring sites but it remains difficult to distinguish global or regional trends from local ones whatever the monitoring site and method are. Hence, only a network of monitoring stations and an overall analysis are able to disentangle the effects of local conditions.

A Dynamic Factor Analysis (DFA) to investigate the different trends in recruitment over Europe

Dynamic Factor Analysis is a technique used to detect common trends in a set of time series. While other dimension-reduction techniques like PCA are usually applied to treat large data sets, these are not able to take account of time in any way. DFA is a dimension-reduction technique designed for time series data and it can be used to model short and nonstationary time series in terms of common trends, effects of explanatory variables and interactions between the response variables (Zuur et al., 2003b; Zuur et al., 2003a).

Mathematics underlying DFA

The DFA model is given by: 𝑌𝑡=𝐶+𝐴𝑍𝑡+𝛽𝑋𝑡+𝜀𝑡 𝜀𝑡~ 𝑁(0,𝜎2𝑉) (1).
Where 𝑌𝑡 is N-by-1 vector containing the values of the N time series at time t, 𝑍𝑡 is M-by-1 vector representing the values of the M common trends at time t, A is the N*M matrix containing factor loadings and determines the exact form of the linear combination of the common trends, 𝜀𝑡 is normally distributed noise (of dimension N-by-1) with expectation 0 and variance 𝜎2𝑉 where V is a covariance matrix. The intercept C is a vector of dimension N-by-1 containing an intercept for each time series. If one wants to include explanatory variables, 𝑋𝑡 would be a vector containing the values of the L explanatory variables at time t and β is a N*L matrix containing regression coefficients. Hence, potential effects of explanatory variables are modeled as in linear regression. If a loading is relatively large and positive, we know that the corresponding time series follows the pattern of the corresponding trend. If a loading is close to zero, we know it doesn’t follow this pattern. A loading that is relatively large and negative indicates that the time series follows the trend opposite pattern. By comparing factor loadings, it can be inferred which group of time-series are related to the same common trends. Moreover, each error component 𝜀𝑡 has a different variance and the error components of different time series can be allowed to covary through an unstructured positive definite error covariance matrix V. The trends, which represent the underlying common patterns over time, follow a random walk, which is mathematically defined by: 𝑧𝑡= 𝑧𝑡−1+ 𝛿𝑡 𝛿𝑡~𝑁(0,𝜎𝛿2) (2) 𝜎𝛿2 is a diagonal error covariance matrix.

Correlations analysis at different time scales

Both short-term and long-term correlations were tested at time lags +1, +2, and +3: recruitment estimates from GEREM 𝑌𝑡+𝜏 were correlated with environmental variables 𝑋𝑡 with 𝜏=1, 2 or 3. We also analyzed correlations between recruitment 𝑌𝑡 and 𝑋𝑡−𝜏∗̅̅̅̅̅̅̅̅̅̅̅̅ where 𝑋𝑡−𝜏∗̅̅̅̅̅̅̅̅̅̅̅̅ accounts for average conditions over the years t- 𝜏 to t. We choose to test correlations between recruitment estimates and environmental data up to a time lag of 3 years to account for processes potentially occurring during the spawning period within the Sargasso Sea and during the subsequent larval migration to the European coasts. Indeed, maximum migration duration is estimated to be not longer than 2-3 years (Bonhommeau et al., 2010). We didn’t take into account the period during which silver eels migrate back to the Sargasso Sea, assuming that this stage was not affected by large-scale oceanic factors due to its capacity to adapt to contrasted environments.

Interannual correlations

We first aimed to analyze correlations between inter-annual variations of recruitment and large scale ocean-atmospheric time series by means of the cross-correlation function (CCF). All environmental data were used in this first analysis except the SST for which we only wanted to test correlations based on its long term evolution. The CCF function (1) enables to investigate time-lagged relationships between environmental time series (NAO, GSI and PP in the Sargasso Sea) and recruitment estimates and should help to identify the nature of the relationship and how they are correlated in time. 𝜌𝑋𝑌(𝜏)=𝐸[(𝑋𝑡−𝜇𝑋)(𝑌𝑡+𝜏−𝜇𝑌)]/(𝜎𝑋𝜎𝑌) (1).
Here, the influential time series called the input time series is 𝑋𝑡, whereas the output time series is 𝑌𝑡. However, short-term correlations between the two time series can be masked by overlaying trends in both time-series, and time series models are required to separate autocorrelations within a single time-series from cross-correlations with another. Here, we removed the autocorrelations from each time series of data by fitting an ARIMA model to the data and by using the resulting residuals (procedure called prewhitening) (Appendix VIII). ARIMA models combine autoregressive (AR) and moving average (MA) models into a single statistical model with the option to difference the time series so that it aims to model the various degrees of autocorrelations within time-series. The modelling steps follow the Box-Jenkins methodology (Box and Jenkins, 1976). The CCF function is then applied to the stationary and prewhitened time-series.

Correlations and long-term fluctuations

Prewhitening the time series may hide coinciding long-term trends in the input and output time-series because removing the autocorrelation is equivalent to removing low-frequency variability from data. Therefore, significant correlations on larger time scales may not be visible after prewhitening. Because the environmental descriptors used here such as the SST, GSI and NAO may reflect slowly changing processes over large temporal scales, it would be useful to investigate this second type of correlation between recruitment index and environmental data. As recommended by Pyper and Peterman (1998), a MA model was used to smooth time series data by removing high frequency variation in order to analyze correlations resulting from low-frequency variability. However, smoothing can increase autocorrelations in the time series and Pyper and Peterman (1998) advised to modify the hypothesis testing procedure by computing a corrected degree of freedom for the sample correlation, following this equation: 1𝑁∗≈1𝑁+2𝑁Σ𝑟𝑥𝑥(𝑗)∗𝑟𝑦𝑦(𝑗)𝑗 (2).

The European eel stock situation through the model GEREM Robustness of the model

Gelman and Rubin diagnostics confirmed that the chain converged after 80000 iterations (40000 burn-in and 40000 samples for inference). 𝑅̅ Statistics were less than 1.03 for all parameters (Table 1).
The posterior distributions of some 𝜏𝐼𝐴𝑖,𝜏𝐼𝐸,𝜏𝑈𝑖,𝜏𝐼𝑃𝑖 (precision of observed time-series) were influenced by their respective prior distributions (Appendix IV). However, since the precision in a lognormal distribution is nearly equal to the inverse of the squared coefficient of variation, we considered that the precision should be greater than 1 (i.e. a coefficient of variation inferior to 1). Greater values leaded to unlikely dissymmetries with median and modes very distinct from the mean, especially when the mean is low (because of the Laurent correction). Moreover, the parameters log (𝑎𝑐) and log (g) (scaling parameters for trap series and total catches in the Somme estuary) were also influenced by their respective prior distributions in some cases (Appendix IV) but we thought that prior knowledge was reliable enough to keep them. The discrepancy between prior and posterior distributions for those parameters is probably due the model assumptions of similar trends and densities of glass eels in a same zone. On the other hand, the model fitted quite well available abundance indices (Appendix V). On 41 time series available, only 8 abundance indices were less well fitted by the model (Appendix V and VI). As such, the model with its data, prior distributions, and modeling assumptions, was able to fit 80% of abundance indices, which suggests that the model assumptions were reasonable. The model fitted less well the SeGEMAC, Tiber, Bresle, Somme, Vaccares, Katwijk, Inagh and Erne time series which received a negative R2 (a negative R2 is sign of a systematic bias). SeGEMAC, Bresle and Somme abundance indices were poorly fitted because they may reflect recruitment densities slightly higher or lower than the average zonal densities within their respective zones. For the other ones, they were less well fitted by the model because they showed a general trend a bit different from the other series of the same zone. Both Tiber time series displayed a sharper decline than other series from the Mediterranean region, which was not visible during the DFA analysis. This was due to differences in data transformations: data were log transformed by removing null data in the model whereas a log (x+1) transformation was performed in the DFA analysis. This only impacted the trend of the Tiber time series which contained very low values at the end of the series, and consequently smoothed the trend in the last years.
Apart from these particular cases, the model was able to fit most abundance indices. This suggests that the model is likely to be trustworthy and reveal average behaviors and processes that govern recruitment of glass eel. Some time series were less well fitted because they might reveal local specificities in terms of density or temporal trend.

Recruitment estimates at different spatial scales

The overall European glass eel recruitment has been estimated and, unsurprisingly shows a clear decline from 1980 onwards. The recruitment was already decreasing in 1960, showing however a rise in the late 1970s but dropping again in 1980 (Fig 11). It seems that the overall recruitment is increasing since 2010 after a minimum level reached in 2009. The recruitment estimates for 2015 seems to be however lower than 2014. These findings are similar to the recruitment analyses performed by ICES (2016, see annex 8).

Table of contents :

I. Introduction
1. The European eel life cycle
2. A declining population at the European scale: an alarming state
3. Which framework to evaluate the status of the European eel stock?
4. A lack in the existing quantitative tools to assess the European eel population
5. A population potentially impacted by large scale climatic factors
II. Materials and Methods
1. Recruitment time series data available at the European level
Glass eel migration process
Relative recruitment indices
Absolute recruitment indices
2. Presentation of the model GEREM
State-space model
Observation model
Prior information and expertise
Modelling assumptions
Bayesian inference
Sensitivity analysis
3. A Dynamic Factor Analysis (DFA) to investigate the different trends in recruitment over Europe Mathematics underlying DFA
Setting up the DFA model
4. Correlation analysis
Data used in the correlation analysis
Correlations analysis at different time scales
III. Results
1. The different trends in recruitment highlighted by the DFA
A widespread decline over the European continent
Is there any spatial pattern in eel recruitment?
2. The European eel stock situation through the model GEREM
Robustness of the model
Recruitment estimates at different spatial scales
Sensitivity analysis
3. Correlation analysis at the European scale
No correlations between recruitment indices and the GSI or the SST
Short term correlations between certain recruitment indices and PP in the Sargasso Sea
The NAO index affect the recruitment at a broad scale
IV. Discussion
1. Modeling methodologies and underlying assumptions
Recruitment modelling
Information included in the model
Recruitment time series used in the model
2. GEREM: a stock assessment model to assess a widely distributed population
Towards an improvement of the trend-based approach
GEREM: a starting point for improving stock-recruitment relationships?
3. Application of GEREM to the European eel population
4. The use of GEREM estimates to analyze correlations between large-scale environmental factors and glass eel recruitment at the European scale
Impact of large scale SST patterns and transport on eel survival
Primary production in the Sargasso Sea affects glass eel recruitment success
Glass eel influenced by the physical and biological structure of the North Atlantic
Are there different migration durations and migratory routes?
Recruitment estimates from GEREM
V. Conclusion
VI. Bibliography

GET THE COMPLETE PROJECT