Assessing changes in size at maturity for the European hake ( Merluccius merluccius ) in Atlantic Iberian waters

Summary: European hake ( Merluccius merluccius ) is a commercially important resource in Iberian Atlantic waters. Despite the recovery plan implemented in 2006 and the multiannual management plan for western waters, fishing mortality is still higher than that corresponding to the maximum sustainable yield for the southern European hake stock. The biological processes underlying the dynamics of this stock and its life history traits are essential for assessing population productivity and resilience, making them basic information for management . We analysed the temporal variability of size at maturity (L 50 ) of this species and the main factors influencing it in Atlantic Iberian waters from 1982 to 2019. The annual variability of L 50 for each sex was modelled with generalized additive models, considering explanatory environmental variables (Atlantic Multidecadal Oscillation, North Atlantic Oscillation and sea surface temperature) and biological variables (biomass, spawning biomass at length and relative condition factor). The results showed that the L 50 of males decreased by a total of 12.9 cm and L 50 of females decreased by a total of 10.9 cm from 1982 to 2019. For females the significant explanatory variables were year, spawning biomass at length, biomass and the North Atlantic Oscillation, while for males only year was an explanatory variable. These results are important for understanding the status of the European hake population, signalling that L 50 is a good indicator for predicting future population dynamics.


INTRODUCTION
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 (Dominguez-Petit et al. 2008). The estimation of the size at which 50% of individuals are mature (L 50 ) can be used to define the minimum legal landing size for exploited stocks. However, in most stock assessments the causes of L 50 variability are not analysed, though there are intrinsic and extrinsic factors affecting maturity that may act at a temporal scale (Lorenzen and Camp 2019). Size and age at maturity are especially sensitive to population abundance (Dominguez-Petit et al. 2008, Lorenzen andCamp 2019) owing to density-dependent effects. When stock size decreases, individuals tend to mature at large sizes and early ages because of compensatory growth (Ali et al. 2003). However, some environmental and anthropogenic drivers, including high rates of fishing mortality, may also impact on size at maturity (Pepin 2015). Moreover, when the effects of these drivers are intense and/or sustained over a long period of time, they could lead to evolutionary changes because of the selection of certain genotypes (Hollins et al. 2018). Subsequently, understanding mechanisms driving maturity variability within species could provide preliminary evidence of the species's adaptability to climate change (Dalgleish et al. 2010), habitat degradation (Öckinger et al. 2010) and overfishing (Hobday et al. 2011). Barot et al. (2004) identified a downward shift in size at maturity in the Georges Bank Atlantic cod (Gadus morhua) that supports the hypothesis that an evolutionary change, likely caused by high fishing mortality, is partially responsible for the observed decline in age and size at maturity in these cod stocks (Olsen et al. 2005). Regardless of the causes of variation in maturation, size and age at maturity can have significant effects on the reproductive performance of the stock by regulating its reproductive potential. For example, fecundity and the quality and viability of eggs and larvae are directly related to the size of spawning females (Trippel et al. 1997, Marteinsdottir andBegg 2002). The extent to which these effects may extend over several generations and affect the adaptive capacity of populations is not yet well understood, although understanding these processes is critical to comprehending stock dynamics.
European hake (Merluccius merluccius) is a resource of great commercial importance in Atlantic Iberian waters. This species is assessed by the International Council for the Exploration of the Sea (ICES) in two units: the northern and the southern stocks (ICES 2021). The southern stock, which is the object of this study, is distributed in the Atlantic Iberian waters that correspond to ICES divisions 27.8.c and 27.9.a. Though this stock has been subjected to a recovery plan since 2006 and the multiannual management plan for Western Waters, fishing mortality is still above that corresponding to the maximum sustainable yield (F msy ) (ICES 2021). For this stock, the main indicator of stock status is the spawning stock biomass, which is calculated as numbers, mean weight and proportion of mature fish. The proportion of mature fish was established for the entire study period  by estimating annual maturity ogives for both sexes combined based on macroscopic observations of hake specimens, although a time-constant female-only maturity ogive was used from 2020 onwards owing to the change of the assessment model. The impact of using different reproductive potential metrics (combined vs only female) was evaluated by Cerviño et al. (2013) for this stock, showing that, although the absolute values change, their impact on the management recommendations is slight because the reference points also change in the same direction.
Several studies have analysed the size at maturity of this stock (Dominguez-Petit et al. 2008, Murua 2010, but generally only for females and never for the two sexes separated. The objective of this study was to provide information on size at maturity for the period 1982-2019 in order to corroborate the decreasing trend of the L 50 in this stock for each sex and to determine which factors explain this trend. Understanding which factors influence the size at maturity of the southern European hake stock can contribute to better assessment and management. To this end, we first estimated the size at maturity of the European hake for each sex in Atlantic Iberian waters from 1982 to 2019. Second, the variability of the L 50 for each sex was modelled with generalized additive models, considering as explanatory variables environmental factors (Atlantic Multidecadal Oscillation, North Atlantic Oscillation and sea surface temperature), and biological variables (biomass, spawning biomass at length and relative factor condition) to test density-dependent effects and the year to assess the interannual variations.

European hake data
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. Cross-referencing 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 Table S1), which makes monthly analysis difficult. To avoid influence of these biases on analyses results, both areas (27.8.c and 27.9.a) were considered together and year was used as a temporal variable.
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 density-dependent 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) (ICES 2019) for assessment purposes.

