Modelling the spatio-temporal distribution of age-1 Bay of Biscay anchovy ( Engraulis encrasicolus ) at spawning time

1 AZTI-Tecnalia, Marine Research Division, Txatxarramendi Ugartea z/g, 48395 Sukarrieta, Spain. E-mail: libaibarriaga@azti.es 2 AZTI-Tecnalia, Marine Research Division, Herrera Kaia Portualdea z/g, 20110 Pasaia, Spain. 3 General Fisheries Commission for the Mediterranean (GFCM), Food and Agriculture Organization of the United Nations (FAO), Palazzo Blumenstihl, Via Vittoria Colonna 1, 00193, Rome, Italy. 4 King Abdullah University of Science and Technology (KAUST), Red Sea Research Center, Thuwal 23955-6900, Saudi Arabia.


INTRODUCTION
Spatial pattern at different life stages is a key factor determining population dynamics.For marine fish populations, processes influencing distribution of early life stages can be reconstructed using data from research surveys that target a specific life stage (Cotano et al. 2008), or by using simulation models, such as individual-based models (Allain et al. 2007).Information about spawner distribution contributes to understanding reproductive strategy and life cycles (Brochier et al. 2009, Mullon et al. 2002).Although not exclusive of each other, hypotheses about factors determining the spawning habitat can be divided into two large lines of thought: environmental factors and homing behaviour (e.g.Jørgensen et al. 2008).Homing behaviour is generally accepted for anadromous fish and shown for fish species spawning in estuaries (Thorrold et al. 2001).For small pelagic fish homing behaviour is not usually considered although it has been shown that it could explain spawning patterns (Mullon et al. 2002).Further, habitat preferences can also be influenced by demographic structure (Planque et al. 2011a).For example, different size or age groups may show different spatial patterns (Werner and Giliam 1984).Social interaction has also been hypothesized to play a role in maintaining migrations to the same suitable spawning areas from one generation to the next (ICES 2007, Planque et al. 2011a).
These considerations have clear implications for management.Obviously, spatial protection measures such as marine protected areas are different for species with homing or roaming behaviours (Dunlop et al. 2009, Thorrold et al. 2001).Also, to help stocks recover from collapse it is necessary to apply measures not only to maintain the overall level of spawning biomass but also to rebuild the extant spatial structure (Dickey-Collas et al. 2010, Petitgas et al. 2010).
Bay of Biscay anchovy (Engraulis encrasicolus, L.) show a distinct spatial structure, which is believed to vary for life stages and between years of contrasting environmental and/or population structure.The adult population is generally dominated by one-year-old individuals, except in the years when recruitment fails and two-year-old individuals become predominant (Uriarte et al. 1996).Data from acoustic (Simmonds and MacLennan 2005) and daily egg production method (DEPM;Lasker 1985, Parker 1980, Stratoudakis et al. 2006) surveys conducted since 1989(ICES 2008b, Santos et al. 2011) have been used to identify the characteristics of the main spawning areas (Irigoien et al. 2008, Motos et al. 1996, Bellier et al. 2007) and provide a comprehensive spatial overview of the adult population.Cort et al. (1976) and Motos et al. (1996) described the adult age structure at spawning.Petitgas et al. (2003) incorporated a linear regression between fish length and bottom depth into the kriging system to examine inter-annual variability in the spatial length distribution of the anchovy.However, none of these studies has provided a spatially explicit age-structured model.
The objective of the paper is to better understand the relation between position, environmental conditions and age structure of the population with a view to future spatial protection measures.We analysed the spatial and temporal variability in the age structure of the anchovy population in the Bay of Biscay at spawning time using data from DEPM surveys conducted in 1990-1992, 1994-1995, 1997-1998 and 2001-2007, which is one of the longest series of DEPM surveys worldwide (Stratoudakis et al. 2006).The proportion of age-1 individuals in the population in relation to geographical variables was modelled using generalized additive models (GAMs, Hastie andTibshirani 1990, Wood 2006).The potential benefits of including environmental covariates, such as sea surface temperature (SST) and sea surface salinity (SSS), were also evaluated.In addition, we tested for differences in distribution of the age-structured population depending on whether the age-1 proportion level was high or low.The final model was used to hindcast the age composition of the stock from 2000 to 2005, identifying the main similarities and differences between these years.

