Spatial and Temporal Variability in the Occurrence and Abundance of European Hake Larvae, Merluccius merluccius, on the Galician Shelf (NE Atlantic)

The European hake (Merluccius merluccius) is represented as one of the most valuable fisheries in the Galician shelf. We analyzed the distribution, abundance, and environmental conditions of the southern-stock European hake larvae from the Galician shelf during the two main spawning peaks, winter-spring and summer, based on the data from three ichthyoplankton surveys (March 2012, March 2017, and June 2017). A total of 395 larvae in March 2012, 121 in March 2017, and 69 in June 2017 were captured. The northeast section of the study area, close to Estaca de Bares, primarily between 100 and 200 m isobaths, had the highest presence of the European hake larvae in all surveys. Generalized additive models (GAMs) indicated that the occurrence of larvae was significantly different between the surveys and was associated negatively with the temperature, while the abundance of larvae was significantly different between sampling years and was the highest at a temperature around 13.36°C and at sea surface heights of about −0.48 m. Studies of the distribution of early life stages and their relation to external conditions are essential to the understanding of the complex process of recruitment, especially in the exploited species and in highly dynamic environments like the Galician shelf.


INTRODUCTION
Larval dynamics play an important role in the regulation of fish populations. Early life history (ELH) stages, i.e., eggs and larvae, are characterized as the most vulnerable phase of the fish life cycle (Hjort, 1914) due to the existence of a high mortality phase when > 99% of larvae usually die (Houde, 1987).
Numerous factors affect the survival of the offspring at ELH stages, both endogenous and exogenous. Maternal attributes affect survivorship through offspring quality (Green, 2008;García-Fernández et al., 2020), while oceanographic conditions affect the survivorship through offspring drift, distribution, physiological rates, or food availability (Agostini et al., 2006;Jansen et al., 2015) collectively impacting recruitment (Landaeta and Castro, 2012). Several environmental factors that operate at different scales (from micro to ocean basin), such as temperature, salinity, dissolved oxygen, water turbulence, and currents, control directly or indirectly the survival of the offspring (Govoni, 2005;Agostini et al., 2006;Houde, 2009;Jansen et al., 2015). So, favorable environmental conditions during reproduction and ELH stages generate larger cohorts of offspring due to higher survival rates (Stenseth et al., 2002). Fluctuations in larval mortality due to variability in conditions are considered one of the major causes of variability in marine fish production regulating year-class size (Houde, 1987;Chambers and Trippel, 1997).
Offspring survival is directly related to recruitment, defined as the number of juveniles in a population that annually survives the egg and larval stages (Miller and Kendall, 2009). It means that poor offspring survival frequently generates weak recruitment while a good survival can lead to strong recruitment (Houde, 2009). Indeed, recruitment success variability is one of the most influential processes in regulating exploited populations, and as such, the understanding of the recruitment dynamics is essential for stock management (Cushing, 1990) including recruitment forecasting (Chambers and Trippel, 1997). The cornerstone in fisheries management is the stock-recruitment relationship which is often impaired. This mismatch is mainly caused by the lack of proportionality between spawning stock biomass and larval survival due to parental effects (stock reproductive potential) and the environmental conditions which strongly modulate recruitment strength (Lasker, 1981;Cury and Roy, 1989;Trippel et al., 1997). Thus, the knowledge of all mechanisms that affect ELH stages and their survival is necessary to understand the interannual recruitment variability (Houde, 1987). Such information highlights the potential role of larval dynamics as a tool in fisheries management, in addition to simply increasing our understanding of fish biology and ecology (Fuiman, 2002).
The European hake (Merluccius merluccius) is one of the most important species in marine demersal ecosystems of temperate regions in the northeast Atlantic Ocean. It is abundant and widely distributed in the regions, from Norway to the Gulf of Guinea, in addition to the Mediterranean and the Black Sea. The European hakes inhabit coastal and oceanic waters at depths between 70 and 370 m with water temperature between 10 and 14 • C (Cohen et al., 1990). They are top predators and are very plastic species, having a high capacity to adapt their reproductive potential to environmental conditions and energy availability (Domínguez-Petit et al., 2008;Domínguez-Petit and Saborido-Rey, 2010). The species also exhibits spatio-temporal variability in the reproductive traits even at the regional level (Korta et al., 2010). The European hake is one of the most economically important species in Europe (Casey and Pereiro, 1995). For assessment purposes, the species is divided into the northern and southern stocks, the latter being the focus of this study (Korta et al., 2015). Both stocks had been exploited outside safe biological limits since the 1990s when their highest landings were recorded. In 2004, it was declared to be in critical status due to overfishing (Korta et al., 2010). Subsequently, in 2005, a recovery plan was implemented [European Council [EC], 2005]. Moreover, both the stocks had a recruitment pulse in the second half of the 2000s. Currently, it is known that the northern stock is at full reproductive capacity while the southern stock is in healthy status but still shows reduced reproductive capacity. Indeed, it is still overexploited (in 2018, Fishing mortality = 0.60), exceeding the maximum rate of fishing mortality (0.25) but not above the precautionary limits (0.75) [International Council for the Exploration of the Sea [ICES], 2019].
Along the NW Iberian Peninsula (the northern limit and the best studied sector of the southern stock), the species has a very protracted spawning season and high population reproductive asynchronicity, with two or occasionally three peaks in spawning (Domínguez-Petit, 2007;Serrat et al., 2019). In general, the primary peak in spawning takes place in late winter and early spring, around March; the secondary one in late spring-summer, around June; and the third one around October (although there is greater variability in both the intensity and timing among years). Nevertheless, it remains unknown to what extent these three peaks correspond (or not) to different spawning stock components (Serrat et al., 2019). Indeed, the spawning seasonality plays a fundamental role in hake larval dynamics due to the presence of different intra-annual recruitment pulses in the target area.
The Galician shelf is the main spawning area for the southern stock of the European hake (Sánchez and Gil, 2000). This region is located off the NW Iberian Peninsula, an area with complex bathymetry and circulation over the shelf and slope regions (Torres and Barton, 2007). It is influenced by the Iberian Poleward Current, an unstable poleward flow that transports warm and salty waters and yields abundant mesoscale features, such as eddies and fronts (Peliz et al., 2003). In this NE Atlantic zone, there are two major water masses: the upper layer (down to 1,000 m depth) is characterized by the Eastern North Atlantic Central Water (ENACW) with temperatures between 10.5 and 12 • C and salinity between 35.45 and 35.75 while waters below 1,000 m are characterized by the Mediterranean Sea run-off water with warmer temperatures (12.7 • C) and higher salinity concentrations (38.4) (Valencia and Franco, 2004). Galician waters are characterized as highly productive, as both the latitude and oceanographic conditions promote large phytoplankton blooms in the spring and in the summer (Casas et al., 1999). In addition, the region is influenced by the Eastern North Atlantic Upwelling System that produces the renewal of water masses and the advection at the coast of cool and enriched waters from deeper areas. This process presents seasonality from winter downwelling to summer upwelling, which begins around April and peaks in mid-summer (Fraga, 1981). Another important supply of nutrients and organic matter is the freshwater input from the continental run-off (Villegas- Ríos et al., 2011).
The aim of this study was to understand the main environmental drivers that define the spatio-temporal variability of the presence and abundance of European hake larvae in the highly dynamic environment of the Galician shelf. The study was based on three ichthyoplankton surveys within two main spawning periods in 2012 and 2017 (winter-spring and summer). Results will help to elucidate the larval European hake distributions in different spawning seasons and principal mechanisms that determine such distributions, as well as the potential impact of implications for hake recruitment dynamics.

