Environmental drivers of spawning and recruitment of anchovy in the Bay of Biscay

Abstract

To prevent future collapse events as observed in the mid-2000s, environmental and biological drivers of anchovy (Engraulis encrasicolus) spawning habitat variability in the Bay of Biscay need to be better understood and monitored. To this end, this study uses an extended time series of egg density observations, new high-resolution Essential Ocean Variables (EOVs), and complementary statistical and mechanistic modelling approaches. The relationship between the spawning habitat and the adult spawning stock biomass (SSB) is also analysed. A delta generalized additive model (delta-GAM) and a mechanistic SEAPODYM-Hs model were applied to evaluate the relative roles of temperature, prey availability, predation, spawning stock biomass (SSB), as well as species biology (seasonal reproduction) and behaviour (coastal attraction). Bottom temperature emerged as a more relevant indicator of spawning habitat suitability than sea surface temperature, possibly reflecting conditions during spawning or egg development. Zooplankton availability and more particularly micronekton predation, essentially by migrant mesopelagic species, were identified as key mechanisms shaping egg mortality and the spatial distribution of successful spawning between coastal and slope areas. However, the seasonality of the reproduction remained essential to account for the monthly spatial distribution, whereas SSB was the primary predictor of recruitment interannual variability. Low SSB combined with persistently unfavourable environmental conditions over the species’ lifespan (~3 years) increased the risk of population collapse as observed in the mid-2000s. The mechanistic model, despite its parsimonious parameterization, showed strong predictive skills comparable to that of the delta-GAM. The resulting index (Hs) provides a robust indicator of the suitability of environmental conditions for anchovy spawning habitat and can be predicted in real-time using the Copernicus marine service variables.

1 Introduction

Climate variability significantly affects marine ecosystems, leading to fluctuations in fish stocks (e.g. Lehodey et al., 2006), especially pelagic fishes (Trenkel et al., 2014). Small pelagic fishes like sardines and anchovies are known to show strong inter-annual and/or multi-decadal biomass variability (Schwartzlose et al., 1999; Ma et al., 2019). These species are highly sensitive to changes in environmental conditions such as sea temperature and primary and zooplankton productivity, which affect their distribution, reproduction, and abundance (Cushing and Harris, 1973; Cushing, 1975; Cushing, 1995; Rothschild, 2000; Peck et al., 2013; Koenigstein et al., 2022). Natural fluctuations in abundance combined with overfishing have resulted in the collapse of several small pelagic fish populations (Corten, 1986; Alheit and Niquen, 2004; Shelton and Mangel, 2011) requiring adaptive management approaches integrating sources of uncertainties to develop recovery strategies (Siple et al., 2021). Anticipating the future trajectories of fish populations requires a good understanding of the mechanisms at play and also a capacity to forecast their key environmental drivers (Tommasi et al., 2017; Ward et al., 2024). The anchovy population in the Bay of Biscay has experienced a stock collapse, which led to the closure of its fishery between 2006 and 2010 (ICES, 2019). This moratorium, and likely more favourable environmental conditions facilitated the recovery of this anchovy population. By 2010, the population had sufficiently recovered to allow a reopening of the fishery under strict management measures including restricting total allowable catch and monitoring programs.

The search for understanding recruitment mechanisms in fish populations has a long history, dating back over more than a century (Hjort, 1914). Scientists have sought to explain how environmental factors and biological processes influence the survival of larvae and then juvenile fish. Many hypotheses, such as mismatches between larval emergence and food availability, predation, shifts in spawning phenology, and changing ocean currents, have been proposed (e.g. Cushing and Harris, 1973; Rothschild, 2000; Irigoien et al., 2007; Marchal et al., 2025). However, these mechanisms are difficult to unravel due to limited data, complex and possibly non-linear dynamics in marine ecosystems. This complexity has made it challenging to establish clear, predictive models of fish recruitment, despite decades of research efforts (Hjort, 1914; Shelton and Mangel, 2011). Nevertheless, most studies dedicated to Bay of Biscay anchovy highlighted a favourable range of temperature for spawning, within which the productivity and food (zooplankton) would control the larval recruitment success and the survival of subsequent life stages. This was shown directly or indirectly using statistical models, with significant relationships emerging between fish abundance at different life stages and proxies of productivity such as upwelling and stratification indices, wind regimes and climate metrics (e.g. Motos et al., 1996; Planque et al., 2007; Irigoien et al., 2007; Borja et al., 2008; Erauskin-Extramania et al., 2019). The role of river flow was also highlighted and often included in these statistical models using sea surface salinity to reflect the effect of estuaries on the dynamics of anchovy (Engraulis sp.). In fact, the variability of river flow modifies coastal water temperature and salinity, from surface to bottom, nutrient inputs, as well as the mixing and productivity in the river plume, the extent of which also depends on wind dynamics (Lazure and Jégou, 1998; Lloret et al., 2004). Variability in the flows of the Gironde and Loire rivers has been positively correlated to mesozooplankton productivity and nutritional conditions of anchovy (Bergeron et al., 2010; Bergeron et al.,2013). Unfortunately, there are multiple correlations between various environmental variables that make the causal interpretation of statistical model results difficult.

In fisheries science, recruitment is generally defined as the first age or size at which fish become vulnerable to exploitation. In population dynamics, recruitment is more conveniently represented by the abundance of the first age group represented in the exploited biomass being modelled. When a model simulates the full life cycle of the species, recruits in the first age cohort may correspond to the larval stage. In all cases, variability in recruitment is strongly associated with spawning success, although density-dependent processes can weaken this relationship. Therefore, spawning biomass, together with environmental changes influence eggs abundance, as well as larval and juvenile recruitment. Bay of Biscay anchovy recruits to the fisheries at age one, while reproductive biomass consists mainly of individuals aged 2 and 3 years (Ibaibarriaga et al., 2008). The anchovy reproductive biomass therefore depends in every year on the fishing mortality exerted on the population in the previous 2 years. Change in spawning biomass and environmental conditions need to be analysed together. For instance, Irigoien et al. (2007) explored the link between mean zooplankton biomass and anchovy recruitment for the period 1998–2006 and found a negative correlation. However, their data indicate also that the lowest mean zooplankton biomass values in this period were observed in 1998 and 1999, just before the abrupt decline of the stock. It could therefore be hypothesised that (at least) two successive unfavourable years of zooplankton production combined with the exploitation of the population led to this abrupt decline. Several years were needed to rebuild the reproductive biomass despite higher zooplankton biomass after 2001, when anchovy reproductive biomass was at its lowest (Ibaibarriaga et al., 2011).

The objective of this study is to explain and analyse the spawning and subsequent recruitment variability of Bay of Biscay anchovy. We used high-resolution environmental variables obtained by combining global models with in situ and satellite observations. These variables include spatially explicit biomass of zooplankton (the main prey of anchovies after the larval stage) and of micronekton (predators of eggs and larvae). Zooplankton biomass was predicted by a model forced by satellite-derived (ocean colour) primary production and by physical variables (temperature and currents), calibrated primarily for global open ocean waters. We note that the same zooplankton product had been used previously by Menu et al. (2023), along with temperature. These authors showed that the zooplankton product explained a significant fraction of the observed regional differences in individual anchovy growth. However, it did not explain the observed temporal trend (decrease between 17 and 31% in the Bay of Biscay over the study period). Therefore, here we will adapt the open ocean zooplankton product to our application in the shelf waters of the Bay of Biscay.

Three types of fish spawning habitats can be considered (e.g. Planque et al., 2007). The potential spawning habitat is the broadest category, referring to areas where environmental conditions are suitable for spawning, whether or not spawning is currently observed and without the influence of fish population dynamics. The realized spawning habitat is where spawning occurs, within the potential spawning habitat, and is shaped by spawning biomass, i.e. the biomass of mature fish contributing to reproduction, and by the seasonal cycle of reproduction, if it exists. Many species, especially in the tropics, are opportunistic spawners (Kahar et al., 2023). Finally, the successful spawning habitat includes areas where spawning results in successful recruitment (i.e. offspring survive and contribute to the exploited biomass). The successful habitat is a subset of the realized habitat, and it can only be identified retrospectively, by using juvenile data or modelling approaches. In the Bay of Biscay, the realized spawning habitat can be inferred from the egg density sampled at sea during PELGAS surveys combined with data on the environmental conditions and on the distribution of spawning abundance and biomass.