MATERIALS AND METHODS
Adult samples from 15 Bay of Biscay anchovy DEPM surveys carried out by AZTI-Tecnalia were analysed: two surveys conducted in 1990 (May and June), 1991, 1992, 1994, 1995, 1997, 1998and 2001-2007(ICES 2008b, Santos et al. 2011, Somarakis et al. 2004).Samples of adults were obtained from several different sources, including the research vessel conducting both egg and adult sampling, a research vessel dedicated only to adult sampling, the French acoustic survey conducted by IFREMER in approximately the same area and over the same time period, and occasionally from commercial vessels operating at the survey time in the study area (ICES 2008b).Immediately after fishing, anchovy were sorted from the catch and a sample consisting of a minimum of 1 kg or 60 anchovies was selected at random, weighed and measured into length classes of 0.5 cm.Whenever possible the otoliths of the fish were extracted and classified at the laboratory into age 1, age 2 and age 3 or older groups.When samples were obtained from the commercial purse seines, a sample of nearly 2 kg was selected from the catch and placed directly into 4% formaldehyde for later processing in the laboratory.When no otoliths were collected, ad hoc, concurrently constructed age-length keys were used to transform the observed length-class frequencies into age groups.The geographical position, bottom depth, fishing depth, and sampling date and time were recorded routinely for each sample.
Environmental data were not always collected during the adult sampling; however, SSS and SST were available from the egg surveys conducted in parallel and covering the same area at the same time.The SSS and SST were modelled for each survey using GAMs with a normal distribution, identity link and bivariate smooth function for the location variables longitude and latitude with a high degree of freedom.These models were used to obtain SST and SSS estimates for each of the adult samples.
The main effects of the available variables, except those related to the geographical position, i.e. bottom depth, SST and SSS, were determined by comparing their distribution with the distribution of the age-1 proportion.The age-1 proportion was modelled using generalized additive models with binomial distribution and logit link (Hastie andTibshirani 1990, Wood 2006): where p is the expected probability of age-1 individuals, s j (x) is a smooth function and x j is a vector of covariates.Different options were tested for the intercept of the model, b 0 .The base case considered a unique common intercept for all the surveys.The alternative case involved changing the intercept depending on the overall age-1 proportion for each of the surveys, so that the average differences between surveys were accounted for in the model.The overall age-1 proportion of each of the surveys was given either as a factor (denoted as P factor ) or as a number between 0 and 1 (denoted as P depm ).The former variable P factor had the classes very low, low, medium, high and very high, derived from the mixture probability distribution of the age-1 proportion from the stock assessment (Ibaibarriaga et al. 2008, ICES 2008a) as is explained below.In this case, the intercept of the model changes for each class of P factor .The latter variable P depm was the age-1 proportion estimated from the DEPM surveys (ICES 2008b, Santos et al. 2011, Uriarte 2001).In this case the model intercept was a smooth function of the variable P depm with the upper limit of the smooth basis dimension set equal to 5.
The initial model was constructed using a bivariate smoothing function that depends on the geographical variables, longitude and latitude.The rest of the variables (SST, SSS and the logarithm of bottom depth to correct for skewness) were first added as univariate smooth functions.We then analysed whether a bivariate smooth function representing the interaction of these variables was statistically more appropriate.Formal comparison of nested models including univariate and bivariate smooth functions of the same covariates can be used to test whether additive or multiplicative structures are more appropriate.In our case, models including covariates as main effects (within univariate smooths) and as interactions (within bivariate smooths) were attempted (though not shown here), but they were considered overfitted for the amount of data available.Therefore, the inclusion of each covariate was restrict-ed either as a main effect or in conjunction with some other covariate as an interaction.
In addition, we examined whether changing the bivariate smooth function of latitude and longitude according to the overall age-1 proportion in the survey (low or high) helped to improve the model.A new categorical variable, P low-high , representing the overall age-1 proportion of each of the surveys, was derived from P factor by collapsing its classes into two levels.This new variable took the value "low" when P factor was very low, low, or medium, and "high" when P factor was high or very high.A bivariate smooth function of latitude and longitude was fitted for each level of the factor P low-high .
All the models were fitted using the package mgcv in the R statistical software (R Development Core Team, 2011; www.R-project.org).This package provides functions for generalized additive modelling with automatic smoothing selection (Wood 2000(Wood , 2004)).As the scale parameter of the binomial distribution is known (Wood 2006), the smoothing parameters for each model term were chosen to minimize the UBRE (Un-Biased Risk Estimator; Craven and Wahba 1979) score.Model selection was based on the procedure given by Wood and Augustin (2002).Consistently with the smoothing parameter selection, it was considered that the lower the UBRE score, the better the model.
The sensitivity of the best model for the surveys used in model fitting was tested by removing each survey in turn and analysing changes in the fitted smooth functions and the corresponding UBRE scores.
The best model was used to predict the spatio-temporal distribution of the age-1 proportion from 2000 to 2005.SST and salinity fields needed for this prediction were derived from the results of the high-resolution 3D prognostic ocean model ROMS (Regional Ocean Model System, Shchepetkin and McWilliams 2005), which was applied to the region and forced by detailed atmospheric, hydrologic and oceanic forcing.The model domain covers the whole Bay of Biscay, extending from the French and Spanish coasts (40.5°N) to the south of the United Kingdom (about 52.5°N) and 13°W.In this present configuration, which is an extension of a configuration limited to the southern Bay of Biscay (Ferrer et al. 2009), the bathymetry was obtained by interpolation following an optimization analysis of the ETOPO2 (2 minute digital Elevation TOPOgraphic model, Smith & Sandwell, 1997), GEBCO (General Bathymetric Chart of the Oceans), and IBCM (International Bathymetric Chart of the Mediterranean) data sets.Bottom depth was obtained from this bathymetry.ROMS for the Bay of Biscay were used to calculate the primitive equations on a 6.6-km grid on the horizontal, and a grid of 32 unequally distributed s-levels (bathymetric levels, with more resolution at the surface to resolve stratification) on the vertical.The initial and boundary conditions for temperature and salinity were interpolated on the grid from the World Ocean Atlas 2005 (WOA05).The water level was specified for the initial condition but also at each time step along the