Size at maturity
The percentage of mature females for every length class was fitted to a logistic equation as described by Ashton (1972): where PL is the probability of an individual being mature at a given length, the intercept and the slope of the curve. These parameters are estimated by iteratively reweighted least squares, assuming a binomial distribution (mature-immature).
The length at which 50% of specimens were mature (L 50, size at maturity) was estimated as where and are the intercept and slope estimates, respectively, estimated using the SizeMat (Torrejón-Magallanes 2020) package of the R software (R Core Team 2021).
In addition, the classical logistic models including and not including sex as a covariable and were compared through a Chisq likelihood-ratio test (P-val-ue<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 imputeTS package (Moritz and Bartz-Beielstein 2017, see details in the supplementary materials).

Le Cren's condition factor
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 (Kn) (Le Cren 1951) was selected. Le Cren's factor is defined as the observed weight of an individual divided by its predicted weight, which is obtained from the linear regression of the length-weight relationship (after linearizing using logarithms).
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 lm function of the R software through the linear regression model log(W)=log(a)+blog(L). Gutted weight was missing for a large percentage of the individuals, so total weight was used instead. However, the presence of missing values of total weight prevented us from estimating the relative condition factor for males between 1998 and 2001 and for females between 1999 and 2001. These missing values were calculated by moving averages with a window of four using the imputeTS package (Moritz and Bartz-Beielstein 2017).

Environmental data
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 large-scale 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 (Nye et al. 2009(Nye et al. , 2014. These indices have disproportionate effects on the early life stages of fishes because of alterations in food sources and temperature (Nye et al. 2009(Nye et al. , 2014. SST was included because it affects the fish community through thermal tolerance and species migration (Marshall and Elliott 1998). The AMO and NAO annual-times series were downloaded from the US National Oceanic and Atmospheric Administration (NOAA, http://www. noaa.gov/). For SST, annual averages covering the entire study period  and sampling area were extracted from the NOAA NCDC platform (https:// www.ncdc.noaa.gov/oisst), which provides a global 1/4° gridded dataset of the advanced very high-resolution radiometer.

Generalized additive models
Generalized additive models (GAMs, Wood 2017) were implemented to test the influences of temporal (year), environmental (SST, AMO and NAO) and biological (relative condition factor, biomass, and spawning biomass at length) variables on the L 50 .
GAMs are semi-parametric extensions of generalized linear models (GLMs) for which the strictly linear predictor g(μ(X)) =β 0 +β 1 X 1 +···+β p X p , where X=(X 1 ,...,X p ) are covariables, μ(X)=E(Y|X) is the conditional exception of the response variable Y, g is the link function (explained below) and β 0 , β 1 ,...,β p are the unknown parameters, is replaced by 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 g is applied to μ(X). In this study, it was assumed that the Gaussian distribution correctly fits the response variable size at maturity (L 50 ). Standardized data exploration techniques were used to identify any outliers and any correlation and collinearity between the explicative variables (Zuur et al. 2010). In particular, correlation among variables was checked using a Pearson's correlation test with the corrplot package (Wei et al. 2017) of the R software. In particular, correlation among variables was checked by performing a Pearson's correlation test with the corrplot package (Wei et al. 2017) of the R software. By calculating the generalized variance inflation factors (GVIF), which are the values of VIF corrected for the number of degrees of freedom of a predictor variable (Fox and Weisberg 2019), collinearity was tested. GVIF was evaluated using the corvif function in R software. In particular, as both the Pearson's correlation (r>0.70) and the GVIF (>3) were high between the year and AMO, the latter was not further considered in models for avoiding collinearity (see more details in supplementary materials).
Covariables were selected with a backward stepwise procedure based on the Akaike information criterion (AIC) and an adjusted R-square. The best (and most parsimonious) model was finally chosen based on the compromise between low AIC values, a high adjusted R-square 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 Ljung-Box test (Wood 2006). Additionally, results of Jarque-Bera and Shapiro-Wilk tests (normality assumption), a Dickye-Fuller test (stationarity assumption) and a Student t-test (zero mean assumption) were also explored.
The GAMs were performed using the mgcv package (Wood 2011) of the R software. The R code for the applied models can be found on the GitHub repository (https://github.com/Lojamodav/PAPER.git).

Size at maturity estimates
Firstly, a single value of the size at maturity of European hake was estimated for the whole period 1982-2019 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 ( Fig. 1). In conclusion, the partial effect of the sex covariable was significant. Secondly, sex-specific size at maturity was also estimated for each year of the time series (Fig. 2). The male L 50 time series showed a decreasing trend over the whole period ( Fig. 2A), while the female L 50 time series showed two trends (Fig. 2B): an increasing trend starting in 1990 and reaching its highest value of 56.2 cm in 1996, followed by a decreasing trend in the L 50 for the remaining years. The male L 50 in 1982 was 36.1 cm and reached the lowest value of 23.2 cm in 2019, which is a decrease of 19.2 cm (40.6% in relative terms) with respect to the mean value for male L 50 . However, the female L 50 started at 49.7 cm in 1982 and declined to 38.8 cm at the end of the series, which represented a decrease of 24.5% in relative terms with respect to the mean value for female L 50 .

Relative condition factor estimates
For the whole period 1982-2019 the relative condition factor of European hake was Kn=0.88 for both sexes combined, 0.61 for males and 1.10 for females. When Kn was estimated by year, it was observed that males showed a decreasing trend from the beginning of the time series until 1993, from which an increase in Kn was observed (Fig. 3B). Likewise, female Kn followed the same pattern as male Kn. In particular, the lowest Kn value for males and females was 0.37 and 0.35, respectively, recorded in 1993. This translates into a decrease in Kn of 39% and 68.2% (with respect to the mean) for males and females, respectively. In contrast, the highest Kn value was 0.87 (2014) and 2.48 (1998) for males and females, respectively, an increase of 42.6% and 125.45% for males and females in 1998, respectively. It should be noted that the high Kn value for females in 1998 was due to a bias in the sample because only 15 mature hakes showed very high weights.

Analysis of temporal variability of L 50
The final GAM for European hake males derived from the backward stepwise procedure included only the year as a significant covariable (Tables 1 and S4).
The model results show that 55.3% of the L 50 variability of European hake males was explained by the temporal covariable year. Figure 4 shows estimated smoother effect of year of the fitted GAM, which showed a negative relationship with the L 50 estimates.
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 (Tables  2 and S5). This model explained 76.8% of the variability of the size at maturity of female hake.
In particular, year showed a dome-shaped relationship with L 50 estimates (Fig. 5). Density-dependent ef-  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 1982-2019 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 trade-off. This is the only way to understand how environmental pressures influence the population dynamics of this species at a regional scale. fects were detected. The results showed a positive relationship of biomass with L 50 estimates and a negative relationship of spawning biomass at length with L 50 estimates (Fig. 5). Finally, the NAO index also showed a slight positive relationship with L 50 (Fig. 5).
In both GAMs all the theoretical assumptions were respected, and the results of the autocorrelation (ACF) and partial autocorrelation (PACF) functions and the Ljung-Box test of the residuals are presented in the supplementary material.

Understanding biological processes underlying southern hake stock dynamics in Atlantic Iberian waters
The results of this study showed that differentiated long-term patterns were observed for the two sexes, with major long-term 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 Kn, both sexes showed the same structure, that is, a tendency to increase from 1993 onwards, in which year both males and females reached the minimum values of 0.37 and 0.35, respectively. This translates into a decrease in Kn of 39% for males and 68.2% for females with respect to the mean. In contrast, the highest Kn value was 0.87 in 2014 for males and 2.48 in 1998 for females, which is an increase of 42.6% and 125.45% percentage points, respectively. It should be noted that the high Kn value for females in 1998 is due to a bias in the sample, as only 15 mature hakes showed very high weights. This can be explained by the fact that when total weight (heavily influenced by gonad size/maturity stage and stomach size/repletion stage) is used to calculate Kn, these weights can fluctuate markedly because of factors unrelated to the actual physiological condition of the hake specimen.
From the above, it can be deduced that though both sexes follow the same pattern for Kn, male and female maturity responded differently to over time, likely linked to external drivers, including fishing pressure.
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 (Godø andHaug 1999, Oven 2004). When fishing pressure removes old and large individuals from the spawning stock, it may favour density-dependent growth, leading to an increase in size at maturity (Olsen et al. 2004); however, if fishing pressure is too intense and maintained for a long time, it can remove the fast-growth or large-size genotypes from the stock, resulting in a decrease of size at maturity and even causing evolutionary changes in the stock (Trippel 1995, Jørgensen et at. 2007, Hidalgo et al. 2012. 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 (Dominguez-Petit 2007, Korta et al. 2010. Currently, the stock is considered to be in a relatively healthy status but still shows reduced reproductive capacity, according to the annual observed recruitment. It is still intensely exploited (in 2018, F=0.60), exceeding the rate of fishing mortality for the maximum sustainable yield (F msy = 0.25) but not above the precautionary limits (F pa =0.75) (ICES 2019). In addition, European hake is the target of mixed fisheries delivered by the industrial and artisanal fleets with different gears: gillnet, longline and different trawl gears (Korta et al. 2015, ICES 2019 targeting medium-large specimens (<40cm). Gillnets and longlines, have accounted for 16% and 13% of catches during the last decade, respectively, targeting hake adults (especially longliners that target large hake). However, Dominguez-Petit (2007) and Korta et. al (2015), reported a decrease in the mean hake size in the landings in the last few decades independently of the gear type used, supporting the perception of overexploitation.
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 (Law 2000), because the size/age truncation of the stock leads to a decrease in reproductive potential (Hixon et al. 2014). However, no proxy for fishing effort was considered in our models (for both males and females) owing to the lack of detailed and reliable fleet information (with the same spatial and temporal coverage) that could support this hypothesis. On the contrary, biological (biomass, spawning biomass at length and Kn) and environmental factors (NAO and SST) were considered as drivers of the variations in the decrease of L 50 for each sex. Several studies have addressed the changes in the maturity ogives in different species (Haug and Tjemsland 1986, Junquera et al. 1999, Marteinsdottir and Begg 2002, Engelhard and Heino 2004. In particular, Dominguez-Petit et al. (2008) studied these changes for European hake, although males were not considered in their study. They observed for the southern stock that the size at maturity increased even when total biomass and spawning biomass continued to decrease, which is contrary to the compensatory theory, probably as a consequence of environmental drivers such as NAO and upwelling. However, their results indicate that the pattern of the size at maturity in the Bay of Biscay (northern stock) was different, as a steady decline was observed throughout the study period.
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 density-dependent 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 Liermann and Hilborn (2001), there are four mechanisms of depensation: i) reduced probability of fertilization, ii) altered group dynamics, iii) environmental influences and iv) predator saturation. These four mechanisms are responsible for the lack of a compensatory response, although population density and intraspecific competition decrease.
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 (Albo-Puigserver et al. 2021). L 50 is expected to be more variable than other traits, particularly in heterogeneous environments (Hidalgo et al. 2014).
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 (Greene andPershing 2000, Köster et al. 2005). Effects of the NAO ripple influence life history traits and population dynamics and do so across trophic levels from primary production to predators (Ottersen et al. 2001). In fact, many examples show that NAO drives the life history and dynamics of many fish populations (Ottersen et al. 2001, Hjermann et al. 2004, Sullivan and Cowen 2005 and in particular of the European hake stocks (Meiners-Mandujano 2007, Dominguez-Petit et al. 2008, Goikoetxea and Irigoien 2013. Similarly, Drinkwater (2005) showed that NAO explains approximately 50% of the variability in the growth increase of northern cod between ages three and five in Newfoundland, which is undoubtedly a determinant of L 50 , as these are the ages at which cod mature in this area. For female hake something similar may occur: positive NAO values imply higher growth, which implies higher L 50 . The NAO index may influence the maturation process of hake indirectly through alterations in ecosystem composition and resource availability (Köster et al. 2005, Kell et al. 2005.

CONCLUSIONS
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 (Dominguez-Petit et al. 2022).

Correlation and multicollinearity analysis
In this section, data exploration techniques were presented to identify possible correlation and collinearity between the explanatory variables (Zuur et al. 2010) used for males and females in L 50 GAM models. In particular, correlation among variables was checked by performing a Pearson correlation test with the corrplot package (Wei et al. 2017) of the R software. Collinearity was tested by computing the GVIFs, which are the VIF values corrected by the number of degrees of freedom of a predictor variable (Fox and Weisberg 2017).

Missing values
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.

Backward procedure of models
Covariables for each final GAM were selected with a backward stepwise procedure based on the AIC and adjusted R-square. The best (and most parsimonious) model was finally chosen based on the compromise between low AIC values, high adjusted R-square and significant predictors. Table S4. -Comparison between the generalized additive models of the L 50 of European hake males during the study period . NAO, North Atlantic Oscillation; SST, sea surface temperature; SSB_length, spawning stock biomass at length; BIO, biomass; Kn, relative condition factor; D2 (%), deviance explained; AIC, Akaike information criterion. With respect to the model, s() indicates smooth function. bs specifies the type of function we are going to approximate, in this case we select "cr" cubic regression; K indicates the degree of the polynomial of the curve; and REML is the restricted maximum likelihood used to estimate these parameters.

Residual analysis for generalized additive models
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. The lower graphs represent the simple and partial correlations of the residuals.