Abstract To date, there are numerous transport simulation studies demonstrating the relevance of the hydrodynamics for the advection, dispersion and recruitment of early stages of marine organisms. However, the lack of data has conditioned the use of realistic locations for the model setup and configuration in transport studies. This work (I) demonstrates the key role played by the use of the realistic initial position of the eggs of small pelagic fishes in the analysis of late-larval recruitment in coastal nursery areas and (II) provides a general solution for deriving future egg positions and abundances from adult biomass obtained from acoustic surveys and available fecundity data. Using European anchovy in the NW Mediterranean as a case study, we first analyzed the impact of the initial location, timing, egg buoyancy and diel vertical migration of larvae on the potential late-larval recruitment to coastal areas. The results suggested that prior knowledge of the initial spawning grounds may substantially affect the estimates of potential recruitment. We then integrated biological and acoustics-derived data (the biomass and size structure, sex ratio, a weight-batch fecundity model, mean weight, number of fish and mean spawning) to build a predictive model for interannual egg production. This model was satisfactorily contrasted with field data for two years obtained with the Daily Egg Production Method (DEPM). We discuss our results in the context of the fluctuations of European anchovy egg abundance from 2003 through 2010 in the NW Mediterranean and in terms of the potential applicability of the acoustics-based spatial predictive egg production model.

Citation: Ospina-Álvarez A, Bernal M, Catalán IA, Roos D, Bigot J-L, Palomera I (2013) Modeling Fish Egg Production and Spatial Distribution from Acoustic Data: A Step Forward into the Analysis of Recruitment. PLoS ONE 8(9): e73687. https://doi.org/10.1371/journal.pone.0073687 Editor: Hans G. Dam, University of Connecticut, United States of America Received: May 17, 2013; Accepted: July 22, 2013; Published: September 16, 2013 Copyright: © 2013 Ospina-Álvarez et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Funding: This work was supported by project PERSEUS (FP7-287600). Egg samples were obtained from the MPOCAT-DEPM surveys financed by Generalitat de Catalunya. A. Ospina-Álvarez benefited from a PhD grant of the JAE program (Spanish Science and Technology Council, CSIC). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist.

