# Monitoring small pelagic fish in the Bay of Biscay ecosystem, using indicators from an integrated survey

Mathieu Doray, Pierre Petitgas, Martin Huret, Erwan Duhamel, Jean Baptiste Romagnan, Matthieu Authier, Christine Dupuy, Jérôme Spitz

► To cite this version:

Mathieu Doray, Pierre Petitgas, Martin Huret, Erwan Duhamel, Jean Baptiste Romagnan, et al. Monitoring small pelagic fish in the Bay of Biscay ecosystem, using indicators from an integrated survey. Progress in Oceanography, 2018, 166, pp.168-188. 10.1016/j-pocean.2017.12.004. hal-03015724

HAL Id: hal-03015724

https://hal.science/hal-03015724v1

Submitted on 19 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Monitoring small pelagic fish in the Bay of Biscay ecosystem, using indicators from an integrated survey

Doray Mathieu<fnref n="1" />,*, Petitgas Pierre<fnref n="1" />, Huret Martin<fnref n="2" />, Duhamel Erwan<fnref n="3" />, Romagnan Jean-Baptiste<fnref n="1" />, Authier Matthieu<fnref n="4" />, Dupuy Christine<fnref n="5" />, Spitz Jérome<fnref n="4" />

## Abstract :

This paper is a practical contribution to two important debates raised by the implementation of marine ecosystem based management: (i) which ecosystem data can be provided by a fisheries survey optimised for ecosystem monitoring; and (ii) how to combine/select potential indicators to derive useful information on marine ecosystem status and dynamics? A suite of 143 potential indicators, including spatial indices, representing small pelagic fish and their biotic and abiotic environments are presented. Indicators were routinely derived from the PELGAS integrated ecosystem survey conducted in spring in the Bay of Biscay (BoB). The general patterns over time in this suite of 5–16 years, non-stationary time series are characterised using a methodology based on min-max autocorrelation factors (MAF), to select the most continuous indicators within, and across, several ecosystem components: hydrology, phytoplankton, mesozooplankton, small pelagic fish and megafauna. Potential interactions between selected indicators and external forcing variables, including climate and fishing, were assessed. The results confirm the importance of river discharges, bottom temperature, chlorophyll-a and mesozooplankton biomass in the dynamics of the BoB pelagic ecosystem. Small pelagic fish species appear to have followed distinct trajectories during the last 15 years. A marked decrease in anchovy and sardine mean weights at ages 1 and 2 over the last 15 years was highlighted; potentially caused by density-dependent competition. The quasi-absence of significant correlation between selected survey indicators and climate and fishing pressure proxies suggests a moderate exploitation rate of small pelagic fish resources, and confirms the so far limited effects of large-scale climate forcing on the BoB pelagic ecosystem. Perspectives for the assessment of marine ecosystem status on the basis of suites of indicators derived from integrated ecosystem surveys are discussed.

<footnote n="1">¹ Unité Écologie et Modèles pour l’Halieutique, Ifremer Nantes, Rue de l’Ille d’Yeu, BP 21105, 44300 Nantes Cedex 3, France</footnote>

<footnote n="2">² Laboratoire de Biologie Halieutique, Centre Ifremer Bretagne, ZI de la Pointe du Diable – CS 10070, 29280 Plouzané, France</footnote>

<footnote n="3">³ Laboratoire de Technologie et Biologie Halieutique, Station Ifremer de Lorient, 8, rue François Toullec, 56100 Lorient, France</footnote>

<footnote n="4">⁴ Observatoire PELAGIS, UMS 3462, Université de La Rochelle/CNRS, La Rochelle, France</footnote>

<footnote n="5">⁵ Université de la Rochelle-CNRS, UMR 7266, Littoral Environnement et Sociétés (LIENSs), 2 rue Olympe de Gouges, 17000 La Rochelle cedex, France</footnote>

<footnote>* Corresponding author : Mathieu Doray, email address : mathieu.doray@ifremer.fr</footnote>
**Highlights**

► A method for selecting relevant ecosystem indicators is applied to PELGAS integrated survey. ►
River plume, bottom temperature, phyto & mesozooplankton selected as hydrobiology indicators. ►
Small pelagic fish species appear to have followed distinct trajectories over the last 15 years. ► A
marked decrease in mean weights of age 1 & 2 anchovy and sardine over the last 15 years is
highlighted. ► Low impact of fishing and no effect of climate forcing on pelagic ecosystem in spring are
confirmed.

**Keywords** : Min/max autocorrelation factors, Time correlation, Survey-based indicators, Marine Strategy Framework Directive, Ecosystem variability, Bay of Biscay, Pelagic ecosystem
Introduction

An indicator-based approach has been endorsed for decades to assess and manage specific human impacts,
such as commercial fishing, on marine ecosystems, (Garcia, 2000; Rice and Rivard, 2007). This sectoral
approach to marine-environment management has shown limitations to address human impacts on other parts
of the system (Kupschus et al., 2016). A more holistic, ecosystem approach to marine-environment
management has been proposed to overcome these limitations (Link et al., 2002, p. 20; Livingston et al.,
2005; Kupschus et al., 2016). The ecosystem approach tends to encompass all human activities impacting
marine ecosystems. It requires extended suites of indicators, to inform on the status of marine environment
and assess the efficiency of management policies. The Marine Strategy Framework Directive (MSFD,
2008/56/EC) aims at achieving Good Environmental Status (GES) in European Union (EU) waters by 2020.
Good Environmental Status of EU waters will be assessed based on indicators collected in all EU regional
seas, within the framework of national monitoring programmes (MPs). For cost-efficiency reasons, MPs
have been developed by extending existing monitoring activities (Shephard et al., 2015).

Fisheries surveys conducted for informing the EU Common Fisheries Policies (CFP) could play a central
role in national MPs (Shephard et al., 2015), because they are the only platforms that routinely sample the
interior of offshore EU waters. However, the potential for integration between fisheries and environmental,
MSFD-related surveys is not clearly assessed (ICES, 2012), because extending sampling activities cannot be
carried out endlessly without compromising data quality (Shephard et al., 2015). Besides questions about the
amount and quality of ecosystem data that will be eventually provided by MPs, significant gaps also remain
on the means of combining the indicators, criteria and descriptors into a holistic assessment of the
environmental status (Borja et al., 2013). This paper is a practical contribution to two important debates
raised by the implementation of marine ecosystem-based management: i) which ecosystem data can be
provided by a fisheries survey, optimised for ecosystem monitoring; and ii) how to combine/select potential
indicators to derive useful information on marine ecosystem status and dynamics? A list of potential pelagic
ecosystem-state indicators, derived from an integrated ecosystem survey is first presented. A methodology is
proposed for selecting a reduced set of relevant indicators to inform on large-scale ecosystem dynamics.

This article is based on data collected during the PELGAS integrated ecosystem survey conducted by
Ifremer, in collaboration with La Rochelle University and CNRS since 2000 in the BoB (Doray et al., 2000; Doray et al., this volume-a). The PELGAS survey's primary objective is to assess the BoB small pelagic fish biomass in springtime, to inform EU CFP. However, the survey has pioneered the move toward ecosystem monitoring, by rapidly extending its sampling scheme to cover the BoB pelagic main ecosystem components since the mid 2000's (Doray et al. this volume-a). The PELGAS survey consortium has moreover been producing numerous potential indicators of the state of the various components of the BoB pelagic ecosystem: hydrology, phytoplankton, meso-zooplankton, small pelagic fish and megafauna (cetaceans, seabirds, and large fish).

The focus of the analysis is concentrated on the small pelagic fish ecosystem component, for which the most detailed information is available, but temporal dynamics of other ecosystem components are also analysed to depict the small pelagic fish within their biotic and abiotic environments. Indicators displaying the most continuous variations in time are selected in each ecosystem compartment, by adapting a time-series analysis methodology developed by the Fisboat project (Woillez et al., 2009a). A selection procedure is applied to indicators collected in different components, and across components, to highlight the most continuous large-scale ecological process and their potential interactions. Selected survey-based indicators series are eventually compared to time series of external forcing variables to assess the potential effects of climate, survey timing, river discharge and fishing on the BoB pelagic ecosystem low-frequency dynamics. The perspectives for the assessment of marine ecosystem status on the basis of suites of indicators derived from integrated ecosystem surveys are eventually discussed.

# Material and methods

## Survey data

Data were collected during the spring acoustic surveys PELGAS (Doray et al., this volume-a), 2000–2015, conducted in May by IFREMER on board RV "Thalassa", over the French shelf of the BoB. The survey design comprised parallel line transects, perpendicular to the isobaths and regularly spaced 12 nautical miles (n.m.) apart, from 43.58N to 48.88N and from the coast (20 m depth) to the shelf break (Figure 1). Diurnal fish acoustic densities were registered at 10 knots along the transects at the 38 kHz frequency, using monobeam echosounders, calibrated using a standard method (Foote et al., 1987). Midwater trawl hauls were adaptively performed to inform the relative species compositions and biological parameters (length, weight
and age) of fish echotraces, using R/V Thalassa pelagic otter trawls with a headline of 76 m and a foot rope of 70 m (or 57 m × 52 m at depths below 50 m). Marine megafauna (marine mammals, marine turtles, large fish, birds), macro-litters (i.e. litter larger than 30 cm) and ships were recorded during the daytime along acoustic transects, by 2 trained observers. During the night, vertical casts of a Seabird SBE 19 conductivity-temperature-depth (CTD) probe fitted with a fluorometer and Niskin bottles, and of three 57 cm diameter, 200 µm mesh WP2 mesozooplankton nets, were performed on a regular grid of stations (Figure 1), to collect information on the hydrology, phytoplankton and mesozooplankton ecosystem components. Details on the PELGAS survey protocols can be found in Doray et al. (this volume-a) and (2014).

## Indicators derivation

Deriving potential indicators of the BoB pelagic ecosystem temporal dynamics involved 4 steps: i) selecting parameters relevant to characterise the temporal dynamics of each ecosystem component; ii) spatial smoothing of the georeferenced parameters on a common grid, to produce spatially comparable standard gridded maps; iii) selecting common spatial zones and time periods for the analysis; and iv) deriving statistics and spatial indicators from gridded maps as potential indicators of the ecosystem component dynamics. Table 1 presents the ecosystem parameters derived from PELGAS survey data that were selected to compute potential indicators of the BoB pelagic ecosystem temporal dynamics, with references to raw data measurements protocols and parameter derivation methodologies.

Hydrological indicators were chosen to characterise the main features of the BoB physical environment: i.e. a stratification indicator (deficit of potential energy), riverine influence (equivalent freshwater height) and thermal forcing (bottom and surface temperature) (Table 1). Vertically integrated chlorophyll-a (Chl-a) concentration derived from Niskin bottles samples was selected to characterise the phytoplankton abundance, including the sub-surface Chl-a maxima commonly observed in the area in spring and summer (Gohin et al., 2005). The total mesozooplankton biomass collected with WP2 nets in the entire water column was used to characterise the mesozooplankton standing stock. The small pelagic fish species that were consistently observed in the BoB during the PELGAS series were selected to characterise the fish component. Selected species include Atlantic chub mackerel (*Scomber colias*), Atlantic mackerel (*Scomber scombrus*), horse mackerel (*Trachurus trachurus*), blue whiting (*Micromesistius poutassou*), European anchovy (*Engraulis encrasicolus*), European sardine, (*Sardina pilchardus*) and sprat (*Sprattus sprattus*). Small pelagic fish were characterised by their biomass per species, estimated in 1852 m long sub-sections of line transects, as well as
in the whole BoB. Available parameters describing the species size distribution (fish length quantiles) and mean weight-at-age of anchovy and sardine at the BoB scale were also retained (Table 1).

Megafauna was represented by the groups of cetacean and seabird species that had been consistently observed in the BoB during PELGAS. Selected groups of megafauna species comprised auks, gannets, and small and big delphinids. These animals were characterised by their densities, estimated following standard line transect methodology for density estimation, with distance-sampling methods (Buckland et al., 2015), in 2000 m long sub-sections of the line transects (Table 1).

Spatially explicit ecosystem parameters are presented in Table 1, and were averaged over a common grid to compensate for spatial disparities in sampling schemes between components. The spatial interpolation was performed using the block-averaging methodology described in Petitgas et al. (2009, 2014). The grid mesh size selected was 0.25° in latitude and 0.25° in longitude, with origin x0 at 43°N and 6°W. This grid was chosen to maximise the number of samples of fish, megafauna and hydrobiological parameters averaged in each block (or cell) on one hand, and the grid spatial resolution on the other hand. The spatial smoothing procedure provided annual gridded maps of ecosystem parameters over the time periods specified in Table 1. An estimation polygon comprising the largest area consistently sampled in all components during the survey series was defined (Figure 1). Time periods during which available parameters had been collected in all of the estimation polygon were identified in each biotic ecosystem component: 2009–2015 for phytoplankton and mesozooplankton, 2000–2015 for fish, 2004–2014 for megafauna.