Sample Collection and Treatment
Three different surveys were carried out off the NW Iberian Peninsula, NE Atlantic (Figure 1) from Finisterre to Estaca de Bares: (i) in the late winter of 2012 (from February 29th to March 12th), (ii) in the late winter of 2017 (from March 14th to March 20th), and (iii) in the summer of 2017 (from June 18th to June 24th). The timing of the sampling was based on the knowledge of the spawning activity of the European hake from previous studies (Domínguez-Petit, 2007). Larval collections were performed on a predefined grid along transects perpendicular to the coast with a transect spacing of 8 nm, and stations along each transect of 4 nm apart (Figure 2). In March 2012, 79 stations were sampled, and 86 stations were sampled in both surveys in 2017. Ichthyoplankton samples were collected with different gears in 2012 and 2017. In 2012, the samples were collected using a multiple opening/closing net MultiNet Midi (0.25 m 2 opening with five nets of 200-µm mesh) while in 2017, a 60 cm diameter Bongo net (0.28 m 2 , 500µm mesh) was used. While the change in the gear type might have led to changes in catchability, we expect this to be minimal for multiple reasons. First, both the gears were towed obliquely covering the entire vertical distribution range of the larvae and the operational procedures (tow speed, towing depth and winch rate) were identical. The volume filtered between the gears was also similar, suggesting that there were minimal clogging differences despite the differences in the mesh size. The difference in the mesh size can also affect the filtration performance and the clogging of the plankton nets, which is an important aspect because it can affect the efficiency of the captures. This can be detected if the filtered volume is significantly lower than expected. This has not been observed in the samples. Mesh size differences will have a significant effect on the plankton selection since mesh width is the most important measurement affecting this. However, it remains possible that there could be differences in catchability due to head wake and the angle of opening which may be important in interpreting inter-survey differences. That is not an issue considering the size of hake larvae dimensions (Palomera et al., 2005) and that 500 µm is considered adequate for sampling fish larvae (Rodriguez et al., 2017). In all surveys, sampling was performed to a maximum depth of 200 or 5 m above the bottom in shallower areas. The volume of the filtered water was measured by an electronic flowmeter located at the mouth of the net and it was similar in all surveys. During the survey, the hake larvae were identified using morphological characteristics, photographed, and preserved for later analyses. In the laboratory, zooplankton abundances (n/m 3 ) were semi-automatically counted using an imaging analyzer system (Bachiller and Fernandes, 2011), including individuals ranging from 0.20 to 2.0 mm based on the net mesh.
Environmental parameters were obtained from in situ CTD data and HYCOM. At each station, a Sea-Bird SBE25 CTD was deployed to record the vertical profiles of temperature and salinity from the surface to 5 m above the bottom (or 200 m at deeper stations). The daily average sea surface height (SSH) for each station in the three surveys was extracted from the Hybrid coordinate ocean model (HYCOM) 1 (1/12 degrees resolution) for examining how oceanographic features, such as cyclonic and anticyclonic eddies influence the distribution of European hake larvae. In addition, estimates of geostrophic meridional and zonal currents (u and v vectors) were calculated from interpolated fields of dynamic height throughout the region. Eddy kinetic energy (m 2 /s 2 ) is a proxy to measure the energy in the surface ocean (Zainuddin et al., 2006) calculated using data from geostrophic velocities as: where EKE is the eddy kinetic energy, u and v are the zonal and meridional velocity, respectively.