In this study, we applied a statistical delta-GAM model to link anchovy egg density data as observed during the annual PELGAS (PELagiques GAScogne) surveys (Doray et al., 2018b) with explanatory variables describing environmental and reproductive biomass variability. The results inform the set-up of functional relationships that describe the spawning habitat in the SEAPODYM (Spatial Ecosystem And Population Dynamics Model) mechanistic model (Lehodey et al., 2008; Hernandez et al., 2014; Senina et al., 2020) using the same predictors. This model incorporates temperature, prey availability, and predation effects by micronekton, the latter being a rarely tested factor (Bailey and Houde, 1989; Pepin, 2015).

Finally, we compare both approaches and evaluate whether the spawning habitat index derived from SEAPODYM can describe past and recent fluctuations in the abundance and distribution of anchovy egg density. We discuss if this variability can be explained by the observed changes in the climatic and oceanographic conditions, as well as by the intrinsic biology and population dynamics of Bay of Biscay anchovy.

2 Materials and methods

2.1 PELGAS eggs density and zooplankton data

Since 2000, PELGAS research cruises have been conducted each spring, covering the Bay of Biscay continental shelf along the French coast, with the exception of the COVID year in 2020. These surveys include sampling of fish eggs and larvae, phytoplankton and zooplankton, as well as measurements of hydrological variables (Doray et al., 2018a, 2018b). The PELGAS data are available in a grid-averaged format with a spatial resolution of 0.25° × 0.25°. A specific processing procedure is applied to avoid averaging artefacts related to the position of the grid origin. The data are freely accessible in Doray et al. (2018a) for the period 2000–2016 for all variables (except zooplankton, which is available for 2003–2016). Egg data for the period 2017–2023 were made available in this work by Martin Huret and were processed here by using the same block-averaging method applied in the previous series.

The mesozooplankton data collected in the Bay of Biscay with the PELGAS annual cruises are available since 2003 (Grandremy et al., 2024). Zooplankton is sampled during the PELGAS survey using a WP2 net with 200 μm mesh size at approximately 80 night-stations, targeting the upper 100m of the water column. Each sample is fractionated into four size classes (200, 500, 1000 and 2000 μm) to determine the dry weight of each class as well as the total dry weight. Because vertically migrating zooplankton occupies the upper layer at night, these data are assumed to represent the total zooplankton biomass in the water column. Details of the methodology can be found in Doray et al. (2018a).

Anchovy egg density is measured by using the Continuous Underway Fish Egg Sampler (CUFES), which continuously pumps sea water from a depth of 3 m and collects approximately 10m³ of surface water over 3 nautical miles per sample. Sampling is performed during daylight hours, and the collected material is analysed to determine egg counts. Since 2015, samples have been analysed by using the ZooCAM Flow imager (Colas et al., 2018), which is an inflow imaging particle and plankton analyser; prior to this, analyses were performed by using a binocular microscope. The processed samples provide estimates of egg density, which do not represent the total number of eggs in the water column but offer a reliable indication of their spatial distribution.

2.2 Predictors

Monthly and annual Spawning Stock Biomass (SSB) estimated by the most recent ICES stock assessment group (ICES, 2024) have been included among the statistical model predictors, as well as six environmental predictors: sea surface temperature (SST), sea bottom temperature (SBT), bathymetry (BATHY), net primary production (NPP), zooplankton biomass (ZOO) and micronekton biomass (MICRO). These environmental predictors were obtained from the Copernicus Marine Service (CMEMS) catalogue and used directly or following a pre-treatment described in Table 1. They have been interpolated to fit the spatial resolution of the egg density dataset (0.25 degrees).

The mean zooplankton abundance observed during PELGAS cruises was compared with SEAPODYM model outputs at a spatial resolution of 0.25 degrees and a temporal resolution of one week. Predictions were substantially higher in coastal areas compared to the observed data, likely due to the parameterization and the model design, which focuses on the pelagic system and ignores the flux of one part of the energy to the benthic system. To address this bias, we applied a bathymetry-based correction (see Supplementary Materials) and validated the result with independent zooplankton observations. For micronekton, no observations of abundance are available in the Bay of Biscay, and we assume that SEAPODYM micronekton biomass is also overestimated in shallow waters. Therefore, the same bathymetry-based correction is applied (see Supplementary Materials.).

A cross-correlation analysis among the predictors was performed to assess if collinearities and redundancies would diminish the reliability of the statistical analysis (see Supplementary Materials.). We found a relatively high correlation score (r = 0.54) between ZOO and MICRO, with more moderate correlations between MICRO and SST (0.36), BATHY and NPP (0.34), and temperature and ZOO (0.26 for SST and 0.33 for SBT), and not surprisingly between the two temperature variables SST and SBT (0.26). All other r scores are very weak. These correlation levels are much below the value of 0.7 proposed by Dormann et al. (2013) to indicate problematic multicollinearity and therefore did not preclude the inclusion of all variables in the GAM.

2.3 Statistical analysis with delta-GAM

Statistical models were fitted to the data to identify the environmental variables best explaining realized annual spawning habitats, by using observed egg densities as dependent variables.

The egg density dataset contains a high proportion of zero observations, which must be properly addressed to account for zero inflation and overdispersion (Jensen et al., 2005; Grüss et al., 2014; Tanaka et al., 2017). One approach involves modelling only the positive observations to simplify the statistical fitting process—commonly using a Generalized Additive Model (GAM) that accommodates multiple predictor variables and their potential nonlinear effects. However, the observed absence of eggs also provides valuable ecological information. To address this, we use a delta-GAM approach, also known as a hurdle model or two-stage GAM, that explicitly separates the modelling of zeros and the modelling of positive values. Delta-GAM models have been widely applied in ecological studies involving zero-inflated continuous data (Welsh et al., 1996; Jensen et al., 2005; Grüss et al., 2014; Tanaka et al., 2017). This two-part approach first uses a binomial logistic GAM to estimate the probability of occurrence (presence/absence) and then applies a GAM to the positive values only to estimate conditional abundance. We selected the Gamma GAM because egg density values are continuous, strictly positive, non-Gaussian, and right-skewed (Wood, 2017). The final predicted abundance is calculated as the product of the predicted probability of presence and the predicted conditional abundance. All models are developed in python with the pygam library (Servén and Brummitt, 2018; see also code availability).

Several combinations of predictors are tested (cf. above) and the best models are ranked according to the significance of the variable and standard metrics: the explained variance (pseudo-R2, an analogous of R2 for the case of models where the outcome is not continuous and normally distributed), and the Corrected Akaike Information Criterion AICc (Bradley, 1997; Cavanaugh, 1997). The AICc is a model selection metric that balances model fit and complexity, especially for small sample sizes. Like AIC, it penalizes models with more parameters to avoid overfitting, but AICc adds a correction for finite sample sizes. Lower AICc values indicate a better model fit (relative to others tested). The significance of each model term is also verified, as well as the residuals.

Ideally, in addition to key environmental drivers that are assumed to strongly influence spawning success, the distribution and abundance of the fish actively participating in the reproduction should also be provided as predictor to expect a good fit to the data. Such detailed information is not available, but we used as a proxy the annual total stock spawning biomass (SSB). A final consideration concerns the representativeness of the data samples. In the Bay of Biscay, egg density sampling is conducted annually in May, based on prior evidence that the spawning season of the Bay of Biscay anchovy occurs predominantly between April and June (Bellier et al., 2007). This temporal window has been confirmed by monitoring the gonadosomatic index (Erauskin et al., 2019). To incorporate this prior knowledge into the statistical model, zero values –constrained to 10% of the annual sample size– are randomly assigned to the months from September to March, based on the existing list of latitudes and longitudes of the observations in the dataset.