Spatial distributions of marine organisms have been extensively studied, because they reveal the effects of internal populations dynamics (e.g. MacCall, 1990; Petitgas et al., 2010), environmental forcing (Perry et al., 2005; Petitgas et al., 2014) or human activities, such as commercial fishing (Greenstreet and Hall, 1996; Fisher and Frank, 2004). Taking advantage of the good spatial resolution of gridded maps, indicators of the spatial distribution of biotic parameters derived from PELGAS surveys were included in the analysis. Spatial indicators proposed by Woillez et al. (2007) and Woillez et al. (2009b), to capture spatial patterns of fish population were computed on biotic indicators gridded maps values, within the estimation polygon, over the above-mentioned time periods. Spatial indicators used in this study included two spatial-location indicators (longitude and latitude of center of gravity, xcg, ycg), three indicators of spatial spread (inertia, I; Isotropy, Iso; positive area, PA) and two indicators of spatial distribution granularity (equivalent area, EA; microstructure index, microS). Spatial indicators were computed with the following parameters for all
ecosystem parameters: no. of discretisation cells to compute areas of influence and microstructure index:
200; maximum distance of influence of samples: 15 nautical miles (NM); mean lag between samples: 10
NM. Mean values and coefficients of variations of annual gridded maps of ecosystem parameters were used
as potential indicators of the average conditions and spatial variability in the hydrology component.

The indicator derivation procedure led to the production of a suite of 143 potential ecosystem indicators from PELGAS survey data, summarised in Table 1.

Min/max autocorrelation factors

The BoB is a highly dynamic system (Koutsikopoulos and Le Cann, 1996) that displays marked interannual variability in springtime (Huret et al., this volume). In an attempt to retain essential information without unnecessary detail (Levin, 1992), we assumed that the essential information needed to inform on ecosystem status and dynamics were associated with the low-frequency components of the indicator time series, whereas higher-frequency components could be treated as noise at the scale of our analysis. Following Woillez et al. (2009a), we used minimum/maximal autocorrelation factors (MAFs) to select the ecosystem indicators that carried the most continuous low frequency components, out of the large suite of potential indicators.

Minimum/maximal autocorrelation factors is a multivariate statistical method, that was first proposed as an
alternative to Principal Component Analysis (PCA), to separate signals from noise in multivariate imagery
analysis (Switzer and Green, 1984). Instead of seeking, as in PCA, a transformation that maximises the
variance of projected data, MAFs aim to maximise the autocorrelation between neighbouring observations.
When applied to multivariate time series, minimum/maximal autocorrelation factors analysis (MAFA)
decomposes the original suite of variables into a series of MAFs, in which autocorrelation decreases from the
first factor to the last (Shapiro and Switzer, 1989; Woillez et al., 2009a). Selecting the first MAFs enables to
extract the time-series components that are the most continuous in time, in a way akin to selecting the first
sinusoid cycles in a Fourier transformation.

Minimum/maximal autocorrelation factors analysis is based on the implementation of 2 successive PCAs:
the first one transforms the multivariate time series into principal components (PCs). The second PCA is
performed on the PCs increments, so as to maximise/minimise their variance, and hence the PCs
autocorrelation at a chosen time lag (Shapiro and Switzer 1989; Woillez et al., 2009a). Minimum/maximal
autocorrelation factors are particular linear combinations of the original variables that are uncorrelated with
each other at all times. They are also uncorrelated with each other at a given time lag (Shapiro and Switzer
1989), here chosen as the time between successive PELGAS surveys (1 year) (Woillez et al., 2009a).
Another practical advantage of MAFs is to be applicable to relatively short, non-stationary, multivariate time
series of indicators (Erzini, 2005).

Indicator selection procedure

We adapted the MAF-based methodology proposed by Woillez et al. (2009a), to select the most relevant
indicators of the BoB pelagic ecosystem large-scale dynamics. The method is summarised in Figure 2 and
involves three steps: i) initial selection of the most continuous indicators; ii) calculation of MAFs on this
restricted set of indicators; and iii) calculation of a continuity index, on which the final indicator selection is
based.

The initial selection of indicators was based on the one-lag variogram of each indicator, scaled to the
indicator variance and ranked (Woillez et al., 2009a). Indicators with a normalised one-lag variogram values
less or equal to 1 were considered as sufficiently continuous to be included in the analysis. However, as
MAFA is based on PCA, the number of MAFs cannot exceed the number of variables included in the
analysis, or the number of years minus 1. Indeed, if the number of variables is equal to or larger than the
number of years, the time series of MAF k (k = 1, 2, ...) will have a constant period equal to (no.years - 1) x
2/k, which is biologically meaningless (Woillez et al. 2009a). Woillez et al. (2009a) proposed a robust
procedure to overcome this MAFs limitation, based on adding random white noise to variables during MAFs
computation. As adding random noise to our data did not overcome the MAF maximum number of variables
limitation, we have adopted a procedure ensuring to select MAFs with biologically meaningful periods.

Let us denote *N* the number of indicators with a normalised one-lag variogram value equal to or lesser than one. If *N* was lesser than the number of years minus one, the selection procedure was conducted using *N* indicators (Figure 2). If *N* exceeded the number of years minus one, a sensitivity analysis was performed to determine the maximum number of indicators to be included in MAFs computation (Figure 2). The sensitivity analysis involved first ranking the *N* indicators in ascending order of their normalised one-lag variogram. A first MAF set was computed with the first two indicators in the list. A total of *N*-2 other MAF sets were therefore computed, while sequentially adding one new indicator from the list at each step, until including *N* indicators. The MAF1&2 periods of each MAF set were plotted against the number of indicators used in the computation. The maximum number of indicators to use in the MAF computation was defined as
the maximum number of indicators for which MAF1&2 periods were not constant (Figure 2).

After MAF computation, a continuity index $c_j$ was derived after Woillez et al. (2009a) (Figure 2):

$$c_j = \sum_{k=1}^{p} \lambda_{k,j}^2 (1 - \mu_k)$$

where $\lambda_{k,j}$ denotes the loading of variable j on MAF k and $\mu_k$ the one-lag variogram value of MAF k. This continuity index was computed based on the loadings and one-lag variogram values of the first 2 MAFs ($p=2$). The threshold used to select the most continuous indicators was defined as the $c_j$ values beyond which no substantial drop in the $c_j$ values ranked in decreasing order was observed. Potential indicators whose $c_j$ values were higher than this threshold were selected as the most representative of the changes occurring along the time series (Figure 2). This method will be referred to as the “MAF-based indicator selection procedure” in the rest of the text.

## Indicator selection strategy

The number of potential ecosystem indicators derived from PELGAS data (143) was far larger than the maximum number of year increments in the longest time series (16). Therefore, applying the MAF-based indicator selection procedure described in the previous section to all potential indicators at once was not possible, owing to the limitation on the maximum number of variables to be included in MAFA. A stepwise indicator selection strategy summarised in Figure 3 was adopted to accommodate the large suite of potential ecosystem indicators derived from the PELGAS survey. This approach involved first applying the MAF-based selection procedure (Figure 2) to the indicators of each ecosystem components separately (Figure 3).

The numbers of indicators and time series lengths in the hydrology, phytoplankton and mesozooplankton components allowed for the direct application of the MAF-based indicator selection procedure (Figure 3). In the small pelagic fish and megafauna components, the number of indicators was too large relative to the time series length, to allow for a direct application of the indicator selection procedure. Indicator selection was thus first performed per species and group of species in the fish and megafauna components, respectively (Figure 3). Indicators selected for each species and group of species were then pooled to define sets of indicators that were the most representative of the fish and megafauna components (Figure 3). Finally, a cross-ecosystem component analysis was conducted for the components with the longest time series: hydrology and fish. The most representative indicators in the hydrology and fish components were pooled and submitted to the MAF-based indicator selection procedure (Figure 3). The results were analysed to
compare the continuity of the processes occurring in both components, and to identify potential hydrological
forcing on fish.

Ecosystem external forcing