Data Analysis and Statistical Approach
A Fisher exact test and a Kruskal-Wallis test were performed to determine if surveys are different in terms of the presence/absence and the abundance of larvae. Afterward, a posterior multiple pairwise comparison between surveys was calculated to identify where pairwise differences occurred.
We used hurdle models to analyze the habitat associations of European hake larvae. These models consist of two parts, first, a binary model [presence-absence (P/A)] followed by a count model for stations where hake larvae were present (the relative abundance) (Cragg, 1971;Ver Hoef and Jansen, 2007). The parameters at each station that were considered in the model development were as follows: temperature ( • C) at 50 m, salinity (PSU) at 50 m, an abundance of zooplankton < 0.55 mm in size (number of individuals/m 3 ), bottom depth (m), SSH (m), and EKE (m 2 /s 2 ). Oceanographic data were considered at 50 m because previous studies suggested that the highest abundance of the European hake larvae is located at this depth stratum (Rodriguez et al., 2015a). Additionally, spatial (latitude and longitude) and temporal factors (month and year of sampling) were included to account for inter-survey differences that may be related to changes in the overall abundance and presence, in addition to the change in gear type. Zooplankton abundance was log-transformed before analyses. In the P/A model, the logtransformed water volume filtered (m 3 ) for each station was used as an offset to standardize and account for the variability in the sampling effort for all surveys. Moreover, the abundance of hake larvae was standardized and expressed as the number of individuals per 10 m 2 at each positive station for the abundance model, thus accounting for the variability in the volume filtered.
A correlation matrix was built using Pearson's correlation coefficients (r) for each set of covariates (Supplementary Figure 1) to determine and ultimately remove collinear predictors. Salinity was highly correlated with temperature (r = −0.67, p < 0.001 for P/A data, Supplementary Figure 1A; r = −0.80, p < 0.001 for the abundance data, Supplementary Figures 1B, 2) and, as such, only temperature was included in the model because temperature resulted in higher explanatory power. We also explored using density as a predictor, though decided to use temperature instead as temperature out-performed density in terms of deviance explained and log-likelihood in our models.
Generalized additive models (GAMs) with the ·· mgcv ·· R package (version 1.8.31; Wood, 2011) were used to explore the influence of environmental, spatial, and temporal conditions on European hake larvae P/A and abundance. GAMs are the commonly used method for modeling fish habitats because they can be used for elucidating complex non-linear relationships between the response and explanatory variables (Hastie and Tibshirani, 1990;Valavanis et al., 2008) that are not observable by using linear models and simple correlations (Kitchens and Rooker, 2014). The P/A models were fitted with a binomial distribution (with a logit link function), while abundance models were fitted with a negative binomial distribution (with a log link function) as abundance data were highly over-dispersed (Poisson overdispersion = 12.50). Full models included physical (temperature), oceanographic (SSH, EKE, depth), biological (zooplankton abundance in P/A models), temporal (month and year of the survey as an interaction for considering the temporal variation), and spatial factors (coordinates of points were included as smoothed terms to account for spatial autocorrelation). To minimize overfitting, four degrees of freedom were allowed for each spline(s) except for the latitude and longitude of the sampling stations (Rooker et al., 2012). A tensor product spline (te) was used for the latitude and longitude to allow for the terms to be penalized anisotropically. The best model for P/A and abundance was selected from the following full models, respectively: Hake larvae presence = offset (log (volume filtered)) + s (temperature at 50 m) + s (log− zooplankton abundance) + s (depth of the station) + s (sea surface height) +s (eddy kinetic energy) + te (longitude, latitude) + month × year.
Hake larvae abundance (n larvae/10 m 2 ) = s (temperature at 50 m) + s (log − zooplankton abundance) The Akaike Information Criterion (AIC) was used for the model selection (Akaike, 1974), with a manual backward stepwise procedure. The best model for P/A was tested with bootstrapping to make a receiver operating characteristic curve (ROC) and measure the area under the curve (AUC) (ROCR package, Sing et al., 2005). A random subsample of the stations (half of the total data) was used as a training data set (1,000 iterations) and the remaining data served as the test data set. With this analysis, the ROC curve of the true positive rate, the false positive rate, and the AUC value are estimated. The AUC is a proxy of the model goodness-of-fit, with values typically ranging from 0.5 to 1 (higher values suggest high discriminatory ability, 0.5 suggests no explanatory power) (Maggini et al., 2006).
Additionally, the within-data predictive capability of the larval abundance model was tested using leave-one-out cross-validation (LOOCV) (Mosteller and Tukey, 1968;Stone, 1974). Values range from 0 to 1 (higher values of r 2 represent better model predictive capability within the dataset). LOOCV iteratively fits the model to all response values except one term and then predicts the leftout response value. R 2 values are then calculated between the predicted and the observed values.
After the predictions of occurrence and abundance were estimated from the best fit models, we multiplied the predicted occurrence probability by the predicted abundance (if present) to estimate the final prediction of the expected value of abundance at each station per survey (using a constant offset based on the mean filtered volume for all stations). All statistical analyses were conducted with R software version 3.6.2 (R Core Team, 2019) considering a significance level of p < 0.05.