RESULTS
The main characteristics of the DEPM adult surveys used in this paper are summarized in Table 1.In most years, the raw average age-1 proportion was above 0.5, except in 1994, 2002 and 2005.The lowest average age-1 proportion (less than 0.2) was observed in 2002.Figure 1 shows the spatial location of all the samples.Even if the space was irregularly covered, the adult samples in the DEPM were taken from where anchovy were detected, and are therefore assumed to represent the population correctly.
The distribution of SST, SSS and bottom depth was compared with the age-1 distribution conditioned by these covariates (Fig. 2).The age-1 proportion increased with SST at temperatures above 14°C and decreased with salinity up to 35, above which it increased again.However, there were few samples with very low salini-ties (less than 32) and most of these corresponded to the Gironde river plume in 2004.When the age-1 proportion was compared according to bottom depth, we found that the largest proportions were within the first 100 m.
The mixture distribution of the age-1 proportion from the ICES stock assessment (Fig. 3) allowed the proportion of age-1 biomass to be separated into five classes with well-defined peaks separated by local minima at 0.37, 0.52, 0.66 and 0.88.Depending on the interval in which the median age-1 proportion fell, each survey was categorized as very low, low, medium, high and very high age-1 proportion (P factor , Table 1).Note that the surveys corresponding to the low age-1 proportion class, i.e. 1987 and 1989, are not investigated in the present paper.The simplified age-1 proportion factor (P low-high ) classified into low and high age-1 proportions (Table 1) was given by the local minimum at 0.66.
The initial model with a bivariate smooth function that depends on the geographical variables, longitude and latitude, with a common intercept for all surveys, explained 25% of the deviance of the model with a UBRE statistic of 22.When the intercept was modelled either by P factor or a smooth function of P depm , the UBRE statistic decreased by up to 10.5 and 9.9 respectively, which significantly improved the model fit (Table 2).The coefficients of P factor increased as the age-1 biomass proportion increased (-1.3 for very low, 1.7 for medium, 2.3 for high and 4.2 for very high), but this relationship was not linear, as shown by the smooth function of P depm in Figure 4 (top row, left panel).
When a univariate smooth function of the logarithm of bottom depth, SST and SSS was added to the initial model with a constant intercept, the UBRE statistic decreased and the percentage of deviance explained increased (to 13.6% and 53%, respectively).Modelling the intercept depending on the age-1 biomass proportion further improved the model by decreasing the UBRE statistic to 9.3 for P factor and 9.1 for the smooth function of P depm (Table 2)., 2008a).The vertical dashed lines define five levels of relative recruitment (very low, low, medium, high and very high) that will be used to define the categorical factor P factor representing the overall age-1 proportion in each survey.Taking into account the interaction between the logarithm of bottom depth and the SST and SSS covariates, the best model was obtained with a bivariate smooth function for longitude and latitude, a bivariate smooth function for the logarithm of bottom depth and SSS and a univariate smooth function for SST (Table 2 and top row in Fig. 4 and left panel in Fig. 5).The pair SSS and the logarithm of bottom depth showed the highest correlation (0.47), which further justified their inclusion in the model as a bivariate smooth function.Again, the best model (UBRE=8.7)was that with the intercept modelled as a function of P depm .
Allowing the bivariate smooth function for longitude and latitude to change depending on whether the age-1 proportion was low or high (P low-high ) reduced the UBRE score to 7.6 and increased the percentage of deviance explained to 73%.The residual diagnostics showed no pattern (not shown here).The univariate smooth functions depending on P depm and SST and the bivariate smooth of the logarithm of bottom depth and SSS remained basically the same (bottom row in Fig. 4).The bivariate smooth functions for longitude and latitude in the models that can change depending on P low-high are shown in Figure 5 (mid-and right panels).For both low and high cases the largest age-1 proportion (above 0.8) was found close to the coast, especially in the Gironde and Adour rives plumes.The proportion then decreased as depth increased, with a Fig. 5. -Bivariate geographical smooth function depending on longitude and latitude when the bivariate geographical smooth function of longitude and latitude is independent (left panel) and dependent (mid and right panels) on P low-high .The mid-and right panels correspond to the low and high age-1 proportion levels, respectively.
higher drop (up to 0.2) in the case of low proportions.
The low case showed an increase in the age-1 proportion on the western Cantabrian coast, whereas the high case showed an increase in the age-1 proportions up to 0.7 in the oceanic regions after the shelf break.The bivariate smooth function for longitude and latitude regardless of P low-high (left panel in Fig. 5) represented the largest age-1 proportion close to the French coast, which decreased sharply towards the shelf edge, and was then followed by a slight increase in more oceanic waters and in the western areas of the Cantabrian coast.
The sensitivity to the input data of the final model (bivariate smooth function for longitude and latitude that changed depending on P low-high , the intercept modelled as a function of P depm , a univariate smooth function for SST and a bivariate smooth function for SSS and the logarithm of bottom depth) was analysed by removing one survey at a time.The UBRE score considerably improved when 1991, 1997 and 2003 were removed (Fig. 6).However, the patterns in the smooth functions were almost unaltered (not shown here), which indicates that the results were quite robust to the individual surveys. The

