European hake (
La merluza europea (
Life history parameters, and in particular reproductive traits such as spawning cycle, size and age at maturity, fecundity and spawning potential are the basis for assessing the productivity and resilience of fish stocks, making them essential information for stock assessment and management (
European hake (
Several studies have analysed the size at maturity of this stock (
Length and weight measurements were taken from historical records collected by the commercial sampling programme of the Spanish Institute of Oceanography (IEO) on a monthly basis between 1982 and 2019 from bottom trawlers that operate in the Atlantic Iberian waters, ICES divisions 27.8.c and 27.9.a. A total of 27692 hake specimens were sampled, from which individuals’ length and weight (total and gutted), maturity stage and sex were recorded. Both maturity and sex, were assessed based on macroscopic observations. Of the 27692 hake specimens, the majority (56%) were females. In addition, the majority of the specimens sampled (54%) were mature. Crossreferencing these variables (sex and maturity stage), the data indicate that 63% of the males and 46% of the females were mature.
In addition, 90% of the hake specimens were caught in division 27.9.a versus 10% in 27.8.c, which indicates that the data were not proportional by area. Also, data are not available in certain months for each year throughout the study period (For more details, see
Owing to its complexity, dealing with a large amount of data coming from different sources of information is not exempt from human error. To avoid this, a thorough review of the matrix was carried out to ensure data consistency. Finally, to test for densitydependent effects, we included data on the total biomass (in tonnes) and spawning stock biomass (SSB, in tonnes) at length of the European hake for the study period, obtained from the data compiled by the ICES Working Group for the Bay of Biscay and the Iberian Waters Ecoregion (WGBIE) (
The percentage of mature females for every length class was fitted to a logistic equation as described by
where PL is the probability of an individual being mature at a given length, β_{0} the intercept and β_{1} the slope of the curve. These parameters are estimated by iteratively reweighted least squares, assuming a binomial distribution (matureimmature).
The length at which 50% of specimens were mature (L_{50,} size at maturity) was estimated as
where
In addition, the classical logistic models including and not including sex as a covariable and were compared through a Chisq likelihoodratio test (Pvalue<0.05) to check whether a partial effect of sex existed.
It should be noted that the lack of observations of immature male individuals for the years 2014 and 2015 owing to a reduction in the sampling prevented the size at maturity of males from being estimated in these years. To address this handicap, missing values were calculated using moving averages with a window of four (eight observations were taken into account; four on the left and four on the right) using the
Body condition is a key indicator of the health status of fish since it is closely related to important fitness variables such as growth, reproduction, behaviour and survival. Consequently, the body condition can partially explain the L_{50} variability. Among the several indices proposed to measure body condition, Le Cren’s relative condition factor (
where W is the total weight (g), L the total length (cm), the constant estimate and the exponent estimate usually between 2.5 and 4. The â and ̂b lengthweight relationship parameters were estimated using the
In order to analyse the effect of the environmental conditions on the L_{50}, we used the sea surface temperature (SST in °C) as a local environmental index and the North Atlantic Oscillation (NAO) and the Atlantic Multidecadal Oscillation (AMO) as largescale climate indices representative of sea level pressure patterns. The AMO and NAO indices were included because largescale climatic events affect not only temperature but also ocean masses and circulation, which affect plankton and zooplankton (
Generalized additive models (GAMs,
GAMs are semiparametric extensions of generalized linear models (GLMs) for which the strictly linear predictor
where
where f_{j} (X_{j}) is the unknown smooth partial effect of X_{j} on the predictor. Hence, GAMs avoid the assumption of a linear relation between the response variable and the covariables, thus providing a more flexible model. Note that GLMs are an extension of linear models for which the distribution of the response variable can be other than Gaussian. For this reason, in the previous models a link function
Standardized data exploration techniques were used to identify any outliers and any correlation and collinearity between the explicative variables (
Covariables were selected with a backward stepwise procedure based on the Akaike information criterion (AIC) and an adjusted Rsquare. The best (and most parsimonious) model was finally chosen based on the compromise between low AIC values, a high adjusted Rsquare and significant predictors.
Cubic regression splines basis were used for explanatory variables, restricting the dimension of the basis (k) to four to allow a high degree of flexibility without overfitting problems. When a spline was not required, the variable was considered a linear effect.
For each final GAM, a diagnosis of the model assumptions was carried out, among which the independence of the residuals was verified using the autocorrelation (ACF) and partial autocorrelation (PACF) functions and the LjungBox test (
The GAMs were performed using the
Firstly, a single value of the size at maturity of European hake was estimated for the whole period 19822019 and for both sexes combined, resulting in L_{50}=37.5 cm. European hake shows sexual dimorphism, hence the was also estimated by sex in this period. The data exploration analysis showed that the estimate was 44.4 cm and 31.8 cm for females and males, respectively (
The grey dots represent the observations, the blue line is the fitted model and the red line represents the size at maturity.
Secondly, sexspecific size at maturity was also estimated for each year of the time series (
For the whole period 19822019 the relative condition factor of European hake was
The final GAM for European hake males derived from the backward stepwise procedure included only the year as a significant covariable (
Parametric coefficients  

Estimate  se  t value  Pr (>t)  

31.0263  0.4684  66.24  <2e16 *** 
Significance of smooth functions  
Edf  Ref. df  F  Pvalue  

1.616  1.969  22.77  1.1e06 *** 
The model results show that 55.3% of the L_{50} variability of European hake males was explained by the temporal covariable year.
The blue circles represent the estimated annual values and the blue shading represents the 95% confidence intervals of the year’s smooth component.
For European hake females, in addition to the temporal year covariable, the final GAM included the biological covariables, biomass and spawning biomass at length, and the NAO environmental covariable (
Parametric coefficients  

Estimate  Std error  t value  Pr (>t)  

44.4441  0.3035  146.4  <2e16 *** 
Significance of smooth functions  
Edf  Ref. df  F  Pvalue  

2.595  3.245  2.541  0.04783 * 

3.198  3.634  4.992  0.00474 ** 

2.249  2.788  7.734  0.00144 ** 

1.000  1.000  8.094  0.00821 ** 
In particular, year showed a domeshaped relationship with L_{50} estimates (
The pink circles represent the estimated annual values and the pink shading represents the 95% confidence intervals of the smooth components.
In both GAMs all the theoretical assumptions were respected, and the results of the autocorrelation (ACF) and partial autocorrelation (PACF) functions and the LjungBox test of the residuals are presented in the supplementary material.
Understanding biological processes underlying southern hake stock dynamics in Atlantic Iberian waters can help to improve their assessment and management. The overall aim of this study was to assess changes in the size at maturity (L_{50}) of both sexes in the period 19822019 and determine the potential drivers of its variation.
It should be mentioned that the data used in this study were limited to the northern part of the southern stock distribution and that, owing to the missing data and the way the data were compiled, these results should be interpreted with caution. However, our results highlight the importance of understanding the factors affecting the balance between the energy invested in reproduction, maintenance and growth and the interannual variability of this tradeoff. This is the only way to understand how environmental pressures influence the population dynamics of this species at a regional scale.
The results of this study showed that differentiated longterm patterns were observed for the two sexes, with major longterm declines for male hakes. For males, the size at maturity decreased by 12.9 cm over the study period, from 36.1 cm in 1982 to 23.2 cm, the lowest observed value, in 2019. However, for females two trends in L_{50} can be distinguished: an increasing trend starting in 1990 and reaching its highest value of 56.2 cm in 1996, followed by a decreasing trend for the rest of the time series, reaching a value of 38.8 cm in 2019.
In relation to
From the above, it can be deduced that though both sexes follow the same pattern for
The reproductive system of fish reacts to any changes in life conditions, so variations in growth and reproduction dynamics of fish populations substantially affects fish production not only quantitatively but also qualitatively (
In Iberian Atlantic waters, European hake has been exploited beyond safe biological limits since the 1990s, with the largest landings being recorded in this period, but in 2004 it was declared in a critical state because of overfishing (
In line with the above, fishing pressure has been widely reported as a selection factor for larger individuals with an effect that is especially problematic when the oldest and largest individuals of the population are removed (
Several studies have addressed the changes in the maturity ogives in different species
For male hake, both environmental and biological factors failed to explain the trend towards decreasing L_{50}for our entire study period. This may well be explained by the fact that fishing pressure has eliminated larger individuals (mostly females), which could lead to a clear alteration of growth rates, causing profound changes in L_{50}, either because fishing intensity has caused a regime shift, which has caused the main drivers of maturation/growth to change, or even because there has been a regime shift in the fishery itself, so fishing pressure has lost influence relative to other drivers. Considering the above, fishing pressure and behaviour of the trawl fleet could explain the decline, but further work is required because males and females do not show the same response.
In the case of female hake, the reduction in stock biomass, usually as a consequence of high fishing exploitation, leads to a reduction in size at maturity. The results showed a negative relationship between biomass and size at maturity for female hake, which is to be expected in a densitydependent relationship. The higher the biomass, the greater the intraspecific competition, the lower the growth and the lower the L_{50}. However, the results showed a positive relationship between biomass at spawning and L_{50}, which is contrary to the compensatory theory. The lack of compensatory response is known as depensation. According to
Major observed declines in L_{50} are usually explained by unfavourable environmental conditions that reduce growth or by high fishing pressure (or a combination of both) as a mechanism to maximize fitness (
The only significant environmental variable included in our model for European hake females was the NAO index, which showed a positive relationship with the L_{50}. The NAO is the most important source of variability in the North Atlantic region. The NAO index acts as an integrator of several environmental factors that can have a synergistic effect on fisheries dynamics. NAO impacts on air temperature, wind and precipitation regime and on large water mass distribution and flow, with the subsequent influence that all these factors have on marine ecosystems (
The size at maturity of the southern European hake stock has declined in both sexes in the last few decades. The present results show that the maturation of southern hake is influenced by a combination of demographic, anthropogenic and environmental factors, leading to a deeper understanding of the factors explaining the annual decrease in L_{50} for both sexes, which was the objective of this study.
However, determining which variables explain the decreasing trends in L_{50}for both sexes requires applying more complex models than an overall linear trend to all the data, because there are much more complex relationships between L_{50}and the factors that may explain the variability in L_{50}.
Total biomass, spawning biomass at length and the NAO partially explain the decline in the L_{50} of females, but for males no relevant factor was found among the analysed ones that could explain this drastic decline. Overall, life history traits that rapidly respond to external changes, such as L_{50}, might be good indicators to anticipate further declines in populations and require close monitoring and evaluation. That is why this study reiterates the need for a posteriori studies to confirm the hypothesis of fishing pressure as a maturity driver in order to have a better understanding of the underlying dynamics of this stock. Additionally, the impact of the reduction in the L_{50} in males on the reproductive potential of the stock should be also investigated as paternal effects have proved to be more important than previously thought in stock dynamics (
This study is a contribution to the IMPRESS (RTI2018099868BI00) project, ERDF, Ministry of Science, Innovation and Universities  State Research Agency. Data were provided by samplings of commercial fleets that were cofunded by the European Union through the European Maritime and Fisheries Fund (EMFF) within the national programme of collection, management and use of data in the fisheries sector and support for scientific advice regarding the Common Fisheries Policy.
1  2  3  4  5  6  7  8  9  10  11  12  

1980  0  0  0  0  177  115  333  0  97  0  0  99 
1981  50  0  0  0  0  0  0  0  0  0  0  0 
1982  81  43  152  61  31  30  33  0  0  0  0  0 
1983  0  93  49  151  0  24  0  0  0  0  0  0 
1984  22  142  302  232  0  0  48  22  0  119  122  0 
1985  167  168  21  50  252  245  328  294  0  277  155  251 
1986  63  181  45  0  0  0  0  103  0  93  0  0 
1987  0  0  0  0  0  5  165  0  0  0  3  0 
1988  0  76  0  0  0  462  240  231  91  0  0  0 
1989  0  29  290  57  345  230  38  198  0  0  0  0 
1990  0  98  130  43  14  66  0  142  0  34  58  0 
1991  58  0  35  77  332  25  41  74  0  0  0  0 
1992  0  26  134  0  0  0  0  0  0  0  0  0 
1993  0  0  0  38  32  135  0  0  0  0  0  33 
1994  0  83  60  0  143  36  0  0  12  20  117  35 
1995  0  81  154  175  0  94  67  32  10  0  80  0 
1996  0  5  227  8  83  152  71  199  201  16  41  18 
1997  12  73  101  8  47  119  0  143  0  110  0  0 
1998  0  0  101  0  162  0  23  0  5  0  0  0 
1999  0  76  83  51  34  0  2  0  0  2  0  3 
2000  31  43  7  50  0  40  22  0  76  0  55  0 
2001  10  14  42  0  22  0  0  0  0  0  0  0 
2002  0  0  0  0  74  177  133  0  35  40  27  12 
2003  95  280  164  159  323  298  358  146  48  103  195  140 
2004  272  162  284  164  200  281  80  165  140  63  318  64 
2005  138  198  260  75  272  3  215  419  0  0  286  0 
2006  248  116  80  348  0  33  224  0  0  0  0  52 
2007  162  39  258  161  116  143  118  78  148  52  0  109 
2008  152  225  25  246  65  75  84  157  0  0  76  153 
2009  140  21  29  100  27  0  0  169  0  0  0  0 
2010  147  106  66  112  48  0  151  0  41  0  11  25 
2011  0  0  286  0  20  0  62  112  0  0  106  0 
2012  44  71  121  0  0  0  0  0  0  0  0  0 
2013  0  30  0  98  47  97  0  0  56  0  0  36 
2014  0  51  0  23  110  0  0  0  0  0  0  18 
2015  0  73  0  0  18  0  0  0  0  0  0  0 
2016  0  76  95  37  79  49  0  0  0  0  0  0 
2017  0  0  0  20  270  0  50  0  0  0  42  38 
2018  0  65  0  0  72  15  0  0  0  0  43  0 
2019  46  46  0  56  119  0  0  0  0  0  0  0 
In this section, a graphical analysis is presented to illustrate the data for size at maturity (for males and females) and the relative condition factor (for males and females) during the study period before imputing missing values.
In this section, data exploration techniques were presented to identify possible correlation and collinearity between the explanatory variables (
Colour intensity (light to dark) and size of circles are proportional to Pearson correlation coefficients; positive correlations are represented in blue and negative correlations in red. Significant values are shown (p<0.05). NAO, North Atlantic Oscillation; AMO, Atlantic Multidecadal Oscillation; SST, sea surface temperature; SSB, spawning stock biomass at length; BIO, biomass; Kn, relative condition factor.
Predictor variables  GVIF 

Year  2.944590 
AMO  3.234100 
NAO  1.135156 
SST  2.432250 
Kn  1.161385 
BIO  2.366113 
SSB_length  1.304750 
Colour intensity (light to dark) and size of circles are proportional to Pearson correlation coefficients; positive correlations positive correlations are represented in blue and negative correlations in red. Significant values are shown (p<0.05). NAO, North Atlantic Oscillation; AMO, Atlantic Multidecadal Oscillation; SST, sea surface temperature; SSB, spawning stock biomass at length; BIO, biomass; Kn, relative condition factor.
Predictors variables  GVIF 

Year  3.278021 
AMO  3.436279 
NAO  1.130783 
SST  2.499691 
Kn  1.494685 
BIO  2.304125 
SSB_length  1.459627 
Covariables for each final GAM were selected with a backward stepwise procedure based on the AIC and adjusted Rsquare. The best (and most parsimonious) model was finally chosen based on the compromise between low AIC values, high adjusted Rsquare and significant predictors.
Models  D2 (%)  AIC 

m<gam(L_{50} ~s(Year,bs=”cr”,k=6)+s(NAO,bs=”cr”,k=6)+s(SST,bs=”cr”,k=6)+s(SSB_length,bs=”cr”,k=6)+s(BIO,bs=”cr”,k=6)+s(Kn,bs=”cr”,k=6),method=”REML”,data=L_{50}mb)  63.8%  195.07 
m<gam(L_{50} ~s(Year,bs=”cr”,k=6)+s(NAO,bs=”cr”,k=6)+s(SSB_length,bs=”cr”,k=6)+s(K_{n},bs=”cr”,k=6),method=”REML”,data=L_{50}mb)  61.6%  191.94 
m<gam(L_{50}m~s(Year,bs=”cr”,k=6)+s(NAO,bs=”cr”,k=6)+s(SSB_length,bs=”cr”,k=6),method=”REML”,data=L_{50}mb)  59.8%  191.71 
m<gam(L_{50}m~s(Year,bs=”cr”,k=7)+s(NAO,bs=”cr”,k=7),method=”REML”,data=L_{50}mb)  57.5%  191.75 
m<gam(L_{50}m~s(Year,bs=”cr”,k=5),method=”REML”,data=L_{50}mb)  55.3%  193.65 
Models  D2 (%)  AIC 

m<gam(L_{50}~s(Year,bs=”cr”,k=6)+s(SSB_length,bs=”cr”,k=6)+s(NAO,bs=”cr”,k=6)+s(SST,bs=”cr”,k=6)+s(BIO,k=6,bs=”cr”)+s(Kn, k=6, bs= “cr”), method = “REML”,data=L_{50})  76.8%  172.31 
m<gam(L_{50}~s(Year, bs = “cr”, k= 6)+s(SSB_length,bs=”cr”,k=6)+s(NAO,bs=”cr”,k=6)+s(BIO,k=6,bs=”cr”)+s(Kn,k=6,bs=”cr”),method=”REML”,data=L_{50})  76.8%  170.49 
m<gam(L_{50}~s(Year,bs=”cr”,k=7)+s(SSB_length, bs=”cr”, k=7)+s(NAO,bs=”cr”,k=7)+s(BIO,k=7,bs=”cr”),method=”REML”,data=L_{50})  76.8%  169.21 
This section presents the contrasts of the basic hypotheses on the residuals of the GAMs addressed in this study to explain the annual variability of the size at maturity for male and female European hakes.
Test  Hypotheses  Pvalue 

JarqueBera  Normality  0.587 
ShapiroWilk  Normality  0.5972 
DickyeFuller  Stationarity  0.01 
Student’s t  Zero mean  1 
The upper left corresponds to the Pvalues of the LjungBox. The upper lefthand side is the pvalue of the LjungBox independence contrast and the upper righthand side is the sign of the homoscedasticity of the model residuals. The lower graphs represent the simple and partial correlations of the residuals.
Test  Hypotheses  Pvalue 

JarqueBera  Normality  0.9266 
ShapiroWilk  Normality  0.9836 
DickyeFuller  Stationarity  0.01 
Student’s  Zero mean  1 
The upper lefthand side is the pvalue of the LjungBox independence contrast and the upper righthand side is the sign of the homoscedasticity of the model residuals. The lower graphs represent the simple and partial correlations of the residuals.