Spatio-Temporal Assessment of the European Hake (Merluccius merluccius) Recruits in the Northern Iberian Peninsula

ORIGINAL RESEARCH published: 28 January 2021 doi: 10.3389/fmars.2021.614675 Frontiers in Marine Science | www.frontiersin.org 1January 2021 | Volume 8 | Article 614675 Edited by: Antonio Di Franco, Stazione Zoologica Anton Dohrn, Italy Reviewed by: Germana Garofalo, Institute for Biological Resources and Marine Biotechnology (CNR), Italy Nerea Goikoetxea, AZTI Foundation, Spain *Correspondence: Francisco Izquierdo [email protected] orcid.org/0000-0002-0781-1354 Specialty section: This article was submitted to Marine Conservation and Sustainability, a section of the journal Frontiers in Marine Science Received: 06 October 2020 Accepted: 04 January 2021 Published: 28 January 2021 Citation: Izquierdo F, Paradinas I, Cerviño S, Conesa D, Alonso-Fernández A, Velasco F, Preciado I, Punzón A, Saborido-Rey F and Pennino MG (2021) Spatio-Temporal Assessment of the European Hake (Merluccius merluccius) Recruits in the Northern Iberian Peninsula. Front. Mar. Sci. 8:614675. doi: 10.3389/fmars.2021.614675 Spatio-Temporal Assessment of the European Hake (Merluccius merluccius) Recruits in the Northern Iberian Peninsula Francisco Izquierdo1,2*, Iosu Paradinas 3, Santiago Cerviño2, David Conesa1, Alexandre Alonso-Fernández4, Francisco Velasco5, Izaskun Preciado5, Antonio Punzón5, Fran Saborido-Rey4and Maria Grazia Pennino2,6 1Departament d’Estadıstica i Investigació Operativa, Universitat de Valéncia, Burjassot, Spain, 2Instituto Español de Oceanografía, Centro Oceanográfico de Vigo, Vigo, Spain, 3Asociación Ipar Perspective, Sopela, Spain, 4Instituto de Investigaciones Marinas (IIM-CSIC), Vigo, Spain, 5Centro Oceanográfico de Santander, Instituto Español de Oceanografía, Promontorio de San Martín, Santander, Spain, 6Fishing Ecology Management and Economics, Departamento de Ecología, Universidade Federal do Rio Grande do Norte, Natal, Brazil The protection of fish nurseries has been recognized as a useful tool to efficiently manage fisheries given that protected areas enhance the recruitment of target species. To identify and locate potential nursery areas, a solid understanding of species-environment relationships and their spatio-temporal dynamics is needed. Within this context, in this study we assess where European hake (Merluccius merluccius) recruits persistently aggregate in the northern continental shelf of the Iberian Peninsula. Hake recruit data collected during scientific trawl surveys between 2005 and 2016 were analyzed using Bayesian hurdle hierarchical spatio-temporal models, considering the environmental variables bathymetry, sea bottom temperature and salinity. Additionally, three different spatio-temporal structures (i.e., persistent, progressive, or opportunistic) were compared to assess the temporal persistence of nurseries over time. Among all the environmental variables analyzed, bathymetry was the most important. The preferential habitat of recruits was found to be within a bathymetric range of 120–200 m. Our findings clearly show that there is a temporally persistent main nursery located along the continental shelf of the Artabrian gulf (off La Coruña) in addition to several areas with high aggregations of hake recruits but with strong inter-annual variability. We argue that the analytical framework applied in this study allowed us to identify European hake nurseries in the northern continental shelf of the Iberian Peninsula, as well as their spatio-temporal fluctuations throughout the study period (2005–2016), and to assess which environmental factors, among bathymetry, sea bottom temperature and salinity, influence the occurrence and abundance of recruits in the study area. Results of our models also produce a new abundance index that could be useful for improving traditional stock assessment models. Keywords: bayesian models, European hake, hurdle-model, recruits, stock assessment, spatial ecology, INLA approach Izquierdo et al. European Hake Recruit’s Aggregations 1. INTRODUCTION Contemporary sustainable fisheries management advisory panels have acknowledged the importance of fisheries restricted areas (FRA) as a complementary tool to conventional stock management policies (FAO, 2016). Increased overfishing can lead to fluctuations in species stock abundance, and therefore, fully accounting for the age-structure and natural growth potential of fish stocks is crucial, being recruitment one of the most vulnerable phases (Diekert et al., 2010). There are some measures to protect this life stage, such as increasing mesh size or using a different net geometry (Carlucci et al., 2009), but other technical measures that involve protecting spatial or spatio-temporal FRAs may be key to accomplishing sustainable management (Colloca et al., 2015; Paradinas et al., 2015) as they are less susceptible to annual assessment strategies or political agreements (Pitchford et al., 2007). European hake (Merluccius merluccius) is one of the most important commercial species in the northeast Atlantic and is assessed by the International Council for the Exploration of the Sea (ICES) in two units: the northern and the southern stocks. Despite the recovery plan for the southern stock implemented in 2006 (EC 2166/2005) and the Multiannual Management Plan for Western Waters (EU, 2019), the fish mortality rate is still above the maximum sustainable yield (Fmsy), even though the spawning stock biomass is above precautionary levels (ICES, 2018a). Diverse spatial studies have identified nursery areas for this species and stock (Casey and Pereiro, 1995). Sánchez and Gil (2000) found persistent nursery areas in the northwest coast of the Iberian Peninsula, while a more recent study by Pennino et al. (2019) identified four areas with different levels of overlap with the past study. This difference in results demonstrates the importance of dynamic spatio-temporal modeling, given that the effect of environmental conditions over time can vary the distribution, shape and abundance of hake nursery grounds. Several studies have assessed whether different physical, environmental and biological factors can affect growth and survival rates of European hake during the early life stages (Alvarez et al., 2001, 2004; Bartolino et al., 2008; Korta et al., 2015). In the northeastern Atlantic, previous studies have determined that recruitment success has a direct relationship with regime shifts expressed through the sea surface temperature and the Ekman’s transport index (Goikoetxea and Irigoien, 2013). Along the continental Galician shelf, random variations in the physical environment, such as North Atlantic Oscillation (NAO) and the upwelling, have been observed to directly affect the growth and size of mature hake females (Domínguez-Petit et al., 2008). Moreover, the interplay of westward larvae transport within anticyclonic eddies that collide with the shelf grounds favors high recruitment to nursery areas at depths of between 90 and 180 meters in the western shelf of the Cantabrian Sea (Sánchez and Gil, 2000). Valuable information regarding ecological hake dynamics is being addressed both through various studies and the stock assessment working group developed by ICES, which relies on information from commercial landings and oceanographic surveys provided by research centers, such as the “Instituto Español de Oceanografía” (IEO) and the “Instituto Português do Mar e da Atmosfera” (IPMA). Since 2014, this data has been used as an input to the GADGET (Globally applicable Area Disaggregated General Ecosystem Toolbox) model (length based), which is currently one of the most complete methods to study the state of a population over time (ICES, 2018b). However, most of these stock assessment studies maintain the same assumption that spatial and environmental variability is implicitly intrinsic in the data, and therefore, it is not explicitly taken into consideration in assessment models. Survey abundance indices are used to calibrate dynamic models assuming some proportionality between survey abundance and modeled abundance. To guarantee that this proportionality is constant on time the sampling design is kept the same every year, allowing a random selection of hauls in each strata. However, changes on fish distribution can happen from year to year affecting the design based abundance estimation. If these abundance indices ignore spatial structures in estimates it can lead to biased results (Shelton et al., 2014; Cao et al., 2017). Studying the spatial distribution of a species like hake over time is a medium-large scale process, which is influenced by changing environmental features. Geostatistical models take into account spatial autocorrelation, which allows spatial dependency between neighboring locations to be incorporated, following the principle that nearby locations are more similar than distant ones (Tobler, 1970; Dormann, 2007). Similarly, when analyzing time series data consecutive measurements of months and/or years are more similar than separated ones, and therefore, display temporal autocorrelation (Tsitsika et al., 2007). Taking explicitly spatial and temporal correlation in models into account is essential since observations are often collected close in time and space and are subject to similar environmental characteristics (Hurlbert, 1984; Hefley et al., 2017). Ignoring spatial and temporal dependence in this type of analysis can generate misleading results (Latimer et al., 2006; Agostini et al., 2008; Fahrmeir and Kneib, 2008). Recently, the use of Bayesian hierarchical spatio-temporal models has demonstrated to be a convenient approach when modeling species habitat preferred conditions (Muñoz et al., 2013) and also to identify species nurseries (Colloca et al., 2009; Paradinas et al., 2015, 2020), given their ability to incorporate the spatio-temporal correlation as well as the associated uncertainties. Moreover, to be useful for fisheries management, the identified nurseries need to be areas that persist over time (Fiorentino et al., 2003; Garofalo et al., 2011; Colloca et al., 2015; Druon et al., 2015). A variety of research approaches have been proposed to assess temporal persistence: Colloca et al. (2009) proposed Bayesian kriging maps by year associated to aggregation curves that identify the level of persistence; Paradinas et al. (2015) proposed a Bayesian approach that compared two spatio-temporal structures, one fitting persistent scenarios (i.e., spatial distributions that do not change with time), and the other accommodating opportunistic scenarios (i.e., spatial distributions that change unrelatedly across subsequent time events). Recently, Paradinas et al. (2017, 2020) enhanced the latter approach by including a third spatiotemporal structure to fit progressive processes (i.e., distributions that are neither totally persistent nor opportunistic but change progressively over time). Frontiers in Marine Science | www.frontiersin.org 2January 2021 | Volume 8 | Article 614675 Izquierdo et al. European Hake Recruit’s Aggregations FIGURE 1 | Sampling locations (black dots) of SP-NSGFS Q4 surveys (2005–2016). Bathymetric lines indicate the 100 and 800 m isobaths. Closure areas to trawl in Galician and Cantabrian Sea Waters. Colored polygons are permanent closed while White areas correspond to the temporal ones (BOE, 2015). Within this context, here we apply a Bayesian hierarchical spatio-temporal hurdle-model (B-HSTHMs) to (1) identify European hake nurseries in the northern continental shelf of the Iberian Peninsula and their spatio-temporal fluctuations over time (2005–2016); (2) assess which environmental factors, among bathymetry, sea bottom temperature and salinity, influence recruit occurrence and abundance in the study area, and (3) produce a new habitat-based abundance index of recruits for traditional stock assessment models. 2. METHODS 2.1. Study Region The region of interest for this study is the northern continental shelf of the Iberian Peninsula (see Figure 1), a narrow area (10– 60 km) of almost 18,000 km2that is characterized by important and marked hydro dynamism (Abad et al., 2019). Over the shelf, currents are driven by regional factors, such as tides and wind (Gil and Sánchez, 2003). In the winter, a warm and saline poleward current moves eastward along the Cantabrian coast and enters the Bay of Biscay (Korta et al., 2015). In addition, the coastal upwelling off the Galician and Portuguese coasts appears during spring and summer which, combined with hydrographic mesoscale activities, has a strong influence on the primary production of the area (Sánchez and Olaso, 2004). In this area, hake is caught by trawls, gillnetters, long-liners, and artisanal fleets (ICES, 2018a). The bottom-trawling fleet operate throughout the year exploiting the fishing grounds, but seasonal and permanent spatial closures for this fleet have been implemented since before 2000 (BOE, 1996, 2015). Specifically, the fishing grounds with seasonal closures for hake recruitment are La Coruña/Cedeira, where trawling is forbidden from October to January, and La Carretera, where trawling is forbidden from September to March (Figure 1). 2.2. Data Collection Data on hake recruits were collected by the IEO for the scientific survey series SP NSGFS Q4, which is carried out annually in autumn (September to October), and the study period includes data from 2005 to 2016. This survey used a stratified sampling design based on depth with three bathymetric strata: 70–120, 121–200, and 201–500 m. Moreover, some samples from additional trawls at greater depths were also included. Sampling stations consisted of 30 min trawling hauls located randomly within each strata at the beginning of the survey design (Figure 1). Approximately 115 hauls divided between the three bathymetric strata were performed every year in this zone, using baka 44/60 gear (ICES, 2017) and following the International Bottom Trawl Survey Working Group protocols (IBTSWG) of ICES. All individuals smaller than 21 centimeters were considered recruits (Murua, 2010). It should be noted that the fishing gear used is well-designed to catch recruit individuals (cod-end mesh size 10 mm), who settle at the bottom with about 10 cm (Sánchez and Gil, 2000). Two different response variables were analyzed to characterize the spatio-temporal behavior of European hake recruits. First, we considered a presence/absence variable to measure the occurrence probability of the species. Second, we used the total number of individuals standardized by 30 min of trawling (i.e., number per unit effort, NPUE) as an indicator of the conditionalto-presence abundance of the species. Frontiers in Marine Science | www.frontiersin.org 3January 2021 | Volume 8 | Article 614675 Izquierdo et al. European Hake Recruit’s Aggregations TABLE 1 | Summary of fitted spatio-temporal models U(s,t). Model Notation Description Opportunistic U(s,t)=WtDifferent and uncorrelated realizations of the spatial field every year. Persistent U(s,t)=W+f(t) A common realization of the spatial field for all years and an additive temporal trend f(t). Progressive U(s,t)=Wt+ρUst−1Spatial realizations change over time through a first order autoregressive model. ρcontrols the level of correlation between subsequent time events. Wrepresents a geostatistical spatial field, f(t)is a temporal trend function and ρis an autoregressive correlation parameter bounded to [0,1]. Due to the demersal nature of this species, only three environmental variables were considered to be potential drivers. Two oceanographic variables, Sea Bottom Temperature (SBT in ◦C) and Sea Bottom Salinity (SBS in PSU), and bathymetry (in meters). SBT and SBS were added to the analysis as they are strongly related to marine system productivity given that they affect nutrient availability and water stratification (Muñoz et al., 2013). SBT and SBS values were collected during the survey with a sounding CTD SeaBird (conductivity, temperature and depth) in different random sampling points of the study area. Average SBT (Supplementary Figure 1) and SBS (Supplementary Figure 2) maps for the months of the survey (September–October) over the entire area were obtained for each year of the period studied with the Radial basis functions (RBF) tool in ArcGIS 10.1. Bathymetry (Supplementary Figure 3) was obtained from the European Marine Observation and Data Network with a spatial resolution of 0.02 ×0.02 decimal degrees. In order to ensure the same spatial resolution, all environmental data were aggregated to the lower spatial resolutions (0.05 ×0.05 decimal degrees) using the raster package (Hijmans, 2018) in the R software (R Core Team, 2017). All covariates were explored for collinearity, outliers, and missing values following the approach of Zuur et al. (2010) before being included in the models. Correlation among variables was tested using the Spearman’s correlation (Supplementary Figure 4), whereas collinearity was computed using Generalized variance-inflation factors (GVIF) (Fox and Weisberg, 2011) (Supplementary Figure 5). To avoid numerical confounding within the model selection, the explanatory variables were standardized (difference from the mean divided by the corresponding standard deviation) (Gelman, 2008). 2.3. Characterizing the Spatio-Temporal Behavior of European Hake Recruits This study used the spatio-temporal geostatistical model structure comparison proposed by Paradinas et al. (2017) to categorize the spatio-temporal behavior of European hake as either opportunistic, persistent or progressive (see Table 1 and Figure 2). These models use a geostatistical term that accounts for spatial autocorrelation in the data and can easily include different spatiotemporal extensions. In particular, opportunistic structures indicate that species change their spatial patterns every year without following any specific pattern, persistent structures imply that species have a spatial distribution that does not change every year, and progressive structures indicate that species change their spatial patterns in a correlated way from 1 year to another. The progressive structure contains an autoregressive ρtparameter (see Table 1) that controls the persistence level of the spatial effect. This ρtparameter is bound to [0, 1], where parameter values close to 0 represent more opportunistic behaviors and parameter values close to 1 represent more persistent distributions over time. The spatial field (Ws) (see Table 1) was modeled as a normal multivariate distribution with a zero mean and a Matérn covariance function that depends on its range (rw) and variance (σw). The temporal trend f(t) could follow any suitable function, for example, a linear effect, a smooth effect, an unstructured random term, etc. 2.4. Modeling European Hake Occurrence and Abundance Distribution Spatio-temporal fisheries abundance data often involves a large proportion of zeros in observations, i.e., zero inflated data. These data are generally tackled using independent two-part models, also known as hurdle models. In these models, the occurrence and conditional-to-presence abundances (NPUE) are modeled independently. However, the variables abundance and detection probability are often related (Kéry et al., 2005), which violates the independence assumption in classic delta models. This study acknowledges that both processes could be related by fitting shared environmental and/or spatio-temporal effects through joint modeling techniques (Hogan and Laird, 1997; Held et al., 2005). We have a semi continuous dataset regarding the presence/absence and abundance data, as both processes are measured in the [0, ∞) interval (Paradinas et al., 2017). Moreover, the proportion of zeros in the data gets close to 30% for some of the years. Within this context, the application of a hurdle model enables that the first process determines the occurrence and the second the intensity when the response is non-zero, what makes the modeling process more informative. Specifically, Y(s,t) and Z(s,t) denote the spatio-temporal distribution occurrence and the conditional-to-presence abundance (NPUE), respectively, where s=1, ... ...,nt is the spatial location and t=1, ...., Tthe temporal index, with i=1, ...,Irepresenting the environmental variable in location s. As is generally done with these kinds of variables, we used a Bernoulli distribution to model the occurrence of Y(s,t). With respect to NPUE Z(s,t), we chose to model it with a gamma distribution, which is a probability distribution that helps to capture any possible overdispersion in the data. The mean of both variables was then related via the usual link functions (logit and log, respectively) to the environmental and Frontiers in Marine Science | www.frontiersin.org 4January 2021 | Volume 8 | Article 614675 Izquierdo et al. European Hake Recruit’s Aggregations FIGURE 2 | Simulated spatio-temporal scenarios extracted from Paradinas et al. (2017). An opportunistic distribution shows big unrelated changes in the spatial structure among neighboring times. A persistent scenario has very similar spatial structures at every time. A progressive scenario has correlated spatial realization changes over time. (A) Opportunistic distribution. (B) Persistent distribution with trend. (C) Progressive spatio-temporal distribution. spatio-temporal effects: Y(s,t)∼Ber(π(s,t)) Z(s,t)∼Gamma(µ(s,t), φ) logit(π(s,t)) =αY+ I X i=1 fi(Xi(s,t)) +UY(s,t) log(µ(s,t)) =αZ+ I X i=1 θifi(Xi(s,t)) +θUUZ(s,t) (1) where π(s,t) represents the probability of occurrence at location sat time t, and µ(s,t) and φrepresent the mean and dispersion of the conditional-to-presence abundance. The linear predictors, which indicate how parameters π(s,t) and µ(s,t) behave in space and time, are formed by: αYand αZ, i.e., the terms representing the intercepts of each variable; f() represents any function applied to our environmental variables Xi, which in our case were either a constant function (fixed effects) or a second order random walk (RW2) function to fit non-linear relationships (Fahrmeir and Lang, 2001); and the final terms, UY(s,t) and UZ(s,t) refer to the spatiotemporal structure of occurrence and the conditional-topresence abundance, respectively, and may follow any of the three spatio-temporal structures described in the previous section. The scale parameters θiand θUpermit differences in scale across both linear predictors in shared effects. It is important to note that sharing occurrence and conditional-to-presence abundance information to fit common environmental effects may not always improve the models, thus we tested all possible model combinations. Vague prior distributions with a zero-mean and a standard deviation of 100 were implemented for all the fixed effects, the abundance process variance, and the scaling parameter (θ) of the shared effects. PC priors (Fuglstad et al., 2018) were used to describe prior knowledge of hyperparameters of the geostatistical terms and the ρparameters of the RW2 functions. These priors were set as follows: if the prior probability of the spatial range was smaller than 0.5, it was set at 0.05, if the probability of the spatial variance was larger than 0.6, it was also set at 0.05, and if the probability that the precision of the RW2 effects was larger than 0.5, it was set at 0.01. A sensitivity analysis of the choice of priors was performed by verifying that the posterior distributions Frontiers in Marine Science | www.frontiersin.org 5January 2021 | Volume 8 | Article 614675 Izquierdo et al. European Hake Recruit’s Aggregations TABLE 2 | Spatio-temporal structures comparison for the conditional-to-presence abundance distribution European hake recruits’ model based on WAIC and LCPO scores. Model WAIC LCPO Time (s) Persistent shared effects 15879.45 2.90 80.91 Persistent not shared effects 16001.28 2.92 118.08 Opportunistic shared effects 16095.17 2.95 59.82 Opportunistic not shared effects 16231.99 2.95 79.56 Progressive shared effects 16774.70 3.05 401.62 Progressive not shared effects 15846.09 3.11 7138.10 Time scores refer only to the estimation process of the model. The best model is highlighted in bold. concentrated well within the support of the priors (Zuur et al., 2017). Model selection was performed by testing all possible combinations among the possible spatio-temporal structures and variables using the Watanabe Akaike Information Criterion (WAIC) (Watanabe, 2010) for goodness of fit and the mean logarithm of the Conditional Predictive Ordinate (LCPO) (Roos et al., 2011) as a predictive quality measure. For both measures, the smaller the score, the better the model. All models were fitted using the integrated nested Laplace approximation (INLA) (Rue et al., 2009) approach in the R-INLA software. The R-INLA code for the applied models can be found on this GitHub repository. From the final spatio-temporal model output we obtained the mean yearly habitat-based recruits abundance index. For comparison purposes, we created a stock assessment strata-based index for the number of hake recruits (<21 cm), standardized for the area of their corresponding bathymetric strata (from 70 to 500 m) (ICES, 2017). Finally, we defined hake nurseries as high recruit aggregation areas which are persistent over the time. Thus, for the average of a given n-years period we are able to identify potential nursery areas with the contour polygon based on the NPUE above 90th percentile, as described in Paradinas et al. (2020). 3. RESULTS 3.1. Spatio-Temporal Structure Selection With respect to the spatio-temporal structures, the shared components showed better goodness of fit in the opportunistic and persistent models (Table 2). However, the best structure was the progressive model with non-shared spatio-temporal effects. This means that the spatial occurrence and abundance patterns of recruits change in a correlated way from 1 year to the next, but that both sub-processes behave considerably different from one another. This structure was selected for the final model because it had the lowest WAIC value. It should be mentioned that the progressive structure had the largest computational time, with a notable difference (Table 2). 3.2. Variable Selection Concerning the environmental variables, bathymetry was the most important predictor to defining the distribution occurrence TABLE 3 | Environmental effects comparison for the conditional-to-presence abundance distribution European hake recruits’ model based on WAIC and LCPO scores. Model PNSE WAIC LCPO Time (s) B shared effects 15659.88 3.02 13667.78 SBS shared effects 15848.98 3.11 7168.39 SBT shared effects 15800.53 3.15 11032.17 B + SBS shared effects 15655.22 3.05 16488.46 B + SBT shared effects 15657.85 3.07 17097.45 SBS + SBT shared effects 15804.95 3.16 11683.53 B not shared effects 15668.76 3.03 10143.00 SBS not shared effects 15852.73 3.11 10662.15 SBT not shared effects 15798.90 3.14 9416.98 B + SBS not shared effects 15672.92 3.03 14104.07 B + SBT not shared effects 15672.60 3.06 15135.95 SBS + SBT not shared effects 15805.43 3.14 11152.92 B, Bathymetry; SBS, Sea bottom salinity; and SBT, Sea bottom temperature. Note all models include the selected Progressive Not Shared Effects spatio-temporal structure (PNSE). Time scores refer only to the estimation process of the model. The best model is highlighted in bold. and abundance of hake recruits in the study area (Table 3). Although, in terms of WAIC, the best models were those that included bathymetry and either SBS or SBT, the difference between these models and the model that only included bathymetry was rather negligible (i.e., lower than 5 units). This means that the variability of the hake recruit abundance is mostly affected by bathymetry and spatio-temporal effects. Therefore, and following a parsimony principle, the selected model included only a shared bathymetric effect as a covariate. 3.3. Final Model Given that the only relevant variable included in the final model was bathymetry, original scale units were used to facilitate interpretation of the partial effect. There is an optimum bathymetric range between 140 and 200 m for both the occurrence and abundance of hake recruits (see Figure 3). This means that hake recruit abundance decreases gradually with bathymetry, as presented in Figure 3. The selection of an autoregressive temporal term in the model suggests the presence of a certain degree of temporal persistence in the spatial distribution of hake recruits in the study area. These results were supported by the high temporal correlation parameters (ρ) of the progressive spatio-temporal structures (0.99 and 0.96 for the occurrence and conditional-to-presence abundances, respectively). With respect to the posterior distribution mean of the occurrence spatial effect, only slight changes were observed between years (Figure 4). This result highlights a clear stable structure for the presence of hake recruits in the study area that remains constant throughout the whole period. However, the mean spatial effect for abundance had high inter-annual variation, thereby reflecting changes to recruit concentration areas over the time (Figure 5). We can observe similar patterns in the standard deviation maps of the posterior distribution of Frontiers in Marine Science | www.frontiersin.org 6January 2021 | Volume 8 | Article 614675 Izquierdo et al. European Hake Recruit’s Aggregations FIGURE 3 | Smoothed bathymetric effect (“RW2”) for the sub-process Binomial (A) and Gamma (B) in the linear predictor scale (logarithmic link). FIGURE 4 | Spatio-temporal occurrence maps for Merluccius merluccius recruits in the northeast Atlantic showing mean spatial effects (log-scale) for the probability of presence. the spatial effect for both the occurrence and abundance subprocesses (see Supplementary Figures 6,7). The median of the posterior predictive distribution for the occurrence revealed that the areas with the presence of hake recruits have a constant spatial pattern over time (Figure 6). However, abundance did not show the same spatio-temporal behavior. Indeed, we found areas with a persistent high concentrations of recruits over time, as well as intermittent zones (Figure 7). Within the study area, the northwest area of the Iberian Peninsula presented the highest concentration of recruits Frontiers in Marine Science | www.frontiersin.org 7January 2021 | Volume 8 | Article 614675 Izquierdo et al. European Hake Recruit’s Aggregations FIGURE 5 | Spatio-temporal abundance (NPUE) maps for Merluccius merluccius recruits in the northeast Atlantic showing mean spatial effects (log-scale). throughout the entire period (Figure 7). Specifically, from west to east three persistent areas were identified: the area in front of Vigo, the Artabrian gulf (off La Coruña) and La Carretera (Figures 7,1). Three intermittent areas were identified: one in front of Santander and the other two corresponding to the Bermeo and Guetaria-Fuenterrabía spatial closure areas (Figures 7,1). The years with the greatest hake recruit abundance were 2005 and 2009, whereas the years with the lowest values were 2008 and 2014 (Supplementary Figure 8). Regarding the variability of the generated predictions, estimates corresponding to the first and third quartiles showed a spatial pattern very similar to the predicted posterior median, thereby suggesting a good degree of accuracy with respect to the uncertainty of the model for both the occurrence and abundance of recruits (Figure 8). The yearly maps that correspond to the three quartiles for occurrence and abundance can be found in Supplementary Figures 8–13. In addition, we identified the main aggregation areas with a NPUE above the 90th percentile to be potential hake nurseries, taking into account the average of all study years (2005–2016) (Figure 8). The overlapping area between the FRA of La Coruña and the main identified nursery area is about the 27.46% (800.250 km2) (see Supplementary Figure 14). The mean predicted posterior abundance trend of European hake recruits was similar to the trend of the designed stratabased index that would be commonly used as an input for stock assessment models. However, the output from our model showed lower mean abundance values within the entire study period, specifically for the years 2009 and 2015, where the strata based index values were high (Figure 9). 4. DISCUSSION In order to implement fisheries restricted areas (FRAs), a spatial characterization of the key life–cycle habitats of commercial stocks, such hake nursery or spawning grounds, is required. Here, we apply a methodological approach to identify nursery grounds based on the spatio–temporal persistence of recruit aggregation areas. Coinciding with Paradinas et al. (2017), the best spatiotemporal structure was the progressive one with non-shared effects. This result suggests that hake recruitment data is generated in two different steps : (1) identifying the probability of observing hake recruits, and if present, (2) identifying their abundance. However, the nature of the process under study leads us to believe that this apparent independence is a consequence of the survey sampling effort, rather than a result of two different processes. The SP-NSGFS Q4 survey trawls a relatively big area, and therefore the probability of observing at least one individual of an abundant fish species, such as hake, is quite high in environmentally not-too-challenging areas. However, if effort had been increased to deeper areas, the detection probability Frontiers in Marine Science | www.frontiersin.org 8January 2021 | Volume 8 | Article 614675 Izquierdo et al. European Hake Recruit’s Aggregations FIGURE 6 | Median occurrence probability maps of the posterior predictive distribution for European hake recruits within the study period. Bathymetric lines correspond to 200 and 800 m isobath. would have decreased proportionally, and thus, generated a lot more zeros in our dataset. Occurrence and abundance have been shown to respond equal and differently depending on the environmental variables, as it has been observed for different elasmobranch species in the Mediterranean, therefore, they are modeled by semi-independent processes (Lauria et al., 2015). This result is highly related to the life habits of a species and their age-size, since the individuals have a wide spatial distribution range where they are present, but they can be more abundant in areas where environmental preferred conditions take place (i.e., trophic niche). In our case of study, due to the benthic habits of European hake recruits, we consider that the occurrence, abundance and environmental variables are modeled by the same generating process. Among all the covariates tested, SBS and SBT were the least relevant at explaining the occurrence and abundance of hake recruits in the study area. Similar results were observed by Paradinas et al. (2015), who found that surface temperature and type of substrate were not relevant to recruiting hake in the Mediterranean. However, different results have been reported in literature. A study of regime shifts carried out by Goikoetxea and Irigoien (2013) found that temperature increase anomalies had a positive effect on hake recruitment from the northeast Atlantic. In the Ligurian Sea northern hake nursery area, Abella et al. (2008) found a negative correlation between the North Atlantic Oscillation (NAO) index and recruitment success. In the Mediterranean, high peaks in water temperatures during summer resulted in a decreased abundance of hake recruits in autumn (Bartolino et al., 2008), therefore demonstrating that the changes produced in a given month have indirect effects on the abundance of individuals in subsequent months. Results from the present study indicate that SBT did not have a big influence on recruitment in the same month, however, it would be interesting to test environmental variables with temporal lags for future studies. Regarding bathymetry, in recent studies it has been observed to be an important variable that contributes to explaining the spatial distribution of both hake recruits (Sánchez and Serrano, 2003; Abella et al., 2005; Colloca et al., 2015; Paradinas et al., 2015; Pennino et al., 2019) and adults (Agostini et al., 2008). The final model predicted the highest occurrence and abundance probabilities of hake recruits to be at a range of between 100 and 250 m deep, with a constant optimum range between 140 and 200 m. The same bathymetric range for hake at age 0 along the shelf break was identified by the authors Colloca et al. (2015) in Mediterranean nurseries. In contrast, our results differ from the study carried out by Sánchez and Gil (2000) for the southern hake stock in the study area, which found a range between 90 and 180 m deep with a maximum peak of abundance at 100 m. However, these authors analyzed a previous period of years and the defined size of the recruits was <17 cm instead of <21 cm. These differences could be Frontiers in Marine Science | www.frontiersin.org 9January 2021 | Volume 8 | Article 614675