External forcing variables were selected to assess the potential effects of climate, survey timing, river discharge and fishing on the PELGAS survey final set of hydrology-fish indicators (Table 2). Climate forcing was represented by two sets of variables. The first set of climate proxies comprised the North Atlantic Oscillation (NAO), East Atlantic (EA) and Scandinavian (SCAND) indices computed by the National Oceans and Atmosphere Administration Climate Prediction Center (http://www.cpc.ncep.noaa.gov/). These indices represent patterns of low-frequency variability of the atmosphere (Wallace and Gutzler, 1981), and are known for their impacts on precipitation and/or temperature in the northeast Atlantic, and for subsequent potential effects on marine ecosystems (e.g. Borja et al., 2008; Reid et al., 2001). Daily climate indices were averaged over the year preceding the PELGAS surveys (June to May), and by trimesters, defined relatively to the survey average timing (in May), so that trimester 1 (T1) was from June to August, trimester 2 (T2) from September to November, trimester 3 (T3) from December to February, and trimester 4 (T4) from March to May (Table 2). The second set of climate proxies considered were the daily occurrences of 4 North Atlantic-Europe weather regimes (NAE-WR) identified by Cassou et al. (2004) and Cassou et al. (2011). Cassou et al. (2011) showed that weather regimes, and associated wind regimes, can be good predictors for surface ocean variability, which in turn can affect marine-organism connectivity patterns (Thomas et al., 2014), and potentially marine ecosystem dynamics. Occurrences of NAE-WR computed over winter months (December to February) before annual PELGAS survey were selected to represent the potential effects of weather regimes on the survey selected hydrology and fish indicators (Table 2).

River discharge is known to have a major influence on the BoB pelagic ecosystem functioning (Koutsikopoulos and Le Cann, 1996). Water discharges of the main rivers (Vilaine, Loire, Garonne, Adour) flowing onto the BoB French platform were averaged per year before the survey (June to May), and per survey-referenced trimesters, to test for the riverine influence on PELGAS selected ecosystem indicators (Table 2).

Huret et al. (this volume) proposed a corrected date for the PELGAS survey, aimed at positioning the survey within the seasonal ecosystem cycle. Annual means of PELGAS corrected dates were selected to test the
effect of survey timing on the selected ecosystem indicator series.

The fishing pressure exerted on small pelagic fish populations sampled in the BoB during PELGAS was represented by the fishing mortality (F) estimates derived from the analytical stock assessment models implemented for anchovy (ICES, 2016a), blue whiting, horse mackerel and Atlantic mackerel (ICES, 2016b). Because no stock assessment model was available for sardine, total landings were also included in the analysis, as a proxy for fishing effort. Stock spawning biomass (SSB) of analytically assessed fish species was also selected, to account for the effects of intrinsic population dynamics, possibly occurring outside of the BoB, in the case of widely distributed species (blue whiting, horse mackerel and Atlantic mackerel) (Table 2).

## Correlations between ecosystem indicators and external forcing

Correlations between time series of selected ecosystem indicators and external forcing variables were analysed, to detect potential ecological interactions. Spearman's rho statistic (Spearman, 1904) was used to compute correlation coefficients, as it is a nonparametric measure of rank correlation, less sensitive to non-normal data, and capable of assessing linear and nonlinear relationships between variables. To avoid elaborating on spurious (i.e. accidental correlations, Haig, 2003), pairs of indicators time series with absolute correlation larger than 0.6 were inspected, to retain only those for which plausible ecological relationships could be hypothesised. A relationship between two indicators was deemed ecologically plausible if: i) a well document ecological process could explain the relationship; and ii) the scales of the processes at play were compatible. Specifically, intra-guild competition for food, egg predation and egg cannibalism in case of shortage of planktonic preys are the main potential biotic interactions between small pelagic fish in the BoB (Bachiller & Irigoien, 2015). In the absence of evidence of planktonic prey shortage, and therefore of presence of intra-guild competition, significant correlations between indicators of small pelagic fish of different species derived on one hand from the survey, and on the other hand from stock assessment, were not retained.

Statistical computations were performed using the R statistical language version 3.3.1 (R Core Team, 2016), with the EchoR package for fish biomass computations and block averaging (Doray et al., 2016), the SpatialIndices package (Woillez, 2013) for spatial indices computation. The R codes used were those proposed by Woillez et al. (2009a) for MAF analysis, and the marmap (Pante and Simon-Bouhet, 2013), rgdal (Bivand et al., 2016), raster (Hijmans, 2016) and rasterVis (Perpinan Lamigueiro and Hijmans, 2016)
packages for mapping.

# Results and Discussion

Indicators selected to represent each BoB pelagic ecosystem component, and the hydrology and fish
components combination, are summarised in Table 3.

## Hydrology

Indicator selection in the hydrology component confirmed that thermal forcing and riverine influence are important drivers of the BoB hydrology (Koutsikopoulos and Le Cann, 1996; Woillez et al., 2010). Normalised one-lag variograms and continuity indices of the 8 hydrology potential indicators are presented in Figure 4a-b. The 5 most continuous indicators, with respect to variogram-at-lag 1 values (Figure 4a), were retained for MAF computation. A total of 3 indicators were selected to represent the hydrology component, based on the continuity index criterion (Figure 4b): equivalent freshwater height average and bottom temperature average and coefficient of variation. The mean bottom temperature was relatively stable during 2000–2015, except for 4 years with below-average bottom temperatures: 2005, 2006, 2009, and 2010 (Figure 4c).

Previous studies have shown that bottom temperature is more relevant than surface temperature for characterising small pelagic fish spawning (Planque et al., 2007) or adult (Petitgas et al., 2014, Doray et al., this volume-b) habitats during spring. Our results confirm that in-situ bottom temperature time series derived from the PELGAS survey are more continuous than surface temperature ones. This is probably due to bottom temperature recorded during PELGAS having less high frequency spatial variability than surface temperature, because it is not affected by seasonal stratification and storm-induced temporary mixing of surface layers, which occurred during the survey (Huret et al., this volume). Long-term analysis of sea temperature trends in the BoB from early 1970s to 2000s have revealed a general increase of 0.1°C per decade in the 0 to 100 m depth layer (Michel et al., 2009), and 0.2°C in surface temperature (deCastro et al., 2009). Such temperature trends were not detected in the in-situ bottom (or surface) temperature series derived from PELGAS data, because they were probably masked by the seasonal variability of sea temperature at the time of the survey (~1°C on average for sea surface temperature, Huret et al. this volume).

The mean equivalent freshwater height has decreased since 2000, until reaching the series minimum in 2011.
It has since then increased to reach in 2014–2015 levels similar to those observed at the beginning of the
series (Figure 4d). The average equivalent freshwater height time series has shown greater variability than
bottom temperature, and its annual spatial variability (gridded map CV) was negatively correlated with the
map average, indicating a higher spatial variability in case of lower riverine influence (Figure 4e).

Phytoplankton

Normalised one-lag variograms and continuity indices of the 8 phytoplankton potential indicators are
presented in Figure 5a-b. The 5 most continuous indicators, with respect to variogram-at-lag 1 values were
retained for MAF computation (Figure 5a). A total of 3 indicators were selected to represent the
phytoplankton component, based on the continuity index criterion (Figure 5b): vertically integrated Chl-a
average (Figure 5d), equivalent area (Figure 5e), and microstructure index (Figure 5f). Mean Chl-a biomass
has rapidly increased from 2003 to 2005, and steadily decreased thereafter, with secondary drops in 2006 and
2013 (Figure 5c). Focussing on the 2009–15 period, for which spatial indices are available, reveals an
increase of the Chl-a biomass microstructure index, and a synchronous decrease of the map average and
equivalent area in 2012–13. This indicates that the Chl-a biomass distribution was less continuous, including
micro-scale patches and isolated higher density values, when the map average Chl-a biomass was lower in
2012–13 (Figure 5d-f). Annual gridded maps confirm that most of the phytoplankton biomass observed in
2013 was concentrated in a localised bloom near the shelfbreak at ~46°N (Figure S1, supplementary
material).

The steady decrease in Chl-a average biomass derived from PELGAS data in springtime from 2007 to 2015
(Figure 5c) conforms with global models forecasting that climate change could cause long-term decreases in
primary production in the North Atlantic, including the BoB, due to increased stratification and nutrient
decrease (Steinacher et al., 2009; Chust et al., 2011). This forecast should be considered with caution,
because the Bay Biscay nutrient concentration also heavily depends on river discharges; the evolution of
which is not properly modelled at regional scales (Chust et al. 2011). Moreover, annual gridded maps of Chl-
a biomass reveal that the greatest concentrations of Chl-a biomass were almost exclusively located near the
shelf break, north of 44°N (Figure S1, supplementary material). This implies that most of Chl-a biomass
observed during the PELGAS survey might be produced by residual post-winter blooms, known to occur in
offshore areas in March and April (Pingree and Garcia-Soto, 2014). Based on survey data only, it is not
possible to ascertain if the observed decrease in the survey Chl-a index resulted from changes in blooms
phenology, survey timing, or a combination of both factors, or from a genuine decrease of the BoB primary
production. PELGAS survey data should be combined with ecosystem model outputs and satellite data to further investigate this issue.

## Mesozooplankton

Normalised one-lag variograms and continuity indices of the 8 mesozooplankton potential indicators are presented in Figure 6a-b. The 5 most continuous indicators, with respect to variogram-at-lag 1 values (Figure 6a), were selected for MAF computation. A total of 3 indicators were finally retained to represent the mesozooplankton component, based on the continuity index criterion (Figure 6b): vertically integrated mesozooplankton biomass equivalent area (Figure 6d), positive area (Figure 6e), and average (Figure 6f). Over the 2006–2015 period, the mean mesozooplankton biomass peaked in 2007, then decreased until 2013 and increased again after 2014 (Figure 6c). Focussing on the 2009–2015 period, for which spatial indices are available, reveals that the mean mesozooplankton biomass decreased since 2010, reaching a minimum in 2013, then increased again to reach 2010 levels in 2015. The positive area followed a similar trend, indicating that the mesozooplankton spatial extension is positively correlated with biomass fluctuations. The equivalent area dropped concurrently to the mean mesozooplankton biomass increase since 2014, suggesting a recent increase of the distribution patchiness, confirmed by gridded maps (Figure S2, supplementary materials). Annual gridded maps reveal that the highest concentrations of mesozooplankton sampled during the PELGAS surveys were located in mid-shelf areas, in the north western part of the bay (north of 46°N).

Previous mid- to long-term studies did not observe any significant trend in mesoplankton biomass in the Southern BoB and Eastearn Cantabrian Sea (Irigoien et al., 2009; 1998–2006), or limited changes in Biscay offshore waters (Pitois and Fox, 2006; 1960-1999) or the western Cantabrian Sea (Valdés et al., 2007; 1992–2006). A longer series of mesozooplankton biomass, ideally separated by size, and/or by main taxonomic groups, would necessary to further interpret the trends observed in the PELGAS survey mesozooplankton indicators.

## Small pelagic fish

Normalised one-lag variograms and continuity indices of the 87 small pelagic fish potential indicators are presented in Figure 7a-b. The 14 most continuous indicators, with respect to variogram-at-lag 1 values (Figure 7a), were retained for MAF computation. Indicators of small pelagic fish species whose core distributions areas are located in the BoB (anchovy, sprat, and small sardine) appeared to be more continuous, and were therefore preferentially selected by the MAF-based procedure. The biomass estimate of
the more cosmopolitan horse mackerel was also retained, as it displayed a very smooth trend (Figure 7). A total of 12 indicators were selected to represent the small pelagic fish component, based on the continuity index criterion (Figure 7b): *Sardina pilchardus* mean weight at age 1 (Figure 7c), mean weight at age 2 (Figure 7d), center of gravity longitude (Figure 7e), and length 25th percentile (Figure 7f); *Trachurus trachurus* biomass estimate (Figure 7g); *Sprattus sprattus* equivalent area (Figure 7h) and isotropy (Figure 7i); *Engraulis encrasicolus* biomass estimate (Figure 7j), isotropy (Figure 7k), mean weight at age 1 (Figure 7l), mean weight at age 2 (Figure 7m), and positive area (Figure 7n).

*Sardina pilchardus* mean weight at age 1 steadily decreased since 2001, in spite of a peak in 2003. Mean weight at age 2 dropped between 2009 and 2012 and remained low ever-since. *Sardina pilchardus* center of gravity longitude has steadily increased since 2015, indicating a sardine distribution shift toward coastal areas. *Trachurus trachurus* biomass markedly decreased from 2006 to 2012, with a peak in 2011, and has consistently increased since 2013. *Sprattus sprattus* equivalent area decreased from 2000 to 2012, and increased after 2013. Its isotropy decreased from 2000 to 2011, with a marked drop in 2005, and steadily increased after 2012. This result indicates that the sprat distribution has redistributed along the coast from 2000 to 2011, and expanded afterwards, with increasing patchiness from 2000 to 2012. Gridded maps inspection reveals that these trends in the sprat spatial indicators were caused by density variations in its core distributional areas in Gironde and Bretagne. *Engraulis encrasicolus* biomass collapsed from 2003 to 2005, leading to a closure of the BoB fishery in 2005. The subsequent anchovy biomass increase from 2006 to 2015 led to re-opening of the fishery in 2010. Positive area and isotropy were positively correlated with biomass, indicating that the anchovy spatial distribution shrunk and elongated along the shore when they the biomass decreased. This result tends to confirm the existence of density-dependent processes influencing the spatial dynamics of this species in the BoB (Petitgas et al., 2014). Anchovy mean weight at age 1 steadily decreased from 2001 to 2015, with high frequency oscillations around its linearly decreasing trend. Mean weight at age 2 dropped in 2012 and remained low thereafter.

## Megafauna

Normalised one-lag variograms and continuity indices of the 32 megafauna potential indicators are presented in Figure 8a-b. The 14 most continuous indicators, with respect to variogram-at-lag 1 values (Figure 8a), were retained for MAF computation. A total of 5 indicators were selected to represent the megafauna component, based on the continuity index criterion (Figure 8b): gannets center of gravity longitude (Figure
8c), inertia (Figure 8d), and equivalent area (Figure 8e); auks positive area (Figure 8f), and small delphinids gridded map microstructure index (Figure 8g). Gannets is the megafauna species group most commonly detected during PELGAS (Authier et al., this volume) in coastal and midshelf areas of the BoB. Its center of gravity longitude decreased from 2005 to 2011, before increasing from 2012 to 2014, with a peak in 2013. It was negatively correlated with equivalent area. This indicates that the gannets distribution has shifted to western offshore areas, while displaying less high density patches from 2005 to 2011, before shifting back toward eastern Biscay coastal areas, and becoming more patchy. Gannets inertia has increased from 2004 to 2014, with 4 years period oscillations (maxima in 2009 and 2013, minima in 2006 and 2010). This suggests that the gannets spatial distribution has consistently spread around its center of gravity over the series. Small delphinids were commonly observed near the shelbreak and in mid-shelf areas during PELGAS. Their microstructure index has increased from 2004 to 2008, decreased until reaching a minimum in 2011, and increased again since 2012, revealing trends in small delphinids microscale patchiness. Auks have been consistently observed during PELGAS in the coastal Gironde area. Their positive area decreased after 2004 until reaching a minimum in 2006. It therefore displayed high values after 2009 (with a maximum in 2011). This reveals that the maximum extension of auks distribution has shrunk at the beginning of the series, before expanding again. Map averages of megafauna species group densities displayed less continuous trends than spatial indices and were not selected.

### Cross ecosystem components indicator selection

The results of cross ecosystem components indicator selection were solely presented for the joint analysis of hydrology and small pelagic fish, for which long enough (16 years) time-series data were available. Normalised one-lag variograms and continuity indices of the 15 hydrology and small pelagic fish potential indicators are presented in Figure 9a-b. The 13 most continuous indicators, with respect to variogram-at-lag 1 values (Figure 9a), were retained for MAF computation. A total of 9 indicators were selected to represent the hydrology and small pelagic fish components, based on the continuity index criterion (Figure 9b): equivalent freshwater height average (Figure 9c) and CV (Figure 9d), *E. encrasicolus* biomass estimate (Figure 9e), mean weight at age 1 (Figure 9f), and isotropy (Figure 9g), *S. pilchardus* mean weight at age 1 (Figure 9h), and length 25th percentile (Figure 9i), *T. trachurus* biomass estimate (Figure 9j), and *S. sprattus* isotropy (Figure 9k).

Figure 10 presents a graphical summary of the first 2 MAFs computed based on the 13 most continuous
indicators representing the hydrology and fish component. MAF1 has displayed positive values from 2000 to 2009, dropped below zero in 2010 and subsequently steadily decreased (Figure 10a). The variogram in Figure 10b characterises the temporal autocorrelation of the MAF1 time series at different temporal lags (x-values). It indicates that the temporal autocorrelation in MAF1 values has a period of 10 years, corresponding to the time lag (or range) at which the variogram reaches a plateau (or sill). Temporal autocorrelation almost appears as a linear trend in MAF1, as its period represents two third of the time series length. Figure 10c presents the contribution of indicators to MAF1. Indicators contributing the most to MAF1 have negative loadings and include sprat and anchovy distributions isotropy, anchovy mean weight at age 1, horse mackerel biomass and average freshwater height. These indicators have displayed higher values since 2010, leading to the overall decrease of MAF1. MAF2 has displayed positive values at the beginning (2000-2004) and end (2014-2015) of the series, and negative ones from 2005 to 2013 (Figure 10d). MAF2 spatial autocorrelation has a period of 5 years (Figure 10e), appearing as periodic oscillations in its time series (Figure 10d). The indicators contributing most to MAF2 had positive loadings and included: average freshwater height, anchovy distributions isotropy, sardine and anchovy mean weights at age 1 and horse mackerel biomass. MAF2 hence accounts for ~5-year period oscillations around the lower period, almost linear, MAF1 trend.

Indicator selection was also performed across all sampled ecosystem components, over the 6-year period (2009–2014) for which consistent time series indicators were available. As the time period was deemed too short, results were not shown and time-series trends were not compared. Nonetheless, a result worth mentioning is that the selected indicators interestingly belonged to 3 different ecosystem components: phytoplankton (Chl-a biomass map average), mesozooplankton (mesozooplankton biomass positive area), and fish (mean sardine weight-at-age 1&2). Longer time series should obviously be collected to draw further conclusions on the BoB pelagic ecosystem dynamics, but these first results suggest that the MAF-based selection methodology proposed in this paper could efficiently select indicators of major ecological processes across ecosystem compartments.

## Correlations between ecosystem indicators

Correlation analysis between selected survey indicators (Table 4) revealed an inverse relationship (rho = -0.74) between equivalent freshwater height (Heq) map average and CV, indicating that the Heq spatial variability increases when the mean riverine influence decreases. No significant correlation was found
between Heq and fish indicators. However, potential ecological interactions between river plumes and small pelagic fish populations occurred at specific scales (Levin, 1992), which were not necessarily those at which indicators have been computed. Rivers indeed deliver nutrients required for early season plankton production (Labry et al., 2001) in small pelagic fish coastal spawning areas. Equivalent freshwater height was moreover described as one of the main hydrological factor structuring the distribution of small pelagic fish communities in the BoB in spring (Doray et al., this volume-b). River plumes probably play an important role in small pelagic fish life cycles, but at smaller scales than the BoB scale at which indicators were computed.

Anchovy distribution isotropy was positively correlated (rho=0.63) with total anchovy biomass, suggesting that the anchovy distribution elongates when its biomass decreases. Positive correlations (rho=0.64) between horse mackerel biomass and anchovy and sardine weights at age 1 were not deemed ecologically meaningful, as those fish do not occupy the same habitats in springtime (Doray et al., this volume-b). Sardine mean weight at age 1 was positively correlated with sardine length 25th percentile, which was expected, as fish length and weight classically covary. The quasi-absence of correlation between selected fish indicators in Table 4 indicates that the BoB small pelagic fish species have followed distinct trajectories from 2000 to 2015. This suggests that the biotic interactions between small pelagic fish species (intra-guild competition for food, egg predation and egg canibalism, cf. Bachiller and Irigoien, 2015 and Bachiller et al., 2015) might be limited. A notable exception is the positive correlation (rho=0.72) between sardine and anchovy mean weights at age 1 (mW.age1) and 2 (mW.age2).

Anchovy and sardine mW.age1 followed strikingly similar decreasing trends from 2001 to 2015. Sardine mW.age2 has dropped 2 years before (2010) anchovy mW.age2 (2012), but both indicators have remained low thereafter. Fish mean lengths followed the same trends. The decreases in mW.age1 of anchovy and sardine and of mW.age2 of anchovy seem to be caused by increases in the abundances (in number of fish) of young (ages 1–2) fish, at a higher rate than concurrent biomass increases. Year 1–2 anchovy and year 1 sardine abundances have steadily increased since 2007, while no particular trend was observed in anchovy age 3+, and sardine age 4+ mean weights, biomasses and abundances (Figure 11). The increase in age 2 sardine abundance is however less marked than for anchovy. These results suggest that the survival rates of anchovy and sardine young stages and age 1 fish might have increased since 2007, causing increases in age1–2 fish abundances.
The carrying capacity of young anchovy and sardine habitats seemed to have been sufficient to support a concurrent increase of young fish biomass, but at a slower rate than the abundance increase, hence yielding a decrease in anchovy *mW.age1* and *mW.age2* and sardine *mW.age1*. As young anchovy, sardine and sprat share the same coastal habitats in the Gironde area from late autumn to summer (Doray et al., this volume-b; ICES, 2010), the almost synchronous trends observed in both species mean weights at ages 1 and 2 might be due to density-dependent competition for food occurring in these areas. Further investigations involving the use of anchovy and sardine life cycle simulation models (e.g. Politikos et al., 2015) are needed to understand the drivers, circumstances and ecological implications of these synchronous decrease of anchovy and sardine mean weights at ages 1 and 2. Contrary to what was recently observed in the Mediterranean Sea (Van Beveren et al., 2014), older and larger anchovy and sardine were still observed in the BoB, with unchanged mean weights (Figure 11).

## External forcing

Most of the correlations between selected survey indicators and potential external forcing variables were considered spurious, as they could not be explained by well documented ecological process, and/or occurred at incompatible scales (Table 5). A negative correlation was found between the PELGAS horse mackerel biomass estimate in the BoB in spring, and the species SSB at the scale of the stock. This correlation was not deemed ecologically plausible, as the survey covers an area too small compared to the overall distribution of the stock, and as the connectivity patterns are too poorly known to safely assume that SSB fluctuations at the stock scale may influence the horse mackerel biomass observed in the BoB in springtime.

A positive correlation (rho=0.72) was found between the horse mackerel biomass estimated in the BoB in spring, and the NAO index averaged over summer months (June to August). Previous studies have suggested that horse mackerel spring feeding migration could be facilitated by an increase in the transport of Atlantic waters into the North Sea (Iversen et al., 2002), potentially correlated with positive winter NAO anomalies (Reid et al., 2001). The majority of horse mackerel sampled by the PELGAS survey in spring were likely age 1 fish (ICES, 2016b), and probably too young to undertake feeding migrations to the North and Norwegian seas (Iversen et al. 2002). Moreover, if positive summer NAO index anomalies are known to generate summer climate extremes in European (Folland et al., 2009), complementary process studies should be conducted to relate their effects to the BoB horse mackerel biomass.

Ecologically plausible relationships included positive correlations between equivalent freshwater height and
river discharges at various scales (Table 5). This result tends to confirm that the freshwater equivalent height derived from PELGAS survey data is a good proxy for riverine influence in the BoB. The PELGAS survey anchovy biomass estimate was positively correlated (rho=0.65) with the BoB anchovy SSB estimated by the CBBM stock assessment model (Table 5). This positive relationship was expected, as the PELGAS index is used as an input in the CBBM model (Ibaibarriaga et al., 2011).

In summary, almost no ecologically meaningful correlation was found between the selected ecosystem indicators derived from PELGAS and proxies of climate change, fishing pressure and seasonal ecological scheduling. If commercial fishing is thought to have played a role in the BoB anchovy stock collapse in early 2000s, the fishery closure and more conservative harvest control rules adopted in 2009 (COM, 2009; ICES, 2016c) seem to have reduced the impact of fishing on the stock, as no significant correlation was found between fishing pressure proxies and anchovy population indicators over the 2000-2015 period. The quasi-absence of significant correlation between selected survey indicators and climate indices confirm that the BoB is a highly variable area, so far not significantly impacted by large scale climate forcing (Trigo et al., 2002; Chust et al., 2011). As in previous studies (Koutsikopoulos and Le Cann, 1996) our results moreover indicate that the system is likely more influenced by local physical drivers, such as river discharge.

## Indicator trend analysis

This study was aimed at: i) characterising the general patterns over time of a large suite of relatively short (5-15 years), non-stationary time series of indicators, to select the most continuous ones; ii) assessing the potential interactions between the selected indicators; and iii) evaluating the potential relationships between selected indicators and external forcing variables. Commonly used time series analysis techniques such as spectral analysis (Priestley, 1981), wavelet analysis (Ogden, 2002), ARIMA and Box-Jenkins models (Box et al., 1994) were not suitable for achieving the first objective, as they require long (>25 y, Erzini et al., 2005) stationary time series (Zuur et al., 2003). The MAF-based indicator selection methodology proposed by Woillez et al. (2009a) was implemented, as it can be used to rank relatively short, non-stationary multiple time series of indicators, based on a continuity criterion. Because the indicators had been derived from data collected in a dynamic ecosystem, during a changing season (Huret et al., this volume), we have assumed that the most informative indicator signals resided in their lower frequency components, namely the pseudo-linear trends carried by MAF1, and the 5-years period oscillations carried by MAF2 (Figure 10). We have computed the continuity index proposed by Woillez et al. (2009a) to combine variogram at lag 1 and
MFA1,2 values in a synthetic metric, that was used to select time series oscillating at mid frequency, around a continuous low frequency trend. This methodology facilitated discarding indicators displaying flat or bumpy time series, and led to the selection of a restricted set of 26 indicators, out of the initial 143 potential ecosystem indicators derived from PELGAS survey data.