Introduction The recruitment success of small pelagic fishes (SPF) is strongly dependent on the vulnerability of their planktonic early life stages. For these species, the origin and trajectories of eggs and larvae linked with the source-sink dynamics in metapopulations might be the cause of the frequently observed large-scale variations in recruitment [1]. Since the advent and widespread use of individual-based models (IBMs) that combine realistic ocean dynamics with transport behavioral models in fish [2]–[5], the study of recruitment has become less of a black box and more of a heuristic exercise. Spatially explicit individual-based models (SEIBMs) allow the understanding of numerous processes affecting population fluctuations in recruitment and have provided insight into past fluctuations as well as into potential changes in a suite of scenarios [2], [6]–[9]. Within SEIBMs, egg and larval transport is defined as the horizontal translocation of a particle in a bidimensional plane from an initial to a final position [9]. In contrast, egg and larval dispersal refers to the spread of eggs and larvae from a spawning source to a nursery site [9]. This definition is consistent with the natal dispersal concept from terrestrial ecology [10], [11]. The use of mathematical models that estimate the trajectories of eggs and larvae using average currents provides an opportunity to investigate general ecological questions as well as to provide qualitative assessments regarding the connectivity of specific regions and populations [12]. Larval dispersal models can be described with a dispersal kernel estimating the probability density function of the number of larvae and the distance from the spawning location (i.e., the dispersal kernel) [13]. A dispersal kernel combines both physical (e.g., ocean currents, temperature, salinity) and biological (e.g., buoyancy, growth, vertical migration) processes [14]. The dispersal of early life stages has implications for the structure and dynamics of populations and for the evolution and distribution of species [14]–[16]. In pelagic species, spawning time and location are important factors conditioning dispersal patterns and larval survival, as are other factors influencing spawning, including adult abundance, the age and condition of spawners and fertilization success. As important factors determining egg and larval dispersal, the release location, spatial abundance and aggregation of particles (eggs or larvae) have been addressed from various points of view. The traditional approach is to release particles from fixed points [17]–[20] or randomly over the total extent of the spawning area [7], [21]–[24]. Other, more realistic models, use as initial fields the particle distribution data from surveys, the long-term mean distribution of eggs or the abundance and seasonal occurrence of eggs and larvae, derived from available information [25]–[28]. However, to our knowledge, there are no existing comparison of methods and no consensus on the most accurate approach. Concurrently with these developments, other methods aimed at providing data for €classical€ stock assessments, such as the Daily Egg Production Method (DEPM) or the analysis of acoustic data. The spatial component and the information on adult biomass are, potentially, very useful for SEIBM exercises. DEPM has several advantages over other methods, such as acoustic assessment, because DEPM provides information about the reproductive condition of females, the distribution of the spawning and reproductive habitats, the egg production and mortality during the peak of spawning [29]–[32]. However, DEPM requires considerable time, both in the field and in the laboratory, and it is therefore relatively expensive. Moreover, in some areas and for some species time series from DEPM surveys are unavailable. In contrast, acoustic techniques have been used to estimate the size of fish populations for sufficiently long-time periods because these techniques require less time and, consequently, are less expensive. In this work, we take the European anchovy Engraulis encrasicolus in the NW Mediterranean as a typical case study of a relatively data-rich species for which further progress may be possible within the framework presented above. Although there is a reasonable understanding of the localities of spawning and nursery grounds in European shelf seas, information on late-larval and juvenile habitats is sporadic, and transport between these locations is less well understood [33], [34]. The spawning areas and the peak dates of high egg production in the NW Mediterranean for European anchovy have been established based on the information from egg abundance samples collected during DEPM surveys in 1993, 1994, 2007 and 2008 [35]–[38]. However, despite the exhaustive research on anchovy in the region, there are no available long time series on the distribution and production of eggs. Such time series are essential for conducting studies of the interannual variability of anchovy recruitment. However, a time series of acoustic surveys conducted by the French Research Institute for Exploration of the Sea (IFREMER) (2003–2010) to assess the biomass of the anchovy population is available for the region. In this context, the objectives of this study were (I) to test if the final position of small pelagic fish larvae in transport simulation studies significantly differs between i) random and ii) realistic initial spatial distribution of eggs. If differences are significant, we aimed (II) to provide a method to estimate the daily egg production (number of stage I eggs per square meter of sea surface) from adult biomass obtained by means of acoustic surveys and available fecundity data. To accomplish our first objective, we used a SEIBM to perform transport simulation experiments addressing the spatial aggregation and distribution of the early eggs of the anchovy. For the second objective, to predict anchovy egg production interannually, we built an integrated model using the spawner biomass and fecundity data series from the acoustic surveys. The initial state of the simulations was determined using acoustic data from the years 2003–2010. The novelty of this work consist in providing a predictive egg production model based on acoustic data and biological information when egg production surveys are not available. Given the high costs of research on fisheries, the scientific community is obliged not only to provide timely and increasingly accurate assessments but also to optimize the return from the funds invested in data collection.

