Front. Mar. Sci.Frontiers in Marine ScienceFront. Mar. Sci.2296-7745Frontiers Media S.A.10.3389/fmars.2021.570354Marine ScienceOriginal ResearchChanges in Life History Traits of Small Pelagic Fish in the Western Mediterranean SeaAlbo-PuigserverMarta12*PenninoMaria Grazia3BellidoJose María4ColmeneroAna Isabel1GiráldezAna5HidalgoManuel6Gabriel RamírezJohn1SteenbeekJeroen7TorresPedro5Cousido-RochaMarta3CollMarta11Renewable Marine Resources Department, Institute of Marine Sciences (ICM-CSIC), Barcelona, Spain2Centro de Ciências do Mar, Universidade do Algarve, Faro, Portugal3Centro Oceanográfico de Vigo, Instituto Español de Oceanografía (IEO-CSIC), Vigo, Spain4Centro Oceanográfico de Murcia, Instituto Español de Oceanografía (IEO-CSIC), San Pedro del Pinatar, Spain5Centro Oceanográfico de Málaga, Instituto Español de Oceanografía (IEO-CSIC), Fuengirola, Spain6Centro Oceanográfico de Baleares, Instituto Español de Oceanografía (IEO-CSIC), Ecosystem Oceanography Group (GRECO), Palma, Spain7Ecopath International Initiative Research Association, Barcelona, Spain
Edited by: Manuel J. Zetina-Rejón, Instituto Politécnico Nacional (IPN), Mexico
Reviewed by: M. Cristina Mangano, Stazione Zoologica Anton Dohrn Napoli, Italy; Claire Saraux, Institut Français de Recherche pour l’Exploitation de la Mer (IFREMER), France; Emigdio Marín-Enríquez, National Council of Science and Technology (CONACYT), Mexico
*Correspondence: Marta Albo-Puigserver, marta.albo.puigserver@gmail.com
This article was submitted to Marine Fisheries, Aquaculture and Living Resources, a section of the journal Frontiers in Marine Science
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Small pelagic fish (SPF) in the western Mediterranean Sea are key elements of the marine food web and are important in terms of biomass and fisheries catches. Significant declines in biomass, landings, and changes in the age/size structure of sardine Sardina pilchardus and anchovy Engraulis encrasicolus have been observed in recent decades, particularly in the northern area of the western Mediterranean Sea. To understand the different patterns observed in SPF populations, we analyzed key life history traits [total length at age, length at maturity, gonadosomatic index (GSI), and body condition (Kn)] of sardine and anchovy collected between 2003 and 2017, from different fishing harbors distributed along a latitudinal gradient from northern to southern Spain. We used Generalized Linear Models (GLM) to estimate the length at maturity and Generalized Additive Models (GAMs) to test the relationship with environmental variables (seawater temperature, water currents, and net primary productivity). The life history traits of both species presented seasonal, interannual and latitudinal differences with a clear decline in length at age, length at first maturity, and body condition, for both species in the northern part of the study area. In the southern part, on the contrary, life history traits did not present a clear temporal trend. The environmental conditions partially explained the long-term changes in life history traits, but the selected variables differed between areas, highlighting the importance of regional oceanographic conditions to understand the dynamics of small pelagic fish. The truncated length-at-age pattern for both species with the disappearance of the larger individuals of the population could have contributed to the poor condition of small pelagic fish populations in the northern part of the western Mediterranean Sea in recent years. In the south area, recent declines in body condition for sardine and anchovy were observed and could be a possible first sign for future population declines. This study highlights the importance of understanding the trade-off between the energy invested in reproduction, maintenance and growth at seasonal and interannual level to advance our knowledge on how environmental and human pressures influence population dynamics of small pelagic fish at local and regional scales.
sardineanchovylength at agematurityreproductionbody conditionCTM2017-88939-RMinisterio de Ciencia, Innovación y Universidades10.13039/100014440Introduction
Marine ecosystems are subjected to different global changes and anthropogenic disturbances. Strong pressures such as changes in the environmental conditions, fishing exploitation or pollution are impacting their functioning (Halpern et al., 2015; Ramírez et al., 2017). As a consequence, fluctuations in life history traits of marine fish affecting their population dynamics have been widely described (Rochet, 1998; De Roos et al., 2003). Changes in the reproductive period or the size at first maturity are usually related to changes in somatic growth and condition and ultimately affect natural mortality and recruitment success (Lloret et al., 2013; Stawitz and Essington, 2018). All of these changes are a consequence of the trade-off between the energy invested in growth, maintenance and reproduction (Stearns, 1989; McBride et al., 2015).
Because fishing is often the main cause of mortality for exploited species, it may induce phenotypic adaptive responses with changes in length and/or age at first maturity and declines in biomass (Nash et al., 2000; Ernande et al., 2004). Life history changes can be reversible if species exhibit phenotypic plasticity, while they can become irreversible if the species evolutionary adapt through selection (Stenseth and Rouyer, 2008). For example, size-selective mortality can trigger a selection for those individuals that have an early maturation (De Roos et al., 2006; Jørgensen et al., 2007). Fishing-induced truncation of size structure has been mostly found in species with protracted demography such as gadoids (Jørgensen et al., 2007), while evidence in species with short-life cycle, such as small pelagic fish, is scarce (but see Walsh et al., 2006; Enberg and Heino, 2007; Sharpe and Hendry, 2009; Dickey-Collas et al., 2010). Besides, populations that suffered from demographic erosion can become more sensitive to environmental fluctuations because of their increased dependence on young age classes and recruitment (Hsieh et al., 2006; Anderson et al., 2008; Planque et al., 2010), a process observed for various species in the western Mediterranean Sea, such as European hake (Merluccius merluccius) (Hidalgo et al., 2011, 2012).
Small pelagic fish have relatively short life-cycles, with fast growth, high mobility and plankton-base feeding. Therefore, they are highly sensitive to fluctuations in environmental conditions, including those related to human-induced climate change, and fishing (Agostini and Bakun, 2002; Checkley et al., 2009). Collapses of small pelagic fish have been observed in different ecosystems, such as the Pacific sardine (Sardinops sagax) in the California current and in the Benguela upwelling region, the Peruvian anchoveta (Engraulis ringens) in southeast Pacific Ocean or more recently the European anchovy (Engraulis encrasicolus) in the Bay of Biscay (Checkley et al., 2009; Taboada and Anadón, 2016). The collapse of these populations has been mainly attributed to a combination of fishing pressure and environmental-dependent recruitment success processes. However, the underlying processes driving population changes are not well understood, with varying environmental conditions, density-dependence processes and fishing pressure put forth as potential candidates of small pelagic fish population dynamics and changes (Silva et al., 2006; Peck et al., 2013; Mangano et al., 2020; Véron et al., 2020a).
Due to their large biomass, fast growth and key trophic position, the fluctuations in small pelagic fish populations can have ultimate ecological and socio-economic consequences. This is the case in the Mediterranean Sea, where sardine (Sardina pilchardus) and anchovy (Engraulis encrasicolus) have a long fishing history and account for almost 40% of the total catch (FAO, 2018). Moreover, both species are key components of the functioning of Mediterranean ecosystems (Palomera et al., 2007; Coll et al., 2008). However, in the last decades, declines in landings and stock biomass of sardine and anchovy have been observed in the northwestern Mediterranean Sea, accompanied by noticeable changes in life history traits, such as declines in body condition and in length at age and length at maturity observed in the Gulf of Lions (northwestern Mediterranean Sea) (Van Beveren et al., 2014; Brosset et al., 2016b, Brosset et al., 2017; Quattrocchi and Maynou, 2017). The main hypothesis proposed for the Gulf of Lions was that these changes could be explained by bottom-up control mechanisms due to a change in plankton (Saraux et al., 2019). Similarly, in other areas of the western Mediterranean Sea, the stock assessment showed that both species are overexploited and a combination of environmental and fishing pressure could be important drivers to explain the depletion of the stocks (Coll and Bellido, 2019; Coll et al., 2019; GFCM, 2019). Interestingly, the reported changes in the life history traits of the northwestern Mediterranean populations have not been observed yet in the southwestern area (Brosset et al., 2017; GFCM, 2019), suggesting that potential regional differences in life history traits of sardine and anchovy may exist across a latitudinal gradient in the western Mediterranean Sea.
This study assesses the changes in life history traits of sardine and anchovy throughout different areas in the western Mediterranean Sea; a better understanding of the small pelagic fish population’s resilience to fishing and climate change may lead to improved stock assessments and ecological analyses. Earlier studies have reported a decline in body condition for anchovy and sardine for the entire Mediterranean, with the exception of the northern Alboran Sea and Sicily channel (Brosset et al., 2017). However, general patterns to explain these changes in body condition were not found and were attributed to local environmental and anthropogenic pressure (Brosset et al., 2017). Based on this previous study, and on the fact that the stocks of sardine and anchovy have not declined homogeneously in the western Mediterranean Sea (Pennino et al., 2020), we expected life history traits of sardine and anchovy to present different interannual patterns among areas. This prediction is based on the fact that environmental and oceanographic conditions in the southwest are different from the northwestern Mediterranean (Viñas et al., 2014; Coll and Bellido, 2019). Accordingly, our main research question was: do the life history traits of sardine and anchovy show different interannual patterns in different latitudinal areas? And if so, can this be explained by differences in the local environment? By using generalized linear and additive models (GLMs and GAMs) we investigated changes in length at age, length at first maturity, reproductive activity and body condition, and explored the synchrony of life history traits with local environmental conditions to identify potential local adaptations to exogenous drivers.
Materials and MethodsStudy Area and Data Collection
The study area is located in the western Mediterranean Sea, covering the Geographical Sub-Areas (GSAs) 01 and 06 (Figure 1). This area is characterized by a narrow continental shelf in the south that broadens to the north, reaching its maximum width in the Ebro delta continental shelf. We studied the area in three separated regions, that either function as different fishing grounds or with significant different potential habitat for sardine and anchovy.
Study area with the Geographical Sub Areas 06 and 01 delimitated. Sampling ports are indicated with a dot (Tarragona-North; Torrevieja-Central; Málaga, Fuengirola-South). The shaded polygon corresponds to the area from where the environmental data were extracted (north, central, and south). The bathymetric line corresponds to the 200 m depth.
The south area, GSA01 (the Alboran Sea, 36.6976° N – −3.4513° E/36.3358° N – −3.4656° E and 36.4039° N – −5.1904° E/36.2214° N – −5.3030 E) is a highly productive area, where the influx of Atlantic Oceanic waters that are rich in nutrients form important hydrographic mesoscale features. The area is delimitated in the east by the Almeria-Oran hydrographic front (Agostini and Bakun, 2002).
The central area is the Levantine region (GSA06; 38.7313° N – 0.2205° E/38.5590° N – 0.5594° E and 37.6550° N – −0.7066° E/37.4749° N – −0.2477° E), which is delineated by the Ibiza Channel (Ramírez-Amaro et al., 2018). This area is considered a more oligotrophic area (D’Ortenzio et al., 2009; Salgado-Hernanz et al., 2019), highly influenced by the dynamics of the Balearic channels (Ibiza and Mallorca channels) (Pinot et al., 2002; Balbín et al., 2014).
The north area, the Catalan coast (GSA06; 41.3492° N – 2.1875° E/41.1284° N – 2.3694° E and 40.2842° N – 0.3520° E/39,9671° N - 1.0044° E), is influenced by a southwestward current that follows the continental shelf and is considered a highly productive area due to the river run-off of the Ebro River with a marked seasonality in the primary production due to the late winter-early spring phytoplankton bloom (Salat et al., 2002; Bosc et al., 2004).
The North and South areas are both important spawning areas for sardine and anchovy in autumn-winter and spring-summer, respectively (García and Palomera, 1996; Agostini and Bakun, 2002; Palomera et al., 2007; Bellido et al., 2008). Monthly samples were obtained from commercial landings (purse seiners) in all three areas from 2003 to 2017 (see Supplementary Table 1). In the north samples were obtained from the harbor of Tarragona, and in the central area samples were from the harbor of Torrevieja. In the southern area samples were obtained from the harbors of Málaga and Fuengirola (Figure 1). Sardines were sampled in all three areas, whereas anchovies were only sampled in the northern and the southern areas, since fishing landings of anchovy in the central area were low and in consequence biological data was not collected in the period studied. Due to fishing closure and adverse weather conditions depending on the year, data in certain months were not available. Total length (cm), total weight (g), sex (male; female; indeterminate), maturity stage by visual examination [1, immature; 2, maturing; 3, mature; 4, spawning; 5, post-spawning; 6, resting; (ICES, 2008)] and gonad weight (g) were recorded for sardine and anchovy. To estimate the length at age, otoliths were collected. The number of samples for each year and area is reported in Supplementary Table 1.
Life History Parameters
Length at age, length at first maturity, gonadosomatic and body condition index were calculated from the data collected in the monthly biological sampling. The temporal variation of these parameters for each area was analyzed using different methodologies that are described below. In the case of length at age and length at first maturity, the relationship with environmental variables was not analyzed since changes in these traits are a consequence of long-term effects experienced by individuals throughout their life (Véron et al., 2020a) and are not driven by current environmental conditions. Instead, in the case of gonadosomatic index and body condition, the temporal relationship with environmental conditions was explored since these are two biological parameters that can respond rapidly to recent changes in environmental conditions (Lloret et al., 2013). The different areas were analyzed separately since each area presents differentiated oceanographic conditions as described above.
Length at Age and Length at First Maturity
Age-length data for sardine and anchovy were obtained from otoliths collected year-round in each area of study and thus, reflect the average status of the population all year-round. An average of 40 otoliths per month were read, covering the size range of each species sampled per month following the protocol of ICES (2017).
For each species and area, the mean length at maturity, defined as 50% of individuals reaching the sexual maturity (L50), was estimated including females, males and indeterminate individuals together. Only maturity data corresponding to the reproduction period of each species was used in the determination of the L50 (October to March in the case of sardine, and April to September in the case of anchovy). Mature individuals were defined as those with maturity stages 2–6, and immature were defined as the individuals with maturity stage 1. Individuals at maturity stage 2 (maturing) of age 0, with otolith reading available, were considered immature, since they have not spawned previously. Generalized Linear Models (GLM) with a binomial error distribution and a logit link were used to estimate the length at first maturity (L50) for each area and species by year. The proportion of mature fish and the length class at 0.5 cm were the dependent and independent variables, respectively. The lack of sufficient immature individuals did not allow us to calculate the L50 in certain years and further statistics to analyze the temporal variability were not applied (Supplementary Table 1).
Gonadosomatic Index and Body Condition
To determine the reproductive activity, only females were considered for the calculation of the gonadosomatic index (GSI), since the offspring production is frequently limited to a greater degree by egg than sperm production (Murua and Saborido-Rey, 2003). GSI was calculated as a proportion of the gonad weight to the total body weight minus the gonad weight [%GSI = 100⋅Gonad weight/(Total weight – Gonad weight)].
The analysis of temporal variability on somatic body condition was based on the calculation of the relative condition factor that has been widely used for small pelagic fish (Le Cren, 1951; Brosset et al., 2015a). The relative condition factor Kn was obtained as the ratio between the weight of the fish (expressed as total weight minus the gonad weight; Wgutted) and the predicted weight (Wp) for a fish of the same length (Le Cren, 1951) [Kn = Wgutted/Wp].
The Wp was estimated applying the linear length-weight relationship log(Wgutted) = log(a) + b × log(TL) to all fish sampled in the three areas together, since we wanted to investigate differences between areas. TL is the total length, and a and b are the regression coefficients and n the sample size (sardine: a = 0.00394, b = 3.2492, n = 35362; anchovy: a = 0.00398, b = 3.1795, n = 21591). This relationship was used to calculate the somatic predicted weight (Wp = a∗TLb). The selection of Wgutted instead of total weight was made to avoid changes related to reproduction and analyze exclusively changes in somatic body condition. By definition, values above 1 represent individuals in better condition than a standard individual, while values below 1 are individuals in worse condition than a standard individual of the studied population.
Environmental Variables
Satellite-derived environmental variables, obtained from models that assimilate satellite data, were used to model the interannual variability in gonadosomatic index and body condition. Four different variables were used: Sea Temperature (ST150 in°C) averaged over the first 150 m, mean Meridional (MC in m/s, current from south to north, positive northward) and Zonal Component of the water current field (ZC in m/s, current from west to east, positive eastward) and Net Primary Productivity (NPP in mol/m3/s) (Supplementary Figure 1). These variables were selected as they have been previously seen to influence sardine and anchovy biological processes (Giannoulaki et al., 2013; Quattrocchi et al., 2016; Pennino et al., 2020; Fernandez-Corredor et al., 2021). These variables are related to food availability and changes in temperature; the two main hypothesized drivers in previous studies to explain changes in Kn and GSI (Brosset et al., 2017; Véron et al., 2020a). Sea temperature averaged over the first 150 m was selected taking into account the bathymetric distribution of both species that can occur at depths of 100 and 200 m for sardine and anchovy, respectively (Palomera, 1992; Tugores et al., 2011). Both MC and ZC are involved in different processes of prey retention (Quattrocchi et al., 2016). The MC in the north and central area is controlled by the southwestward current that follows the continental slope, which would correspond to negative values of MC (Salat, 1996). In both areas, freshwater inputs from river discharge increase the ZC. However, the local mesoscale activity produced by local wind stress and vertical mixing can have a major influence on the circulation patterns (Salat, 1996; Agostini and Bakun, 2002). In the case of the south area, the dominance of the eastward Jet Atlantic current favor positive values of ZC (Ruiz et al., 2013). Monthly averages of oceanographic variables covering the entire study period (2003–2017) were extracted from the Copernicus Mediterranean Sea biogeochemistry reanalysis (Ref: MEDSEA_REANALYSIS_BIO_006_008; Date: 27/03/2020; Teruzzi et al., 2019; Bolzon et al., 2020) and the Copernicus Mediterranean physical reanalysis (Ref: MEDSEA_ REANALYSIS_PHY_006_004; Date: 27/03/2020; Simoncelli et al., 2019).
Statistical Analysis
Temporal changes in length at age for each area were investigated by the mean of Generalized Additive Models (GAMs; Wood, 2006). Different models were tested following the procedure described in Véron et al. (2020a). First, length was modeled as a function of age class (introduced as categorical effect). Then a smooth function of the temporal co-variable year (defined as a continuous variable) was added to the previous model. Finally, the third model explored the partial effect of year in the length for each age, for that we used the “by” argument of the “gam()” function to fit a separate smooth function for each age group [i.e., s(year,by = age)], instead of one smooth function of year. In this way, we allowed the relationship between year and length to vary across age. As mentioned, age was included as a factor and year as a continuous variable through a plate regression spline. The number of knots was restricted to 4 to avoid over-fitting. GAMs were applied for the individual area, including the same predictors in each area. Each model was fitted using a Gaussian distribution validating on residuals the theoretical assumptions of normality, homoscedasticity and independence. The independence assumption of the residuals was checked plotting the residuals against the month time variable (Supplementary Figure 4), and no pattern suggesting lack of independence was identified. The partial autocorrelation at monthly scale was not used in this model since the response variable was a multiple time series, in the sense that, several observations were obtained at the same time. The best model was selected according to the lowest Akaike Information Criterion (AIC) and the highest adjusted R-square. The statistical analyses were performed using the mgcv (Wood, 2011) package, and graphs were plotted using visreg (Breheny and Burchett, 2017) ggplot2 (Wickham, 2016) package of R 4.0.3 (R Core Team, 2020).
To investigate the potential environmental effects on long-term and seasonal changes in GSI and Kn, the time series of each area were decomposed in three components: (1) the trend component that was used to analyze the long-term interannual relationship between environmental variables and Kn and GSI trends, (2) the seasonal component that was used to compare the phenology of Kn and GSI in each area in relation to the seasonal behavior of the environmental variables and (3) the residual component that was not used further in the analysis. The function “decompose” (R Core Team, 2020) that uses the classical additive seasonal decomposition by moving averages (Kendall and Stuart, 1983), was implemented to decompose each response variable (GSI and Kn) and the environmental variables (NPP, ST150m, ZC, and MC) into Trend + Seasonality + Residuals. Since the time series of GSI and Kn had missing values in certain months of specific years (Supplementary Table 1), before performing the decomposition, the package imputeTS (Moritz and Bartz-Beielstein, 2017) was used to complete the time series applying Kalman smoother imputation method. After the decomposition, imputed values were removed from the time series for further analyses. Long-term de-seasonalized time series of environmental variables are presented in Supplementary Figure 1.
An exploratory analysis of seasonal pattern was carried out, while the de-seasonalized time series of GSI and Kn were modeled with GAMs for each species and area against the de-seasonalized environmental variables (ST150, NPP, MC, and ZC). Prior to performing the analyses, standard techniques were used to identify possible correlation and collinearity between the explanatory variables (Zuur et al., 2010). In particular, correlation among variables was checked by performing a Pearson’s correlation test with the corrplot package (Wei and Simko, 2017). Collinearity was tested by computing the generalized variance-inflation factors (GVIF), which are the corrected VIF values by the number of degrees of freedom of a predictor variable (Fox and Weisberg, 2011). The GVIF was assessed using the “corvif” function in R software. Pairs of variables with high correlation values (Pearson’s correlation above 0.7) and GVIF higher than 3, were identified. In the north area, the ZC and the MC presented high correlation and GVIF higher than 3. In the central area also ZC and MC presented high correlation and GVIF higher than 3, but also ST150 and NPP (Supplementary Figure 5). In those cases, only one of the variables within correlated pairs were included in the modeling process. In the GAMs, plate regression splines basis was used for explanatory variables, restricting the dimension of the basis (k) to 4, to allow a high degree of flexibility without overfitting problems. From the GAM including all the covariates, the explanatory variables included in the final model were selected with backward stepwise procedure based on Akaike Information Criterion (AIC) and adjusted R-square (R2-adj.) (Supplementary Tables 2, 5, 8). Considering the distribution of each response variable, GSI and Kn, a Gaussian distribution and identity link function were used. Temporal correlation in the residuals was also ultimately checked using the partial ACF (Wood, 2006). For each final GAM, a diagnosis of the model assumptions was carried out, among which the independence of the residuals was verified. In some cases, a clear temporal dependence in the residuals was detected. For solving the violation of the independence assumption the following two analyses were performed (Table 1). Firstly, a new temporal explanatory variable was included through a smooth function in the final GAM, such variable was defined as a sequence from 1 to the total number of observations of the response variable with step 1, representing the month in which observation was recorded. This variable, named “Time” (T), modeled the remaining time dependence of the response time series variable avoiding correlated residuals. If the diagnosis of the residuals of the previous model was corrected, a final model was achieved. In other cases, time series models could be useful to solve the temporal dependence in the residuals. More precisely, the best autoregressive–moving-average model (ARMA), see Brockwell and Davis (2002), according to the AIC criteria, must be identified for the residuals time series. The selection of the best ARMA model using AIC criteria was done through “auto.arima” function of forecast package (Hyndman and Khandakar, 2008; Hyndman et al., 2020). Once, the values of the optimal orders of the ARMA were identified, as described previously, the argument correlation of “gamm()” function of mgcv package (Wood, 2011) allowed to adjust both models simultaneously, the GAM model with ARMA terms, this approach has been widely used (Yang et al., 2012; Holmes et al., 2021). The best (and most parsimonious) model was finally chosen based on the compromise between low AIC values, high R2-adj., and significant predictors. The code used for the analysis carried out can be found in the GitHub repository1.
Comparison between the Generalized Additive Model selected (GAM) of GSI and Kn, the GAM combined with the temporal explanatory variable (time variable) (GAM + Time) and the GAM combined with an ARMA model (corARMA), for each species (sardine and anchovy), trait (GSI and Kn) and for each area.
Area
Model type
Parameters included
R2-adj.
AIC
Sardine GSI
North
GAM
MC + ST150 + NPP
0.526
23.1918
GAM + Time
MC + ST150 + NPP + T
0.675
−11.3313
GAM + Time + ARMA
ST150
0.343
−437.9223
Central
GAM
ZC + NPP
0.496
−3.60679
GAM + Time
ZC + NPP + T
0.748
−49.31529
GAM + Time + ARMA
ZC
0.0735
−142.75485
South
GAM
ZC + ST150 + NPP
0.368
69.7838
GAM + Time
ZC + ST150 + NPP + T
0.398
73.7103
GAM + ARMA
ZC
−0.108
−354.4148
Anchovy GSI
North
GAM
MC + ST150 + NPP
0.492
179.35468
GAM + Time
MC + ST150 + NPP + T
0.922
−44.81075
GAM + Time + ARMA
−
0.724
−307.31535
South
GAM
MC + ZC
0.388
87.00971
GAM + Time
MC + ZC + T
0.447
84.70779
GAM + ARMA
MC + ZC
−0.419
−236.5413
Sardine Kn
North
GAM
MC + ST150 + NPP
0.397
−691.2924
GAM + Time
MC + ST150 + NPP + T
0.653
−765.8567
GAM + Time + ARMA
T
0.279
−1115.3969
Central
GAM
MC + NPP
0.530
−665.8308
GAM + Time
MC + NPP + T
0.691
−692.0069
GAM + Time + ARMA
MC + NPP + T
0.621
−881.1502
South
GAM
MC + ST150 + NPP
0.27
−619.5305
GAM + Time
MC + ST150 + NPP + T
0.683
−717.0292
GAM + Time + ARMA
MC + ST150 + T
0.255
−1182.2717
Anchovy Kn
North
GAM
ZC + ST150 + NPP
0.523
−763.9239
GAM + Time
ZC + ST150 + NPP + T
0.830
−847.8216
GAM + Time + ARMA
ZC + NPP
−0.0254
−1067.5924
South
GAM
MC + ZC + ST150 + NPP
0.745
−715.4204
GAM + Time
MC + ZC + ST150 + NPP + T
0.816
−740.2001
GAM + Time + ARMA
MC + ST150
0.269
−1015.6587
The parameters that were statistical significant (p-value < 0.05) are reported. Statistics acronyms are: R2-adj., coefficient of determination adjusted; AIC, Akaike Information Criterion. The model selected is highlighted in bold. The letters describing the parameters included in the model stands for: meridional current (MC), zonal currents (ZC), average sea temperature of 150 m (ST150), net primary production (NPP), and time (T).ResultsChanges in Length at AgeSardine
Sardine individuals of age classes from 0 to 8 were found, with few individuals belonging to the oldest categories (276 out of 23,269 individuals sampled, had an age of 6–8 years). Most of the older individuals were sampled in the south area (Figure 2 and Supplementary Figure 2). In the north area, there was a sharp decline in length at age during the period 2012–2017 (Figure 2 and Supplementary Figure 2). Moreover, individuals of older ages (4, 5, and 6) disappeared from the population in the north and central area.
Summary figure with annual mean and standard deviation values of total length of age 1 (see Supplementary Figures 2, 3 for a plot of all ages) (A,E), mean and standard error of length at maturity (L50) (B,F), mean and standard deviation values of the gonadosomatic index trend (values obtained after decomposition; GSI trend) (C,G) and trend of the relative condition factor (values obtained after decomposition; Kn trend) (D,H). All values are reported by species, for sardine on the left and for anchovy on the right side. Each color represents an area (North, Central, and South).
Length modeled as a function of age and a partial effect of year for each age was the best model selected (Supplementary Table 2). The GAM fitted to the sardine length-at age showed no apparent violation of the normality of the residuals or of the homogeneity of variances. Overall, comparing the models in the three areas, length at age followed a latitudinal gradient with larger fish found in the south than in the central area and the smaller fish found in the north for all ages. GAM of the north area showed that the decline in length at age was consistent through the age classes 0–4, and the model that include the covariates age and year had an adjusted R-square (R2-adj.) of 0.653 (Figure 3 and Supplementary Table 2). In the central area, an increase in length at age 0 was observed from 2013 until 2017, while in ages 1–4 there was a significant decline in length at age. The fitted model including age and year covariates had a R2-adj. of 0.669. In the south area, a continuous increase of length at age for age 0 and an increase for ages 1–4 until 2014 were observed, with a latter decline in length at age for the period 2015–2017 for ages 1–4. In this case, the model had a R2-adj. of 0.803 (Figure 3 and Supplementary Tables 2, 3).
Long-term changes in total body length at age for sardine (Sardina pilchardus) and anchovy (Engraulis encrasicolus). Partial effect of the GAMs model are represented for each species and area; sardine North (A), Central (B), South (C); anchovy North (D), south (E). Lines correspond to the fit of the GAM model. The shaded areas indicate the 95% confidence interval. Only age classes with significant effects are represented.
Anchovy
The ages of the majority of individuals of anchovy were comprised between 0 and 2 years of age and individuals of age 3 and 4 were only recorded before 2013 in the two areas sampled (Figure 2 and Supplementary Figure 3). Opposite to sardine, in anchovy, latitudinal differences in length at age were not observed.
Length modeled as a function of age and a partial effect of year for each age was the best model selected (Supplementary Table 2). The GAM fitted to the anchovy length-at age showed no apparent violation of the normality of the residuals or of the homogeneity of variances. In the north area, the age class and the smooth function of year by age class had a R2-adj. of 0.497, and a decline in length at age was observed in age 0 and 1. An increase in length at age 2 was observed from 2013 onward (Figure 3 and Supplementary Table 3). In the south, the model had a R2-adj. of 0.328. The dynamics in the south were different from the north area. Ages from 0 to 2 showed an increase in length at age. For age 3 a slight increase was observed followed by a decline after 2008 (Figure 3). In general, the fitted GAM models of anchovy had lower R2-adj. values compared to sardine (Supplementary Tables 2, 4).
Changes in Length at First Maturity (L<sub>50</sub>)Sardine
Different temporal patterns were observed in sardine L50 between areas (Figure 2). Although in the three areas there were yearly fluctuations in L50, a clear temporal decline was only observed in the north area. Consequently, in the period 2012–2017 the north area presented lower L50 than the other two areas. The average L50 in the north area in the period 2012–2017 was 10.45 ± 0.09 cm (mean ± standard error). Instead, the L50 in the central and southern areas showed larger L50 (12.41 ± 0.12 cm, and 12.39 ± 0.13 cm, respectively; Figure 2).
In the north area, for the period 2004–2011, it was only possible to calculate the L50 in 2007 and 2010, due to the lack of immature individuals in the other years. Despite the missing values, it was observed a progressive decline in L50 with the highest value of L50 observed in 2007 (12.58 ± 0.17 cm) and the lowest values observed in 2016 (9.51 ± 0.52 cm). In the central area, the L50 presented fluctuations with maximum L50 observed in 2014 (13.01 ± 0.18 cm) and minimum values in 2010 and 2012 (10.51 ± 1.19 cm and 11.22 ± 1.00 cm, respectively). Similar to the central area, the south also presented fluctuations in the L50. The highest L50 values were observed in 2007 and 2013 (13.68 ± 0.26 cm and 13.55 ± 0.20 cm, respectively) and the lowest in 2005 (10.85 ± 0.71 cm).
Anchovy
Length at first maturity of anchovy estimated per year and area showed high interannual variability in both areas. In the period 2008–2011 the L50 in the north area was clearly higher than in the south. In the last period of this study (2012–2017), the average L50 was similar between both areas (9.56 ± 0.06 cm and 9.39 ± 0.1 cm, north and south areas, respectively).
In the north area, the maximum L50 values were observed in 2011 (11.56 ± 0.08 cm), followed by three consecutive years of L50 decrease. The minimum L50 was observed in 2014 (8.66 ± 0.44 cm). In the south area, the maximum L50 was observed in 2003 and 2006 (10.72 ± 0.30 and 10.92 ± 00.16 cm, respectively) and declined, with fluctuations, until 2012, when the lowest L50 value was observed (8.22 ± 1.06 cm) (Figure 2).
Changes in Gonadosomatic Index (GSI)Sardine
The seasonal component of the decomposed gonadosomatic index (GSI) of sardine females showed maximum values from November to February in the three areas (Figure 4A). In the north the seasonal component had a lower amplitude, while the central area and the south presented similar seasonal patterns. Regarding interannual trends (Figure 2C), lower values of GSI were observed in the north compared to the central and south areas for 2004–2017. In the north area, the GSI was highest in 2005 and declined until 2011. Then, from 2011 to 2014, the GSI presented interannual fluctuations without a clear long-term trend. Instead, the other two areas, central and south area, showed fluctuations for the entire study period. The central area had a minimum value of GSI in 2011 coinciding with the minimum GSI of the north area and a maximum GSI value in 2014. In the south, the lowest GSI were observed in 2006 and 2010 and the highest in 2012 and 2017 (Figure 2C).
Seasonal component of the detrended time series of the gonadosomatic index (GSI seasonal) and relative condition factor (Kn seasonal) for sardine (A,C) and for anchovy (B,D). Environmental seasonal components for net primary production (NPP; E), Sea Temperature at 150 m (ST150m; F) and the zonal current (ZC; G) and meridional current (MC; H). Each color represents an Area (North, Central, and South).
The interannual trend of NPP and ST150 in the central area, and MC and ZC in the north and central area, presented high correlation (Pearson’s correlation above 0.7; Supplementary Figure 5), and hence only one of both explanatory variables, in particular MC variable in the north and ZC and NPP in the central area, was included in the GAM model. Therefore, all variables used in the models had a GVIF lower than 3. The environmental variables retained in the models, by the backward selection procedure, were MC, ST150, and NPP in the north and south area, and ZC and NPP in the central area (Supplementary Table 5). Results of the GAM analysis of the interannual trend of GSI showed marked differences between areas in the form of the partial effect of selected environmental variables. The GAM fitted to the sardine gonadosomatic index showed no apparent violation of the normality of the residuals or of the homogeneity of variances. Given the high autocorrelation of long-term trends time series, some partial autocorrelation at monthly scale once the seasonal cycles were removed, persisted in the residuals at the first and second order. Hence, as mentioned in section “Statistical Analysis,” the temporal explanatory variable (month variable) was included in the model for avoiding correlated residuals. The resulting model showed that for the north and central area the R2-adj., AIC and partial autocorrelation function of the residuals improved (GAM + Time; Table 1). Despite this improvement, autocorrelation in the residuals was not removed completely (Supplementary Figure 6). Hence, we carried out the second analysis, described in section “Statistical Analysis,” which models the dependence in the residuals through an ARMA model. The GAM combined with the ARMA model was able to remove the temporal dependency in the residuals and the obtained AIC was lower. However, the models were not meaningful and the R2-adj. values were very low in the three areas (Table 1). Therefore, the selected models were the GAM including the temporal variable (month variable; GAM + Time) in the north and central area. In the south the selected model was the GAM, since the addition of a temporal variable had a small improvement of the R2-adj. and the AIC increased (Table 1).
In the north area, the selected model had a R2-adj. of 0.675 (Supplementary Tables 5, 6). The GSI displayed a negative relationship with NPP and a dome shaped relationship with MC. Maximum GSI values were observed when the meridional current was close to 0. The sea temperature positively affected the GSI at sea temperatures below 15.4°C. Above this value significant negative effects were found (Figure 5A). The final model of the central area had a R2-adj. of 0.748 and included the ZC and NPP variables (Supplementary Tables 5, 6). The GSI of sardine had a positive relationship with ZC (eastward current) (Figure 5B). Unlike the north area, the GSI in the central area displayed a positive relationship with NPP. In the south area, the model selected had a R2-adj. of 0.368 lower than in the other two areas, because the temporal variable was not included in this model (Table 1 and Supplementary Table 5). In the south, GSI displayed a negative relationship with southward flowing current (corresponding to negative values of MC) and a positive relationship with northward flowing current (positive values of MC) (Supplementary Table 6). The ST150 and NPP had a positive effect on GSI of sardine (Figure 5C).
Partial GAM plots for gonadosomatic index (GSI) of sardine females (Sardina pilchardus) for each area; (A) North (purple), (B) Central (green), (C) South (yellow). Significant partial effects of meridional current (MC), zonal currents (ZC), average sea temperature of 150 m (ST-150) and net primary production (NPP) are represented. The shaded areas indicate the 95% confidence interval.
Anchovy
The seasonal component of the decomposed GSI of anchovy females showed high values from May to August in the north area and from April to August in the south area (Figures 4B,E–H). The south area presented higher maximum values than the north. The interannual trend in GSI (Figure 2F) showed similar values between areas in the period 2006–2011, after which a marked decline in the GSI of the north area was observed, while the GSI of the south was maintained. Consequently, the GSI values of the north were lower than the GSI values of the south in the period 2012–2017 (Figure 2F).
As mentioned above for sardine, MC and ZC in the north presented high correlation (Pearson’s correlation above 0.7; Supplementary Figure 5), and hence only one of both explanatory variables, in particular MC variable, was included in the GAM model. Therefore, all variables used in the models had a GVIF lower than 3. The environmental variables retained in the models, by the backward selection procedure, were MC, ST150 and NPP in the north, and MC and ZC in the south (Supplementary Table 5). The GAM fitted to the anchovy gonadosomatic index, in the north and south area, showed no apparent violation of the normality of the residuals or of the homogeneity of variances. Similar to sardine and following the same procedure described in the previous section, partial autocorrelation at monthly scale once the seasonal cycles were removed, persisted in the residuals at the first and second order. In the north area, the GAM including the temporal variable (GAM + Time) was selected. In the south, the GAM without the temporal variable was maintained as the best model, since the slight improvement in terms of the R2-adj. and the AIC criterion, when the temporal variable was included, did not compensate for the increase in the model complexity (Table 1 and Supplementary Figure 6).
The final model of the north area had a R2-adj. of 0.922 (Supplementary Tables 5, 7). High values of GSI in the north were related to MC close to 0 or positive (south to north current) and the ST150 had a positive effect. The NPP positively affected the GSI, but for high values of NPP the effect on GSI was negative (Figure 6A). The final model of the south area had a R2-adj. of 0.368, this value was lower than in the north since the temporal variable was not included in the model (Supplementary Table 7). In the final model the two components of the current were selected. The MC, south to north current, had a positive effect with increasing GSI. The ZC also had a slight positive effect for low positive values, with maximum values of GSI between 0.4 and 0.6 m/s of the west to east current (Figure 6B).
Partial GAM plots for gonadosomatic index (GSI) of anchovy female (Engraulis encrasicolus) for each area; (A) North (purple) and (B) South (yellow). Significant partial effects of meridional current (MC), zonal currents (ZC), average sea temperature of 150 m (ST-150) and net primary production (NPP) are represented. The shaded areas indicate the 95% confidence interval.
Changes in Body Condition (Kn)Sardine
The seasonal component of the decomposed body condition index (Kn) of sardine had different peaks between areas (Figure 4C). Maximum values of seasonal Kn in the north were observed from April to August, while in the central and south areas the maximum values of sardine seasonal Kn were observed from May to October. When looking at the interannual trend (Figure 2D), lower values of Kn were observed in the north and central area compared to south for the entire study period. In the north, the Kn interannual trend declined from 2006 to 2009, after which it increased from 2010 until 2013. Kn then sharply declined, followed by a recovery on the period 2015–2017. The central area showed Kn pattern similar to the north, with a decline of Kn from 2007 to 2011 and a moderate increase from 2012 to 2017 (Figure 2D). In the south area the maximum values of Kn were observed in the period 2005–2008, whereas the period with lower Kn values was from 2009 onward.
After the exploration of correlation between environmental variables explained in the section “Changes in Gonadosomatic Index (GSI),” the analysis of the relationship between interannual trend of sardine body condition and the environmental variables showed that the retained variables in the models differed between areas. In the north the environmental variables selected were MC, ST150 and NPP, in the central area MC and NPP and in the south MC and ST150 (Supplementary Table 8). The GAM fitted to the sardine body condition showed no apparent violation of the normality of the residuals or of the homogeneity of variances. Given the high autocorrelation of long-term trends time series some partial autocorrelation in the residuals persisted at the first and second order. Hence, as mentioned in sections “Statistical Analysis” and “Changes in Gonadosomatic Index (GSI),” two alternative models were applied to improve the temporal autocorrelation in the residuals. The best model attained in the north was the GAM model with the temporal variable (GAM + Time; Table 1) since the inclusion of an ARMA model on the residuals leads to a model where all environmental variables were not significant, and consequently to a substantial decrease in the R2-adj value. In the central and south area the selected models were the GAM model combined with ARMA model on the residuals (GAM + Time + ARMA; Table 1 and Supplementary Figure 7). This model was selected in the south area in spite of its lower R2-adj. value since previous models showed a violation of the independence assumption on the residuals whereas the GAM + Time + ARMA attained a valid model removing such residuals pattern.
In the north area, from the correlated variables MC and ZC, only MC was included in the GAM model (Supplementary Table 8). The selected model had a R2-adj. of 0.653 (Supplementary Table 9). The MC had a negative effect on Kn at negative values (negative values correspond to north to south current) and positive effect for low current (values close to 0). The ST150 and NPP had a positive effect on Kn (Figure 7A). In the central area from the pair correlated environmental variables, MC and NPP were included. The final GAM model combined with ARMA model had a R2-adj. of 0.621 (Supplementary Table 9). Kn of sardine in the central area decreased with higher values of northward current (MC) and the NPP had a clear positive effect on Kn (Figure 7B). In the south area the final GAM with a ARMA model had a R2-adj. of 0.255. The MC had a positive relationship with Kn, and the ST150 had a negative effect on Kn (Figure 7C).
Partial GAM plots for body condition index (Kn) of sardine (Sardina pilchardus) for each area; (A) North (purple), (B) Central (green), and (C) South (yellow). Significant partial effects of meridional current (MC), zonal currents (ZC), average sea temperature of 150 m (ST150) and net primary production (NPP), are represented. The shaded areas indicate the 95% confidence interval.
Anchovy
When comparing the north and south area, the seasonal component of the decomposed body condition index (Kn) presented different seasonal pattern (Figure 4D). While in the north there was a clear peak in Kn of anchovy in April with a continued decline of Kn from May to December, in the south area the Kn of anchovy was maintained with similar Kn values from March until October and with monthly variability. When looking at the interannual trend of anchovy Kn (Figure 2G), lower values of Kn were observed in the north compared to the south area for the entire period studied (2004–2017), with the exception of 2015 that presented extremely low Kn values and was similar to the values observed in the north. The Kn in the north declined until 2012 and in the period 2012 to 2017 values of anchovy Kn in the north were lower than in the period 2004–2011.
In the interannual analysis of the Kn trends the four environmental variables were retained in the south and in the north since MC and ZC had a high correlation, only ZC was included in the GAM model. The GAMs fitted to the anchovy body condition, in the north and south areas, showed no apparent violation of the normality of the residuals or of the homogeneity of variances. Partial autocorrelation at monthly scale once the seasonal cycles were removed, persisted in the residuals at the first and second order. After adding to the GAM the explanatory month variable (GAM + Time), in the north and south area the R2-adj., AIC and partial autocorrelation of the residuals improved (Table 1 and Supplementary Figure 7). The GAM combined with the ARMA model was able to remove the temporal dependency in the residuals and the obtained AIC was lower. However, the models were not meaningful and the R2-adj. values were very low in the two areas compared to the GAM and GAM + Time models (Table 1). Therefore, the selected models were the GAM including the temporal variable (GAM + Time) in both areas (Table 1 and Supplementary Tables 9, 10).
The final GAM model of Kn for anchovy in the north had a R2-adj. of 0.83 (Supplementary Table 10). The ZC had a positive effect at negative values of current (east to west current), and a negative effect at positive values of ZC (west to east current). The ST150 and NPP had a positive effect on Kn. In the case of sea temperature, the highest values of Kn of anchovy were observed between sea temperatures of 15.8 and 16.0°C (Figure 8A). In the south area, the final model of anchovy Kn in the south had a R2-adj. of 0.816 (Supplementary Table 10). The MC had a positive effect until 0.010 m/s, above which the effect on Kn was negative. Similarly, the ZC had a positive effect on Kn for west to east currents up to 0.08 m/s. Both, ST150m and NPP had a clear positive effect on Kn of anchovy in the south (Figure 8B).
Partial GAM plots for body condition index (Kn) of anchovy (Engraulis encrasicolus) for each area; (A) North (purple) and (B) South (yellow). Significant partial effects of meridional current (MC), zonal currents (ZC), average sea temperature of 150 m (ST-150) and net primary production (NPP), are represented. The shaded areas indicate the 95% confidence interval.
Discussion
The overall aim of this study was to assess changes in life history traits of sardine and anchovy in three areas of the western Mediterranean Sea. Differentiated long-term patterns were observed between northern and southern areas in the traits studied in both species (i.e., length at age, length at maturity, GSI, and Kn) with major long-term declines observed in the north (Table 2). The environmental conditions partially explained the changes in GSI and Kn, but the selected variables differed between areas, highlighting the importance of regional oceanographic conditions to understand the dynamics of small pelagic fish (Table 3). Although the overall seasonal phenology of the gonadosomatic index and body condition was similar to those described in previous studies for sardine and anchovy, monthly differences between areas were observed in both parameters, demonstrating that latitudinal differences exist at seasonal and interannual level. Thus, management plans of sardine and anchovy should take into consideration both large scale and regional dynamics.
Summary table of the main trends observed for each life history trait studied; length at age, length at maturity (L50), gonadosomatic index (GSI) and body condition index (Kn), for sardine (Sardina pilchardus) and anchovy (Engraulis encrasicolus) in the western Mediterranean Sea.
For each species, the latitudinal differences between north “N,” central “C,” and south “S” areas, and the temporal pattern for each area separately, are described. Decreasing trends are indicated with red and the combination of increasing and decreasing trends or fluctuations are indicated in yellow.
Summary table of the relationship of gonadosomatic index (GSI) and body condition index (Kn) of sardine (Sardina pilchardus) and anchovy (Engraulis encrasicolus) with environmental variables (NPP, net primary production; ST150, average sea temperature 150 m; MC, meridional current; ZC, zonal current).
Decrease or negative (−) relationship is indicated with red, increase or positive relationship (+) is indicated in green and (∩) dome shaped and (∪) U shaped relationships are indicated in yellow.
Looking at the spatial differences in life history traits, the north and south area presented clear differences for both species, with higher GSI and Kn values in the south area and in the case of sardines, higher length at age and L50 as well. Spatial differences in L50 within the same species are common (Stahl and Kruse, 2008). For example, in the case of sardine, Silva et al. (2006) described lower size at first maturity in the Gulf of Cadiz than the northeastern Atlantic. Also, previous works described a higher growth for anchovy populations and a higher length at age for sardine in the Alboran Sea compared to the Ebro delta area (Alemany and Álvarez, 1993; Ventero et al., 2017). These differences between areas were linked to the existence of two genetically differentiated stocks that are separated by the Almeria-Oran front (Viñas et al., 2014; Pascual et al., 2017; Coll and Bellido, 2019). In the northwestern Mediterranean larval dispersal between the Gulf of Lions and the Ebro Delta has been described for anchovy. The central area of our study has received less attention and studies on larval dispersion are not available, but certain connectivity between north and central area of study could be expected due to the effect of the Northern Current (Ospina-Alvarez et al., 2015). In other studies, better feeding conditions (quality and quantity) in anchovy and sardine larvae were proposed as the driver to explain the faster growth observed in the Alboran Sea compared to the northwestern Mediterranean Sea (García et al., 2003; Mercado et al., 2007; Quintanilla et al., 2015). As the south area presented higher net primary production levels than the central and north area, apart from genetic differences between populations, differences in length at age for adult stages of sardine and anchovy could be influenced by the early life stages growth and condition related to plankton dynamics. Similarly, the higher Kn and GSI values observed in the south could also be explained by the higher NPP observed in the south area.
Differences between areas in food (energy) availability or timing may influence the energy acquisition and allocation patterns and ultimately affect the trade-off between maintenance, reproduction and growth of sardine and anchovy. The seasonal cycles of life history traits play a key role in understanding these trade-offs. When exploring the relationship between seasonal patterns in GSI, body condition and environmental variables, we observed that for sardine and anchovy in the north, the body condition peaked and declined earlier than in the central and south area. Similar to our results, Giráldez and Abad (1995) also described different seasonal peaks of body condition in the south for anchovy. These differences in the seasonal pattern of body condition might play a key role in the dynamics of the reproduction of sardine, as this species accumulates energy during spring and summer for reproduction in autumn and winter (Albo-Puigserver et al., 2020). Low body condition before the start of the reproduction in October could have affected the quality of the spawning and ultimately the recruitment success in the northern area, as has been suggested for sardine off the coast of Portugal (Garrido et al., 2007). Instead, the delayed peak in body condition of the central area could explain the maintenance of higher GSI values. Surprisingly, the seasonal Kn pattern of anchovy was highly different between areas, with a clear decline of the body condition after the peak in April in the north. Anchovy reproduces in spring and summer and obtains energy from feeding that is directly used for reproduction (Ganias, 2009; Giráldez, 2009). Seasonal variability in the body condition of anchovy tends to be less marked than in sardine (Albo-Puigserver et al., 2020), but the marked seasonality of the body condition of anchovy observed in the north suggest potential food limitations and the use of reserves for reproduction. Seasonal environmental patterns could partially explain differences in the phenology of life history traits. In the north area NPP depends mainly on nutrients from river discharge and local wind stress (Lloret et al., 2004; Martín et al., 2008). In recent decades, the discharges of the main rivers, such as the Ebro River, in the western Mediterranean basin, have decreased which could have been a driver of the altered NPP and plankton dynamics in the north, reducing energy availability for higher trophic levels (Lloret et al., 2004; Ludwig et al., 2009). This could explain the lower body conditions observed in the north before the reproduction period in the case of sardine (Brosset et al., 2016b; Albo-Puigserver et al., 2020). In contrast, higher NPP values were observed in the south than in the north and central areas. Higher NPP could thus favor a more energetic and abundant zooplankton production, and in consequence, a higher accumulation of fat reserves in sardine and anchovy during spring that would allow maintaining better body conditions during summer (Dessier et al., 2018). These differences in seasonal environmental conditions, in combination with genetic differences, could explain different phenotypic adaptations in each area.
Temporal declines in the life history traits of sardine and anchovy were observed in the studied period. However, the declines were not homogenously observed in all the areas and a latitudinal gradient was apparent. In the north, there was a synchronized temporal decline in GSI and body condition for both species that was not observed in the south area. In the case of sardine, the decline in the north started in 2006, while in anchovy the major changes in GSI and Kn appeared in 2012 coinciding with the decline in L50. Also, the decline in length at age was more pronounced in the north. In the case of sardine, older age classes disappeared and individuals of up to 7 and 8 years of age were only showed in the south in the entire study period. Previous studies showed that sardine of 7 and 8 years of age were present in the 1980’s in the northwestern Mediterranean (Morales-Nin and Pertierra, 1990). Therefore, older age classes of sardine in the north must have disappeared before 2004. These changes in the population age structure may have impacted other life history parameters; for instance, reproduction output has been directly related to the age and size structure of the population (Barneche et al., 2018). In our study, the gonadosomatic index of sardine and anchovy in the north presented a temporal decline, probably linked to the absence of older age classes.
Similar to the decline in length at age in sardines from the north, a decline in length at maturity was also observed. While in 2007 the L50 of sardine was 12.58 cm, individuals of the period 2012–2017 presented a L50 more than 2 cm smaller than the L50 in the central and southernmost areas. The decline in L50 to levels below 10 cm and the disappearance of older individuals and the slower growth has also been reported in the Gulf of Lions (northwestern Mediterranean Sea) for sardine and anchovy (Van Beveren et al., 2014; Brosset et al., 2016b). A similar pattern was also apparent in the central area of our study but was less pronounced and without an important decline in L50. Therefore, the Gulf of Lions (GSA07) and the northern area of our study (north-GSA06) presented more similar changes in L50 and length at age than the central area for sardine. Important observed declines in L50 are usually explained by unfavorable environmental conditions that reduce the growth or by a high fishing pressure (or a combination of both) as a mechanism to maximize the fitness (Basilone et al., 2017). These mechanisms have been previously reported in short-lived species with high plasticity, like small pelagic fish (Silva et al., 2006; Neuheimer and Grønkjær, 2012; Hunter et al., 2015; Brosset et al., 2016b). Somatic growth is a more plastic trait compared to maturation schedules, and then easily used to adapt L50, particularly in heterogeneous environments (Hidalgo et al., 2014; Véron et al., 2020b). Fishing pressure has also been widely reported as a driver for the selection of larger individuals, with an effect resulting in spawners of smaller sizes that would affect the trade-offs in the energy allocation between reproduction and growth (Jørgensen et al., 2007; Hidalgo et al., 2011). In the western Mediterranean Sea, there is a long history of small pelagic fish exploitation (Pertierra and Lleonart, 1996; Van Beveren et al., 2016), and high fishing effort has been maintained during the last decades (Coll and Bellido, 2019). Moreover, the minimum landing size of 9 cm for anchovy and 11 cm for sardine were below the historical lengths at first maturity for both species (Palomera et al., 2007). In a recent study, the combination of climate impacts and fishing pressure on sardine and anchovy populations have been described to be more severe in the northwestern Mediterranean Sea (corresponding to the north area of our study and the Gulf of Lions; Ramírez et al., 2021) than to the central and southern areas of our study. It is thus plausible that the high fishing pressure in the north area has been an important driver behind the observed age structure truncation. This result is of high relevance as most of the previous studies mainly reported age-truncation in species with broad demography but there are few studies related to small pelagic fish (but see, Van Beveren et al., 2014; Brosset et al., 2016b).
When analyzing the overall long-term patterns with environmental variables, the net primary production (NPP) had a significant effect on the GSI and Kn for sardine and anchovy except for GSI of anchovy and Kn of sardine in the south. The positive effect of NPP with sardine Kn in the north and central area and anchovy in the north and south is probably related to higher food availability that could improve the energy reserves and the reproduction output. Specifically, for 2015 in the south, we observed that the lowest Kn values of the time series for sardine and anchovy coincided with an abrupt decline in NPP. Similar to our results, Basilone et al. (2006) found a positive relationship of Chl-a (as a proxy of primary production) with the body condition and GSI of anchovy in the Strait of Sicily, which was confirmed by a mechanistic energetic model (Mangano et al., 2020). Accordingly, previous studies highlighted that the lower quality or smaller food size could be driving the declines in sardine and anchovy in the northwestern Mediterranean Sea (Brosset et al., 2016a; Queiros et al., 2019; Feuilloley et al., 2020) and latitudinal trend in the diet of sardine and anchovy in the western Mediterranean Sea, with lower energetic prey present in the northern latitudes, have been recently reported (Bachiller et al., 2020). Surprisingly, in the north, the NPP had a negative effect on GSI of sardine that could be explained by more complex interactions, including time-lag effects, changes in quality of food or changes in the energy allocation. In the Gulf of Lions, already Brosset et al. (2016b) suggested that small pelagic fish under poor environmental conditions might modify the allocation trade-off between reproduction and maintenance.
In the case of sardine, sea surface temperature (SST) has been described as an important driver limiting the spawning activity (Ganias, 2009) and different authors have proposed the hypothesis that with the increase of SST the spawning window of sardine would be retracted (Palomera et al., 2007). At long-term, the effect of sea temperature of 150 m depth averaged (ST150) presented contrasting effects between areas and species. Since sardine spawns in winter, with a preference for colder waters (Palomera et al., 2007; Ramírez et al., 2021), we expected a negative relationship between ST150 and GSI or Kn. However, opposite to what we expected, for sardine in the north the ST150 did not present a strong negative effect on GSI, and had a positive effect on Kn of sardine. Warmer waters have been related to more stratified and less enriched waters while years with cooler sea temperatures can be indicative of other processes, such as wind mixing, river runoff and upwelling that could favor appropriate conditions for sardine and anchovy (Bellido et al., 2008; Quattrocchi et al., 2016; Fernandez-Corredor et al., 2021). Accordingly, we observed that the hydrological conditions of each area also had an effect in the GSI and Kn of both species. For anchovy, ST150 had a positive effect on GSI in the north, that could be related with an increasing number of days with optimal temperatures for their spawning, as previously described in the Bay of Biscay (Erauskin-Extramiana et al., 2019).
The meridional and zonal current are parameters involved in the prey retention and transport (Sabatés et al., 2007) and in the generation of mesoscale anticyclonic eddies that increase food availability (Bakun, 2006; Quattrocchi and Maynou, 2017). Depending on the species and trait the effect of the components of the current differed, probably due to the major complexity in the interactions between Kn and GSI and the eddy kinetic processes since the north and central area had local different process compared with the south area (Macías et al., 2014; Lorente et al., 2015). The Kn of sardine in the central area was negatively related to MC. Sardine in the central area had higher body conditions during periods with a north to south current, which could be related with increased inflow of main current in the Levantine area (Northern Current) that flows from north to south (Font, 1990; Font et al., 1990), instead a positive relationship was observed for sardine landings per unit effort in the north (Quattrocchi and Maynou, 2017). In the south, the long-term change in GSI of anchovy was exclusively explained by the MC and ZC. The south area has highly differentiated oceanographic conditions from the other areas. In the south area, there is an enrichment processes where the Atlantic Jet current increases the kinetic energy in the area creating upwelling events with an increase in primary productivity during spring and summer (Ruiz et al., 2013; Macías et al., 2014). Accordingly, the meridional and zonal currents partially explain the variability in GSI and Kn of both sardine and anchovy in the south.
Future Directions and Implications for Management of Small Pelagic Fish
The life history traits of sardine and anchovy investigated here showed temporal and latitudinal differences within our study area. Overall, we observed that the changes in life history traits in the north, and to a lesser extent, in the central area were partially similar to the changes described for the same period and both species in the Gulf of Lions (Van Beveren et al., 2014; Brosset et al., 2015b, 2016b; Saraux et al., 2019). However, the patterns were different in the southernmost area, the Alboran Sea. Differences in life history traits in the south were expected since the Alboran Sea is a productive area with well-differentiated conditions from the central and north areas (Ruiz et al., 2013; García-Martínez et al., 2019). Although global change also affects the southern area, a plausible explanation could be the genetic differentiation between populations (Coll and Bellido, 2019) with a higher environmental heterogeneity and variability, and consequently more plastic traits in the south than in the north, as previously suggested by Viñas et al. (2014) for anchovy populations.
For sardine in the northern area of this study, the observed declines in length at age, length at first maturity and age structure are critical. The implications of the disappearance of older individuals and truncation of the age-structure, in terms of the capacity of populations to buffer environmental events and in the reproductive capacity, are mechanisms that have been well studied and defined in previous studies (Hsieh et al., 2006; Barneche et al., 2018). Therefore, there is an urgent need to manage this resource proactively and precautionary to promote the recovery of the age-structure of the sardine population in the north and the central areas of the Western Mediterranean Sea and ultimately increase the resilience of the population to environmental extreme events and further changes in environmental conditions. Interestingly, a recent study has identified regions within these areas as potential “future climate refugees” for sardine and anchovy (Pennino et al., 2020) that are candidate areas to consider for future management interventions (Ramírez et al., 2021).
Our study showed clear spatially distinct changes in the life history traits of anchovy and sardine populations in the western Mediterranean Sea. A formal consideration of these differences could alter the current perception of the stock status from the stock assessment arena. Consequently, definition of the stock status, ecological analyses and management advice, in general, could benefit from taking into account these heterogeneously distributed differences. Currently, stock assessment evaluations of sardine and anchovy do not consider the north and central areas (GSA-06) differently. Although both areas are genetically very similar, we observed how input parameters that are important in stock assessment, such as length at age and length at first maturity of sardine, differed between areas. Therefore, the life history traits of sardine described in this study should be taken into consideration when evaluating the stock status but further research is needed to identify if these differences between areas will persist in the future. Life history traits that rapidly respond to environmental changes, such as body condition, might be good indicators to anticipate further declines in populations as has been previously highlighted by other studies (Lloret et al., 2012; Véron et al., 2020a). Most notably, Vargas-Yáñez et al. (2020) found that in the Alboran Sea the body condition of spawning sardine, which is closely linked to food availability in a preceding year, is the main success factor for recruitment. Therefore, the recent decline in body condition that has been observed for sardine and anchovy in the Alboran Sea should be heeded as a possible first sign for future declines of these species in the southernmost area of this study, with a future risk to show similar declines as observed in the north. This issue requires close monitoring and evaluation.
Data Availability Statement
The datasets presented in this article are not readily available because the dataset is available by request to the Spanish Institute of Oceanography. Requests to access the datasets should be directed to webmaster@ieo.es.
Ethics Statement
Ethical review and approval was not required for the animal study because the biological samples were obtained from landings of commercial fishing vessels.
Author Contributions
MA-P, MP, MH, JB, and MC designed the study. AC, AG, and PT performed the data collection and created the database. JS performed the environmental data analysis. MA-P, MP, JG, MH, MC-R, and MC performed and analyzed the data. MA-P wrote the first version of the manuscript assisted by MC and MP. All authors contributed significantly to revisions and improvements of the final manuscript.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Funding. This study was carried out within the Spanish Research project PELWEB (CTM2017-88939-R) funded by Spanish Ministry of Science, Innovation and Universities and the European Research Contract SPELMED (SC NR 02-TENDER EASME/EMFF/2016/032XXX) funded by EC EASME. Fisheries data collection has been co-funded by the EU through the European Maritime and Fisheries Fund (EMFF) within the National Program of collection, management and use of data in the fisheries sector and support for scientific advice regarding the Common Fisheries Policy (Regulation, EU 2017/1004).
We wish to thank all the participants in the fisheries data collection through all these years in the different harbors along the Spanish Mediterranean coast. We are also grateful to Marc Farré and Federico Quattrocchi for their useful advice on the analysis of data, and reviewers for their constructive comments and suggestions. This study has been conducted using E.U. Copernicus Marine Service Information. MP thank the IMPRESS Project (RTI2018-099868-B-I00), ERDF, Ministry of Science, Innovation and Universities – State Research Agency.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.570354/full#supplementary-material
ReferencesAgostiniV. N.BakunA. (2002). ‘Oceantriads’in the Mediterranean Sea: physical mechanisms potentially structuring reproductive habitat suitability (with example application to European anchovy, Engrauli sencrasicolus).Fish. Oceanogr.11129–142. 10.1046/j.1365-2419.2002.00201.xAlbo-PuigserverM.SánchezS.CollM.BernalM.Sáez-LianteR.NavarroJ. (2020). Year-round energy dynamics of sardine and anchovy in the north-western Mediterranean Sea.Mar. Environ. Res.159:105021. 10.1016/j.marenvres.2020.10502132662445AlemanyF.ÁlvarezF. (1993). Growth differences among sardine (Sardina pilchardus Walb.) populations in Western Mediterranean.Sci. Mar.57229–234.AndersonC. N. K.HsiehC. H.SandinS. A.HewittR.HollowedA.BeddingtonJ. (2008). Why fishing magnifies fluctuations in fish abundance.Nature452835–839. 10.1038/nature0685118421346BachillerE.Albo-PuigserverM.GiménezJ.PenninoM. G.Marí-MenaN.EstebanA. (2020). A trophic latitudinal gradient revealed in anchovy and sardine from the Northwestern Mediterranean Sea using a multi-proxy approach.Sci. Rep.10:17598. 10.1038/s41598-020-74602-y33077761BakunA. (2006). Wasp-waist populations and marine ecosystem dynamics: navigating the “predator pit” topographies.Prog. Oceanogr.68271–288. 10.1016/J.POCEAN.2006.02.004BalbínR.López-JuradoJ. L.FlexasM. M.RegleroP.Vélez-VelchíP.González-PolaC. (2014). Interannual variability of the early summer circulation around the Balearic Islands: driving factors and potential effects on the marine ecosystem.J. Mar. Syst.13870–81. 10.1016/j.jmarsys.2013.07.004BarnecheD. R.WhiteC. R.MarshallD. J. (2018). Fish reproductive-energy output increases disproportionately with body size.Science360642–645.BasiloneG.GuisandeC.PattiB.MazzolaS.CuttittaA.BonannoA. (2006). Effect of habitat conditions on reproduction of the European anchovy (Engraulis encrasicolus) in the Strait of Sicily.Fish. Oceanogr.15271–280. 10.1111/j.1365-2419.2005.00391.xBasiloneG.ManganoS.PulizziM.FontanaI.GiacaloneG.FerreriR. (2017). European anchovy (Engraulis encrasicolus) age structure and growth rate in two contrasted areas of the Mediterranean Sea: the paradox of faster growth in oligotrophic seas.Mediterr. Mar. Sci.15739–752. 10.12681/mms.2059BellidoJ. M.BrownA. M.ValavanisV. D.GiráldezA.PierceG. J.IglesiasM. (2008). Identifying essential fish habitat for small pelagic species in Spanish Mediterranean waters.Hydrobiologia612171–184. 10.1007/s10750-008-9481-2BolzonG.CossariniG.LazzariP.SalonS.TeruzziA.FeudaleL. (2020). “Mediterranean sea biogeochemical analysis and forecast (CMEMS MED-biogeochemistry 2018-2021),” in Copernicus Monitoring Environment Marine Service (CMEMS).10.25423/CMCC/MEDSEA_ANALYSIS_FORECAST_BIO_006_014_MEDBFM3(accessed July 25, 2018).BoscE.BricaudA.AntoineD. (2004). Seasonal and interannual variability in algal biomass and primary production in the Mediterranean Sea, as derived from 4 years of SeaWiFS observations.Global Biogeochem. Cycl.18:GB1005. 10.1029/2003gb002034BrehenyP.BurchettW. (2017). Visualization of Regression Models Using Visreg.R J.956–71.BrockwellP. J.DavisR. A. (2002). Introduction to Time Series and Forecasting, 2 Edn.Berlin: Springer.BrossetP.FromentinJ.-M.MénardF.PernetF.BourdeixJ.-H.BigotJ. L. (2015a). Measurement and analysis of small pelagic fish condition: a suitable method for rapid evaluation in the field.J. Exp. Mar. Bio. Ecol.46290–97. 10.1016/j.jembe.2014.10.016BrossetP.MénardF.FromentinJ.BonhommeauS.UlsesC.BourdeixJ. (2015b). Influence of environmental variability and age on the body condition of small pelagic fish in the Gulf of Lions.Mar. Ecol. Prog. Ser.529219–231. 10.3354/meps11275BrossetP.FromentinJ. M.Van BeverenE.LloretJ.MarquesV.BasiloneG. (2017). Spatio-temporal patterns and environmental controls of small pelagic fish body condition from contrasted Mediterranean areas.Prog. Oceanogr.151149–162. 10.1016/j.pocean.2016.12.002BrossetP.Le BourgB.CostalagoD.BǎnaruD.Van BeverenE.BourdeixJ. H. (2016a). Linking small pelagic dietary shifts with ecosystem changes in the Gulf of Lions.Mar. Ecol. Prog. Ser.554157–171. 10.3354/meps11796BrossetP.LloretJ.MuñozM.FauvelC.BeverenE.Van (2016b). Body reserves mediate trade - offs between life - history traits: new insights from small pelagic fish reproduction.R. Soc. Open Sci.3:160202. 10.1098/rsos.16020227853538CheckleyD.AlheitJ.OozekiY.RoyC. (eds) (2009). Climate Change and Small Pelagic Fish.Cambridge: Cambridge University Press.CollM.Albo-PuigserverM.NavarroJ.PalomeraI.DambacherJ. M. (2019). Who is to blame? Plausible pressures on small pelagic fish population changes in the northwestern Mediterranean Sea.Mar. Ecol. Prog. Ser.617277–294. 10.3354/meps12591CollM.BellidoJ. M. (2019). Evaluation of the Population Status and Specific Management Alternatives for the Small Pelagic Fish Stocks in the Northwestern Mediterranean Sea (SPELMED) - Final Report SC NR 02. TENDER EASME/EMFF/2016/32 – SPELMED.Luxembourg: Publications Office of the European Union.CollM.PalomeraI.TudelaS.DowdM. (2008). Food-web dynamics in the South Catalan Sea ecosystem (NW Mediterranean) for 1978–2003.Ecol. Model.21795–116. 10.1016/J.ECOLMODEL.2008.06.013De RoosA. M.BoukalD. S.PerssonL. (2006). Evolutionary regime shifts in age and size at maturation of exploited fish stocks.Proc. R. Soc. B Biol. Sci.2731873–1880. 10.1098/rspb.2006.351816822746De RoosA. M.PerssonL.McCauleyE. (2003). The influence of size-dependent life-history traits on the structure and dynamics of populations and communities.Ecol. Lett.6473–487. 10.1046/j.1461-0248.2003.00458.xDessierA.DupuyC.KerricA.MornetF.AuthierM.BustamanteP. (2018). Variability of energy density among mesozooplankton community: new insights in functional diversity to forage fish.Prog. Oceanogr.166121–128. 10.1016/j.pocean.2017.10.009Dickey-CollasM.NashR. D. M.BrunelT.van DammeC. J. G.MarshallC. T.PayneM. R. (2010). Lessons learned from stock collapse and recovery of North Sea herring: a review.ICES J. Mar. Sci.671875–1886. 10.1093/icesjms/fsq033D’OrtenzioF.AlcalM. R.BiologicaO.DohrnS. Z. A.ComunaleV. (2009). On the trophic regimes of the Mediterranean Sea: a satellite analysis.Biogeosci. Discuss52959–2983. 10.5194/bgd-5-2959-2008EnbergK.HeinoM. (2007). “Fisheries-induced Life History Changes in Herring (Clupea harengus),” in ICES CM Documents 2007 - ICES Annual Science Conference.Laxenburg: International Institute for Applied Systems Analysis.Erauskin-ExtramianaM.AlvarezP.ArrizabalagaH.IbaibarriagaL.UriarteA.CotanoU. (2019). Historical trends and future distribution of anchovy spawning in the Bay of Biscay.Deep. Res. Part II Top. Stud. Oceanogr.159169–182. 10.1016/j.dsr2.2018.07.007ErnandeB.DieckmannU.HeinoM. (2004). Adaptive changes in harvested populations: plasticity and evolution of age and size at maturation.Proc. R. Soc. B Biol. Sci.271415–423. 10.1098/rspb.2003.251915101701FAO. (2018). FAO Yearbook: Fisheries and Aquaculture Statistics 2016.Rome: FAO.Fernandez-CorredorE.Albo-PuigserverM.PenninoM. G.BellidoJ. M.CollM. (2021). Influence of environmental factors on different life stages of European anchovy (Engraulis encrasicolus) and European sardine (Sardina pilchardus) from the Mediterranean Sea: a literature review.Estuar. Coast. Shelf Sci.41:101606. 10.1016/j.rsma.2020.101606FeuilloleyG.FromentinJ.-M.StemmannL.DemarcqH.EstournelC.SarauxC. (2020). Concomitant changes in the Environment and small pelagic fish community of the Gulf of Lions.Prog. Oceanogr.186:102375. 10.1016/j.pocean.2020.102375FontJ. (1990). A comparison of seasonal winds with currents on the continental slope of the Catalan Sea (northwestern Mediterranean).J. Geophys. Res.951537–1545. 10.1029/JC095iC02p01537FontJ.SalatJ.JuliàA. (1990). Marine circulation along the Ebro continental margin.Mar. Geol.95165–177. 10.1016/0025-3227(90)90114-YFoxJ.WeisbergS. (2011). An R Companion to Applied Regression.Thousands Oaks: Sage Publications.GaniasK. (2009). Linking sardine spawning dynamics to environmental variability.Estuar. Coast. Shelf Sci.84402–408. 10.1016/j.ecss.2009.07.004GarcíaA.CortésD.RamírezT.GiráldezA.CarpenaÁ (2003). Contribution of larval growth rate variability to the recruitment of the Bay Málaga anchovy (SW Mediterranean) during the 2000-2001 spawning seasons.Sci. Mar.67477–490. 10.3989/scimar.2003.67n447733311142GarcíaA.PalomeraI. (1996). Anchovy early life history and its relation to its surrounding environment in the western Mediterranean basin.Sci. Mar.60155–166.García-MartínezM.delC.Vargas-YáñezM.MoyaF.SantiagoR.MuñozM. (2019). Average nutrient and chlorophyll distributions in the Western Mediterranean: RADMED project.Oceanologia61143–169. 10.1016/j.oceano.2018.08.003GarridoS.RosaR.Ben-HamadouR.CunhaM. E.ChícharoM. A.van der LingenC. D. (2007). Effect of maternal fat reserves on the fatty acid composition of sardine (Sardina pilchardus) oocytes.Comp. Biochem. Physiol. B. Biochem. Mol. Biol.148398–409. 10.1016/j.cbpb.2007.07.00817720577GFCM. (2019). Working Group on Stock Assessment of Small Pelagic Species (WGSASP). Scientific Advisory Committee on Fisheries (SAC). 9–14 December 2019.Rome: FAO.GiannoulakiM.IglesiasM.TugoresM. P.BonannoA.PattiB.De FeliceA. (2013). Characterizing the potential habitat of European anchovy Engraulis encrasicolus in the Mediterranean Sea, at different life stages.Fish. Oceanogr.2269–89. 10.1111/fog.12005GiráldezA. (2009). Estudio de la variabilidad temporal de los parámetros reproductivos del boquerón (Engraulis encrasicolus L.) en el mar de Alborán.Spain: Instituto Español de oceanografía.GiráldezA.AbadR. (1995). Aspects on the reproductive biology of the Western Mediterranean anchovy from the coasts of Málaga (Alborán Sea).Sci. Mar.5915–23.HalpernB. S.FrazierM.PotapenkoJ.CaseyK. S.KoenigK.LongoC. (2015). Spatial and temporal changes in cumulative human impacts on the world’s ocean.Nat. Commun.6:7615. 10.1038/ncomms861526172980HidalgoM.OlsenE. M.OhlbergerJ.Saborido-ReyF.MuruaH.PiñeiroC. (2014). Contrasting evolutionary demography induced by fishing: the role of adaptive phenotypic plasticity.Ecol. Appl.241101–1114. 10.1890/12-1777.1HidalgoM.RouyerT.BartolinoV.CerviñoS.CiannelliL.MassutíE. (2012). Context-dependent interplays between truncated demographies and climate variation shape the population growth rate of a harvested species.Ecography35637–649. 10.1111/j.1600-0587.2011.07314.xHidalgoM.RouyerT.MolineroJ. C.MassutíE.MorantaJ.GuijarroB. (2011). Synergistic effects of fishing-induced demographic changes and climate variation on fish population dynamics.Mar. Ecol. Prog. Ser.4261–12. 10.3354/meps09077HolmesE. E.ScheuerellM. D.WardE. J. (2021). Applied Time Series Analysis for Fisheries and Environmental Data.Maryland: Northwest Fisheries Science Center(NOAA).HsiehC. H.ReissC. S.HunterJ. R.BeddingtonJ. R.MayR. M.SugiharaG. (2006). Fishing elevates variability in the abundance of exploited species.Nature443859–862. 10.1038/nature0523217051218HunterA.SpeirsD. C.HeathM. R. (2015). Fishery-induced changes to age and length dependent maturation schedules of three demersal fish species in the Firth of Clyde.Fish. Res.17014–23. 10.1016/j.fishres.2015.05.004HyndmanR.AthanasopoulosG.BergmeirC.CaceresG.ChhayL.O’Hara-WildM. (2020). Forecast: Forecasting Functions for Time Series and Linear Models. R Package Version 8.13.HyndmanR. J.KhandakarY. (2008). “Automatic time series forecasting: the forecast package for R.J. Stat. Softw.261–22.ICES. (2008). Report of the Workshop on Small Pelagics (Sardina pilchardus, Engraulis encrasicolus) Maturity Stages (WKSPMAT), 10-14 November 2008.Italy: Mazara del Vallo.ICES. (2017). Report of the Workshop on Age estimation of European anchovy (Engraulis encrasicolus). WKARA2 2016 Report 28 November - 2 December 2016.Pasaia: ICES.JørgensenC.EnbergK.DunlopE. S.ArlinghausR.BoukalD. S.BranderK. (2007). Ecology: managing evolving fish stocks.Science3181247–1248. 10.1126/science.114808918033868KendallM.StuartA. (1983). The Advanced Theory of Statistics, Vol.3.London: Griffin.Le CrenE. D. (1951). The Length-Weight Relationship and Seasonal Cycle in Gonad Weight and Condition in the Perch (Perca fluviatilis).J. Anim. Ecol.20201–219. 10.2307/1540LloretJ.FaliexE.ShulmanG. E.RagaJ. A.SasalP.MuñozM. (2012). Fish health and fisheries, implications for stock assessment and management: the mediterranean example. Rev. Fish. Sci.20, 165–180. 10.1080/10641262.2012.695817LloretJ.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 Ebre (Ebro) River delta (north-western Mediterranean).Fish. Oceanogr.13102–110.LloretJ.ShulmanG.LoveR. M. (2013). Condition and Health Indicators of Exploited Marine Fishes.Hoboken: John Wiley & Sons.LorenteP.PiedracobaS.Soto-NavarroJ.Alvarez-FanjulE. (2015). Evaluating the surface circulation in the Ebro delta (northeastern Spain) with quality-controlled high-frequency radar measurements.Ocean Sci.11921–935. 10.5194/os-11-921-2015LudwigW.DumontE.MeybeckM.HeussnerS. (2009). River discharges of water and nutrients to the Mediterranean and Black Sea: major drivers for ecosystem changes during past and future decades?Prog. Oceanogr.80199–217. 10.1016/j.pocean.2009.02.001MacíasD.Castilla-EspinoD.García-del-HoyoJ. J.NavarroG.CatalánI. A.RenaultL. (2014). Consequences of a future climatic scenario for the anchovy fishery in the Alboran Sea (SW Mediterranean): a modeling study.J. Mar. Syst.135150–159. 10.1016/j.jmarsys.2013.04.014ManganoM. C.MieszkowskaN.HelmuthB.DomingosT.SousaT.BaiamonteG. (2020). Moving Toward a Strategy for Addressing Climate Displacement of Marine Resources: a Proof-of-Concept.Front. Mar. Sci.7:408. 10.3389/fmars.2020.00408MartínP.BahamonN.SabatésA.MaynouF.SánchezP.DemestreM. (2008). European anchovy (Engraulis encrasicolus) landings and environmental conditions on the Catalan Coast (NW Mediterranean) during 2000-2005.Hydrobiologia612185–199. 10.1007/s10750-008-9482-1McBrideR. S.SomarakisS.FitzhughG. R.AlbertA.YaraginaN. A.WuenschelM. J. (2015). Energy acquisition and allocation to egg production in relation to fish reproductive strategies.Fish Fish.1623–57. 10.1111/faf.12043MercadoJ. M.CortésD.GarcíaA.RamírezT. (2007). Seasonal and inter-annual changes in the planktonic communities of the northwest Alboran Sea (Mediterranean Sea).Prog. Oceanogr.74273–293. 10.1016/j.pocean.2007.04.013Morales-NinB.PertierraP. J. (1990). Growth rates of the anchovy Engraulis ecransicolus and the sardine Sardina pilchardus in the Northwestern Mediterranean Sea.Mar. Biol.107349–356.MoritzS.Bartz-BeielsteinT. (2017). Imputets: time Series Missing Value Imputation in R.R J.9207–218. 10.32614/RJ-2017-009MuruaH.Saborido-ReyF. (2003). Female reproductive strategies of marine fish species of the North Atlantic.J. Northwest Atl. Fish. Sci.3323–31. 10.2960/J.v33.a212778214NashR. D.WitthamesP.PawsonM.AlesworthE. (2000). Regional variability in the dynamics of reproduction and growth of Irish Sea plaice, Pleuronectes platessa L.J. Sea Res.4455–64. 10.1016/S1385-1101(00)00046-0NeuheimerA. B.GrønkjærP. (2012). Climate effects on size-at-age: growth in warming waters compensates for earlier maturity in an exploited marine fish.Glob. Chang. Biol.181812–1822. 10.1111/j.1365-2486.2012.02673.xOspina-AlvarezA.CatalánI.BernalM.RoosD.PalomeraI. (2015). From egg production to recruits: connectivity and inter-annual variability in the recruitment patterns of European anchovy in the northwestern Mediterranean.Prog. Oceanogr.138431–447. 10.1016/j.pocean.2015.01.011PalomeraI. (1992). Spawning of anchovy Engraulis encrasicolus in the northwestern Mediterranean relative to hydrographic features in the region.Mar. Ecol. Prog. Ser.79215–223. 10.3354/meps079215PalomeraI.OlivarM. P.SalatJ.SabatésA.CollM.GarcíaA. (2007). Small pelagic fish in the NW Mediterranean Sea: an ecological review.Prog. Oceanogr.74377–396. 10.1016/j.pocean.2007.04.012PascualM.RivesB.SchunterC.MacphersonE. (2017). Impact of life history traits on gene flow: a multispecies systematic review across oceanographic barriers in the Mediterranean Sea.PLoS One12:e0176419. 10.1371/journal.pone.017641928489878PeckM. A.RegleroP.TakahashiM.CatalánI. A. (2013). Life cycle ecophysiology of small pelagic fish and climate-driven changes in populations.Prog. Oceanogr.116220–245. 10.1016/J.POCEAN.2013.05.012PenninoM. G.CollM.Albo-PuigserverM.Fernandez-CorredorE.SteenbeekJ.GiráldezA. (2020). Current and future influence of environmental factors on small pelagic fish distributions in the Northwestern Mediterranean Sea.Front. Mar. Sci.7:622. 10.3389/fmars.2020.00622PertierraJ. P.LleonartJ. (1996). NW Mediterranean anchovy fisheries.Sci. Mar.60257–267.PinotJ. M.López-JuradoJ. L.RieraM. (2002). The CANALES experiment (1996-1998). Interannual, seasonal, and mesoscale variability of the circulation in the Balearic Channels.Prog. Oceanogr.55335–370. 10.1016/S0079-6611(02)00139-8PlanqueB.FromentinJ. M.CuryP.DrinkwaterK. F.JenningsS.PerryR. I. (2010). How does fishing alter marine populations and ecosystems sensitivity to climate?J. Mar. Syst.79403–417. 10.1016/j.jmarsys.2008.12.018QuattrocchiF.MamouridisV.MaynouF. (2016). Occurrence of adult anchovy in Catalonia (NW Mediterranean) in relation to sea surface conditions.Sci. Mar.80457–66. 10.3989/scimar.04413.24a33311142QuattrocchiF.MaynouF. (2017). Environmental drivers of sardine (Sardina pilchardus) in the Catalan Sea (NW Mediterranean Sea).Mar. Biol. Res.131003–1014. 10.1080/17451000.2017.1331039QueirosQ.FromentinJ. M.GassetE.DuttoG.HuibanC.MetralL. (2019). Food in the sea: size also matters for pelagic fish.Front. Mar. Sci.6:385. 10.3389/fmars.2019.00385QuintanillaJ. M.Laiz-CarriónR.UriarteA.GarcíaA. (2015). Influence of trophic pathways on daily growth patterns of western Mediterranean anchovy Engraulis encrasicolus larvae.Mar. Ecol. Prog. Ser.531263–275. 10.3354/meps11312R Core Team. (2020). R: A Language and Environment for Statistical Computing.Vienna: R Foundation for Statistical Computing.RamírezF.AfánI.DavisL. S.ChiaradiaA. (2017). Climate impacts on global hot spots of marine biodiversity.Sci. Adv.3:e1601198.RamírezF.PenninoM. G.Albo-PuigserverM.SteenbeekJ.BellidoJ. M.CollM. (2021). SOS small pelagics: a Safe Operating Space for small pelagic fish in the Western Mediterranean Sea.Sci. Total Environ.756:144002. 10.1016/j.scitotenv.2020.14400233277012Ramírez-AmaroS.PicornellA.ArenasM.CastroJ. A.MassutíE.RamonM. M. (2018). Contrasting evolutionary patterns in populations of demersal sharks throughout the western Mediterranean.Mar. Biol.1651–16. 10.1007/s00227-017-3254-2RochetM. (1998). Short-term effects of fishing on life history traits of fishes.ICES J. Mar. Sci.55371–391. 10.1006/jmsc.1997.0324RuizJ.MacíasD.RincónM. M.PascualA.CatalánI. A.NavarroG. (2013). Recruiting at the Edge: kinetic Energy Inhibits Anchovy Populations in the Western Mediterranean.PLoS One8:e55523. 10.1371/journal.pone.005552323451027SabatésA.SalatJ.PalomeraI.EmelianovM.De PuellesM. L. F.OlivarM. P. (2007). Advection of anchovy (Engraulis encrasicolus) larvae along the Catalan continental slope (NW Mediterranean).Fish. Oceanogr.16130–141. 10.1111/j.1365-2419.2006.00416.xSalatJ. (1996). Review of hydrographic environmental factors that may influence anchovy habitats in northwestern Mediterranean.Sci. Mar.6021–32.SalatJ.GarciaM. A.CruzadoA.PalanquesA.ArínL.GomisD. (2002). Seasonal changes of water mass structure and shelf slope exchanges at the Ebro shelf (NW Mediterranean).Cont. Shelf Res.22327–348. 10.1016/S0278-4343(01)00031-0Salgado-HernanzP. M.RacaultM. F.Font-MuñozJ. S.BasterretxeaG. (2019). Trends in phytoplankton phenology in the Mediterranean Sea based on ocean-colour remote sensing.Remote Sens. Environ.22150–64. 10.1016/j.rse.2018.10.036SarauxC.Van BeverenE.BrossetP.QueirosQ.BourdeixJ.-H.DuttoG. (2019). Small pelagic fish dynamics: a review of mechanisms in the Gulf of Lions.Deep Sea Res. Part II Top. Stud. Oceanogr.15952–61. 10.1016/J.DSR2.2018.02.010SharpeD. M. T.HendryA. P. (2009). Life history change in commercially exploited fish stocks: an analysis of trends across studies.Evol. Appl.2260–275. 10.1111/j.1752-4571.2009.00080.x25567879SilvaA.SantosM. B.CanecoB.PestanaG.PorteiroC.CarreraP. (2006). Temporal and geographic variability of sardine maturity at length in the northeastern Atlantic and the western Mediterranean.ICES J. Mar. Sci.63663–676. 10.1016/j.icesjms.2006.01.005SimoncelliS.FratianniC.PinardiN.GrandiA.DrudiM.OddoP. (2019). “Mediterranean sea physical reanalysis (CMEMS MED-physics) [data set],” in Copernicus Monitoring Environment Marine Service (CMEMS).10.25423/MEDSEA_REANALYSIS_PHYS_006_004StahlJ. P.KruseG. H. (2008). Spatial and Temporal Variability in Size at Maturity of Walleye Pollock in the Eastern Bering Sea.Trans. Am. Fish. Soc.1371543–1557. 10.1577/T07-099.1StawitzC. C.EssingtonT. E. (2018). Somatic growth contributes to population variation in marine fishes.J. Anim. Ecol.88315–329. 10.1111/1365-2656.1292130381829StearnsS. C. (1989). Trade-Offs in Life-History Evolution.Funct. Ecol.3259–268. 10.2307/2389364StensethN. C.RouyerT. (2008). Destabilized fish stocks.Nature452825–826. 10.1038/452825a18421343TaboadaF. G.AnadónR. (2016). Determining the causes behind the collapse of a small pelagic fishery using Bayesian population modeling.Ecol. Appl.26886–898. 10.1890/15-0006/suppinfoTeruzziA.BolzonG.CossariniG.LazzariP.SalonS.CriseA. (2019). “Mediterranean sea biogeochemical reanalysis (CMEMS MED-biogeochemistry) [data set],” in Copernicus Monitoring Environment Marine Service (CMEMS).10.25423/MEDSEA_REANALYSIS_BIO_006_008TugoresP.GiannoulakiM.IglesiasM.BonannoA.TičinaV.LeonoriI. (2011). Habitat suitability modelling for sardine Sardina pilchardus in a highly diverse ecosystem: the Mediterranean Sea.Mar. Ecol. Prog. Ser.443181–205. 10.3354/meps09366Van BeverenE.BonhommeauS.FromentinJ.-M.BigotJ.-L.BourdeixJ.-H.BrossetP. (2014). Rapid changes in growth, condition, size and age of small pelagic fish in the Mediterranean.Mar. Biol.1611809–1822. 10.1007/s00227-014-2463-1Van BeverenE.FromentinJ. M.RouyerT.BonhommeauS.BrossetP.SarauxC. (2016). The fisheries history of small pelagics in the Northern Mediterranean.ICES J. Mar. Sci.731474–1484. 10.1093/icesjms/fsw023Vargas-YáñezM.GiráldezA.TorresP.GonzálezM.García-MartínezM.delC. (2020). Variability of oceanographic and meteorological conditions in the northern Alboran Sea at seasonal, inter-annual and long-term time scales and their influence on sardine (Sardina pilchardus Walbaum 1792) landings.Fish. Oceanogr.291–14. 10.1111/fog.12477VenteroA.IglesiasM.VillamorB. (2017). Anchovy (Engraulis encrasicolus) otoliths reveal growth differences between two areas of the Spanish Mediterranean Sea.Sci. Mar.81327–3. 10.3989/scimar.04615.21a33311142VéronM.DuhamelE.BertignacM.PawlowskiL.HuretM. (2020a). Major changes in sardine growth and body condition in the Bay of Biscay between 2003 and 2016: temporal trends and drivers.Prog. Oceanogr.182:102274. 10.1016/j.pocean.2020.102274VéronM.DuhamelE.BertignacM.PawlowskiL.HuretM.BaulierL. (2020b). Determinism of Temporal Variability in Size at Maturation of Sardine Sardina pilchardus in the Bay of Biscay.Front. Mar. Sci.7:567841. 10.3389/fmars.2020.567841ViñasJ.SanzN.PeñarrubiaL.AraguasR.-M.García-MarínJ.-L.RoldánM.-I. (2014). Genetic population structure of European anchovy in the Mediterranean Sea and the Northeast Atlantic Ocean using sequence analysis of the mitochondrial DNA control region.ICES J. Mar. Sci.71391–397. 10.1093/icesjms/fst132WalshM. R.MunchS. B.ChibaS.ConoverD. O. (2006). Maladaptive changes in multiple traits caused by fishing: impediments to population recovery.Ecol. Lett.9142–148. 10.1111/j.1461-0248.2005.00858.x16958879WeiT.SimkoV. (2017). R Package “corrplot”: Visualization of a Correlation Matrix (Version 0.84). Available Online at: https://github.com/taiyun/corrplot.WickhamH. (2016). ggplot2: Elegant Graphics for Data Analysis.New York: Springer-Verlag:.WoodS. N. (2006). Generalized Additive Models: An Introduction With R. Texts in Statistical Science. Boca Raton, FL: Chapman & Hall/CRC, 410.WoodS. N. (2011). Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models.J. R. Stat. Soc. Ser. B Stat. Methodol.733–36. 10.1111/j.1467-9868.2010.00749.xYangL.QinG.ZhaoN.WangC.SongG. (2012). Using a generalized additive model with autoregressive terms to study the effects of daily temperature on mortality.BMC Med. Res. Methodol.12:165. 10.1186/1471-2288-12-16523110601ZuurA. F.IenoE. N.ElphickC. S. (2010). A protocol for data exploration to avoid common statistical problems.Methods Ecol. Evol.13–14.