Despite Woillez et al.'s (2009a) robustification procedure, the maximum number of indicators that can be safely included in MAFA remains close to the time series length. That limitation was overcome by implementing sensitivity analysis on the number of indicators to be included (Figure 2), within a stepwise selection procedure (Figure 3). This allowed to analyse a number of indicators compatible with the time series length considered at each step (Figure 2&3).

Since the PELGAS survey initial objective was to monitor small pelagic fish populations in the BoB (Doray et al., this volume-a), more data and indicators were available to describe the fish component. The inclusion of the abundant suite of fish indicators in the analysis de facto set the focus on the fish component. A larger number of indicators (12) was indeed selected to describe this component. However, as indicators were sequentially included in MAF computation based on a continuity criterion (their one-lag variogram values), if the fish indicators had been less continuous than hydrology ones, they would not have been retained in the joint hydrology/fish indicator selection. The fact that a majority of fish indicators were selected in the joint hydrology-fish analysis indicates that their time series were more continuous than those of hydrology indicators, except freshwater height.

Spatial indicators had been so far essentially applied to the study of marine fish population. The joint analysis of indicators representing fish spatial distributions and populations status has enabled to finely characterise exploited fish populations dynamics in space and time (e.g. Fisher&Franck, 2004; Woillez et al. 2007; Gastauer et al., 2016), and to signal thresholds in marine fish stock biomass (Reuchlin-Hugenholtz et al., 2015). Our MAF-based indicator selection suggest that the spatial indicators of biotic ecosystem components can be equally (phyto-, mesozoo-plankton and fish) or more (megafauna) continuous than biological indicators, such as biomass gridded map averages. Our results indicate that spatial indicators proposed by Woillez et al. (2007) could be used to characterise the spatial dynamics of various biotic components of pelagic ecosystems, including phyto and mesozooplankton and megafauna. However, spatial indicators tend to summarise the main features of the spatial distributions of organisms, within a user-defined estimation polygon. Therefore, their relevancy to capture spatial patterns that are relevant at the scale of the
study largely depends on: i) the spatial polygon definition: the polygon must comprise a significant part of
the organisms distribution range, ii) the original map resolution.

The adoption of a common spatial polygon encompassing the core area consistently sampled during the PELGAS survey, the selection of species or group of species whose distributional range was adequately sampled during the survey, and the spatial averaging of all ecosystem parameter on the same grid, ensured that the spatial indicators derived in this paper were comparable in all ecosystem components. It is moreover advised to check the reality of trends revealed by spatial indicators by inspecting the original maps from which they were derived. Correlations between indicators should also be checked by verifying that the spatial distributions of correlated indicators overlap on the original maps.

Because selected hydrology and fish components indicators mainly involved short-lived small pelagic fish species indicators, some of them related to age 1 fish, we considered that potential ecological relationships between components and species could not occur at time lags larger than a year. Cross correlations were computed at lag 0 and 1. Only correlations at lag 0 were reported, as they showed higher coefficients than correlations at lag 1.

Dynamic factor analysis (DFA) is a multivariate time series method commonly used to deal with short, non-stationary time series (Zuur et al., 2003a; Zuur et al., 2003b; Zuur and Pierce, 2004; Erzini, 2005; Erzini et al., 2005). Dynamic factor analysis can be used to identify underlying common patterns in a multivariate time series data set, evaluate interactions between the response variables and determine the effects of explanatory variables on response variables (Zuur et al., 2003a). Application of MAFA and DFA on the same fisheries and environmental datasets provided consistent results (Erzini 2005; Erzini et al. 2005). To our knowledge, the DFA framework does not provide a variable selection procedure similar to the MAF-based variable selection methodology proposed by Woillez et al. (2009a). Our first objective being the selection of relevant ecosystem indicators out of a large suite of potential metrics, MAFA was preferred to DFA in this study.

***Implications for ecosystem management***

In this paper, we have shown that a fisheries survey optimised for ecosystem monitoring can routinely
provide a large suite of potential indicators on all the major component of a pelagic ecosystem. The MAF-
based indicator selection methodology adapted from Woillez et al. (2009a) proved to be efficient to
characterise trends in a tractable sets of indicators, selected from a large suite of potential metrics, collected
in the main components of a dynamic and diverse pelagic ecosystem.