DISCUSSION
In this study we modelled the spatial distribution of the age-1 proportion at spawning of Bay of Biscay anchovy as a function of geographical and environmental variables.We have shown that the largest age-1 proportions are especially associated with the Gironde and Adour river plumes.Lower proportions of age-1 individuals are found on the shelf break and in oceanic waters, where older individuals are usually prevalent, as already described by Cendrero et al. (1981), Cort et al. (1979), Motos et al. (1996) and Vaz and Petitgas (2002).The role of geographical position has also been found to be determinant for the zooplankton distribution in the Bay of Biscay (Irigoien et al. 2011).Permanent hydrographical features are related to the topography (e.g.fronts related to river plumes and shelf break or eddies to the canyons; Albaina and Irigoien 2007).Therefore there seems to exist a "landscape" where the geographical position guarantees the repetition of certain conditions (enhanced food abundance, enhanced retention) rather than the particular value of the environmental factors themselves, such as temperature or salinity.The existence of such permanent landscape structures even in the water column has an evolutionary value in terms of homing behaviour: the survival related to spawning at natal date and position is not aleatory, but related to the existence of some permanent mechanisms that are repeated in time.
This size-dependent spatial pattern of anchovy, modelled in Petitgas et al. (2003) and observed yearly in the spring research surveys (ICES 2008b), has also been observed for other marine fish (Gordoa andDuarte 1991, Loots et al. 2011).This observation does not fit well with a pure homing behaviour: a sizedependent habitat selection (Planque et al. 2011a) would also be required.According to Werner and Gilliam (1984) body size affects resource-use capabilities (such as food and habitat use) and vital rates (growth and mortality) and individuals would tend to occupy the most suitable habitat according to their size and energy requirements.However, at least for anchovy juveniles in the Bay of Biscay, there are no differences between individuals growing close to the coast and those occurring off the shelf break (Irigoien et al. 2007, Aldanondo et al. 2010).The major coastal prevalence of one-year-olds in the following spring shown here could be related to both the homing behaviour at the end of autumn (Boyra et al. in press, Uriarte et al. 2001) and a sufficient productivity during winter at the coastal, estuarine and river mouth-influenced areas to facilitate winter survival of the juveniles prior to their first spawning.Thus, the homing behaviour of juveniles towards coastal areas in late autumn and winter will be much enhanced by the productivity in these areas given the high energy demand for growth and survival in winter.On the other hand, homing behaviour would be relaxed as the energy demands decrease as the fish gets older.This would facilitate the differential spatial distributions observed between one-and two-year-olds or older, which may find sufficient food supply while spawning in the shelf break.
Modelling the intercept as a function of the actual relative level of recruitment significantly improves the predictions and indicates that the degree of prevalence of age-1 individuals over older individuals (2+) in space depends greatly on the overall relative strength of recruitment in the population.Therefore, the age-1 proportion is above 50% in all areas in years of good relative recruitment, while in years of very low recruitment (as in 2002 and 2005), older individuals (age 2+) predominate in most of the surveyed areas.In years with a high age-1 proportion, this age class also becomes dominant in oceanic waters, whereas in the case of a low proportion the presence of age-1 individuals decreases as the water becomes deeper (Motos et al. 1996, Santiago andSanz 1992).The main spatial similarities between the model predictions from 2000 to 2005 could be attributed to the different geographical smoothing depending on the low or high age-1 proportions.Space occupation has also been shown to be density-dependent for this and other anchovies (Barange et al. 2005, Somarakis et al. 2004).More space is occupied when biomass is high, a finding which could also be attributed to the space occupied by the recruiting year class.When preferred habitats are occupied, age-1 individuals move to spawn to other areas.Again, this suggests that, if existing, homing behaviour is not very strong and could be modulated in relation not only to age but also to density.
In terms of environmental conditions, the age-1 proportion increases in SST above 16ºC and in low-salinity areas (around 34) in oceanic waters (around 500 m depth).These patterns represent deviations from the average spatial pattern (represented by the geographic functions) that could be attributed to specific environmental conditions.For instance, the low-salinity areas might be related to larger areas of influence of the river plumes in particular years, as occurred in 1997, when low-salinity waters were found in the oceanic region where one-year-old anchovies were also predominant.
A possible source of bias to our analysis is that the age structure can affect the spatial and also the temporal spawning pattern.Motos et al. (1996) and references therein suggest that older individuals spawn earlier in coastal areas and then move to the shelf break and oceanic waters as the spawning season progresses, whereas spawning of younger individuals starts later (after the spring warming up period, for temperatures above 14°C) mainly in shallow waters near the river plumes.A temporal variable is not included in the model in this study.The DEPM surveys are assumed to provide information for the spawning peak, and therefore the data available do not provide enough contrast for this analysis.Collection and analysis of age structure data obtained routinely from April to July could provide further insight into temporal patterns related to the age structure at spawning.Catch data collected throughout the year could help to obtain an overall view (Uriarte et al. 1996), but good spatial coverage would also be necessary for modelling purposes.Probably genetic and otolith microchemistry studies will be needed to really evaluate the homing behaviour hypothesis.
Robustness of the results was tested by studying the influence of each survey on the model.The model was almost unchanged when the surveys were left out one at a time.This confirms the major model inferences with the data at hand (15 surveys).However, additional caution should be taken when using the model for making projections outside the range of observed years and conditions (Dormann 2007, Planque et al. 2011b).
There is considerable previous research analysing ichthyoplankton samples for characterization of spawning habitat in order to better understand the ecological context of spawning (Bernal et al. 2007, Planque et al. 2007, Stratoudakis et al. 2003).However, the adult samples in routine DEPM applications are usually only used for estimating the daily fecundity in order to obtain the spawning stock biomass.Uriarte (2001) extended the traditional DEPM to obtain population at age estimates and their corresponding variances for Bay of Biscay anchovy.This paper goes a step further and models the age structure in space, using some of the environmental data available from the simultaneous ichthyoplankton sampling.A similar exercise was conducted by Begg and Marteinsdottir (2002) for mature cod in Icelandic waters using data from spring groundfish surveys.The analysis here was restricted to the age-1 class that generally forms the larger fraction of the population, but could be expanded to other age classes as done by Bailey et al. (1998) for herring acoustic surveys.
The results obtained on spawner distribution can help to better understand various elements of the life cycle (Brochier et al. 2009).Egg production can be attributed to the first or repeated spawners and it can be determined whether different age structures at spawning affect recruitment survival (Marteinsdottir and Thorarinsson 1998).Older and larger females might have higher per capita fecundity (Kennedy et al. 2007, Morgan et al. 2009) and are often able to produce higher-quality eggs (Brooks et al. 1997), which results in higher survival (Marteinsdottir and Steinarsson 1998).
The information on the spatial dynamics at spawning according to age classes can be applied directly to designing spatial measures for managing the stock (Dickey-Collas et al. 2010, Petitgas et al. 2010).If the management goal is to maximize the first spawning, it could be useful to close coastal areas in spring (Lehuta et al. 2010).However, if there are low age-1 population levels, and especially when the newly recruited individuals form only a small fraction of the population, it is essential to also protect the areas where the older individuals are more abundant in order to maintain the overall level of the spawning biomass and also rebuild the contingent structure of the population.Alternatively, the implications of different spawning site fidelity assumptions for the design of marine protected areas are discussed in Thorrold et al. (2001) and Dunlop et al. (2009).