Discussion To date, there are numerous transport simulation studies demonstrating the relevance of the hydrodynamics for the advection, dispersion and recruitment of early stages of marine organisms [7], [17]–[28]. However, the lack of data has conditioned the use of realistic locations for the model setup and configuration in transport studies. In this work, (I) we found that the initial position of eggs conditions the localization and density of larvae at recruitment time, and (II) we proposed a method to derive the relative abundance of eggs at spawning time from acoustic surveys to explore population dynamics with SEIBM applications. The proposed method is not only useful for the setup and configuration of SEIBM’s but also for a broad variety of case studies. For example, to calibrate the quantitative predictions of acoustic surveys and continuous under-water fish egg sampler (CUFES) methods for Small Pelagic Fish stock assessments. Currently, there are many small pelagic fisheries with a long temporal series data from acoustic surveys that could benefit from our findings, e.g.: Anchovy Engraulis capensis, sardine Sardinops sagax, and round herring Etrumeus whiteheadi in the southern Benguela upwelling region [79], [80]; european anchovy in the Bay of Biscay [81]; anchoveta Engraulis ringens and sardine in the Humboldt Current ecosystem [60], [82], [83]; anchovy Engraulis japonicus in the East China Sea and the Yellow Sea [84]; anchovy Engraulis anchoita in the Brazilian shelf [85]; sardine and northern anchovy Engraulis mordax in the California Current System [86], [87]. In this study, the SEIBM experiments performed with random and real distributions of anchovy eggs at spawning time showed that the spawning location strongly influenced the larval transport (trajectory), retention and potential destination because of interactions with local hydrographic features. Physical processes act on the transport retention of larvae during their pelagic phase (e.g., wind-driven circulation, eddies, seawater density, fronts). For example, the spawning time and location affect cod larval transport [88]; consequently, the synchronicity of spawning, age, size and condition of spawners and fertilization success are processes and factors that may be associated with late fish larval recruitment [9]. This work confirms that these processes and factors determined the larval recruitment success of anchovy in 2007 and 2008. Also, larvae from the spawning areas GoL East, GoL West, Ebro Delta and GV had higher probability of reaching nursery areas, comparing to larvae from Palamós and Barcelona. However, the principal result of the inclusion of the real distribution of eggs at spawning time in our SEIBM was the retention of 14-mm larvae in the GoL (Fig. 5, 6). The inter-annual variability in the retention of particles in GoL might be linked with mesoscale processes related with the formation of sub-meso and mesoscale eddies [28]. In this work, the eggs and larvae transport and distributions were strongly related to the hydrodynamic structures on the shelf and the offshore circulations associated with the North Current for 2007 and 2008. This work strongly evidences that long-time series of accurate egg distribution data are necessary, as model inputs, to make valuable inferences rather than starting the simulations with a random distribution of eggs over the spawning area. Additionally, the results suggest that the anchovy larvae from the GoL never supplied the southern late larval recruitment in the GV. Although anchovy larvae reached the Balearic Islands, this event was not common and occurred under particular hydroclimatic conditions (e.g., a strong North Current). In the other hand, the combination of spatial and temporal scales, together with the biological characteristics of individuals, may have ecologically meaningful effects on population connectivity [5]. For example, a study on the identification of pelagic fish subpopulations using amino acid (AA) composition differentiated the northern and southern anchovy subpopulations in the Western Mediterranean Sea [89]. However, DNA and AA studies need to be linked with physical dynamics, particularly at intermediate and small scales (i.e., at scales of kilometers or less) to better understand the dynamics of population connectivity. Moreover, it is necessary to develop an improved understanding of larval behavior and of environmental gradients and cues. It is clear that further studies in population connectivity are necessary to determine the importance of these findings. The spatial structure of the spawning distributions of anchovy allowed GAMs to be fitted that explained between 60% and 89% of the deviance in the survey station values of the anchovy biomass and gave significant estimates of spatial and temporal trends for data collected in eight acoustic surveys from 2003 to 2010. This analysis allowed the estimation of egg production (Figs. 10, 11) through the integration of the acoustic surveys with the available fecundity data (Tables 2, 4). There are several advantages of GAMs over stratified sample survey methods. For example, the application of GAMs to survey data from other fish stocks showed a substantial reduction in coefficients of variation of egg abundance, while increasing estimation precision [61], [64]. In this work, the principal benefits of GAMs were that they provided an objective method for data interpolation into unsampled areas; these models allowed density predictions within the survey area and provided maps that allowed a visual inspection of the spawning distribution. In addition, the increase in precision and the spatially explicit results provided insight into the trends of egg production by area. The visual inspection of egg production estimates from the GAMs was consistent with survey station values of eggs per square meter from DEPM surveys in 2007 and 2008 (Fig. 3, 10). The integration of acoustic surveys and fecundity data potentially allows the modeling of changes of egg abundance over time as well as variation in space. The estimates of egg abundance obtained represent a desirable option for filling the gaps in the time series of egg production estimates for anchovy. In the study area, the visual inspection of habitat adequacy maps indicated that the occurrence of suitable sites was subject to temporal and spatial variability. The spawning intensity and location of the spawning grounds have a strong seasonality. The occurrence of anchovy is usually associated with point sources of nutrients that enhance productivity locally, e.g., river runoff or local upwelling [66]. The most persistent suitable locations for anchovy spawning were identified in the northeastern part of the GoL. In this area, the influence of the Rhone River plume and an existing upwelling along the coast yield a local increase in productivity (i.e., 2007 and 2008, Fig. 7). The GoL, receiving discharges from the Rhone River, has much higher nutrient and phytoplankton concentrations than the adjacent open NW Mediterranean [90], [91]. In the GoL, ocean color studies using remote sensors (e.g., CZCS, MERIS) have revealed strong spatial variation in the distribution of chlorophyll [92]–[94]. In this area, the coastal eddies transport organic matter from the coastal zone to the offshore domain [95]. This effect may be the reason that the coastal waters off Languedoc-Roussillon, influenced hydrographically by lagoons and shallow waters, were identified as occasional hotspots for anchovy reproduction. In contrast, certain locations were quite persistent, i.e., areas east of the meridian and the deepest areas with oceanic water influence. However, the location and abundance of anchovy spawners in the GoL could not be directly predicted from ocean productivity; the percentage of variation explained by the model relating anchovy spawners and Chl a concentration was less than 4%. In conclusion, the geographic distribution of Chl a values were only indicative for potential spawning habitat rather than realized spawning habitat. Additionally, similar tests were performed to predict the spawning habitat of anchovy from hydrodynamical conditions (i.e., sea surface temperature, salinity and water stratification); these test were not significant and results were not included in this work. Overall, these findings lead to the hypotheses that the abundance and distribution of the anchovy spawners is a consequence of density-dependent factors rather than the environmental conditions. Changes in the intensity or location of the principal currents, river runoff, upwelling intensity, climate variability and food availability are all known to affect the ecosystem in different ways and could be responsible for changes in the spatial distribution of spawning [66], [67], [96]–[101]. During the sampling time, the natural or anthropogenic environmental effects could explain the differences in the spatial distribution of spawning between years. Also, changes in the stock size and age structure of the parents could also produce these differences. In the GoL from 2002 to 2010, adult biomass estimates were highest in 2003–2004 followed by a marked decrease in 2005. The extension of the spawning area in 2003 was continual and reached a higher level than in subsequent years. The values for 2006–2010 did not reach the biomass estimates previous to 2005. However, 2003 and 2004 did not show evident egg production hotspots. The current study shows that in the GoL, as in large upwelling ecosystems, the spawning areas of anchovy are primarily situated in inshore, coastal waters. In the 2000s (2007 and 2008, see Fig. 3), the eggs showed a more coastal pattern. However, egg abundance in the 1990s (1992, 1993 and 1994 (1992, 1993 and 1994 [42], [96]) showed a more widespread distribution, reaching relatively great depths. Our data for the GoL support the hypothesis that the egg production of anchovy shows density-dependent effects. The stock fluctuations coincided with an expansion of the area of distribution in years of high abundance. An increase in the spatial extent of spawning areas was associated with an increase in egg production. Previous egg surveys in the NW Mediterranean have indicated that the locations of most intensive spawning of anchovy are reasonably consistent throughout the spawning season and are predictable at scales down to ca. 50 km or even less, with additional variability from smaller-scale patchiness and egg dispersal [66], [96]. The choice of oceanographic constraints in spawning locations may well be an evolutionary trait serving to maximize larval recruitment. In the SW Mediterranean, the limited coastal area providing shelter from strong currents has been proposed as a mechanistic driver to explain the spatially restricted anchovy spawning grounds [101]. We only can speculate about the ability of adult anchovy to return to their natal locations, although there are some indications that this may occur [6], [102], [103]. Salmonids have long been recognized for their long-distance homing migrations based on olfactory and other cues [104]. Similarly, we propose that chemosensory imprinting by juvenile could occur in anchovy while they reside in natal waters, as juveniles may spend 2 to 3 months in these nursery areas before migrating from these waters during the fall of their first year of life. Alternatively, the formation of large aggregations of fish during migration to feeding and spawning areas may allow younger individuals to learn migration routes from older population members [105], [106]. However, we do not intend to review or synthesize the vast field of natal or return homing here; for a detailed discussion of homing in related species, see [103], [105]. Furthermore, historical contingencies can shape the ecology and behavior of organisms [58]. Several studies have followed patches of marine eggs and larvae, but these efforts are not true measures of larval dispersal because the spawning locations and the ultimate destinations of the eggs and larvae were inferred [24], [107]–[109]. The current management of small pelagic fishes in the Mediterranean basin ignores the importance of recruitment levels to yearly yields; the sizes of the stocks of short-lived pelagic species depend primarily on recruitment, which is driven primarily by local oceanographic and climatic regimes. Based on the General Fisheries Commission for the Mediterranean (GFCM) recommendations and on current perceptions of proper management choices for stocks of small pelagic fishes [110], the development of methods for the adaptive management of these stocks is required in the region. The recent advances in methodology and standardization of acoustic surveys in the Mediterranean will provide useful data, particularly as a result of the use of such surveys as a means of estimating spawning areas and egg production to furnish a basis for SEIBMs that can predict recruitment areas and recruitment success. The results from this study could provide useful advances related to the set-up and design of this SEIBMs experiments. Moreover, we argue that the integration of the real distribution of eggs at spawning time with coupled hydrodynamic and biological models is a key step forward to improve our understanding of larval survival, growth and dispersal. The following step in SEIBMs for anchovy, and other small pelagic fishes in the NW Mediterranean, is to perform long-time series transport simulation experiments (e.g., the last decade), and to include behavioural rules to analyze relevant emerging properties (i.e. mortality).

Supporting Information Supporting information S1. The compressed file include: • SOURCE_FILE.R is a commented routine with equations, models and the possibility to generate plots and maps presented in this work. • Four functions, including findlimits.fun.R and standardise.fun.R, to assure the reproducibility of methods if geofun and shachar packages are not available from http://sourceforge.net/apps/trac/ichthyoanalysis/wiki. More details in the comments of the script file. • coast.dat is a file with the NW Mediterranean coastline necessary to generate maps. • RAW_2003_2010 PELMED DATA.csv is a comma separated value (CSV) file containing the acoustic PELMED surveys data for anchovy from 2003 to 2010. • RAW_FecundityParamAnchovy.csv is a comma separated value (CSV) file with the data presented in Table 1: Adult reproductive parameters of anchovy obtained with the application of the DEPM in different regions (Somarakis et al., 2004, and MPOCAT 2007 and 2008 data). https://doi.org/10.1371/journal.pone.0073687.s001 (ZIP)

Acknowledgments We are grateful to P. Garreau and A. Nicolle for their support with MARS–3D hydrodynamic model outputs, as well as to P. Verley for his assistance in ICHTHYOP code. We are grateful to S. de Juan for discussions and suggestions that greatly improved previous versions of this manuscript. Finally, the authors gratefully acknowledge the assistance of the captain and crew of the R/V García del Cid and R V LEurope for their help during the cruises.

Author Contributions Conceived and designed the experiments: AO MB IAC DR IP. Performed the experiments: AO DR JB IP. Analyzed the data: AO MB DR IP. Wrote the paper: AO. Designed the software used in analysis: MB AO. Wrote the majority of this paper with input from each listed author: AO. Performed major editorial support: IAC IP MB.