Survey Data
Overall, a total of 251 stations were sampled and 585 European hake larvae were collected during the three surveys in 2012 and 2017 ( Table 1). The number of larvae collected in each survey differed greatly with the lowest catch total occurring in June 2017 (n = 69), followed by the total number of larvae encounter in March 2017 (n = 121), while the largest number of hake larvae were caught in March 2012 (n = 395). There were notable differences in the larval occurrence and abundance among and within the surveys (Table 1 and  (Supplementary Figure 3). Figure 2 shows the horizontal distribution and abundance of hake larvae by station and survey.
Larval presence and abundance were negatively correlated with temperature and salinity (r = −0.67, p < 0.001 for P/A data; r = −0.80, p < 0.001 for abundance data), indicating an association with specific oceanographic features of water masses in each survey (with different oceanographic features). Thus, seawater in March 2012, which had the highest larval abundances, was generally cooler and more saline than the other two surveys (Supplementary Figure 2).

Environmental Data
Oceanographic conditions were highly variable among the surveys. The mean and SD of the environmental variables for each survey are listed in Table 1 and represented in Supplementary Figure 4. In March 2012, temperatures at 50 m exhibited a northsouth gradient with higher temperatures in the southwestern area and cooler in the northern area, but the temperature range was narrow (from 12.25 to 12.75 • C). Zooplankton abundances were nearly uniform, around 500 n/m 3 in the whole area except for one patch close to Estaca de Bares, where the abundance reached more than 3,000 n/m 3 . There was a flow of surface water along the shelf toward the northeastern region. Geostrophic velocities in the region were low (∼10 cm/s) and there was a gradient that resulted in an area of low dynamic heights at the inner shelf (<100 m) and higher currents offshore. Moreover, between Coruña and Estaca de Bares, there was a well-defined cyclonic structure, having the highest geostrophic velocities. During the March 2017 survey, the mean temperature at 50 m was higher than in the same period in 2012, but it was spatially more variable (from 12.75 to 13.4 • C). Warmer waters were located in the coastal areas, mainly in the northeastern area of the Galician shelf, while cooler water masses were observed offshore of that area and along the southwestern coast. Zooplankton abundances showed similar distribution patterns to those in March 2012 but with higher values, with a mean of 2,130 n/m 3 . There was a similar patch close to Estaca de Bares but further west than the one found in 2012, with more than 15,000 n/m 3 . In addition, the SSH showed higher values southward and currents presented a predominantly southwestern flow with low velocities (<10 cm/s). Lastly, in June 2017, temperatures were high in some specific offshore areas closer to the edge of the Galician shelf, and the whole area showed higher temperature variability (from 12.8 to 16 • C). Mean zooplankton abundance was similar mean as in March 2017 (2,380 n/m 3 ), but the highest abundances were along the coast. SSH during this survey was lower than in the rest of the surveys, and the water flow was not evident as in both March surveys. Even that, the water flow presented higher geostrophic velocities (<40 cm/s).

Presence/Absence Model
The best fit GAM for the European hake larval occurrence was composed of five variables: the temperature at 50 m, the interaction between latitude and longitude, and the interaction between month and year. The final model, AIC, and deviance explained (DE) were 255.81 and 41.20%, respectively ( Table 2). The variables in the order of the greatest AIC were latitude and longitude interaction ( AIC = 51.47), temperature ( AIC = 6.27), and the interaction between month and year ( AIC = 5.42).
Response plots indicate a negative relationship between hake larvae presence and temperature and higher occurrence in the north (between 8 • W and 7.5 • W of longitude and between 43.8 • N and 44 • N of latitude) and less occurrence in the southward (Figure 3). The interaction between the month and year underlined the existence of significant differences between all surveys. Finally, when the (probability of) occurrence in the sampling area was predicted for each station by using the best GAM model (1,000 iterations), the average AUC was 0.79 (median = 0.79 and S.D. = 0.04), i.e., the ROC curve for the best model showed that it has a moderately good predictive capability within the dataset (Supplementary Figure 5).

Abundance Model
The best abundance GAM for the European hake included five variables: the temperature at 50 m, SSH, year as a factor, and the interaction between the latitude and longitude. The final model AIC and DE were 964.78 and 58.60%, respectively ( Table 2). The variables in the order of the greatest AIC were longitude and latitude ( AIC = 23.61), year ( AIC = 14.78), SSH Parametric coefficients present estimate, standard error (Std. error), z value, and p-value; while smooth terms present estimated degrees of freedom (edf), Std. error, z value, and p-value. Moreover, Akaike Information Criterion ( AIC) was calculated for the environmental, temporal, and spatial variables and total AIC of the best P/A and Abundance models.
FIGURE 3 | Response plots for abiotic variables on the partial effects of the P/A of hake larvae from final generalized additive models (GAMs). The plot includes the additive effect of temperature at 50 m and spatial (latitude-longitude interaction) and temporal effect (year-month interaction) on hake larvae occurrence. Solid lines denote smoothed values, and shaded areas represent 95% CI.
( AIC = 11.73), and temperature ( AIC = 3.25). Temperature and SSH had a significant influence on the abundance of hake larvae with a dome-shaped relationship (Figure 4). Hence, when present, the highest hake larval abundances were present at a temperature of 13.36 • C and an SSH of -0.48 m. Moreover, temporal (year) and spatial variables were highly associated with abundances, with higher abundances in 2012 than in 2017, in the areas between 8.5 • W and 7 • W of longitude and around 44 • N of latitude (northeast section of the study area), and low hake larvae abundance southward. As with the occurrence model, we generated the prediction of the mean abundance when the hake larvae were present using the best abundance model output. The r 2 of the abundance model was 0.42 and the LOOCV test presented a mean r 2 of 0.26 (Supplementary Figure 6), indicating that parameter estimation was sensitive to data selection. Lastly, the estimated final prediction based on P/A and abundance (full hurdle model) predictions presented a total mean larval abundance of 14.

DISCUSSION
Previous studies on the European hake larval ecology focused either on the northern stock (Alvarez et al., 2004) or descriptive analysis of the larval fish community in the area (Rodriguez et al., 2015b). This is the first effort to examine the environmental drivers defining the spatio-temporal variability of hake larvae distribution and the abundance of the southern stock. Surveys were conducted in only two different years; so the study is both spatial and temporally limited. However, results contribute to the understanding of the larval dynamics of a complex but important species in a highly dynamic environment. Although the sampling in this study was restricted to the shelf off the NW coast of the Iberian Peninsula, this has been reported as the main spawning area of the stock and one of the most important nursery areas (Izquierdo et al., 2021), and thus a key for understanding the recruitment processes.
Our results confirm that this is an important spawning area of the European hake with larvae widely distributed over the study area but displaying a clear spatial structure, i.e., a decreasing abundance gradient from the northeast to the south-western areas. Hake larvae were mainly found offshore, on the shelf between the 100 and 200 m isobaths, close to Estaca de Bares.  The European hake larvae concentrations were also located close to the shelf edge in the NW Mediterranean (Olivar et al., 2003), the Bay of Biscay (Alvarez et al., 2004), or the Celtic Sea (Kacher and Amara, 2005). Mean larval abundance in this study (29.10 larvae/10 m 2 ) was higher than in the NW Mediterranean (4.70 larvae/10 m 2 ) June values (Olivar et al., 2003) and similar to the March values of the Northern stock (with 35 and 25 larvae/10 m 2 in 1995 and 1998, respectively). Conversely, these abundance values are lower than those reported for the Northern stock in April (1995: 45 larvae/10 m 2 ; 1998: 33 larvae/10 m 2 ) (Alvarez et al., 2004). Even so, it is necessary to remark that the larval abundance values indicate larval density, i.e., the concentration of larval patches. Thus, the data suggest that larvae from our study area form higher concentrations than the hake larvae in the Mediterranean and are similar to those in the northern stock, except for April 1998. We observed notable interannual differences in the overall larval occurrence and abundance among the surveys. The abundance in 2012 was almost two times that of 2017 (70.90 and 38.70 larvae/m 2 respectively). This year's effect may be associated with a potential gear effect since two different gears were used that may yield differences in catchability (Johnson and Fogarty, 2013). This is a complex issue as catchability depends on a number of factors (clogging, extrusion, and avoidance) and usually is species-and stagespecific. The differences in many of these aspects are minimized by monitoring the volume filtered. There were no differences in the volume filtered among the years. In a wide comparison of gears performance, the mesh size of the net and towing speed were the major factors affecting catchability (Skjoldal et al., 2013). While gear differences may have played a role in the interannual differences in the abundance observed in this study, we believe the differences in environmental conditions and stock sizes between surveys also contributed to our observations and likely exceeded the differences due to mesh size. In our opinion, the gear differences have also an impact on the observed differences in the abundance between the years. However, environmental and population factors may also explain these differences.

Environmental Variables
Fish larval habitat selection, distribution, and abundance are strongly influenced by environmental conditions (Heath, 1992;Govoni, 2005;Houde, 2009). In concordance, our results showed an association between oceanographic parameters (mainly temperature at 50 m and SSH) and the distribution and abundance of European hake larvae on the Galician shelf. Temperature appears to affect the habitat selection for this life stage. It is widely demonstrated that temperature and salinity are some of the main drivers of distribution and abundance of marine fish larvae, regulating spawning areas and timing and hence, contributing to the variability in the ELH stage dynamics (Houde, 2009(Houde, , 2016Miller and Kendall, 2009). Hake is not an exception and several studies have shown temperature as a key factor for hake larval survival, through its influence on larval growth and mortality (Palomera et al., 2005;Bartolino et al., 2008).
In the Galician shelf region, the occurrence of hake larvae has a negative relationship with the increasing water temperature. A negative relationship was also observed in the Mediterranean stock between peaks in water temperature during summer and the abundance of hake recruits in autumn, possibly triggered by a reduction in offspring survival rates (Bartolino et al., 2008). However, Goikoetxea and Irigoien (2013) observed a positive relation in the northern stock between the European hake recruitment and an increase in temperature anomalies. Conversely, our study shows that larval abundance has a relationship with temperature in a dome-shaped fashion, with higher abundances occurring between 13 and 13.5 • C. In the northern stock, the optimal temperature for hake spawning was expectedly lower, between 10 and 12.5 • C in the Bay of Biscay (Alvarez et al., 2001(Alvarez et al., , 2004, and between 10.5 and 12 • C in the British Isles (Coombs and Mitchell, 1982). Interestingly, the lowest hake egg mortality under laboratory conditions for the southern stock occurred at a slightly lower temperature, 12.7 • C ( Guevara-Fletcher et al., 2016).
The higher salinity observed in March 2012 may be the consequence of unusually low river run-off. In 2017, precipitation (and subsequently river run-off) in both sampled months were close to the mean of the time series (calculated since 1981); however, in March 2012, precipitation was 59% less than the accumulated rainfall compared to the mean of the time series. The convergence of oceanic waters and freshwater discharge is a favorable habitat for ELH stages facilitating larval growth and survival (Grimes and Finucane, 1991). It produces strong vertical temperature gradients (Koutsikopoulos and Le Cann, 1996) with a distinct thermocline in which larvae and their food can concentrate due to lower dispersion and higher prey concentration, as reported in other hake species (Bailey, 1980;Álvarez-Colombo et al., 2011).
The circulation pattern also drives larval transport and retention mechanisms toward optimum habitat, influencing offspring distribution (Houde, 2016). It also determines the mixing level of different water masses, i.e., it indirectly affects temperature and salinity, and hence larval survival. A significant and dome-shaped relationship was observed between the European hake larvae abundance and SSH, with the highest abundances of hake detected at SSH between −0.50 and −0.45 m. The SSH is used as an indicator of upper ocean circulation, and the negative values of SSH correspond to the presence of cyclonic eddies (Häkkinen, 2001). Cyclonic and anticyclonic eddies are common throughout the study area (Ríos et al., 1992). Cyclonic eddies are characterized by high productivity due to the upwelling of cool and nutrient-rich waters into the euphotic zone that can benefit fish larvae through eventual increased food availability (Rooker et al., 2012). As eggs and early larvae are passive life-stages, the currents are the main factor determining the retention in the favorable habitat or offshore drift (Alvarez et al., 2004). Offspring needs to be transported toward suitable nursery grounds, often located in the coastal areas, where the conditions for survival (such as food availability, refuge structures, settlement habitat, and predator abundances) are better (Bradbury et al., 2008). In the case of the European hake, spawning takes place offshore near the continental slope, but juvenile recruitment occurs inshore, over the continental shelf. It indicates the importance of the retention and transport of eggs and larvae from the spawning to the nursery grounds (Alvarez et al., 2004), such as the Artabrian Gulf (off La Coruña), which is a persistent hake recruit aggregation area, favored by the mesoscale eddies produced by the Iberian Poleward Current (Sánchez and Gil, 2000;Izquierdo et al., 2021). These critical processes for larval hake directly impact recruitment: larval dispersion toward the open ocean produces a negative effect (high larval mortality) while the transport to and retention in the nursery grounds produce a positive effect in the annual recruitment, as was demonstrated also for Merluccius productus (Hollowed and Bailey, 1989).

Population Factors
Marine spawning strategies have been shown to be coupled to specific oceanographic conditions that benefit larval retention, survival, and recruitment (Muhling et al., 2013). Additionally, species that inhabit upwelling areas, like the Galician shelf, are adapted to a variable environment and its potentially adverse effects, developing several opportunistic tactics to cope with these conditions, including changes in behavior, variations in spawning timing and duration or existence of multiple annual spawning events (Grote et al., 2007). The spawning fraction, i.e., the percentage of females spawning in a given period of time, was estimated from commercial samplings for the same months of the three different surveys following the methodology of Murua et al. (2003). The spawning fraction for March 2012, March 2017, and June 2017 were 0.27, 0.21, and 0.01, respectively (unpublished results). This spawning activity seems to be in agreement with observed differences in larval abundance between years, with higher spawning activity in March 2012 and consequently, higher larval production.
Inter-annual differences in hake productivity on the Galician shelf were also observed in the autumn bottom-trawl survey indexes [DEMERSALES survey, International Council for the Exploration of the Sea [ICES], 2019]. A clear decrease of recruit biomass and recruitment abundance (<20 cm) between 2012 and 2017 supported our results, with 22% of less biomass (from 8.44 to 6.58 kg/30 min) and 50% of the recruits (from 280 to 140 recruits/30 min). Lower biomass and lower spawning activity produced fewer larvae and hence, fewer recruits.
However, the European hake shows a highly plastic reproductive behavior with a very protracted spawning season that has several spawning peaks that vary among the years in both timing and intensity (Serrat et al., 2019;García-Fernández et al., 2020). The low larval abundance in 2017 may be the result of a timing mismatch between surveys and spawning, or a gear effect as discussed above. It is also worth noting that the predictive capabilities of the models were sensitive to survey, and were particularly poor for June. This may indicate different spawning conditions for the second spawning peak that we were unable to effectively elucidate. In addition, annual differences may also emerge from other potential causes, such as competition or predation where density acts as a controlling factor causing size-dependent mortality or size-dependent growth (Bailey and Houde, 1989;Stige et al., 2019). During the March 2017 survey, an unusually high number of salps (Salpa sp.) were collected in the bongo net. These organisms exhibit high feeding activity and are capable of capturing particles in an extensive size range with high efficiency (Mullin, 1983), potentially affecting the larval survival due to predation.
Due to the absence of previous studies on the ELH hake abundance in the study area, it is not possible to estimate spawning success in this study. However, combining the results of this study with adult population features (García-Fernández et al., 2020), it is likely that spawning success was higher in 2012 compared to 2017.
In this study, we assumed that the location and timing of larvae collections in the Galician shelf can be used to infer spawning in space and time because pelagic eggs and small larvae for most marine species have limited swimming abilities (Houde, 2009). Information about the distribution and abundance of larvae can be combined with environmental data to determine the timing and location of spawning (Govoni, 2005;Suca et al., 2018) as well as to determine when and where suitable nursery habitats occur (Kitchens and Rooker, 2014). Using the larval abundance estimates from early life surveys, we can also estimate other population dynamic parameters, such as potential recruitment (Ingram et al., 2010).
However, the contribution of each hake spawning peak to total annual recruitment is still poorly understood. The observed annual differences in the timing and intensity of each spawning peak (García-Fernández et al., 2020) had even led to the hypothesis that there are different spawning stock components and not a single spawning stock (Serrat et al., 2019). The existence of maternal effects (comprised of the Stock Reproductive Potential concept) in egg production and quality, as well as in the larval performance (García-Fernández et al., 2020) produce an even more complex scenario, especially due to the interaction between maternal effects and environmental conditions during the spawning. Understanding this interaction is essential to know the fate of the offspring and estimate yearclass strength and recruitment. The implications of this could be substantial when considering the global increase in ocean temperatures (Pecl et al., 2017). Specifically, along the NW Iberian Peninsula, there has been a significant increasing trend in sea surface temperatures with a mean rate of 0.24 • C per decade in recent decades (Gómez-Gesteira et al., 2011). Thus, with climate change and warming waters, the location of the habitat for hake spawning and larval development would potentially shift. A multitude of other climate-driven impacts, including mismatches with zooplankton prey, would also be possible (Llopiz et al., 2014), which could jeopardize the production capacity of the population in the long term.

Concluding Remarks
Our results suggest that reproductive phenology and environmental conditions interact to determine the occurrence and abundance of hake larvae. In this study, we have identified that temperature and SSH play important roles in the spatiotemporal distribution of the European hake larvae in the main spawning area of the southern stock. Moreover, we provide valuable information about the larval habitat of a commercially important species and shed light on factors related to larval survival, a critical process for understanding the impact of annual recruitment success on population productivity and sustainability for harvested and diminished stocks (Chambers and Trippel, 1997). Lastly, it is necessary to perform new studies of European hake in the southern stock to increase the information and improve the prediction of models, extending the years and sampling periods. Additionally, the integration of environmental conditions and maternal effects on larval recruitment should be addressed to understand the complex dynamics of the European hake population.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation. Requests to access these datasets should be directed to CG-F, cgarcia@iim.csic.es.

AUTHOR CONTRIBUTIONS
CG-F compiled the database. CG-F, JJS, and JKL designed the statistical analysis. CG-F and JJS performed the statistical analysis, RD-P validated the statistical analysis. FS-R and PÁ made the conceptualization, experimental design, and survey methodology. FS-R administrated and coordinated the project and the funding acquisition. CG-F wrote the original draft of the manuscript. JJS, RD-P, and FS-R helped with the manuscript writing. All authors contributed to manuscript revision, read, and approved the submitted version.

FUNDING
This study was carried out with financial support from the CRAMER (CTM2010-21856-218567-CO3-02) and DREAMER projects (CTM2015-66676-C2-1-R) funded by the Spanish National Research Program. CG-F was supported by a Predoctoral Fellowship from the Fundación Tatiana Pérez de Guzmán el Bueno, as well as the support by Do * Mar (Ph.D. program). JJS was supported by the National Science Foundation Graduate Research Fellowship. JKL was supported by the Academic Programs Office at Woods Hole Oceanographic Institution.