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. 2008Dominguez-Petit R., Korta M., Saborido-Rey F., et al. 2008. Changes in size at maturity of European hake Atlantic populations in relation with stock structure and environmental regimes. J. Mar. Syst. 71: 260-278. https://doi.org/10.1016/j.jmarsys.2007.04.004 ). The estimation of the size at which 50% of individuals are mature (L50) can be used to define the minimum legal landing size for exploited stocks. However, in most stock assessments the causes of L50 variability are not analysed, though there are intrinsic and extrinsic factors affecting maturity that may act at a temporal scale (Lorenzen and Camp 2019Lorenzen K., Camp E.V. 2019. Density-dependence in the life history of fishes: when is a fish recruited?. Fish. Res. 217: 5-10. https://doi.org/10.1016/j.fishres.2018.09.024 ). Size and age at maturity are especially sensitive to population abundance (Dominguez-Petit et al. 2008Dominguez-Petit R., Korta M., Saborido-Rey F., et al. 2008. Changes in size at maturity of European hake Atlantic populations in relation with stock structure and environmental regimes. J. Mar. Syst. 71: 260-278. https://doi.org/10.1016/j.jmarsys.2007.04.004 , Lorenzen and Camp 2019Lorenzen K., Camp E.V. 2019. Density-dependence in the life history of fishes: when is a fish recruited?. Fish. Res. 217: 5-10. https://doi.org/10.1016/j.fishres.2018.09.024 ) 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. 2003Ali M., Nicieza A., Wootton R.J. 2003. Compensatory growth in fishes: a response to growth depression. Fish. Fish. 4: 147-190. https://doi.org/10.1046/j.1467-2979.2003.00120.x ). However, some environmental and anthropogenic drivers, including high rates of fishing mortality, may also impact on size at maturity (Pepin 2015Pepin P. 2015. Reconsidering the impossible - linking environmental drivers to growth, mortality, and recruitment of fish. Can. J. Fish. Aquat. Sci. 73: 205-215. https://doi.org/10.1139/cjfas-2015-0091 ). 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. 2018Hollins J., Thambithurai D., Koeck B., et al. 2018. A physiological perspective on fisheries-induced evolution. Evol. Appl. 11: 561-576. https://doi.org/10.1111/eva.12597 ). Subsequently, understanding mechanisms driving maturity variability within species could provide preliminary evidence of the species’s adaptability to climate change (Dalgleish et al. 2010Dalgleish H.J., Koons D.N., Adler P.B. 2010. Can life-history traits predict the response of populations to changes in climate variability? J. Ecol. 98: 209-217. https://doi.org/10.1111/j.1365-2745.2009.01585.x ), habitat degradation (Öckinger et al. 2010Öckinger E., Schweiger O., Crist T.O., et al. 2010 Life-history traits predict species responses to habitat area and isolation: a cross-continental synthesis. Ecol. Lett. 13: 969-979. https://doi.org/10.1111/j.1461-0248.2010.01487.x ) and overfishing (Hobday et al. 2011Hobday A.J., Smith A.D.M., Stobutzki I.C., et al. 2011. Ecological risk assessment for the effects of fishing. Fish Res. 108: 372-384. https://doi.org/10.1016/j.fishres.2011.01.013 ).
Barot et al. (2004)Barot S, Heino M, O’Brien L, Dieckmann U. 2004. Estimating reaction norms for age and size at maturation when age at first reproduction is unknown. Evol. Ecol. Res. 6: 659-678. 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. 2005Olsen EM, Lilly GR, Heino M, et al. 2005. Assessing changes in age and size at maturation in collapsing populations of Atlantic cod (Gadus morhua). Can. J. Fish. Aquat. Sci. 62:811-823. https://doi.org/10.1139/f05-065 ). 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. 1997Trippel E.A., Kjesbu O.S., Solemdal P. 1997. Effects of adult age and size structure on reproductive output in marine fishes. In Chambers R.C., Trippel E.A. (eds), Early life history and recruitment in fish populations. Chapman and Hall, London, U.K, pp. 31-62. https://doi.org/10.1007/978-94-009-1439-1_2 , Marteinsdottir and Begg 2002Marteinsdottir G., Begg G.A. 2002. Essential relationships incorporating the influence of age, size and condition on variables required for estimation of reproductive potential in Atlantic cod Gadus morhua. Mar. Ecol. Prog. Ser. 235:235-256. https://doi.org/10.3354/meps235235 ). 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 2021ICES. 2021. Working Group for the Bay of Biscay and the Iberian Waters Ecoregion (WGBIE).ICES Sci. Rep. 3:48.). 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 (Fmsy) (ICES 2021ICES. 2021. Working Group for the Bay of Biscay and the Iberian Waters Ecoregion (WGBIE).ICES Sci. Rep. 3:48.). 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 (1982-2019) 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)Cerviño S., Domínguez R., Jardim E., et al. 2013. Impact of egg production and stock structure on MSY reference points. Implications for Southern hake management. Fish. Res. 138: 168-178. https://doi.org/10.1016/j.fishres.2012.07.016 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. 2008Dominguez-Petit R., Korta M., Saborido-Rey F., et al. 2008. Changes in size at maturity of European hake Atlantic populations in relation with stock structure and environmental regimes. J. Mar. Syst. 71: 260-278. https://doi.org/10.1016/j.jmarsys.2007.04.004 , Murua 2010Murua H. 2010. The biology and fisheries of European hake, Merluccius merluccius, in the north-east Atlantic. Adv. Mar. Biol. 58: 97-154. https://doi.org/10.1016/B978-0-12-381015-1.00002-2 ), 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 L50 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 L50 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.
MATERIALS AND METHODS
⌅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 2019ICES. 2019. Working Group for the Bay of Biscay and the Iberian Waters Ecoregion (WGBIE). ICES Sci. Rep. 1:31.) 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)Ashton W.D. 1972. The logit transformation with special reference to its uses in bioassay. Haffner Publishing Co., INC., New York, 88 pp.:
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 (mature-immature).
The length at which 50% of specimens were mature (L50, size at maturity) was estimated as
where and are the intercept and slope estimates, respectively, estimated using the SizeMat (Torrejón-Magallanes 2020Torrejón-Magallanes E.J.2020. sizeMat: Estimate Size at Sexual Maturity. R package version 1.1.2.) package of the R software (R Core Team 2021R Core Team. 2021. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Version 4.1.2. Vienna, Austria.).
In addition, the classical logistic models including and not including sex as a covariable and were compared through a Chisq likelihood-ratio test (P-value<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 2017Moritz S., Bartz-Beielstein T. 2017. Imputets: Time Series Missing Value Imputation in R. R J. 9: 207-218. https://doi.org/10.32614/RJ-2017-009 , 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 L50 variability. Among the several indices proposed to measure body condition, Le Cren’s relative condition factor (Kn) (Le Cren 1951Le Cren E.D. 1951. The length-weight relationship and seasonal cycle in gonad weight and condition in the perch. J. Anim. Ecol. 2: 201-219. https://doi.org/10.2307/1540 ) 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 length-weight relationship parameters were estimated using the lm function of the R software through the linear regression model . 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 2017Moritz S., Bartz-Beielstein T. 2017. Imputets: Time Series Missing Value Imputation in R. R J. 9: 207-218. https://doi.org/10.32614/RJ-2017-009 ).
Environmental data
⌅In order to analyse the effect of the environmental conditions on the L50, 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 large-scale climatic events affect not only temperature but also ocean masses and circulation, which affect plankton and zooplankton (Nye et al. 2009Nye J.A., Link J.S., Hare J.A., Overholtz W.J. 2009. Changing spatial distribution of fish stocks in relation to climate and population size on the Northeast United States continental shelf. Mar. Ecol. Prog. Ser. 393, 111-129. https://doi.org/10.3354/meps08220 , 2014Nye J. A., Baker M.R., Bell R., et al. 2014. Ecosystem effects of the Atlantic Multidecadal Oscillation. J. Mar. Syst. 133, 103-116. https://doi.org/10.1016/j.jmarsys.2013.02.006 ). These indices have disproportionate effects on the early life stages of fishes because of alterations in food sources and temperature (Nye et al. 2009Nye J.A., Link J.S., Hare J.A., Overholtz W.J. 2009. Changing spatial distribution of fish stocks in relation to climate and population size on the Northeast United States continental shelf. Mar. Ecol. Prog. Ser. 393, 111-129. https://doi.org/10.3354/meps08220 , 2014Nye J. A., Baker M.R., Bell R., et al. 2014. Ecosystem effects of the Atlantic Multidecadal Oscillation. J. Mar. Syst. 133, 103-116. https://doi.org/10.1016/j.jmarsys.2013.02.006 ). SST was included because it affects the fish community through thermal tolerance and species migration (Marshall and Elliott 1998Marshall S., Elliott M. 1998. Environmental influences on the fish assemblage of the Humber estuary. U.K. Estuar. Coast. Shelf Sci. 46:175-184. https://doi.org/10.1006/ecss.1997.0268 ). 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 (1982-2019) 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 2017Wood S.N. 2017. Generalized Additive Models: An Introduction with R (2nd edition). Chapman and Hall/CRC, New York, 496 pp.) 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 L50.
GAMs are semi-parametric extensions of generalized linear models (GLMs) for which the strictly linear predictor
where X=(X1,...,Xp) are covariables, 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 fj (Xj) is the unknown smooth partial effect of Xj 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 (L50).
Standardized data exploration techniques were used to identify any outliers and any correlation and collinearity between the explicative variables (Zuur et al. 2010Zuur A.F., Ieno E.N., Elphick C.S. 2010. A protocol for data exploration to avoid common statistical problems. Methods. Ecol. Evol.1: 3-14. https://doi.org/10.1111/j.2041-210X.2009.00001.x ). In particular, correlation among variables was checked using a Pearson’s correlation test with the corrplot package (Wei et al. 2017Wei T., Simko V., Levy M., et al. 2017. Package ‘corrplot’. J. Am. Stat. 56: 316-324.) 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. 2017Wei T., Simko V., Levy M., et al. 2017. Package ‘corrplot’. J. Am. Stat. 56: 316-324.) 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 2019Fox J., Weisberg S. 2019. An R Companion to Applied Regression (Third). SAGE Publications Inc, pp. 608.), 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 2006Wood S.N. 2006. Generalized Additive Models: An Introduction with R. Chapman & Hall/CRC. Series: Chapman & Hall/CRC. Texts in Statistical Science.). 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 2011Wood S.N. 2011. Mgcv: GAMs with GCV/AIC/REML smoothness estimation and GAMMs by REML/PQL.) of the R software. The R code for the applied models can be found on the GitHub repository (https://github.com/Lojamodav/PAPER.git).
RESULTS
⌅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 L50=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 L50 time series showed a decreasing trend over the whole period (Fig. 2A), while the female L50 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 L50 for the remaining years. The male L50 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 L50. However, the female L50 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 L50.
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 L50
⌅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).
| Parametric coefficients | ||||
|---|---|---|---|---|
| Estimate | se | t value | Pr (>|t|) | |
| 31.0263 | 0.4684 | 66.24 | <2e-16 *** | |
| Significance of smooth functions | ||||
| Edf | Ref. df | F | P-value | |
| (year) | 1.616 | 1.969 | 22.77 | 1.1e-06 *** |
The model results show that 55.3% of the L50 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 L50 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.
| Parametric coefficients | ||||
|---|---|---|---|---|
| Estimate | Std error | t value | Pr (>|t|) | |
| 44.4441 | 0.3035 | 146.4 | <2e-16 *** | |
| Significance of smooth functions | ||||
| Edf | Ref. df | F | P-value | |
| (year) | 2.595 | 3.245 | 2.541 | 0.04783 * |
| (biomass) | 3.198 | 3.634 | 4.992 | 0.00474 ** |
| (Spawning biomass at length) | 2.249 | 2.788 | 7.734 | 0.00144 ** |
| (NAO) | 1.000 | 1.000 | 8.094 | 0.00821 ** |
In particular, year showed a dome-shaped relationship with L50 estimates (Fig. 5). Density-dependent effects were detected. The results showed a positive relationship of biomass with L50 estimates and a negative relationship of spawning biomass at length with L50 estimates (Fig. 5). Finally, the NAO index also showed a slight positive relationship with L50 (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.
DISCUSSION
⌅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 (L50) 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.
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 L50 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ø and Haug 1999Godø O.R., Haug T. 1999. Growth rate and sexual maturity in cod (Gadus morhua) and Atlantic halibut (Hippoglosus hippoglossus). J. Northwest Atl. Fish. Sci. 25: 115-123. https://doi.org/10.2960/J.v25.a10 , Oven 2004Oven L.S. 2004. Resorption of Vitellogenous Oocytes as an Indicator of the State of the Black Sea Fish Populations and Their Environment. Journal of Ichthyology, 44: 115-119.). 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. 2004Olsen E.M., Heino M., Lilly G.R., et al. 2004. Maturation trends indicative of rapid evolution preceded the collapse of northern cod. Nature. 428: 932-935. https://doi.org/10.1038/nature02430 ); 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 1995Trippel E.A. 1995. Age at maturity as a stress indicator in fisheries. Bioscience. 45:759-771. https://doi.org/10.2307/1312628 , Jørgensen et at. 2007Jørgensen C., Enberg K., Dunlop E.S., et al. 2007. Ecology: Managing Evolving Fish Stocks. Sci. 318: 1247-1248. https://doi.org/10.1126/science.1148089 , Hidalgo et al. 2012Hidalgo M., Rouyer T., Bartolino V., et al. 2012. Context-dependent interplays between truncated demographies and climate variation shape the population growth rate of a harvested species. Ecography 35: 637-649. https://doi.org/10.1111/j.1600-0587.2011.07314.x ).
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 2007Dominguez-Petit R. 2007. Study of reproductive potential of Merluccius merluccius in the Galician shelf. Doctoral Thesis. University of Vigo, Spain., Korta et al. 2010Korta M., Domínguez-Petit R., Murua H., Saborido-Rey F. 2010. Regional variability in reproductive traits of European hake Merluccius merluccius L. populations. Fish. Res. 104: 64-72. https://doi.org/10.1016/j.fishres.2009.03.007 ). 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 (Fmsy= 0.25) but not above the precautionary limits (Fpa=0.75) (ICES 2019ICES. 2021. Working Group for the Bay of Biscay and the Iberian Waters Ecoregion (WGBIE).ICES Sci. Rep. 3:48.). 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. 2015Korta M., García, D., Santurtún M., et al. 2015. “European Hake (Merluccius merluccius) in the North-east Atlantic,” in Hakes: biology and Explotation, ed. H. Arancibia (Hoboken: John Wiley & Sons, Ltd), 1-37. https://doi.org/10.1002/9781118568262.ch1 , ICES 2019ICES. 2021. Working Group for the Bay of Biscay and the Iberian Waters Ecoregion (WGBIE).ICES Sci. Rep. 3:48.) 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)Dominguez-Petit R. 2007. Study of reproductive potential of Merluccius merluccius in the Galician shelf. Doctoral Thesis. University of Vigo, Spain. and Korta et. al (2015)Korta M., García, D., Santurtún M., et al. 2015. “European Hake (Merluccius merluccius) in the North-east Atlantic,” in Hakes: biology and Explotation, ed. H. Arancibia (Hoboken: John Wiley & Sons, Ltd), 1-37. https://doi.org/10.1002/9781118568262.ch1 , 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 2000Law R. 2000. Fishing, selection, and phenotypic evolution. ICES J. Mar. Sci. 57: 659-668. https://doi.org/10.1006/jmsc.2000.0731 ), because the size/age truncation of the stock leads to a decrease in reproductive potential (Hixon et al. 2014Hixon M.A., Johnson D.W., Sogard S.M. 2014. BOFFFFs: on the importance of conserving old-growth age structure in fishery populations. ICES J. Mar. Sci. 71: 2171-2185. https://doi.org/10.1093/icesjms/fst200 ). 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 L50 for each sex.
Several studies have addressed the changes in the maturity ogives in different species (Haug and Tjemsland 1986Haug T., Tjemsland T. 1986. Changes in size and age distribution and age at sexual maturity in Atlantic Halibut, Hippoglossus hippoglossus, caught in North Norwegian waters. Fish. Res. 4: 145-155. https://doi.org/10.1016/0165-7836(86)90039-1 , Junquera et al. 1999Junquera S., Roman E., Paz X., Ramilo G. 1999. Changes in Greenland halibut growth, condition and fecundity in the Northwest Atlantic (Flemish Pass, Flemish Cap and southern Grand Banks). Variations in maturation, growth, condition and spawning stock biomass production in groundfish. J. Northwest Atl. Fish. Sci. 25: 17-28. https://doi.org/10.2960/J.v25.a2 , Marteinsdottir and Begg 2002Marteinsdottir G., Begg G.A. 2002. Essential relationships incorporating the influence of age, size and condition on variables required for estimation of reproductive potential in Atlantic cod Gadus morhua. Mar. Ecol. Prog. Ser. 235:235-256. https://doi.org/10.3354/meps235235 , Engelhard and Heino 2004Engelhard G.H., Heino M. 2004. Maturity changes in Norwegian spring spawning herring Clupea harengus: compensatory or evolutionary responses? Mar. Ecol. Prog. Ser. 272: 245-256. https://doi.org/10.3354/meps272245 ). In particular, Dominguez-Petit et al. (2008)Dominguez-Petit R., Korta M., Saborido-Rey F., et al. 2008. Changes in size at maturity of European hake Atlantic populations in relation with stock structure and environmental regimes. J. Mar. Syst. 71: 260-278. https://doi.org/10.1016/j.jmarsys.2007.04.004 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 L50for 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 L50, 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 L50. However, the results showed a positive relationship between biomass at spawning and L50, which is contrary to the compensatory theory. The lack of compensatory response is known as depensation. According to Liermann and Hilborn (2001)Liermann M., Hilborn R. 2001. Depensation: evidence, models and implications. Fish. Fish. 2: 33-58. https://doi.org/10.1046/j.1467-2979.2001.00029.x , 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 L50 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. 2021Albo-Puigserver M., Pennino M.G., Bellido J.M., et al. 2021. Changes in life history traits of small pelagic fish in the western Mediterranean Sea. Front. Mar. Sci. 8: 1197. https://doi.org/10.3389/fmars.2021.570354 ). L50 is expected to be more variable than other traits, particularly in heterogeneous environments (Hidalgo et al. 2014Hidalgo M., Rouyer T., Molinero J.C., et al. 2014. Contrasting evolutionary demography induced by fishing: The role of adaptive phenotypic plasticity. Ecol. App. 24:1101-1114. https://doi.org/10.1890/12-1777.1 ).
The only significant environmental variable included in our model for European hake females was the NAO index, which showed a positive relationship with the L50. 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 and Pershing 2000Greene C.H., Pershing A.J. 2000. The response of Calanus finmarchicus populations to climate variability in the Northwest Atlantic: basin-scale forcing associated with the North Atlantic Oscillation. ICES J. Mar. Sci. 57: 1536-1544. https://doi.org/10.1006/jmsc.2000.0966 ,Köster et al. 2005Köster F.W., Möllmann C. Hinrichsen H.H., et al. 2005. Baltic cod Recruitment - the impact of climate variability on key processes. ICES J. Mar. Sci. 62: 1408-1425. https://doi.org/10.1016/j.icesjms.2005.05.004 ). 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. 2001Ottersen G., Planque B., Belgrano A., et al. 2001. Ecological effects of the North Atlantic Oscillation. Oecologia. 128: 1-14. https://doi.org/10.1007/s004420100655 ). In fact, many examples show that NAO drives the life history and dynamics of many fish populations (Ottersen et al. 2001Ottersen G., Planque B., Belgrano A., et al. 2001. Ecological effects of the North Atlantic Oscillation. Oecologia. 128: 1-14. https://doi.org/10.1007/s004420100655 , Hjermann et al. 2004Hjermann D.Ø., Stenseth N.C., Ottersen G. 2004. Indirect climatic forcing of the Barents Sea capelin: a cohort effect. Mar. Ecol. Prog. Ser. 273: 229-238. https://doi.org/10.3354/meps273229 , Sullivan and Cowen 2005Sullivan M.C., Cowen R.K., Steves B.P. 2005. Evidence for atmosphere-ocean forcing of yellowtail flounder (Limanda ferruginea) recruitment in the Middle Atlantic Bight. Fish Oceanogr. 14: 386-399. https://doi.org/10.1111/j.1365-2419.2005.00343.x ) and in particular of the European hake stocks (Meiners-Mandujano 2007Meiners-Mandujano C.G. 2007. Importancia de la variabilidad climática en las pesquerías y biología de la merluza europea Merluccius merluccius (Linnaeus, 1758) de la costa Noroccidental Africana. PH.D. Thesis. Universitat Politècnica de Catalunya (UPC) (Spain). 207 pp., Dominguez-Petit et al. 2008Dominguez-Petit R., Korta M., Saborido-Rey F., et al. 2008. Changes in size at maturity of European hake Atlantic populations in relation with stock structure and environmental regimes. J. Mar. Syst. 71: 260-278. https://doi.org/10.1016/j.jmarsys.2007.04.004 , Goikoetxea and Irigoien 2013Goikoetxea N., Irigoien X. 2013. Links between the recruitment success of northern European hake (Merluccius merluccius L.) and a regime shift on the NE Atlantic continental shelf. Fish. Oceanogr. 22: 459-476. https://doi.org/10.1111/fog.12033 ). Similarly, Drinkwater (2005)Drinkwater K.F. 2005. The response of Atlantic cod (Gadus morhua) to future climate change. ICES J. Mar. Sci. 62: 1327-1337. https://doi.org/10.1016/j.icesjms.2005.05.015 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 L50, 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 L50. The NAO index may influence the maturation process of hake indirectly through alterations in ecosystem composition and resource availability (Köster et al. 2005Köster F.W., Möllmann C. Hinrichsen H.H., et al. 2005. Baltic cod Recruitment - the impact of climate variability on key processes. ICES J. Mar. Sci. 62: 1408-1425. https://doi.org/10.1016/j.icesjms.2005.05.004 , Kell et al. 2005Kell L.T., Pilling G.M., O’Brien C.M. 2005. Implications of the climate change for the management of North Sea cod (Gadus morhua). ICES J. Mar. Sci. 62: 1483-149. https://doi.org/10.1016/j.icesjms.2005.05.006 ).
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 L50 for both sexes, which was the objective of this study.
However, determining which variables explain the decreasing trends in L50for both sexes requires applying more complex models than an overall linear trend to all the data, because there are much more complex relationships between L50and the factors that may explain the variability in L50.
Total biomass, spawning biomass at length and the NAO partially explain the decline in the L50 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 L50, 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 L50 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. 2022Dominguez-Petit R., García-Fernandez C., Leonarduzzi E., et al. 2022. Parental effects and reproductive potential of fish and marine invertebrates: Cross-generational impact of environmental experiences. In: Domínguez-Petit R. (ed). Impact of Environmental Stress on Reproductive Processes of Aquatic Animals. Fishes 7: 188. https://doi.org/10.3390/fishes7040188 ).