Time series of indicators selected in this paper displayed clear trends that could be interpreted without much ambiguity. This might not be always the case, namely when the time series will get longer. The next step toward ecosystem status assessment would be to define thresholds and reference points for the selected indicators, to correctly interpret the variability in the indicator time series, and detect changes with a defined probability of errors (Petitgas, 2009). A statistical process control scheme CUMSUM (Mesnil and Petitgas, 2009) could e.g. be applied to selected ecosystem indicators derived from PELGAS survey data, to highlight alerts that occur when indicators cross reference limit values. CUMSUM results could feed e.g. a traffic light table, to inform managers on ecosystem changes in a user-friendly way (Petitgas, 2003). Another option would consist in combining the series of relevant indicators of each ecosystem component into synthetic indicators, before defining thresholds and reference points, to limit the number of traffic lights in the table.

MAFs could be used as ecosystem synthetic indicators, as exemplified in Figure 10. While studying low-frequency components is useful to identify long-term trends and conduct first order ecosystem status assessment, high-frequency components can capture potentially important discontinuities, ruptures and singularities in the indicator time series. Computing Shewhart individuals control charts (Shewhart, 1931; Woillez et al., 2010) on suites of survey indicators could e.g. help detecting extremely atypical years, to eventually analyse their occurrence and their potential influence on the system dynamics.

Our indicator-based analysis has revealed that river discharges, bottom temperature, chlorophyll-a and mesozooplankton biomass have displayed large scale trends in the BoB in springtime, since the beginning of the century. No relationship was found between these indicators, due to the shortness of plankton time series, and/or to incompatibilities between the scales at which indicators were computed. Nonetheless, the presence of large scale trends in the time series of river plume strength and lower trophic levels indicators are in line with modelling results suggesting that bottom-up processes under riverine influence play a significant role in the global structuring of the BoB marine ecosystem (Lassalle et al., 2011). On the management side, the input of nutrients by rivers is known to be a relatively important pressure in BoB coastal areas, particularly off areas of intensive agriculture and certain industries (ICES, 2016c). Our findings confirm that nutrient and organic enrichment by rivers, and the subsequent impacts on pelagic food webs, should be closely monitored to achieve effective management of the BoB pelagic ecosystem.

Small pelagic fish species appeared to have followed distinct trajectories in the BoB over the last 15 years.
This result suggests that intra-guild biotic interactions between small pelagic fish in the BoB might be
limited, except for small clupeiforms. A marked decrease in anchovy and sardine mean weights at ages 1 and
2 over the last 15 years was indeed highlighted. It might be caused by an increase of the younger fish
survival rate, leading to density-dependent competition for food in spawning habitats. Implementing spatial
management policies to protect the springtime habitats of young clupeiforms in the BoB (cf. Villamor et al.,
2007; Doray et al, this volume-b) might be an effective way to reduce the overall pressure on these
commercially important species. The inclusion of spatial indices in the portfolio of indicators used for
ecosystem management would be in this case essential to assess the efficiency of spatial management
policies (Woillez et al., 2007).

The quasi-absence of significant correlations between selected survey indicators and fishing pressure proxies
confirms that the actual exploitation rate of small pelagic fish resources present in the BoB in springtime is
likely moderate, close to the level that provides maximum sustainable yield (MSY; ICES, 2016c). This
results from the close monitoring of and strict management regulations applied to pelagic stocks in the BoB
(ICES, 2016c), which should be maintained in the wider framework of an ecosystem management.

The absence of significant correlation between survey indicators and climate proxies suggests that the effects
of large-scale climate forcing on the BoB pelagic ecosystem in springtime are so far probably limited.
PELGAS survey time series are however relatively short (phyto and mesozooplankton: 8 years; megafauna:
12 years; hydrology and fish: 16 years) to extract significant long-term trends in all components. The
standardised ecosystem data collection should be pursued during PELGAS, to enable the selection of
indicators across all components, and contribute to the assessment of the effects of global change on the BoB
pelagic ecosystem.

Conclusion

To conclude, we reiterate the importance of integrating information provided by various sources at various
scales, for achieving an effective indicator-based monitoring of marine ecosystems. The addition of new
standard data collection protocols on existing platforms, as envisaged (e.g. in the MSFD framework), might
indeed lead to the collection of data at sub-optimal time or spatial scales, as survey timing or spatial
resolution will not be adapted to study all ecological process on the same platform (Kupschus et al., 2016).
A way forward would consist in combining indicators derived from survey data with indicators provided by other sampling platforms, operational oceanography or ecosystem models (see e.g. Manderson et al., 2011). This would facilitate interpreting survey indicators in a larger, seasonal and/or inter-annual context, provided by remote sensing or model outputs on one hand (e.g. Huret et al., this volume), and to validate operational oceanography products with survey data on the other hand (e.g. Perrot et al., this volume). Integrating indicators from all available sources would bridge sampling gaps without extra sampling effort, to better interpret trends in indicators series, and ultimately improve ecosystem management.

## Acknowledgements

This study was supported by Ifremer. PELGAS surveys have been funded by the European Common Fishery Policy Data Collection Framework, the Marine Strategy Framework Directive, and Ifremer. We want to acknowledge the joint efforts of R/V Thalassa crews and PELGAS scientists that made the survey possible. We are indebted to Dr. Matt K Broadhurst for editing the English text.
References

Aminot, A., Kérouel, R., 2005. Hydrologie des écosystèmes marins. Paramètres et analyses, Ifremer. ed, Méthodes d'analyses du milieu marin.

Authier, M., Dorémus, G., Van Canneyt, O., Boubert, J.-J., Gautier, G., Doray, M., Duhamel, E., Massé, J., Petitgas, P., Ridoux, V., Spitz, J., this volume. Exploring Change in the Relative Abundance of Marine Megafauna in the Bay of Biscay, 2004-2016. Prog. Oceanogr.

Bachiller, E., Cotano, U., Ibaibarriaga, L., Santos, M., Irigoien, X., 2015a. Intraguild predation between small pelagic fish in the Bay of Biscay: impact on anchovy (*Engraulis encrasicolus*) L.) egg mortality. Mar. Biol. 162, 1351-1369. doi:10.1007/s00227-015-2674-0

Bachiller, E., Irigoien, X., 2015b. Trophodynamics and diet overlap of small pelagic fish species in the Bay of Biscay. Mar. Ecol. Prog. Ser. 534, 179-198. doi:10.3354/meps11375

Bivand, R., Keitt, T., Rowlingson, B., 2016. rgdal: Bindings for the Geospatial Data Abstraction Library, R package.

Borja, A., Elliott, M., Andersen, J.H., Cardoso, A.C., Carstensen, J., Ferreira, J.G., Heiskanen, A.-S., Marques, J.C., Neto, J.M., Teixeira, H., Uusitalo, L., Uyarra, M.C., Zampoukas, N., 2013. Good Environmental Status of marine ecosystems: What is it and how do we know when we have attained it? Mar. Pollut. Bull. 76, 16-27.
doi:10.1016/j.marpolbul.2013.08.042

Box, G., Jenkins, G., Reinsel, G., 1994. Time series analysis : forecasting and control, 3rd ed.
Prentice-Hall International, Upper Saddle River N.J.

Buckland, S.T., Rexstad, E.A., Marques, T.A., Oedekoven, C.S., 2015. Distance Sampling: Methods and Applications, Methods in Statistical Ecology. Springer International Publishing, Cham. doi:10.1007/978-3-319-19219-2

Cassou, C., Minvielle, M., Terray, L., Périgaud, C., 2011. A statistical-dynamical scheme for reconstructing ocean forcing in the Atlantic. Part I: weather regimes as predictors for ocean surface variables. Clim. Dyn. 36, 19-39. doi:10.1007/s00382-010-0781-7
Cassou, C., Terray, L., Hurrell, J.W., Deser, C., 2004. North Atlantic Winter Climate Regimes: Spatial Asymmetry, Stationarity with Time, and Oceanic Forcing. J. Clim. 17, 1055–1068.
doi:10.1175/1520-0442(2004)017<1055:NAWCRS>2.0.CO;2

Chust, G., Borja, A., Caballero, A., Irigoien, X., Sáenz, J., Moncho, R., Marcos, M., Liria, P., Hidalgo, J., Valle, M., Valencia, V., 2011. Climate change impacts on coastal and pelagic environments in the southeastern Bay of Biscay. Clim. Res. 48, 307–332.
doi:10.3354/cr00914

COM, 399 final, 2009. Proposal for a COUNCIL REGULATION establishing a long term plan for the anchovy stock in the Bay of Biscay and the fisheries exploiting that stock (SEC(2009) No. 1076 final).

de Castro, M., Gómez-Gesteira, M., Alvarez, I., Gesteira, J.L.G., 2009. Present warming within the context of cooling–warming cycles observed since 1854 in the Bay of Biscay. Cont. Shelf Res. 29, 1053–1059. doi:10.1016/j.csr.2008.11.016

Dessier, A., Bustamante, P., Chouvelon, T., Huret, M., Pagano, M., Marquis, E., Rousseaux, F., Pignon-Mussaud, C., Mornet, F., Bréret, M., Dupuy, C., this volume. The spring mesozooplankton variability and its relationship with hydrobiological structure over year-to-year changes (2003–2013) in the southern Bay of Biscay (North-East Atlantic). Prog. Oceanogr.

Doray, M., Badts, V., Massé, J., Duhamel, E., Huret, M., Doremus, G., Petitgas, P., 2014. Manual of fisheries survey protocols. PELGAS surveys (PELagiques GAScogne). Manuel des protocoles de campagne halieutique Ifremer No. 30259.

Doray, M., Duhamel, E., Huret, M., Petitgas, P., Massé, J., 2000. PELGAS, http://dx.doi.org/10.18142/18

Doray, M., Petitgas, P., Romagnan, J.B., Huret, M., Duhamel, E., Dupuy, C., Spitz, J., Authier, M., Sanchez, F., Berger, L., Dorémus, G., Bourriau, P., Grellier, P., Massé, J., this volume-a. The PELGAS survey: ship-based integrated monitoring of the Bay of Biscay pelagic
ecosystem. Prog. Oceanogr.

Doray, M., Petitgas, P., Huret, M., this volume-b. Spring habitats of small pelagic fish communities
in the Bay of Biscay. Prog. Oceanogr.

Doray, M., Petitgas, P., Saraux, C., Cornou, A.S., 2016. EchoR: R package for computing indices of the state of fish population and communities, based on fisheries acoustic data, R package.
https://forge.ifremer.fr/plugins/mediawiki/wiki/echor/index.php/Accueil

Erzini, K., 2005. Trends in NE Atlantic landings (southern Portugal): identifying the relative importance of fisheries and environmental variables. Fish. Oceanogr. 14, 195-209.
doi:10.1111/j.1365-2419.2005.00332.x

Erzini, K., Inejih, C., Stobberup, K., 2005. An application of two techniques for the analysis of short, multivariate non-stationary time-series of Mauritanian trawl survey data. ICES J. Mar. Sci. 62, 353-359. doi:10.1016/j.icesjms.2004.12.009

Fisher, J., Frank, K., 2004. Abundance-distribution relationships and conservation of exploited marine fishes. Mar. Ecol. Prog. Ser. 279, 201-213. doi:10.3354/meps279201

Folland, C.K., Knight, J., Linderholm, H.W., Fereday, D., Ineson, S., Hurrell, J.W., 2009. The Summer North Atlantic Oscillation: Past, Present, and Future. J. Clim. 22, 1082-1103.
doi:10.1175/2008JCLI2459.1

Foote, K.G., Knudsen, H.P., Vestnes, G., MacLennan, D.N., Simmonds, E.J., 1987. Calibration of acoustic instruments for fish density estimation: a practical guide. ICES Coop Res Rep 144, 57.

Garcia, S., 2000. Sustainability reference systems and indicators for responsible marine capture fisheries: a review of concepts and elements for a set of guidelines. Mar. Freshw. Res. 51, 385-426.

Gastauer, S., Fässler, S.M.M., O'Donnell, C., Høines, Å., Jakobsen, J.A., Krysov, A.I., Smith, L., Tangen, Ø., Anthonypillai, V., Mortensen, E., Armstrong, E., Schaber, M., Scoulding, B., 2016. The distribution of blue whiting west of the British Isles and Ireland. Fish. Res. 183,
32-43. doi:10.1016/j.fishres.2016.05.012

Gohin, F., Loyer, S., Lunven, M., Labry, C., Froidefond, J.-M., Delmas, D., Huret, M., Herbland, A., 2005. Satellite-derived parameters for biological modelling in coastal waters: Illustration over the eastern continental shelf of the Bay of Biscay. Remote Sens. Environ. 95, 29-46.
doi:10.1016/j.rse.2004.11.007

Greenstreet, S., Hall, S., 1996. Fishing and the ground-fish assemblage structure in the north-western North Sea: an analysis of long-term and spatial trends. J. Anim. Ecol. 65, 577-598.

Haig, B.D., 2003. What Is a Spurious Correlation? Underst. Stat. 2, 125-132.
doi:10.1207/S15328031US0202_03

Hijmans, R.J., 2016. raster: Geographic Data Analysis and Modeling, R package.

