Changes in Life History Traits of Small Pelagic Fish in the Western Mediterranean Sea

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.

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 INTRODUCTION 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(Hidalgo et al., , 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 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., 2017Quattrocchi 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;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 , 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.

Study 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.
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 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 (L 50 ), 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 L 50 (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 (L 50 ) 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 L 50 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; W gutted ) and the predicted weight (Wp) for a fish of the same length (Le Cren, 1951) The Wp was estimated applying the linear length-weight relationship log(W gutted ) = 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 * TL b ). The selection of W gutted 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/m 3 /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)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)

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 longterm 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 deseasonalized 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 varianceinflation 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 (R 2 -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 R 2 -adj., and significant predictors. The code used for the analysis carried out can be found in the GitHub repository 1 . The parameters that were statistical significant (p-value < 0.05) are reported. Statistics acronyms are: R 2 -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).

Changes in Length at Age
Sardine 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.
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 lengthat 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 (R 2 -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 R 2 -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 R 2 -adj. of 0.803 (Figure 3 and Supplementary Tables 2, 3).

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 R 2 -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 R 2 -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 R 2 -adj. values compared to sardine (Supplementary Tables 2, 4).

Changes in Length at First Maturity (L 50 )
Sardine Different temporal patterns were observed in sardine L 50 between areas (Figure 2). Although in the three areas there were yearly fluctuations in L 50 , a clear temporal decline was only observed in the north area. Consequently, in the period 2012-2017 the north area presented lower L 50 than the other two areas. The

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 L 50 in the north area was clearly higher than in the south. In the last period of this study (2012)(2013)(2014)(2015)(2016)(2017), the average L 50 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 L 50 values were observed in 2011 (11.56 ± 0.08 cm), followed by three consecutive years of L 50 decrease. The minimum L 50 was observed in 2014 (8.66 ± 0.44 cm). In the south area, the maximum L 50 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 L 50 value was observed (8.22 ± 1.06 cm) (Figure 2).

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) Figure 2C).
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 R 2 -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 R 2 -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 R 2 -adj. and the AIC increased (Table 1).
In the north area, the selected model had a R 2 -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 R 2 -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 R 2 -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).

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 R 2 -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 R 2 -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 R 2 -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).

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. 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 R 2 -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 R 2 -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 R 2 -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 R 2 -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 R 2 -adj. of 0.255. The MC had a positive relationship with Kn, and the ST150 had a negative effect on Kn (Figure 7C).

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)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(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 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 R 2 -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 R 2 -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 R 2 -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 R 2 -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).

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 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. 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.
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 L 50 as well. Spatial differences in L 50 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 . 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 , 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 L 50 . 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 L 50 of sardine was 12.58 cm, individuals of the period 2012-2017 presented a L 50 more than 2 cm smaller than the L 50 in the central and southernmost areas. The decline in L 50 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 L 50 . Therefore, the Gulf of Lions (GSA07) and the northern area of our study (north-GSA06) presented more similar changes in L 50 and length at age than the central area for sardine. Important observed declines in L 50 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 . These mechanisms have been previously reported in short-lived species with high plasticity, like small pelagic fish (Silva et al., 2006;Neuheimer and Grønkjaer, 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 L 50 , 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 tradeoffs 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 . 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 . 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 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  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., 2015bBrosset et al., , 2016bSaraux 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  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.

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).