Fig. 1 .
Fig. 1. -Location of the daily egg production method adult samples from 1990 to 2007.The lines represent the coastline and the 100 and 200 m depth contour lines.

Fig. 2 .
Fig. 2. -Number of stations (vertical line) and distribution of age-1 individuals (lines) in relation to sea surface temperature, sea surface salinity and the logarithm of bottom depth.The solid and dashed lines represent the median age-1 proportion and the 90% probability intervals, respectively.The depth classes are defined by the 0, 50, 100, 200, 500, 1000 and 3000 break points.

Fig
Fig.4.-Smooth functions of the best model in terms of un-biased risk estimator score when the bivariate geographical smooth function of longitude and latitude is independent (top row) and dependent (bottom row) on whether the overall age-1 proportion is low or high (P low-high ).From left to right the univariate smooth function of P depm , the univariate smooth function of sea surface temperature (SST) and the contour lines of the bivariate smooth of sea surface salinity (SSS) and depth are shown.

Fig
Fig. 6. -Sensitivity to surveys.Un-biased risk estimator (UBRE) scores of the final model when each survey (x-axis) is excluded from the model fitting.The horizontal dashed line is the UBRE score of the final model when all the surveys are used.

Fig. 7 .
Fig. 7. -Maps of the age-1 proportion from 2000 to 2005 predicted according to the final model.

Table 1 .
-Summary of the Bay of Biscay daily egg production method (DEPM) adult surveys from 1990 to 2007.The age-1 proportion corresponds to the average proportion of all the samples per survey.P depm , overall age-1 proportion estimated from the DEPM surveys; P factor , a factor of five classes (very low, low, medium, high and very high) indicating the overall age-1 proportion as derived from the stock assessment; P low-high , a factor of two classes (low and high) indicating the overall age-1 proportion as derived from the stock assessment.

Table 2 .
-Summary of the generalized additive models without allowing spatial changes that depend on P low-high .The first two columns describe the model (the intercept and the smooth functions of the covariates).The last two columns give the un-biased risk estimator (UBRE) score and the % of deviance explained by each model.