Huret, M., Bourriau, P., Doray, M., Gohin, F., Petitgas, P., this volume. Survey timing vs. ecosystem scheduling: degree-days to underpin observed variability with application to anchovy and sardine spawning in the Bay of Biscay. Prog. Oceanogr.

Ibaibarriaga, L., Fernandez, C., Uriarte, A., 2011. Gaining information from commercial catch for a Bayesian two-stage biomass dynamic model: application to Bay of Biscay anchovy. ICES J. Mar. Sci. 68, 1435-1446. doi:10.1093/icesjms/fsr094

ICES, 2016a. ICES Advice on fishing opportunities, catch, and effort. Bay of Biscay and the Iberian Coast and the Oceanic Northeast Atlantic ecoregions.

ICES, 2016b. Report of the Working Group on Widely Distributed Stocks (WGWIDE), 31 August-6 September 2016, ICES HQ, Copenhagen, Denmark. ICES CM 500.

ICES, 2016c. Bay of Biscay and the Iberian Coast Ecoregion – Ecosystem overview. ICES Ecosyst. Overv.

ICES, 2012. ICES scientific and advisory services of relevance to the EU Marine Strategy Framework Directive. ICES position paper.

ICES, 2010. Life cycle spatial patterns of small pelagic fish in the Northeast Atlantic. ICES Coop. Res. Rep. 93.
Irigoien, X., Fernandes, J.A., Grosjean, P., Denis, K., Albaina, A., Santos, M., 2009. Spring zooplankton distribution in the Bay of Biscay from 1998 to 2006 in relation with anchovy recruitment. J. Plankton Res. 31, 1–17. doi:10.1093/plankt/fbn096

Iversen, S.A., D. Skogen, M., Svendsen, E., 2002. Availability of horse mackerel (*Trachurus trachurus*) in the north-eastern North Sea, predicted by the transport of Atlantic water. Fish. Oceanogr. 11, 245–250. doi:10.1046/j.1365-2419.2002.00204.x

Koutsikopoulos, C., Le Cann, B., 1996. Physical processes and hydrological structures related to the Bay of Biscay anchovy. Sci. Mar. 60, 9–19.

Kupschus, S., Schratzberger, M., Righton, D., 2016. Practical implementation of ecosystem monitoring for the ecosystem approach to management. J. Appl. Ecol. 53, 1236–1247.
doi:10.1111/1365-2664.12648

Labry, C., Herbland, A., Delmas, D., Laborde, P., Lazure, P., Froidefond, J.M., Jegou, A.M., Sautour, B., 2001. Initiation of winter phytoplankton blooms within the Gironde plume waters in the Bay of Biscay. Mar. Ecol. Prog. Ser. 212, 117–130. doi:10.3354/meps212117

Lassalle, G., Lobry, J., Le Loc’h, F., Bustamante, P., Certain, G., Delmas, D., Dupuy, C., Hily, C., Labry, C., Le Pape, O., Marquis, E., Petitgas, P., Pusineri, C., Ridoux, V., Spitz, J., Niquil, N., 2011. Lower trophic levels and detrital biomass control the Bay of Biscay continental shelf food web: Implications for ecosystem management. Prog. Oceanogr. 91, 561–575.
doi:10.1016/j-pocean.2011.09.002

Levin, S., 1992. The problem of pattern and scale in ecology. Ecology 73, 1943–1967.

Link, J.S., Brodziak, J.K., Edwards, S.F., Overholtz, W.J., Mountain, D., Jossi, J.W., Smith, T.D., Fogarty, M.J., 2002. Marine ecosystem assessment in a fisheries management context. Can. J. Fish. Aquat. Sci. 59, 1429–1440. doi:10.1139/f02-115

Livingston, P., Aydin, K., Boldt, J., Ianelli, J., Juradomolina, J., 2005. A framework for ecosystem impacts assessment using an indicator approach. ICES J. Mar. Sci. 62, 592–597.
doi:10.1016/j.icesjms.2004.12.016
MacCall, A.D., 1990. Dynamic geography of marine fish populations. University of Washington Press, Seattle.

Manderson, J., Palamara, L., Kohut, J., Oliver, M., 2011. Ocean observatory data are useful for regional habitat modeling of species with different vertical habitat preferences. Mar. Ecol. Prog. Ser. 438, 1–17. doi:10.3354/meps09308

Mesnil, B., Petitgas, P., 2009. Detection of changes in time-series of indicators using CUSUM control charts. Aquat Living Resour 22, 187–192.

Michel, S., Treguier, A.-M., Vandermeirsch, F., 2009. Temperature variability in the Bay of Biscay during the past 40 years, from an in situ analysis and a 3D global simulation. Cont. Shelf Res. 29, 1070–1087. doi:10.1016/j.csr.2008.11.019

Ogden, T., 2002. Wavelet Methods for Time Series Analysis. J. Am. Stat. Assoc. 97, 362–363. doi:10.1198/jasa.2002.s460

Pante, E., Simon-Bouhet, B., 2013. marmap: A Package for Importing, Plotting and Analyzing Bathymetric and Topographic Data in R. PLoS ONE 8, e73051. doi:10.1371/journal.pone.0073051

Perpinan Lamigueiro, O., Hijmans, R.J., 2016. rasterVis, R package version. meteoForecast.

Perrot, L., Lampert, L., Ruiz-Pino, D., Gohin, F., Dessier, A., Bourriau, P., this volume. Spatial and vertical variability of coccolithophore blooms in the Bay of Biscay by satellite SPM and in situ turbidity from PELGAS cruises. Prog. Oceanogr.

Perry, A., Low, P., Ellis, J., Reynolds, J., 2005. Climate change and distribution shifts in marine fisheries. Science 308, 1912–1915.

Petitgas, P., Doray, M., Huret, M., Massé, J., Woillez, M., 2014. Modelling the variability in fish spatial distributions over time with empirical orthogonal functions: anchovy in the Bay of Biscay. ICES J. Mar. Sci. 71, 2379–2389. doi:10.1093/icesjms/fsu111

Petitgas, P., Goarant, A., Massé, J., Bourriau, P., 2009. Combining acoustic and CUFES data for the quality control of fish-stock survey estimates. ICES J Mar Sci 66, 1384–1390.
Petitgas, P., Secor, D.H., McQuinn, I., Huse, G., Lo, N., 2010. Stock collapses and their recovery: mechanisms that establish and maintain life-cycle closure in space and time. ICES J. Mar. Sci. 67, 1841–1848. doi:10.1093/icesjms/fsq082

Pingree, R.D., Garcia-Soto, C., 2014. Plankton blooms, ocean circulation and the European slope current: Response to weather and climate in the Bay of Biscay and W English Channel (NE Atlantic). Deep Sea Res. Part II Top. Stud. Oceanogr. 106, 5–22.
doi:10.1016/j.dsr2.2014.07.008

Pitois, S., Fox, C., 2006. Long-term changes in zooplankton biomass concentration and mean size over the Northwest European shelf inferred from Continuous Plankton Recorder data. ICES J. Mar. Sci. 63, 785–798. doi:10.1016/j.icesjms.2006.03.009

Planque, B., Bellier, E., Lazure, P., 2007. Modelling potential spawning habitat of sardine (*Sardina pilchardus*) and anchovy (*Engraulis encrasicolus*) in the Bay of Biscay. Fish Ocean. 16, 16–30.

Politikos, D.V., Huret, M., Petitgas, P., 2015. A coupled movement and bioenergetics model to explore the spawning migration of anchovy in the Bay of Biscay. Ecol. Model. 313, 212–222. doi:10.1016/j.ecolmodel.2015.06.036

Priestley, M., 1981. Spectral Analysis and Time Series Analysis. Academic Press, London.

R Core Team, 2016. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria.

Reuchlin-Hugenholtz, E., Shackell, N.L., Hutchings, J.A., 2015. The Potential for Spatial Distribution Indices to Signal Thresholds in Marine Fish Biomass. PLOS ONE 10, e0120500. doi:10.1371/journal.pone.0120500

Rice, J.C., Rivard, D., 2007. The dual role of indicators in optimal fisheries management strategies. ICES J. Mar. Sci. 64, 775–778. doi:10.1093/icesjms/fsm033

Shapiro, D.E., Switzer, P., 1989. Extracting time trends from multiple monitoring sites. Dep. Stat. Stanf. Univ. Tech Rep.
Shephard, S., van Hal, R., de Boois, I., Birchenough, S.N.R., Foden, J., O'Connor, J., Geelhoed, S.C.V., Van Hoey, G., Marco-Rius, F., Reid, D.G., Schaber, M., 2015. Making progress towards integration of existing sampling activities to establish Joint Monitoring Programmes in support of the MSFD. Mar. Policy 59, 105–111. doi:10.1016/j.marpol.2015.06.004

Shewhart, W.A., 1931. Economic control of quality of manufactured product. D. Van Nostrand Company, Inc., New York.

Spearman, C., 1904. The Proof and Measurement of Association between Two Things. Am. J. Psychol. 15, 72. doi:10.2307/1412159

Steinacher, M., Joos, F., Frölicher, T.L., Bopp, L., Cadule, P., Doney, S.C., Gehlen, M., Schneider, B., Segschneider, J., 2009. Projected 21st century decrease in marine productivity: a multi-model analysis. Biogeosciences Discuss. 6, 7933–7981. doi:10.5194/bgd-6-7933-2009

Switzer, P., Green, A.A., 1984. Min/Max Autocorrelation Factors for Multivariate Spatial Imagery. Dep. Stat. Stanf. Univ. Tech Rep 16.

Thomas, Y., Dumas, F., Andréfouët, S., 2014. Larval Dispersal Modeling of Pearl Oyster Pinctada margaritifera following Realistic Environmental and Biological Forcing in Ahe Atoll Lagoon. PLoS ONE 9, e95050. doi:10.1371/journal.pone.0095050

Trigo, R., Osborn, T., Corte-Real, J., 2002. The North Atlantic Oscillation influence on Europe: climate impacts and associated physical mechanisms. Clim. Res. 20, 9–17. doi:10.3354/cr020009

Valdés, L., López-Urrutia, A., Cabal, J., Alvarez-Ossorio, M., Bode, A., Miranda, A., Cabanas, M., Huskin, I., Anadón, R., Alvarez-Marqués, F., Llope, M., Rodríguez, N., 2007. A decade of sampling in the Bay of Biscay: What are the zooplankton time series telling us? Prog. Oceanogr. 74, 98–114. doi:10.1016/j-pocean.2007.04.016

Van Beveren, E., Bonhommeau, S., Fromentin, J.-M., Bigot, J.-L., Bourdeix, J.-H., Brosset, P., Roos, D., Saraux, C., 2014. Rapid changes in growth, condition, size and age of small pelagic fish in the Mediterranean. Mar. Biol. 161, 1809–1822. doi:10.1007/s00227-014-
Villamor, B., Abaunza, P., Petitgas, P., Massé, J., Uriarte, A., 2007. Bay of Biscay anchovy: a resource in crisis. ICES Insight 9–10.

Wallace, J.M., Gutzler, D.S., 1981. Teleconnections in the Geopotential Height Field during the Northern Hemisphere Winter. Mon. Weather Rev. 109, 784–812. doi:10.1175/1520-0493(1981)109<0784:TITGHF>2.0.CO;2

Woillez, M., 2013. Spatial indices package.

Woillez, M., Rivoirard, J., Petitgas, P., 2009a. Using min/max autocorrelation factors of survey-based indicators to follow the evolution of fish stocks in time. Aquat Living Resour 22, 193–200.

Woillez, M., Petitgas, P., Huret, M., Struski, C., Léger, F., 2010. Statistical monitoring of spatial patterns of environmental indices for integrated ecosystem assessment: Application to the Bay of Biscay pelagic zone. Prog. Oceanogr. 87, 83–93.

Woillez, M., Poulard, J.C., Rivoirard, J., Petitgas, P., Bez, N., 2007. Indices for capturing spatial patterns and their evolution in time, with application to European hake (Merluccius merluccius) in the Bay of Biscay. ICES J Mar Sci 64, 537–550.

Woillez, M., Rivoirard, J., Petitgas, P., 2009b. Notes on survey-based spatial indicators for monitoring fish populations. Aquat Living Resour 22, 155–164.

Zuur, A.F., Fryer, R.J., Jolliffe, I.T., Dekker, R., Beukema, J.J., 2003a. Estimating common trends in multivariate time series using dynamic factor analysis. Environmetrics 14, 665–685. doi:10.1002/env.611

Zuur, A.F., Pierce, G.J., 2004. Common trends in northeast Atlantic squid time series. J. Sea Res. 52, 57–72. doi:10.1016/j.seares.2003.08.008

Zuur, A.F., Tuck, I.D., Bailey, N., 2003b. Dynamic factor analysis to estimate common trends in fisheries time series. Can. J. Fish. Aquat. Sci. 60, 542–552. doi:10.1139/f03-030
Tables