2.4 SEAPODYM spawning habitat model

The Spatial Ecosystem And Population Dynamics Model (SEAPODYM) is a numerical framework for studying fish population dynamics under environmental and fishing pressures. Using ocean model outputs (temperature, currents), satellite data, or biogeochemical models (primary production, phytoplankton), SEAPODYM simulates zooplankton and micronekton in the epipelagic and mesopelagic layers (Lehodey et al., 2010). It then models commercial species and fisheries (Lehodey et al., 2008; Senina et al., 2020), including movement and life stages governed by behavioural rules linked to spawning and feeding habitats, shaped by environmental conditions and species’ preferences. Population dynamics are described with advection–diffusion–reaction equations, which account for movement, natural and fishing mortality, growth, and recruitment. This formulation reduces the number of parameters and facilitates their estimation through optimization methods (Senina et al., 2008; Senina et al., 2020). The SEAPODYM user manual and code are publicly available (https://github.com/PacificCommunity/seapodym-codebase).

A critical aspect of the model is the parameterization of spawning habitat, which is often challenging due to limited observations of early life stages. The SEAPODYM potential spawning habitat model integrates the effects of temperature (with the temperature field specified by the user), prey availability (zooplankton or primary production used as a proxy), and predation (micronekton, including vertically migrating mesopelagic organisms) to define suitable spawning habitats and predict the spatial distribution of early life stages. These effects are represented by a Gaussian function, a Holling type-III function, and a log-normal function, respectively (Table 2).

A MONTH effect is also tested here as in the delta-GAM using a seasonal Gaussian function allowing to identify the timing of the spawning peak and its standard deviation as well as the baseline and amplitude of the cycle, accounting for cyclicity. The product of these functions is in turn multiplied by the result of the local recruitment index derived from the stock-recruit Beverton-Holt function, thus providing the eventual larvae recruitment in the first cohort of the modelled population structure.

To compute the index in practice, one needs to test and calibrate the functions with observations. In our application, anchovy egg density surveys were used, along with the same predictor datasets used in the GAM approach (Section 2.3). The GAM SSB predictor was represented in SEAPODYM-Hs by an equivalent “year-effect” parameter (YEAR), to be also estimated in the range 0-1, and representing the interannual variability.

For the estimation of the SEAPODYM parameter values, we used a two-stage optimization approach that combines a global and local optimization method to achieve the best fit between the model and observed data. This is computed by minimizing the mean squared error between predicted potential spawning habitat suitability and observed egg densities. In the first stage, the global optimization approach uses a differential evolution (i.e. genetic) algorithm (Storn and Price, 1997) to find the best region in parameter space, avoiding local minima. In the second stage, the solution is improved further by means of a local optimization approach, by using the L-BFGS-B gradient method to quickly converge to the nearest minimum, while respecting parameter bounds. Both stages were implemented via the SciPy library in Python (Virtanen et al., 2020). This two-stage approach is flexible and adapted to fit complex ecological models where parameter interactions are non-linear and multiple local optima may exist.

The goodness of fit of SEAPODYM was assessed by comparing the predicted potential spawning habitat probabilities with the observed egg densities and the predictions obtained with the statistical GAM approach. We used standard metrics, such as the root mean square error (RMSE), the mean absolute error (MAE) and the correlation by year. The parameterization of both the GAM and SEAPODYM models was conducted for the period 2010–2023. The SEAPODYM model, using the same parameterization, was then applied retrospectively to the period 1998–2003 to assess whether it could provide insights into the collapse and subsequent recovery of the anchovy stock in the Bay of Biscay during the 2000s.

3 Results

3.1 Variability in observed egg density distributions

The distribution and abundance of eggs showed important changes between 2000 and 2023 (Figure 1). Over that period, the spawning area expanded by approximatively a factor of three and mean eggs density doubled, accompanied by a marked northward spatial shift in spawning distribution between the beginning and the end of the series. The core spawning area common to all years is located in the southern Bay of Biscay, between the coast and the shelf break, at latitudes 43°-46°N (Figures 1A, C). The 2000–2010 period corresponds to the collapse of the population, during which egg presence is largely restricted to this core area (Figure 1A).

Figure 1

Bellier et al. (2007), spawning was considered to occur in each grid cell when eggs density was above a minimum threshold of 5 eggs 10 m^-3. (A, C) show maps of the mean presence probability for the periods 2000–2010 and 2010-2023, respectively, with the corresponding interannual standard deviation shown in (B, D). Blue dashed lines on the maps show isobaths and the two regions used for the time series are delineated by blue (northern zone) and red (southern zone) rectangles. (E) shows the temporal evolution of the total surveyed area (black dashed line), the egg occurrence area (black solid line), and the mean egg density (blue line).

There is interannual variability, with relatively favourable years in 2001, 2007, 2011 and 2022 (Figure 1E). Overall, the population recovered during the period 2010–2023. This recovery was accompanied by a marked increase in egg density and an expansion of the spawning distribution (Figure 1C), which occurred both within the core southern region and across the northern plateau and offshore areas.

3.2 Oceanographic conditions in the Bay of Biscay

The Bay of Biscay’s continental shelf exhibits strong seasonality (see Supplementary Materials.). Sea surface temperature (SST) ranges from 11 °C in winter to 21 °C in summer, with rapid warming during the anchovy spawning period in May, while bottom temperatures peak later in autumn and show less variability. Primary production and zooplankton biomass both peak in late spring (May–June), with zooplankton cycles more pronounced in the north, whereas micronekton biomass peaks later, in July–August.

The inter-annual variability of the environmental variables used in this study is analysed using the z-score (see Supplementary Materials) and illustrated in Figure 2 for contrasting years with negative (2002) and positive (2020) anomalies. Z-scores exhibit pronounced interannual variability with few consecutive years of stable conditions for all variables.

Figure 2

The period 2002-2006, corresponding to the collapse and subsequent low level of the Bay of Biscay anchovy population, is characterized by predominantly negative anomalies in temperature, primary production and zooplankton and in some cases higher micronekton (see Supplementary Materials.). In particular, the year 2002 combines strong negative anomalies of temperature and zooplankton, together with a high positive anomaly in MICRO in the northern region (Figure 2). The year 2003 has the strongest negative anomaly in ZOO in the southern region, that is the historical core habitat for anchovy spawning. Taken together, these effects could potentially act in concert to strongly impair spawning success.

In the period 2010-2023, which coincides with the return to normal population levels, there is a succession of years dominated by negative anomalies (2010, 2012 and 2013) or positive anomalies (2011, 2020, 2022 and 2023) (see Supplementary Materials). The year 2020 is particularly noteworthy because, in association with high positive temperature anomalies in both the northern and southern regions, it also exhibits the highest positive anomalies in the whole 1998–2023 series for NPP (z = 3.43) and zooplankton (z = 2.25) in the southern region, as well as the second-highest values in the northern region (z = 2.17 for NPP and z = 1.23 for zooplankton). Unfortunately, no research cruise was conducted in 2020 due to the COVID-19 pandemic.

3.3 Exploratory analysis using delta-GAM

All possible model combinations were tested and their performance evaluated based on pseudo-R² and AICc criteria. Results for the best selected model are described below.

3.3.1 Logistic GAM (Presence-absence)

The best fitting model retained all predictors except NPP and ZOO (see Supplementary Materials). MONTH is a key predictor, which by itself already gets a high pseudo-R2 of 0.67 compared to the best model (0.82). This is not unexpected given the structure of the dataset, with egg presence recorded only in May during research cruises and absence imposed during September-March. Despite it is not included in the final best model, the variable ZOO is the second strongest explanatory variable for the presence–absence dataset (pseudo-R2 = 0.40) before SSB (0.38) and SST (0.33).

The partial effects for the best model are shown in Figure 3. They reveal an optimum SST around 16 °C; a positive response to SBT above 11 °C; a rapid decline in habitat favourability with increasing bathymetry, with the steepest decline occurring within the first 200 m; and a negative impact of micronekton biomass (MICRO) above 3 g WW m^-2. The SSB effect exhibits an increasing trend above a minimum value of approximately 80, 000 tonnes.

Figure 3

Somewhat counterintuitively, the best model is achieved without ZOO, despite it being the second-strongest single predictor. This is most likely because ZOO encompasses spatio-temporal information already captured by the other retained variables. The cross-correlation matrix supports this interpretation, highlighting in particular a moderate correlation between ZOO and MICRO (r = 0.54).

3.3.1 Gamma GAM (abundance)

The best fitting model for positive values only included all predictors (see Supplementary Materials), with a pseudo-R2 of 0.23. The partial effects of the predictors included in the selected Gamma GAM model (see Supplementary Materials.) resemble those obtained with the logistic GAM model for SST, SBT, BATHY, SSB and MICRO. Temperature, followed by BATHY and MICRO, shows the strongest contributions (as indicated by the range on the y-axes). Bottom temperature has a strong positive effect up to 12.5 °C; SST shows maximum positive effect between 14 and 17 °C; and BATHY has a negative effect, with favourability decreasing rapidly as depth increases from 0 to 200 m and at a slower rate thereafter. The functional responses of ZOO, MICRO and NPP each show a dome-shaped peak in the lower part of their range of values, around ~600 mg C m^-2 for NPP, 2.5 g WW m^-2 for MICRO and 0.6 g C m^-2 for ZOO. At higher values, the effect stabilises for NPP and increases for both ZOO and MICRO. However, as data are sparse at the upper end of these ranges, the reliability of the model responses in these regions is limited. SSB has a more moderate effect, with an overall increasing trend but large fluctuations between low and high values.

3.3.2 Delta GAM (abundance)

The delta GAM is the product of the logistic presence-absence model with the Gamma abundance model. The performance of this delta model is assessed in the intercomparison with the SEAPODYM-Hs models in section 3.5.

3.4 SEAPODYM spawning Habitat model

We used the same dataset as employed for the delta-GAM to estimate the parameters of the SEAPODYM-H_s spawning habitat model. Sea surface temperature (SST) or sea bottom temperature (SBT), and net primary production (NPP) or zooplankton biomass (ZOO), were tested alternatively in the temperature (Gaussian) and prey (Holling type III) functions, with or without the additional effect of bathymetry (equations in Table 2). To account for the effect associated with adult spawning biomass, a YEAR effect is included and estimated jointly with the other parameters.

The two-stage optimization procedure successfully converged, providing parameter estimates that minimized the mean squared error between predicted habitat suitability and observed egg density. The best performance was obtained when including bathymetry and ZOO as the food source, with a pseudo-R² of 0. 387 (Pearson’s r = 0.62) for the full dataset (including winter months pseudo-absence data; N = 3896). The RMSE and MAE are 114.61 and 56.00 (to be compared to 162.7 and 66.6, for the delta-GAM above). The resulting functions and their parameter estimates are shown in Figure 4.

Figure 4

Consistent with the GAM analysis, SBT proved to be a better temperature predictor than SST, with an estimated optimum of 13.69 °C (standard deviation: 1.46 °C). The bathymetry function imposes a rapid decline in the habitat suitability index with increasing depth. There is an increasing effect of predation by micronekton predicted above a threshold of 4.36 g WW m^-2. There is strong negative effect for zooplankton biomass values decreasing below 0.12 g C m^-2.

The year-effect values are compared to the SSB, highlighting the year 2022 as one outlier with very high value approximately four times greater than the mean value observed across the other years. Once this year is excluded, a Beverton-Holt type function, , can be fitted to the year effects, where S represents the spawning stock biomass (see Supplementary Materials). This relationship suggests a biologically meaningful nonlinear relationship indicative of density-dependent processes between the estimated year effect and mean annual SSB, analogous to a stock–recruitment relationship. Consequently, the estimated year effects exhibit behaviour consistent with ecological theory and biological data, supporting the credibility of the year-effect parameterization. It also highlights the anomalous mean value of the egg density index for the year 2022.

Finally, the parametrised SEAPODYM-Hs model is tested on the validation dataset extending back to 1998 and including the population collapse in the 2000s. The year effect is derived from the Beverton-Holt relationship with SSB estimated above (see Supplementary Materials). The spawning stock biomass (SSB) estimated by ICES is compared to the total annual habitat index Hs for the spawning season April-June, computed with and without the year effect in the study domain, i.e. south of 48°N and east of 6°W (Figure 5). Using May only does not alter the time trend, except for the range of Hs values (not shown). Hs series exhibit three consecutive very low values in 2004–06 and two more in 2009–10 with two much more favourable years in between (2007-08). A recovery is predicted to occurs after 2010 with no more very unfavourable years, even if there is still a strong interannual variability.

Figure 5

A strong correlation (r = 0.74) is observed between SSB and Hs when the year effect is included, partly reflecting the relationship with SSB used to construct the year effect. When the year effect is excluded, the habitat index represents only the influence of environmental conditions, and the correlation decreases to a moderate level (r = 0.42). The largest discrepancies occur in a few specific years. In 1998, 2003, and especially 2007, environmental conditions were estimated to be favourable while SSB was estimated to be low or very low; conversely, in 2015 and 2018, SSB was high whereas environmental conditions for Hs were predicted to be unfavourable.

A decline on Hs is predicted after 2003, associated with three consecutive years (2004–2006) of unfavourable environmental conditions combined with low spawning biomass (or year effect). The subsequent improvement in environmental conditions during 2007–2008 is followed by two additional years of poor environmental conditions (2009–2010) and is insufficient to allow a rapid recovery. Recovery of Hs occurs only after 2010, coinciding with a substantial increase in spawning biomass (or high year effects) and generally improved environmental conditions, although poor conditions may still occur sporadically in some particular years (e.g. 2015, 2018, and 2022).

3.5 Intercomparison of models

The Gamma-GAM, Delta-GAM and SEAPODYM-Hs models were compared using the same dataset 2010–2021 with Pearson correlation (r), RMSE and MAE (Table 3). MAE and RMSE provide complementary information on model performance. MAE represents the typical magnitude of the errors, while RMSE is more sensitive to large deviations and therefore highlights the influence of extreme errors.

When each year is considered separately, the Delta model performs best overall, with the highest yearly correlation scores on average. This represents the within-year correlation that measures how well the model predicts the spatial patterns within each year separately. Hs has slightly lower correlation than the delta-GAM, but has the lowest RMSE and MAE, meaning it predicted most accurately in absolute terms. The gamma GAM is the least performing, with lower correlation and larger errors.

When all years are considered together, the scores obtained with the SEAPODYM-Hs model are substantially better than those of the GAM models, indicating that Hs more effectively captures interannual variability. This is likely due to the explicit estimation of year effects in SEAPODYM-Hs, whereas in the GAMs the SSB is imposed and therefore provides less flexibility for explaining interannual differences. Consequently, although the delta-GAM performs well in reproducing spatial patterns, the Hs model offers superior predictive skill with respect to interannual variability. Predicted maps are compared with observations for each year in the Supplementary Materials.

4 Discussion

This study examines the drivers of anchovy spawning habitat variability using an extended dataset of egg density distributions, high-resolution environmental variables, and both statistical and mechanistic modelling approaches. In addition to temperature, which is classically used in such studies, we evaluate the influence of zooplankton and micronekton abundance, which are hypothesized to mediate feeding and predation and are key mechanisms proposed to explain variability in eggs density and subsequent larval and juvenile recruitment successes. In addition to environmental variability, we incorporate factors accounting for changes in adult spawning biomass (SSB; year effect), which is strongly influenced by fishing pressure, and the seasonal cycle (month effect), which is driven by species biology—particularly the reproductive cycle synchronized with seasonal environmental changes, e.g. through photoperiod (Bromage et al., 2001; Migaud et al., 2010)—to ensure optimal conditions for larval survival. Finally, bathymetry is also tested as a potential driver, as it may act as a natural cue for this species, which is known to exhibit a predominantly coastal distribution.

4.1 Drivers of variability in eggs density distribution and abundance

Seasonal reproductive timing is likely an evolved adaptation that optimizes the high energy demand required for reproduction to coincide with the annual window of favourable conditions, combining sufficient food availability and energy reserves for gonad development with favourable environmental conditions for the survival of eggs and rapid growth of larvae (Lowerre-Barbieri et al., 2011). The seasonal cycle of gonad development and maturation in Bay of Biscay anchovy is well documented, with gonadosomatic indices and histological stages showing marked seasonal progression (Erauskin-Extramiana et al., 2018; Coll and Bellido, 2019). However, there is a degree of plasticity that allows adaptation to changing environmental conditions. For example, Erauskin-Extramaniana et al. (2018) showed that the spawning peak of Biscay anchovy advanced at a rate of 5.5 days per decade between 1987 and 2015.

Regarding environmental factors, temperature is one of the main drivers explaining the suitability of environmental conditions for anchovy spawning, with bottom temperature playing a key role, as already highlighted in other studies (Planque et al., 2007). In addition to driving the timing of spawning with other possible cues (e.g. photoperiod), temperature also exerts strong influences on early life-stage dynamics. Experimental and field evidence indicates that incubation duration, larval growth rate, and survival are strongly temperature-dependent, such that temperatures outside optimal ranges (both lower and higher) can lead to elevated egg and larval mortality and modify the duration of vulnerable stages, thus affecting recruitment success (Bernal et al., 2012).

Bottom temperature is possibly a better indicator of habitat quality than SST in these retention zones, because it reflects the actual environmental conditions experienced by the eggs and larvae. Eggs of anchovy in the Bay of Biscay are distributed in the first 50 m of the water column (Boyra et al., 2003). Therefore, the bottom temperature might be more relevant for egg survival and development (Motos et al., 1996), given that surface waters warm faster, and thus could lead to a mismatch with the actual conditions experienced by eggs. Nevertheless, there is also a clear signal detected for an optimal SST around 16 °C.

The use of newly available essential ocean variables for zooplankton and micronekton, delivered through the Copernicus Marine Service catalogue, proves useful for incorporating food-resource availability and predation effects. However, these variables still require bias correction in shallow coastal regions, where values are strongly overestimated, maybe due to the lack of representation of benthic systems that capture part of the energy flux derived from phytoplankton, or the overestimation of coastal primary production based on ocean colour. This bias should be accounted for in future versions of these products.

The effect of micronekton appears particularly important and it is consistently selected across all GAM models. Predation has long been recognized as having the potential to control or regulate recruitment levels (Bailey and Houde, 1989), but it remains difficult to monitor, and explicit hypotheses are rarely tested (Houde, 2008). Krautz et al. (2007) investigated euphausiid predation on anchoveta (Engraulis ringens) eggs in the Humboldt Current and quantified predation-attributable mortality rates accounting for 24–27% of total egg mortality. When combined with mortality due to cannibalism by adult anchoveta, their results suggest that approximately half of the daily anchoveta egg production dies each day. Even higher contributions of predation to total mortality (0–98%) were reported by Voss et al. (2011) in their study of sprat egg predation in the central Baltic Sea. They showed that predation intensity was highest in spring, coinciding with high predator abundance (adult sprat) and strong vertical overlap between adults and eggs. These findings highlight potential sources of the extreme variability observed in egg density estimates obtained during research cruises.

The predation effect associated with micronekton in the GAM models shows a response pattern broadly similar to that described by the log-normal function in the mechanistic H_s model. It is well identified as an explanatory variable over the full range of values predicted within the study domain. The micronekton variable is dominated by the biomass of migrant mesopelagic species, which is particularly high over the continental slope. Predation pressure is therefore likely to be strongest in this area and may contribute to the concentration of successful spawning between the coastal zone and the slope.

The response of the models to the zooplankton variable is ambiguous. In the logistic GAM, ZOO is the strongest single predictor, yet when considered in combination with other variables it is excluded, likely because the other retained predictors collectively capture the variance it would otherwise explain. In the SEAPODYM Hs model, the functional response to ZOO is positive by design, and suggesting that habitat suitability saturates quickly at relatively low zooplankton biomass values.

We also note that the years 2002–2006, when the stock was at very low levels, coincide with large negative anomalies in zooplankton biomass. The importance of prey availability is clear, both for larval survival and for the energy accumulation adult spawners require to produce gametes. However, zooplankton biomass during the spawning month of May may generally be sufficiently abundant so as not to constitute a limiting factor. This is consistent with the fact that the seasonal peak in zooplankton biomass typically occurs in May.

Alternatively, the zooplankton product may lack sufficient accuracy despite the bias correction applied, as suggested by the residual discrepancies we identified. Zooplankton sampling during PELGAS is conducted at night and therefore captures both resident epipelagic and deeper migrant zooplankton, whereas anchovy egg sampling is conducted during the day. This diel sampling discrepancy is a potential source of bias arising from diel vertical migration but cannot be fully resolved with the currently available dataset and the present version of the zooplankton model. Improving the zooplankton model to explicitly represent the resident and migratory fractions would help determine whether this accounts for the systematic offshore discrepancies observed. More broadly, this highlights the need for a deeper evaluation and revision of the zooplankton model to improve its predictive skill in both offshore and coastal shallow-water environments.

Finally, it is possible that additional environmental variables not considered here also influence spawning success and egg density distributions, either directly (e.g. currents, wind mixing or salinity) or indirectly through their effects on larvae and spawners (Coombs et al., 2003; Goarant et al., 2007; Džoić et al., 2022). For instance, dissolved oxygen may influence adult behaviour and larvae mortality (Kreiner et al., 2009; Bertrand et al., 2011). It should be noted that the mixed-layer thickness (from the Copernicus Marine service IBI product) was tested by replacing the mean density index with the product of this index and the mixed-layer depth. This analysis was not shown, as it did not alter the results, confirming that the mean value of the product is representative of the average distribution within the mixed layer. Other potentially relevant variables could be tested in future revisions of this study using the same modelling framework.

4.2 Modelling approach

In our modelling approach, we address the issue of observational absences in the dataset using a delta generalized additive model (delta-GAM) and compare this approach with a more mechanistic model (SEAPODYM-Hs). Incorporating absence data in the delta-GAM framework improves model performance relative to a gamma model based solely on positive observations. This approach also allows the inclusion of a well-known seasonal signal through the addition of pseudo-absences during winter months. The delta-GAM, which explicitly separates the presence–absence and abundance components, performs slightly better than Hs in explaining the observed spatial patterns. The mechanistic SEAPODYM model, which integrates temperature, prey availability, and predation processes, also shows strong predictive skill, with smaller errors, and better prediction of interannual variability thanks to the estimated year effects.

Spawning stock biomass (SSB), or its equivalent year-effect in the Hs model, together with month (representing the seasonal cycle), emerge as the primary drivers of variability, highlighting the importance of stock–recruitment dynamics and the biological adaptation of the species to a strongly seasonal environment. As observed elsewhere (e.g. Anderson et al., 2008; Cury et al., 2014; Szuwalski et al., 2015), this study underscores the need for careful management of fishing pressure on small pelagic species to avoid the harmful combination of overexploitation and consecutive years of unfavourable environmental conditions. We demonstrate that three successive years of unfavourable conditions, corresponding approximately to the mean lifespan of the species, combined with low spawning stock biomass (SSB), are likely to result in the collapse of both the stock and the fishery.

Another key achievement of this study is the demonstration that a mechanistic model, based on appropriate biological data and meaningful functional relationships, performs as well as a more complex GAM with more parameters. Although other sources of variability not explicitly included in the model may influence interannual fluctuations and be confounded with the year effect, this approach offers an independent mean to estimate variability that can be directly compared with the effect of spawning stock biomass (SSB). The spawning habitat index (H_s) could therefore be proposed as a routinely updated indicator. Indeed, while the index without year effect could be produced using seasonal forecasts, the index with the estimated yearly effect could be computed after each new sampling cruise to provide an independent indicator to add to those achieved through the daily egg production method (Hüssy and Røttingen, 1982; Ward et al., 2021) and stock assessment model fitting (ICES, 2024). If the predictive skills of the SEAPODYM-Hs model are confirmed, and eventually improved, management could adopt a more dynamic approach. This would involve closely monitoring SSB levels in combination with consecutive unfavourable environmental years, including the previous two and the newly forecasted year, and proactively reducing catch quotas to prevent another collapse. The same approach can be applied to other small pelagic species (e.g. sardine) in the Bay of Biscay or in other coastal systems.

Climate change, which is now superimposing its effects on natural variability, could also be investigated using this framework by relying on climate projections and decadal forecasts of environmental variables derived from coupled ocean circulation and biogeochemical models. However, for small pelagic fishes in coastal regions such as the Bay of Biscay, achieving sufficiently high spatial resolution in these projections and decadal forecasts is essential. This remains challenging because of the associated computational costs. Nevertheless, downscaling techniques are available, and several ongoing projects are addressing this issue with the aim of delivering such high-resolution products in the coming years (Shapiro et al., 2021; Soares et al., 2024).

This need is becoming increasingly urgent, as climate-driven changes may already be underway. For example, van der Kooij et al. (2024) showed that anchovy spawning activity in the Bay of Biscay has expanded both spatially and temporally, resulting in increased larval transport and survival into the English Channel during the autumns of 2019 and 2020, well beyond the nearest previously known spawning grounds.

References

• 1

  AlheitJ.NiquenM. (2004). Regime shifts in the Humboldt Current ecosystem. Prog. Oceanogr.60, 201–222. doi: 10.1016/j.pocean.2004.02.006

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  AndersonC. N. K.HsiehC.-H.SandinS. A.HewittR.HollowedA.BeddingtonJ.et al. (2008). Why fishing magnifies fluctuations in fish abundance. Nature452, 835–839. doi: 10.1038/nature06851

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 3

  BaileyK. M.HoudeE. D. (1989). Predation on eggs and larvae of marine fishes and the recruitment problem. Adv. Mar. Biol.25, 1–83. doi: 10.1016/S0065-2881(08)60187-X

  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  BehrenfeldM. J.FalkowskiP. G. (1997). Photosynthetic rates derived from satellite-based chlorophyll concentration. Limnol. Oceanogr.42, 1–20. doi: 10.4319/lo.1997.42.1.0001

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  BellierE.PlanqueB.PetitgasP. (2007). Historical fluctuations in spawning location of anchovy (Engraulis encrasicolus) and sardine (Sardina pilchardus) in the Bay of Biscay during 1967–1973 and 2000–2004. Fish. Oceanogr.16, 1–15. doi: 10.1111/j.1365-2419.2006.00410.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 6

  BergeronJ.-P.DelmasD.KouetaN. (2010). Do river discharge rates drive the overall functioning of the pelagic ecosystem over the continental shelf of the Bay of Biscay (NE Atlantic)? A comparison of two contrasting years with special reference to anchovy (Engraulis encrasicolus) nutritional state. J. Oceanogr.66, 621–631. doi: 10.1007/s10872-010-0051-7

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  BergeronJ.-P.KouetaN.MasséJ. (2013). Interannual fluctuations in spring pelagic ecosystem productivity in the Bay of Biscay measured by mesozooplankton aspartate transcarbamylase activity and relationships with anchovy population dynamics. Fish. Res.143, 184–190. doi: 10.1016/j.fishres.2013.02.006

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  BernalM.JiménezM. P.DuarteJ. (2012). Anchovy (Engraulis encrasicolus) egg development in the Gulf of Cádiz and its comparison with development rates in the Bay of Biscay. Fish. Res.108, 240–248. doi: 10.1016/j.fishres.2011.04.010

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  BertrandA.ChaigneauA.PeraltillaS.LedesmaJ.GracoM.MonettiF.et al. (2011). Oxygen: A fundamental property regulating pelagic ecosystem structure in the coastal southeastern tropical Pacific. PloS One6 (12), e29558. doi: 10.1371/journal.pone.0029558

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  BorjaA.FontánA.SáenzJ.ValenciaV. (2008). Climate, oceanography, and recruitment: the case of the Bay of Biscay anchovy (Engraulis encrasicolus). Fish. Oceanogr.17, 477–493. doi: 10.1111/j.1365-2419.2008.00494.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  BoyraG.RuedaL.CoombsS. H.SundbyS.AdlandsviskB.SantosM.et al. (2003). Modelling the vertical distribution of eggs of anchovy (Engraulis encrasicolus) and sardine (Sardina pilchardus). Fish. Oceanogr.12, 381–395. doi: 10.1046/j.1365-2419.2003.00260.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  BradleyA. P. (1997). The use of the area under the ROC curve in the evaluation of machine learning algorithms. Pattern Recognit.30, 1145–1159. doi: 10.1016/S0031-3203(96)00142-2

  • CrossRef
  • Google Scholar

• 13

  BromageN. R.PorterM.RandallC. F. (2001). The environmental regulation of maturation in farmed finfish with special reference to the role of photoperiod and melatonin. Aquaculture197, 63–98. doi: 10.1016/S0044-8486(01)00583-X

  • CrossRef
  • Google Scholar

• 14

  CavanaughJ. E. (1997). Unifying the derivations for the Akaike and corrected Akaike information criteria. Stat. Probab. Lett.33, 201–208. doi: 10.1016/S0167-7152(96)00105-8

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  ColasF.TardivelM.PerchocJ.LunvenM.ForestB.GuyaderG.et al. (2018). The ZooCAM, a new in-flow imaging system for fast onboard counting, sizing and classification of fish eggs and metazooplankton. Prog. Oceanogr.166, 54–65. doi: 10.1016/j.pocean.2018.07.004

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  CollM.BellidoJ. M. (2019). Reproductive biology of European anchovy (Engraulis encrasicolus) in the Bay of Biscay: seasonal dynamics and environmental influence. ICES J. Mar. Sci.76, 2106–2118. doi: 10.1093/icesjms/fsz129

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  CoombsS. H.GiovanardiO.HallidayN. C.FranceschiniG.ConwayD. V. P.ManzuetoL.et al. (2003). Wind mixing, food availability and mortality of anchovy larvae Engraulis encrasicolus in the northern Adriatic Sea. Mar. Ecol. Prog. Ser.248, 221–235. doi: 10.3354/meps248221

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  CortenA. (1986). On the causes of the recruitment failure of herring in the central and northern North Sea in the years 1972–1978. J. Cons.42, 281–294. doi: 10.1093/icesjms/42.3.281

  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  CuryP. M.FromentinJ.-M.FiguetS.BonhommeauS. (2014). Resolving Hjort’s Dilemma: How is recruitment related to spawning stock biomass in marine fish? Oceanography27, 42–47. doi: 10.5670/oceanog.2014.85

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  CushingD. H. (1975). Marine Ecology and Fisheries (Cambridge: Cambridge University Press), 1–278.

  • Google Scholar
  • View reference in article

• 21

  CushingD. H. (1995). The long-term relationship between zooplankton and fish: IV. Spatial/temporal variability and prediction. ICES J. Mar. Sci.52, 611–626. doi: 10.1016/1054-3139(95)80076-X

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  CushingD. H.HarrisJ. G. K. (1973). Stock and recruitment and the problem of density-dependence. Rapp. P.-v. Réun. Cons. Int. Explor. Mer164, 142–155.

  • Google Scholar
  • View reference in article

• 23

  DorayM.HuretM.AuthierM.DuhamelE.RomagnanJ.-B.DupuyC.et al. (2018a). Gridded maps of pelagic ecosystem parameters collected in the Bay of Biscay during the PELGAS integrated survey. SEANOE dataset. doi: 10.17882/53389

  • CrossRef
  • Google Scholar

• 24

  DorayM.PetitgasP.RomagnanJ.-B.HuretM.DuhamelE.DupuyC.et al. (2018b). The PELGAS survey: ship-based integrated monitoring of the Bay of Biscay pelagic ecosystem. Prog. Oceanogr.166, 15–29. doi: 10.1016/j.pocean.2017.09.015

  • CrossRef
  • Google Scholar

• 25

  DormannC. F.ElithJ.BacherS.BuchmannC.CarlG.CarréG.et al. (2013). Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. Ecography36, 27–46. doi: 10.1111/j.1600-0587.2012.07348.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  DžoićT.ZoricaB.MatićF.SestanovićM.Cikes KećV. (2022). Cataloguing environmental influences on spatiotemporal variability of Adriatic anchovy early life stages using neural network analysis. Front. Mar. Sci.9, 997937. doi: 10.3389/fmars.2022.997937

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  Erauskin-ExtramianaM.AlvarezP.ArrizzabalagaH.IbaibarriagaL.UriarteA.CotanoU.et al. (2019). Historical trends and future distribution of anchovy spawning in the Bay of Biscay. Deep-Sea Res. II159, 169–182. doi: 10.1016/j.dsr2.2018.07.007

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  Erauskin-ExtramianaM.ArrizabalagaH.HobdayA. J.CabréA.IbaibarriagaL.ChustG. (2018). Large-scale distribution of tuna species in a warming ocean. Glob. Change Biol.25, 204–216. doi: 10.1111/gcb.14487

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  GoarantA.PetitgasP.BourriauP. (2007). Anchovy (Engraulis encrasicolus) egg density variation with sea surface salinity in the Bay of Biscay. Mar. Biol.151, 1677–1688. doi: 10.1007/s00227-007-0624-1

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  GrandremyN.BourriauP.DachéE.DanielouM-M.DorayM.DupuyC.et al. (2024). Metazoan zooplankton in the Bay of Biscay: a 16-year record of individual sizes and abundances obtained using the ZooScan and ZooCAM imaging systems. Earth Syst. Sci. Data16, 1265–1282. doi: 10.5194/essd-16-1265-2024

  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  GrüssA.DrexlerM.AinsworthC. H. (2014). Using delta generalized additive models to produce distribution maps for spatially explicit ecosystem models. Fish. Res.159, 11–24. doi: 10.1016/j.fishres.2014.02.010

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 32

  HernandezO.LehodeyP.SeninaI.EchevinV.AyonP.BertrandA.et al. (2014). Understanding mechanisms that control fish spawning and larval recruitment: Parameter optimization of an Eulerian model (SEAPODYM-SP) with Peruvian anchovy and sardine eggs and larvae data. Prog. Oceanogr.123, 105–122. doi: 10.1016/j.pocean.2014.03.001

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  HjortJ. (1914). Fluctuations in the great fisheries of northern Europe viewed in the light of biological research. Rapp. P.-v. Réun. Cons. Perm. Int. Explor. Mer19, 1–228.

  • Google Scholar
  • View reference in article

• 34

  HoudeE. D. (2008). Emerging from hjort’s shadow. J. Northwest Atl. Fish. Sci.41, 53–70. doi: 10.2960/J.v41.m634

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 35

  HüssyK.RøttingenI. (1982). A review of the daily egg production method for the estimation of spawning-stock biomass of pelagic fish ( ICES C.M).

  • Google Scholar
  • View reference in article

• 36

  IbaibarriagaL.FernándezC.UriarteA. (2008). A two-stage biomass dynamic model for Bay of Biscay anchovy: a Bayesian approach. ICES J. Mar. Sci.65, 191–205. doi: 10.1093/icesjms/fsn002

  • CrossRef
  • Google Scholar
  • View reference in article

• 37

  IbaibarriagaL.FernándezC.UriarteA. (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

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 38

  ICES (2019). Working group on southern horse mackerel, anchovy and sardine (WGHANSA). ICES Sci. Rep.1, 1–653. doi: 10.17895/ices.pub.4983

  • CrossRef
  • Google Scholar

• 39

  ICES (2024). Anchovy (Engraulis encrasicolus) in subarea 8 (Bay of biscay). ICES Advice: Recurrent Advice. doi: 10.17895/ices.advice.25019084.v1

  • CrossRef
  • Google Scholar

• 40

  IrigoienX.FiksenØ.CotanoU.UriarteA.AlvarezP.ArrizabalagaH.et al. (2007). Could Bay of Biscay anchovy recruit through a spatial loophole? Prog. Oceanogr.74, 132–148. doi: 10.1016/j.pocean.2007.04.011

  • CrossRef
  • Google Scholar
  • View reference in article

• 41

  JensenO. P.SeppeltR.MillerT. J.BauerL. J. (2005). Winter distribution of blue crab (Callinectes sapidus) in Chesapeake Bay: application and cross-validation of a two-stage generalized additive model. Mar. Ecol. Prog. Ser.299, 239–255. doi: 10.3354/meps299239

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 42

  KaharR.AhmadN.AraiT. (2023). Year-round spawning of three tropical Cypriniformes fishes in Southeast Asia. Sci. Rep.13, 8971. doi: 10.1038/s41598-023-36065-9

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 43

  KoenigsteinS.JacoxM. G.Pozo BuilM.FiechterJ.MuhlingB. A.BrodieS.et al. (2022). Population projections of Pacific sardine driven by ocean warming and changing food availability in the California Current. ICES J. Mar. Sci.79, 2510–2523. doi: 10.1093/icesjms/fsac191

  • CrossRef
  • Google Scholar
  • View reference in article

• 44

  KrautzM. C.MuckP.RosalesC. (2007). Euphausiid predation on anchoveta (Engraulis ringens) eggs off Peru. J. Plankton Res.29, 415–425. doi: 10.1093/plankt/fbm030

  • CrossRef
  • Google Scholar
  • View reference in article

• 45

  KreinerA.StenevikE. K.EkauW. (2009). Sardine (Sardinops sagax) and anchovy (Engraulis encrasicolus) larvae avoid regions with low dissolved oxygen concentration in the northern Benguela Current system. J. Fish Biol.74, 270–277. doi: 10.1111/j.1095-8649.2008.02124.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 46

  LazureP.JégouA.-M. (1998). 3d modelling of seasonal evolution of Loire and Gironde plumes on the Bay of Biscay continental shelf. Oceanol. Acta21, 165–177. doi: 10.1016/S0399-1784(98)80006-6

  • CrossRef
  • Google Scholar
  • View reference in article

• 47

  LehodeyP.AlheitJ.BarangeM.BaumgartnerT.BeaugrandG.DrinkwaterK.et al. (2006). Climate variability, fish and fisheries. J. Climate19(20), 5009–5030. doi: 10.1175/JCLI3898.1

  • CrossRef
  • Google Scholar
  • View reference in article

• 48

  LehodeyP.SeninaI.MurtuguddeR. (2008). A spatial ecosystem and populations dynamics model (SEAPODYM) - modelling of tuna and tuna-like populations. Prog. Oceanog78, 304–318. doi: 10.1016/j.pocean.2008.06.004

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 49

  LehodeyP.MurtuguddeR.SeninaI. (2010). Bridging the gap from ocean models to population dynamics of large marine predators: a model of mid-trophic functional groups. Prog. Oceanogr.84, 69–84. doi: 10.1016/j.pocean.2009.09.008

  • CrossRef
  • Google Scholar

• 50

  LloretJ.PalomeraI.SalatJ.SoléI. (2004). Impact of freshwater input and wind on landings of anchovy (Engraulis encrasicolus) and sardine (Sardina pilchardus) in shelf waters surrounding the Ebro River delta. Fish. Oceanogr.13, 102–110. doi: 10.1046/j.1365-2419.2003.00279.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 51

  Lowerre-BarbieriS. K.GaniasK.Saborido-ReyF.MuruaH.HunterJ. R. (2011). Reproductive timing in marine fishes: variability, temporal scales, and methods. Mar. Coast. Fish.3, 71–91. doi: 10.1080/19425120.2011.556932

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 52

  MaS.ChengJ.LiJ.LiuY.WanR.TianY. (2019). Interannual to decadal variability in the catches of small pelagic fishes from China Seas and its responses to climatic regime shifts. Deep-Sea Res. II159, 112–129. doi: 10.1016/j.dsr2.2018.10.005

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 53

  MarchalP.GiraldoC.JohnsD.LefebvreS.LootsC.ToomeyL. (2025). Effects of zooplankton abundance on the spawning phenology of winter-spawning Downs herring (Clupea harengus). PLoS ONE20(2), e0310388. doi: 10.1371/journal.pone.0310388

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 54

  MenuC.PecquerieL.BacherC.DorayM.HattabT.van der KooijJ.et al. (2023). Testing the bottom-up hypothesis for the decline in size of anchovy and sardine across European waters through a bioenergetic modeling approach. doi: 10.1016/j.pocean.2022.102943

  • CrossRef
  • Google Scholar
  • View reference in article

• 55

  MigaudH.DavieA.TaylorJ. F. (2010). Current knowledge on the photoneuroendocrine regulation of reproduction in temperate fish species. J. Fish Biol.76, 27–68. doi: 10.1111/j.1095-8649.2009.02500.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 56

  MotosL.UriarteA.ValenciaV. (1996). The spawning environment of the Bay of Biscay anchovy (Engraulis encrasicolus). Sci. Mar.60, 117–140.

  • Google Scholar
  • View reference in article

• 57

  PeckM. A.RegleroP.TakahashiM.CatalánI. A. (2013). Life cycle ecophysiology of small pelagic fish and climate-driven changes in populations. Prog. Oceanogr.116, 220–245. doi: 10.1016/j.pocean.2013.05.012

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 58

  PepinP. (2015). Death from near and far: alternate perspectives on size-dependent mortality in larval fish. ICES J. Mar. Sci.73, 196–203. doi: 10.1093/icesjms/fsv160

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 59

  PlanqueB.BellierE.LazureP. (2007). Modelling potential spawning habitat of sardine (Sardina pilchardus) and anchovy (Engraulis encrasicolus) in the Bay of Biscay. Fish. Oceanogr.16, 16–30. doi: 10.1111/j.1365-2419.2006.00411.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 60

  RothschildB. J. (2000). Fish stocks and recruitment: the past thirty years. ICES J. Mar. Sci.57, 191–201. doi: 10.1006/jmsc.2000.0645

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 61

  SchwartzloseR. A.AlheitJ.BakunA.BaumgartnerT. R.CloeteR.CrawfordR. J. M.et al. (1999). Worldwide large-scale fluctuations of sardine and anchovy populations. S. Afr. J. Mar. Sci.21, 289–347. doi: 10.2989/025776199784125962

  • CrossRef
  • Google Scholar
  • View reference in article

• 62

  SeninaI.SibertJ.LehodeyP. (2008). Parameter estimation for basin-scale ecosystem-linked population models of large pelagic predators: application to skipjack tuna. Prog. Oceanogr.78, 319–335. doi: 10.1016/j.pocean.2008.06.003

  • CrossRef
  • Google Scholar

• 63

  SeninaI.LehodeyP.SibertJ.HamptonJ. (2020). Integrating tagging and fisheries data into a spatial population dynamics model to improve its predictive skills. Can. J. Fish. Aquat.Sci.77, 576–593. doi: 10.1139/cjfas-2018-0470

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 64

  ServénD.BrummittC. (2018). pyGAM: generalized additive models in Python. Zenodo. doi: 10.5281/zenodo.1208723

  • CrossRef
  • Google Scholar

• 65

  ShapiroG. I.Gonzalez-OndinaJ. M.BelokopytovV. N. (2021). High-resolution stochastic downscaling method for ocean forecasting models and its application to the Red Sea dynamics. Ocean Sci.17, 891–907. doi: 10.5194/os-17-891-2021

  • CrossRef
  • Google Scholar

• 66

  SheltonA. O.MangelM. (2011). Fluctuations of fish populations and the magnifying effects of fishing. Proc. Natl. Acad. Sci. U.S.A.108, 7075–7080. doi: 10.1073/pnas.1100334108

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 67

  SipleM. C.KoehnL. E.JohnsonK. F.PuntA. E.CanalesT. M.CarpiP.et al (2021). Considerations for management strategy evaluation for small pelagic fishes. Fish Fish.22, 1167–1186. doi: 10.1111/faf.12579

  • CrossRef
  • Google Scholar

• 68

  SoaresP. M. M.JohannsenF.LimaD. C. A.LemosG.BentoV. A.BushenkovaA. (2024). High-resolution downscaling of CMIP6 Earth system and global climate models using deep learning for Iberia. Geosci. Model Dev.17, 229–259. doi: 10.5194/gmd-17-229-2024

  • CrossRef
  • Google Scholar

• 69

  StornR.PriceK. (1997). Differential evolution – a simple and efficient heuristic for global optimization over continuous spaces. J. Glob. Optim.11, 341–359. doi: 10.1023/A:1008202821328

  • CrossRef
  • Google Scholar

• 70

  SzuwalskiC. S.Vert-PreK. A.PuntA. E.BranchT. A.HilbornR. (2015). Examining common assumptions about recruitment: A meta-analysis of recruitment dynamics for worldwide marine fisheries. Fish Fish.16, 633–648. doi: 10.1111/faf.12083

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 71

  TanakaK. R.BelknapS. L.HomolaJ. J.ChenY. (2017). A statistical model for lobster shell disease in Long Island Sound. PloS One12, e0172123. doi: 10.1371/journal.pone.0172123

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 72

  TommasiD.StockC. A.PegionK.VecchiG. A.MethotR. D.AlexanderM. A.et al. (2017). Improved management of small pelagic fisheries through seasonal climate prediction. Ecol. Appl.27, 378–388. doi: 10.1002/eap.1458

  • CrossRef
  • Google Scholar

• 73

  TrenkelV. M.HuseG.MacKenzieB.AlvarezP.ArrizabalagaH.CastonguayM.et al. (2014). Comparative ecology of widely-distributed pelagic fish species in the North Atlantic: implications for modelling climate and fisheries impacts. Prog. Oceanogr.129, 219–243. doi: 10.1016/j.pocean.2014.04.030

  • CrossRef
  • Google Scholar

• 74

  van der KooijJ.McKeownN.CampanellaF.BoyraG.DorayM.Santos MocoroaM.et al. (2024). Northward range expansion of Bay of Biscay anchovy into the English Channel. Mar. Ecol. Prog. Ser.741, 217–236. doi: 10.3354/meps14603

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 75

  VirtanenP.GommersR.OliphantT. E.HaberlandM.ReddyT.CournapeauD.et al (2020). SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat. Methods17, 261–272. doi: 10.1038/s41592-019-0686-2

  • CrossRef
  • Google Scholar

• 76

  VossR.HinrichsenH.-H.StepputtisD.BernreutherM.HuwerB.NeumannV.et al. (2011). Egg mortality: predation and hydrography in the central Baltic Sea. ICES J. Mar. Sci.68, 1379–1390. doi: 10.1093/icesjms/fsr061

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 77

  WardT. M.GrammerG. L.IveyA. R.SmartJ. J.McGarveyR. (2021). Increasing the precision of the daily egg production method: 2020’s remix of a 1980’s classic. ICES J. Mar. Sci.78, 1177–1195. doi: 10.1093/icesjms/fsab015

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 78

  WardE. J.HunsickerM. E.MarshallK. N.OkenK. L.SemmensB. X.FieldJ. C.et al. (2024). Leveraging ecological indicators to improve short-term forecasts of fish recruitment. Fish Fish.25, 895–909. doi: 10.1111/faf.12850

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 79

  WelshA. H.CunninghamR. B.DonnellyC. F.LindenmayerD. B. (1996). Modelling the abundance of rare species: statistical models for counts with extra zeros. Ecol. Modell.88, 297–321. doi: 10.1016/0304-3800(95)00113-1

  • CrossRef
  • Google Scholar

• 80

  WoodS. N. (2017). Generalized Additive Models: An Introduction with R (2nd ed.). (Boca Raton: Chapman & Hall/CRC). doi: 10.1201/9781315370279

  • CrossRef
  • Google Scholar