<table><caption>Table1. Parameters and potential indicators derived from the PELGAS survey data to assess the temporal dynamics of the Bay of Biscay pelagic</caption><thead><tr><th>Sampling tool</th><th>Sampling time coverage</th><th>Sampling spatial coverage</th><th>Sampling protocol reference</th><th>Parameter (unit)</th><th>Parameter derivation reference</th><th>Indicator (unit or range)</th><th>Code</th><th>Indicator time coverage</th></tr></thead><tbody><tr><td rowspan="5">Conductivity-Temperature-Depth probe</td><td rowspan="5">2000-2015</td><td rowspan="5">Hydrobiology station</td><td rowspan="5">Doray et al., 2014</td><td>Bottom temperature (°C)</td><td rowspan="5">Huret et al., 2013</td><td>Annual gridded map average (°C)</td><td>BotTemp.mean</td><td rowspan="5">2000-2015</td></tr><tr><td>Surface temperature (°C)</td><td>Annual gridded map coefficient of variation (%)</td><td>BotTemp.CV</td></tr><tr><td>Equivalent freshwater height (EFE, m)</td><td>Annual gridded map average (°C)</td><td>SurfTemp.mean</td></tr><tr><td>Deficit of potential energy (DEP, kg·m<sup>-1</sup>·s<sup>-2</sup>)</td><td>Annual gridded map coefficient of variation (%)</td><td>SurfTemp.CV</td></tr><tr><td>Annual gridded map average (m)</td><td>Annual gridded map coefficient of variation (%)</td><td>Heq.mean</td></tr><tr><td rowspan="5">Niskin bottle</td><td rowspan="5">2003-2015</td><td rowspan="5">Hydrobiology station</td><td rowspan="5">Aminot and Kérouel, 2005</td><td rowspan="5">Vertically integrated Chl-a biomass (mg.m<sup>-2</sup>)</td><td rowspan="5">Dessier et al., this volume</td><td>Annual gridded map average (kg·m<sup>-1</sup>·s<sup>-2</sup>)</td><td>DefEpot.mean</td><td rowspan="5">2003-2015</td></tr><tr><td>Annual gridded map coefficient of variation (%)</td><td>DefEpot.CV</td></tr><tr><td>Annual gridded map average (mg.m<sup>-2</sup>)</td><td>Chloroint.mean</td></tr><tr><td>Annual gridded map center of gravity coordinates (decimal degrees)</td><td>Chloroint.xcg/ycg</td></tr><tr><td>Annual gridded map inertia (square nautical miles)</td><td>Chloroint.I</td></tr><tr><td rowspan="6">WP2 net</td><td rowspan="6">2006-2015</td><td rowspan="6">Hydrobiology station</td><td rowspan="6">Doray et al., 2014</td><td rowspan="6">Vertically integrated mesozooplankton biomass (mg.m<sup>-2</sup>)</td><td rowspan="6">Dessier et al., this volume</td><td>Annual gridded map isotropy ([0,1])</td><td>Chloroint.Iso</td><td rowspan="6">2009-2015</td></tr><tr><td>Annual gridded map positive area (square nautical miles)</td><td>Chloroint.PA</td></tr><tr><td>Annual gridded map equivalent area (square nautical miles)</td><td>Chloroint.EA</td></tr><tr><td>Annual gridded map microstructure index ([0,1])</td><td>Chloroint.microS</td></tr><tr><td>Annual gridded map average (mg.m<sup>-2</sup>)</td><td>Pmesotot.mean</td></tr><tr><td>Annual gridded map center of gravity coordinates (decimal degrees)</td><td>Pmesotot.xcg/ycg</td></tr><tr><td rowspan="6">Fisheries acoustics and midwater trawling</td><td rowspan="6">2000-2015</td><td rowspan="6">1852 m transect segment</td><td rowspan="6">Doray et al., 2014</td><td rowspan="6">Biomass (metric ton) per species (Scomber colias, Scomber scombrus, Trachurus trachurus, Micromesistius poutassou, Engraulis encrasicolus, Sardina pilchardus, Sprattus sprattus)</td><td rowspan="6">Doray et al., 2010</td><td>Annual gridded map inertia (square nautical miles)</td><td>Pmesotot.I</td><td rowspan="6">2006-2015</td></tr><tr><td>Annual gridded map isotropy ([0,1])</td><td>Pmesotot.Iso</td></tr><tr><td>Annual gridded map positive area (square nautical miles)</td><td>Pmesotot.PA</td></tr><tr><td>Annual gridded map equivalent area (square nautical miles)</td><td>Pmesotot.EA</td></tr><tr><td>Annual gridded map microstructure index ([0,1])</td><td>Pmesotot.microS</td></tr><tr><td>Annual total estimate (tons)</td><td>spCode.wbiom</td></tr><tr><td rowspan="6">Sighting</td><td rowspan="6">2004-2014</td><td rowspan="6">2000 m transect segment</td><td rowspan="6">Doray et al., 2014</td><td rowspan="6">Density (No. of individual per km<sup>2</sup>) per species groups (auks, gannets, small delphinids, big delphinids)</td><td rowspan="6">Authier et al., this volume</td><td>Annual gridded map center of gravity coordinates (decimal degrees)</td><td>spCode.xcg/ycg</td><td rowspan="6">2000-2015</td></tr><tr><td>Annual gridded map inertia (square nautical miles)</td><td>spCode.I</td></tr><tr><td>Annual gridded map isotropy ([0,1])</td><td>spCode.Iso</td></tr><tr><td>Annual gridded map positive area (square nautical miles)</td><td>spCode.PA</td></tr><tr><td>Annual gridded map equivalent area (square nautical miles)</td><td>spCode.EA</td></tr><tr><td>Annual gridded map microstructure index ([0,1])</td><td>spCode.microS</td></tr><tr><td></td><td></td><td></td><td></td><td>Length percentiles (cm) per species (Scomber colias, Scomber scombrus, Trachurus trachurus, Micromesistius poutassou, Engraulis encrasicolus, Sardina pilchardus, Sprattus sprattus)</td><td></td><td>Length 25<sup>th</sup>, 50<sup>th</sup>, 75<sup>th</sup> percentiles (cm)</td><td>spCode.Lo.0.25,0.5,0.75</td><td rowspan="2">2004-2014</td></tr><tr><td></td><td></td><td></td><td></td><td>Mean weight-at-age (grams) per species (Engraulis encrasicolus, Sardina pilchardus)</td><td></td><td>Mean weight (g) at age 1 to 5 (g)</td><td>spCode.mW.age1-5</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>Annual gridded map average (No. of individual per km<sup>2</sup>)</td><td>gspCode.mean</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>Annual gridded map center of gravity coordinates (decimal degrees)</td><td>gspCode.xcg/ycg</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>Annual gridded map inertia (square nautical miles)</td><td>gspCode.I</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>Annual gridded map isotropy ([0,1])</td><td>gspCode.Iso</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>Annual gridded map positive area (square nautical miles)</td><td>gspCode.PA</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>Annual gridded map equivalent area (square nautical miles)</td><td>gspCode.EA</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>Annual gridded map microstructure index ([0,1])</td><td>gspCode.microS</td><td></td></tr></tbody></table>
ecosystem. spCode: small pelagic fish species codes: SCOM-JAP: chub mackerel (*Scomber colias*); MICR-POU: blue whiting (*Micromesistius poutassou*); ENGR-ENC: European anchovy (*Engraulis encrasicolus*); SARD-PIL: European sardine, (*Sardina pilchardus*); SCOM-SCO: Atlantic mackerel (*Scomber scombrus*); SPRA-SPR: sprat (*Sprattus sprattus*); TRAC-TRU: horse mackerel (*Trachurus trachurus*). gspCode: megafauna species groups codes: AUKS: auks, GANNETS: gannets, SML.DEL: small delphinids, BIG.DEL: big delphinids.
<table><caption>Table 2. External forcing variables selected to investigate the potential effects of climate, survey timing, river discharge and fishing on the PELGAS survey indicators. Months were grouped into trimesters relative to the PELGAS survey average timing (May) so that trimester 1 (T1) goes from June to August, trimester 2 (T2) from September to November, trimester 3 (T3) from December to February, trimester 4 (T4) from March to May.</caption><thead><tr><th>Forcing</th><th>Variable</th><th>Data source</th></tr></thead><tbody><tr><td rowspan="6">Climate</td><td>NAO+ North Atlantic Winter Climate Regime occurrence (No. of days) during winter months (Dec. To Feb.) preceding PELGAS surveys</td><td rowspan="3">C. Cassou, CERFACS, Christophe.Cassou@cerfacs.fr</td></tr><tr><td>NAO- North Atlantic Winter Climate Regime occurrence (No. of days) during winter months (Dec. To Feb.) preceding PELGAS surveys</td></tr><tr><td>Ridge North Atlantic Winter Climate Regime occurrence (No. of days) during winter months (Dec. To Feb.) preceding PELGAS surveys</td></tr><tr><td>Scandinavian blocking North Atlantic Winter Climate Regime occurrence (No. of days) during winter months (Dec. To Feb.) preceding PELGAS surveys</td><td rowspan="3">National Oceans and Atmosphere Administration Climate Prediction Center, http://www.cpc.ncep.noaa.gov/</td></tr><tr><td>North Atlantic Oscillation index (NAO, unitless): june to may average and average per trimester</td></tr><tr><td>East Atlantic index (EA, unitless): june to may average and average per trimester</td></tr><tr><td>Scandinavia index (SCAND, unitless): june to may average and average per trimester</td><td></td></tr><tr><td rowspan="4">River discharge</td><td>Loire river discharge (m³.s⁻¹): average from june to may and by trimester, estimated at Nantes</td><td rowspan="4">Ministère de l'Écologie, du Développement Durable et de l'Energie, http://www.hydro.eaufrance.fr/</td></tr><tr><td>Vilaine river discharge (m³.s⁻¹): average from june to may and by trimester, estimated at Guichen</td></tr><tr><td>Garonne river discharge (m³.s⁻¹): average from june to may and by trimester, estimated at Tonneins</td></tr><tr><td>Adour river discharge (m³.s⁻¹): average from june to may and by trimester, estimated at Saint-Vincent-de-Paul</td></tr><tr><td></td><td>Loire, Vilaine, Garonne, Adour rivers total discharge (m³.s⁻¹): average from june to may and by trimester</td><td></td></tr><tr><td rowspan="10">Commercial fishing</td><td>Engraulis encrasicolus in ICES subarea 8 Stock Spawning Biomass (SSB, in tons)</td><td rowspan="5">ICES WGHANSA</td></tr><tr><td>Engraulis encrasicolus in ICES subarea 8 fishing mortality (F, unitless)</td></tr><tr><td>Engraulis encrasicolus in ICES subarea 8 total commercial catches (tons) in ICES</td></tr><tr><td>Sardina pilchardus in ICES subareas 7 and 8a,b,c total commercial catches (tons)</td></tr><tr><td>North-East Atlantic Scomber scombrus Stock Spawning Biomass (SSB, in tons)</td></tr><tr><td>North-East Atlantic Scomber scombrus fishing mortality (F, unitless)</td><td rowspan="6">ICES WGWIDE</td></tr><tr><td>North-East Atlantic Scomber scombrus total commercial catches (tons)</td></tr><tr><td>North-East Atlantic Trachurus trachurus Stock Spawning Biomass (SSB, in tons)</td></tr><tr><td>North-East Atlantic Trachurus trachurus fishing mortality (F, unitless)</td></tr><tr><td>North-East Atlantic Trachurus trachurus total commercial catches (tons)</td></tr><tr><td>Micromesistius poutassou in ICES subareas 1–9, 12, and 14 Spawning Biomass (SSB, in tons)</td></tr><tr><td></td><td></td></tr><tr><td></td><td>Micromesistius poutassou in ICES subareas 1–9, 12, and 14 fishing mortality (F, unitless)</td><td></td></tr><tr><td></td><td>Micromesistius poutassou in ICES subareas 1–9, 12, and 14 total commercial catches (tons)</td><td></td></tr><tr><td></td><td>Sardina pilchardus in ICES subareas 7 and 8a,b,c total commercial catches (tons)</td><td></td></tr><tr><td>Ecosystem seasonality</td><td>PELGAS survey corrected starting date</td><td>M. Huret, Ifremer, martin.huret@ifremer.fr</td></tr></tbody></table>
<table><caption>Table 3. PELGAS survey-derived Indicators selected to represent the Bay of Biscay pelagic ecosystem components. Indicators selected to represent the hydrology and fish components combination are shaded in grey.</caption><thead><tr><th>Compartment</th><th>Indicator (unit or range)</th><th>Time period</th></tr></thead><tbody><tr><td rowspan="3">Hydrology</td><td>Bottom temperature gridded map average (°C)</td><td rowspan="3">2000-2015</td></tr><tr><td>Equivalent freshwater height coefficient of variation (CV, %)</td></tr><tr><td>Equivalent freshwater height gridded map average (Heq, m)</td></tr><tr><td rowspan="3">Phytoplankton</td><td>Vertically integrated Chl-a biomass gridded map average (mg.m-2)</td><td>2003-2015</td></tr><tr><td>Vertically integrated Chl-a biomass gridded map equivalent area (square nautical miles)</td><td rowspan="2">2009-2015</td></tr><tr><td>Vertically integrated Chl-a biomass gridded map microstructure index ([0,1])</td></tr><tr><td rowspan="3">Mesozooplankton</td><td>Vertically integrated mesozooplankton biomass gridded map average (mg.m-2)</td><td>2006-2015</td></tr><tr><td>Vertically integrated mesozooplankton biomass gridded map positive area (square nautical miles)</td><td rowspan="2">2009-2015</td></tr><tr><td>Vertically integrated mesozooplankton biomass gridded map equivalent area (square nautical miles)</td></tr><tr><td rowspan="10">Small pelagic fish</td><td>Sardina pilchardus mean weight (g) at age 1</td><td rowspan="10">2000-2015</td></tr><tr><td>Sardina pilchardus mean weight (g) at age 2</td></tr><tr><td>Sardina pilchardus gridded map center of gravity longitude (decimal degrees)</td></tr><tr><td>Sardina pilchardus length 25th percentile (cm)</td></tr><tr><td>Trachurus trachurus biomass estimate (tons)</td></tr><tr><td>Sprattus sprattus gridded map equivalent area (square nautical miles)</td></tr><tr><td>Sprattus sprattus gridded map isotropy ([0,1])</td></tr><tr><td>Engraulis encrasicolus biomass estimate (tons)</td></tr><tr><td>Engraulis encrasicolus gridded map isotropy ([0,1])</td></tr><tr><td>Engraulis encrasicolus mean weight (g) at age 1</td></tr><tr><td>Engraulis encrasicolus mean weight (g) at age 2</td></tr><tr><td>Engraulis encrasicolus gridded map positive area (square nautical miles)</td></tr><tr><td rowspan="5">Megafauna</td><td>Gannets gridded map center of gravity longitude (decimal degrees)</td><td rowspan="5">2004-2014</td></tr><tr><td>Gannets gridded map inertia (square nautical miles)</td></tr><tr><td>Gannets gridded map equivalent area (square nautical miles)</td></tr><tr><td>Auks gridded map positive area (square nautical miles)</td></tr><tr><td>Small delphinids gridded map microstructure index ([0,1])</td></tr></tbody></table>
<table><caption>Table 4. Pairs of selected PELGAS survey indicators with absolute Spearman's rho correlation coefficient higher or equal to 0.6. Relationships considered as ecologically plausible are shaded in grey. Indicators codes can be found in Table 1.</caption><thead><tr><td>Indicator 1</td><td>Indicator 2</td><td>Spearman's rho correlation coefficient</td></tr></thead><tbody><tr><td>Heq.meanZ</td><td>Heq.CV</td><td>-0.74</td></tr><tr><td>ENGR.ENC.Iso</td><td>ENGR.ENC.wbiom</td><td>0.63</td></tr><tr><td>TRAC.TRU.wbiom</td><td>SARD.PIL.mW.age1</td><td>0.64</td></tr><tr><td>ENGR.ENC.mW.age1</td><td>SARD.PIL.mW.age1</td><td>0.72</td></tr><tr><td>SARD.PIL.mW.age1</td><td>SARD.PIL.Lo.0.25</td><td>0.69</td></tr></tbody></table>
<table><caption>Table 5. Pairs of selected PELGAS survey indicators and external forcing variables with absolute Spearman's rho correlation coefficient higher or equal to 0.6. Relationships considered as ecologically plausible are shaded in grey. Indicators codes can be found in Table 1.</caption><thead><tr><th>Indicator 1</th><th>Indicator 2</th><th>Spearman's rho correlation coefficient</th></tr></thead><tbody><tr><td>Heq.meanZ</td><td>Adour.discharge</td><td>0.62</td></tr><tr><td>Heq.meanZ</td><td>All.discharge</td><td>0.77</td></tr><tr><td>Heq.meanZ</td><td>Garonne.discharge</td><td>0.66</td></tr><tr><td>Heq.meanZ</td><td>Loire.discharge</td><td>0.77</td></tr><tr><td>Heq.meanZ</td><td>LoireT2.discharge</td><td>0.60</td></tr><tr><td>Heq.meanZ</td><td>VilaineT3.discharge</td><td>0.63</td></tr><tr><td>Heq.meanZ</td><td>AdourT4.discharge</td><td>0.75</td></tr><tr><td>Heq.meanZ</td><td>AllT4.discharge</td><td>0.61</td></tr><tr><td>Heq.meanZ</td><td>GaronneT4.discharge</td><td>0.70</td></tr><tr><td>ENGR.ENC.Iso</td><td>TRAC.TRU.TotalCatch</td><td>0.76</td></tr><tr><td>ENGR.ENC.Iso</td><td>SCOM.SCO.TotalCatch</td><td>0.69</td></tr><tr><td>TRAC.TRU.wbiom</td><td>naoT1.mean.index</td><td>0.72</td></tr><tr><td>TRAC.TRU.wbiom</td><td>TRAC.TRU.SSB</td><td>-0.62</td></tr><tr><td>TRAC.TRU.wbiom</td><td>SCOM.SCO.SSB</td><td>-0.75</td></tr><tr><td>TRAC.TRU.wbiom</td><td>MICR.POU.SSB</td><td>0.70</td></tr><tr><td>TRAC.TRU.wbiom</td><td>ENGR.ENC.F</td><td>0.70</td></tr><tr><td>TRAC.TRU.wbiom</td><td>SCOM.SCO.F</td><td>0.83</td></tr><tr><td>TRAC.TRU.wbiom</td><td>MICR.POU.F</td><td>0.74</td></tr><tr><td>TRAC.TRU.wbiom</td><td>MICR.POU.TotalCatch</td><td>0.67</td></tr><tr><td>ENGR.ENC.mW.age1</td><td>SCOM.SCO.SSB</td><td>-0.60</td></tr><tr><td>ENGR.ENC.mW.age1</td><td>MICR.POU.F</td><td>0.61</td></tr><tr><td>SARD.PIL.mW.age1</td><td>SCOM.SCO.SSB</td><td>-0.92</td></tr><tr><td>SARD.PIL.mW.age1</td><td>MICR.POU.SSB</td><td>0.65</td></tr><tr><td>SARD.PIL.mW.age1</td><td>SCOM.SCO.F</td><td>0.78</td></tr><tr><td>SARD.PIL.mW.age1</td><td>MICR.POU.F</td><td>0.80</td></tr><tr><td>SARD.PIL.mW.age1</td><td>MICR.POU.TotalCatch</td><td>0.69</td></tr><tr><td>SPRA.SPR.Iso</td><td>SCOM.SCO.SSB</td><td>-0.61</td></tr><tr><td>SPRA.SPR.Iso</td><td>SCOM.SCO.TotalCatch</td><td>-0.65</td></tr><tr><td>Heq.CV</td><td>eaT4.mean.index</td><td>-0.72</td></tr><tr><td>Heq.CV</td><td>AllT2.discharge</td><td>-0.69</td></tr><tr><td>Heq.CV</td><td>LoireT2.discharge</td><td>-0.78</td></tr><tr><td>Heq.CV</td><td>VilaineT2.discharge</td><td>-0.67</td></tr><tr><td>ENGR.ENC.wbiom</td><td>ENGR.ENC.SSB</td><td>0.65</td></tr><tr><td>ENGR.ENC.wbiom</td><td>MICR.POU.SSB</td><td>-0.70</td></tr><tr><td>ENGR.ENC.wbiom</td><td>MICR.POU.F</td><td>-0.70</td></tr><tr><td>ENGR.ENC.wbiom</td><td>SCOM.SCO.TotalCatch</td><td>0.69</td></tr><tr><td>ENGR.ENC.wbiom</td><td>MICR.POU.TotalCatch</td><td>-0.69</td></tr></tbody></table>
Figures

<figure><img src="image_2.png" /><figcaption>Figure 1. PELGAS survey sampling scheme. Solid lines: systematic line transects, red dots: hydrobiology stations. Light grey grid: block averaging grid within survey polygon. Light grey lines: 100, 200, 300, 400, 500 m isobaths..</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2. Schematic outline of the min/max autocorrelation factor (MAF)-based indicator selection procedure.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3. Schematic outline of the min/max autocorrelation factor (MAF)-based indicator selection strategy used to select the most representative indicators of the Bay of Biscay pelagic ecosystem dynamics, based on PELGAS survey data.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4. Hydrology indicators derived from the PELGAS survey: a) normalised one-lag variograms, b) continuity indices. Indicators on the left of the vertical dashed line on figure 3b were selected to represent the hydrology component: c) equivalent freshwater height gridded map average, d) bottom temperature gridded map average, e) bottom temperature gridded map coefficient of variation. Indicators codes can be found in Table 1.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5. Phytoplankton indicators derived from the PELGAS survey: a) normalised one-lag variograms, b) continuity indices, c) vertically integrated chlorophyl-a biomass gridded map average, 2003-2015. Indicators on the left of the vertical dashed line on Figure 5b were selected to represent the phytoplankton component over the 2009-2015 period: vertically integrated Chl-a biomass gridded map average (d), equivalent area (e), and microstructure index (f). Indicators codes can be found in Table 1.</figcaption></figure>
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Figure 6. Zooplankton indicators derived from the PELGAS survey: a) normalised one-lag variograms, b) continuity indices, c) vertically integrated mesozooplankton biomass gridded map average, 2006-2015. Indicators on the left of the vertical dashed line on Figure 6b were selected to represent the mesozooplankton component over the 2009-2015 period: vertically integrated mesozooplankton biomass gridded map equivalent area (d), positive area (e), and average (f). Indicators codes can be found in Table 1.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 7. Small pelagic fish indicators derived from the PELGAS survey: a) normalised one-lag variograms, b) continuity indices. Indicators on the left of the vertical dashed line on Figure 7b were selected to represent the small pelagic fish component over the 2000-2015 period: c) *Sardina pilchardus* mean weight at age 1, d) *S. pilchardus* mean weight at age 2, e) *S. pilchardus* gridded map center of gravity longitude, f) *S. pilchardus*</figcaption></figure>
length 25th percentile, g) *Trachurus trachurus* biomass estimate, h) *Sprattus sprattus* gridded map equivalent area, i) *S. sprattus* gridded map isotropy, j) *Engraulis encrasicolus* biomass estimate, k) *E. encrasicolus* gridded map isotropy, l) *E. encrasicolus* mean weight at age 1, m) *E. encrasicolus* mean weight at age 2, n) *E. encrasicolus* gridded map positive area. Indicators codes can be found in Table 1.
<figure><img src="image_14.png" /><figcaption>Figure 8. Megafauna indicators derived from the PELGAS survey: a) normalised one-lag variograms, b) continuity indices. Indicators on the left of the vertical dashed line on Figure 8b were selected to represent the megafauna component over the 2004-2014 period: c) Gannets gridded map center of gravity longitude, d) Gannets gridded map inertia, e) Gannets gridded map equivalent area, f) Auks gridded map positive area, g) Small delphinids gridded map microstructure index. Indicators codes can be found in Table 1.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 9. Hydrology and small pelagic fish indicators derived from the PELGAS survey: a) normalised one-lag variograms, b) continuity indices. Indicators on the left of the vertical dashed line on Figure 9b were selected to represent the hydrology and small pelagic fish components over the 2000-2015 period: equivalent freshwater height average (c) and CV (d), Engraulis encrasicolus biomass estimate (e), Engraulis encrasicolus mean weight at age 1 (f), Engraulis encrasicolus isotropy (g), Sardina pilchardus mean weight at age 1 (h), Sardina pilchardus length 25th percentile (i), Trachurus trachurus biomass estimate (j), Sprattus sprattus isotropy (k). Indicators codes can be found in Table 1.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 10. Summary of the first 2 min/max autocorrelation factors (MAFs) computed based on 13 PELGAS survey-derived indicators representing the hydrology and fish component: a) MAF1 time series, b) MAF1 variogram, c) indicators loadings on MAF1, d) MAF2 time series, e) MAF2 variogram, f) indicator loadings on MAF2. Indicators codes can be found in Table 1.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 11. Abundance (log(No. of fish+1)) at ages 1,2,3 and 4+ (sum of ages equal to or higher than 4) of anchovy (ENGR-ENC) and sardine (SARD-PIL) derived from the PELGAS survey data. A loess smooth (solid bold line) is shown for each species:age combination.</figcaption></figure>
Supplementary material

<figure><img src="image_18.png" /><figcaption>Figure S1. Annual, mean and standard deviation gridded maps of vertically integrated chlorophyll-a mean biomass (mg.m⁻2) derived from the PELGAS survey data. Grid cell size: 0.25°x0.25°.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure S2. Annual, mean and standard deviation gridded maps of vertically integrated mesozooplankton mean biomass (mg.m-2) derived from the PELGAS survey data. Grid cell size: 0.25°x0.25°.</figcaption></figure>