write a summary of articles

Contents lists available at ScienceDirect

Deep-Sea Research Part II

journal homepage: www.elsevier.com/locate/dsr2

Temporal changes in mesoscale aggregations and spatial distribution scenarios of the Peruvian anchovy (Engraulis ringens)

Giancarlo Morona,b,⁎, Paola Gallosoc, Dimitri Gutierreza,d, Josymar Torrejon-Magallanesa

a Instituto del Mar del Perú, Esquina Gamarra y General Valle s/n Chucuito, Callao, Peru b Universidad Nacional Mayor de San Marcos, Programa de Maestría en Matemática Aplicada, Lima, Peru c Institut de Recherche pour le Développement, Lima, Peru d Universidad Peruana Cayetano Heredia, Programa de Maestría en Ciencias del Mar, Lima, Peru

A R T I C L E I N F O

Keywords: Spatial variations Temporal variations Organism aggregations Recurrent areas Peruvian anchovy

A B S T R A C T

The Peruvian anchovy (Engraulis ringens) is the most important small pelagic of the Humboldt Current System (HCS), supporting the largest mono-specific fishery in the world. The spatial behavior of this species tends to be very dynamic at different spatial scales, influenced mostly by its biomass level and environmental factors. The aim of this study was to analyze temporal and spatial fluctuations in anchovy spatial distribution off Peru, based and modeled on acoustic data, focusing on large- and meso-scale spatial structures. We employed data from 41 scientific surveys (1994–2016) and Bayesian hierarchical spatial models to obtain the anchovy's spatial dis- tribution, allowing us to identify spatial structures at specific scales. Our results showed similar temporal trends in the number of mesoscale structures, their areas and anchovy density, exhibiting altogether two breakpoints in the time series: ~1999 and ~2013. The last period (2013–2016) was similar to the earlier one (1994–1999), in terms of low values of mesoscale structure indicators. On the other hand, we identified four spatial scenarios differentiated by the aggregative behavior, which were highly influenced by seasons and El Niño events. Each scenario had recurrent, or fidelity, areas placed in different locations. For instance, for the ‘El Niño scenario’ a particularly coastal recurrent area was identified, which might be a refuge zone for this species during these unfavorable events. Finally, we assessed differences in biomass estimates of each scenario. The highest biomass values were estimated for the ‘Summer favorable scenario’ and the lowest ones for the ‘El Niño scenario’, which supports the MacCall's basin hypothesis for this species. This study expands the current knowledge of the Peruvian anchovy and it is a first step to understand the effects on this species of the last El Niño events (2014–2016) that occurred in the northern Humboldt Current System.

1. Introduction

The spatial distribution of fish populations tends to be very dynamic at different temporal and spatial scales, even more if we refer to small pelagic fish (Bellier and Planque, 2007; Gutierrez et al., 2007; Saraux et al., 2014). This variability is driven by intra-specific needs such as reproductive success, feeding, predator avoidance (Hixon et al., 2012) and response to the environment, where the last one is commonly ac- cepted as the main cause of these fluctuations (Petitgas et al., 2012). The two features commonly studied to quantify changes in spatial distribution are presence area and population density, which translate in variations of the pattern of spatial aggregations, an essential feature of a fish population that could affect its sensitivity to fishing pressure and natural predation (Gutierrez et al., 2007; Perry et al., 2002).

The Humboldt Current System (HCS) and the Peruvian anchovy

(Engraulis ringens) offer an intriguing study case to research changes in spatial aggregations and factors affecting them. The HCS is classified as one of the four major Eastern Boundary Upwelling Ecosystems (EBUE), a group of ecosystems mainly characterized by a strong variability at different spatial and temporal scales, and a high productivity of plankton and fish, especially pelagic fish (Fréon et al., 2009; Gutierrez et al., 2016). Regarding this last feature, the HCS is the system where the highest number of fish per unit area is produced in the world's ocean (Chavez et al., 2008), essentially as a result of the Peruvian anchovy. The biological productivity of the HCS exhibits strong interannual variability, driven mostly by the El-Niño Southern Oscillation (ENSO) cycle (Bertrand et al., 2004; Chavez et al., 2003). This cycle has a period of 2–7 years and oscillates between a warm condition (El Niño: higher sea surface temperature or SST, weaker upwelling productivity, and deeper oxycline), and a cold condition (La Niña: lower SST,

https://doi.org/10.1016/j.dsr2.2018.11.009

⁎ Corresponding author. E-mail address: [email protected] (G. Moron).

Deep-Sea Research Part II xxx (xxxx) xxx–xxx

0967-0645/ © 2018 Elsevier Ltd. All rights reserved.

Please cite this article as: Moron, G., Deep-Sea Research Part II, https://doi.org/10.1016/j.dsr2.2018.11.009

stronger upwelling, and shallower oxycline), affecting all components of the ecosystem (Bertrand et al., 2008b; Fiedler, 2002).

The Peruvian anchovy sustains the most important mono-specific fishery in the world, and this activity yields an important income to the Peruvian economy (Chavez et al., 2008). This small pelagic species inhabits the upwelling cold coastal waters (CCW) and features fast growth rates, short lifespan, and schooling behavior, which allow it to display a fast response to environmental variability (Bertrand et al., 2008a; Bertrand et al., 2004; Swartzman et al., 2008). Since the Per- uvian anchovy is a key component of the HCS's trophic web, it has motivated numerous research efforts, many of them related to its spatial dynamics (Bertrand et al., 2008a; Bertrand et al., 2004; Bertrand et al., 2008b; Gutiérrez et al., 2007; Joo et al., 2014).

In general, small pelagic fish present an aggregative behavior at various scales, ranging from individuals within schools to the spatial pattern of school patches or school clusters (Barange et al., 2005; Barra et al., 2015; Bertrand et al., 2008a). Bertrand et al. (2008a) identified different spatial scales of organization for the Peruvian anchovy: fine scale (100 m – 1 km), sub-mesoscale (1–20 km), mesoscale (10 s km) and large-scale (100 s km). Fine to mesoscale structures are key for life- cycle processes since they are related to the area of retention, food production, and food aggregation (Bakun, 1996; Bertrand et al., 2008a). The oxycline, the depth range which delimits the top of oxygen- deficient waters, forms a sharp barrier for living organisms intolerant to hypoxia and has notable effects on the formation of these spatial or- ganizations (Bertrand et al., 2011; Ekau et al., 2010), therefore its temporal and spatial variability can trigger changes in the way that anchovies aggregate.

Grados et al. (2016) investigated the variations of the physical structures shaped by the ocean near-surface turbulent physical forcing, using the oxycline as a proxy, from fine-scale internal waves (IW) to mesoscale eddies. They found seasonal differences in the number and dimensions of fine-scale and sub-mesoscale structures between summer and spring, where a less IW activity was detected for summer seasons. Furthermore, the activity of these structures was greater along the shelf break, attributed to the interactions between tidal flows, stratification, and the continental slope. The effects of the spatiotemporal variability of these physical structures not only have a momentous role modifying anchovy aggregations but also in shaping the seascape from zoo- plankton to seabirds (Bertrand et al., 2014).

The present study focuses on the variability of the Peruvian anchovy aggregative behavior. Our first goal was to explore temporal changes in the Peruvian anchovy aggregation patterns from 1994 to 2016. The main tool employed thus far regarding this topic has been spatial in- dicators using raw acoustic data. However, this approach does not provide information about different spatial scales of organization and does not offer insights on the locations where aggregation areas take place (Bertrand et al., 2008b; Gutierrez et al., 2007; Joo et al., 2014). For these reasons, according to the dimensions proposed by Bertrand et al. (2008a), we modeled the total area of distribution of this species and identified two spatial scales of aggregation: large-scale and me- soscale structures. The second goal was to find spatial scenarios with similar aggregative behavior and recurrent, or fidelity, areas for each of them. Since Grados et al. (2016) found seasonal differences in the physical structures that influence anchovy spatial behavior, we also expected to find similar variations for these aggregations. Moreover, a recurrent area might be a proxy of preferred locations for anchovies under a determined scenario. Finally, our third goal was to explore possible influences of the biomass level on these scenarios.

2. Materials and methods

2.1. Data

We used 41 scientific surveys of pelagic resources carried out from 1994 to 2016, covering the area between 5°S and 16°S of the northern

HCS. Each acoustic survey was completed by the Marine Institute of Peru (IMARPE) and consisted of parallel cross-shore transects of ~100 nm long and with ~10–15 nm inter-transect spacing (Castillo et al., 2009) (Fig. S1). The sampled area covered the distribution do- main of the most important stock in terms of biomass: the north-center stock off Peru (Palomares et al., 1987). Simrad scientific echo-sounders were used to record the nautical-area-back-scattering coefficient (NASC, m2 nm−2) at each geo-referenced elementary distance sampling unit (EDSU, 1 nm), using 120 kHz for abundance determination. NASC values are normally used as a proxy of fish abundance in an EDSU (Simmonds and MacLennan, 2005) and we used them to model the anchovy spatial distribution. Extensive mid-water trawl sampling ac- companied the acoustic surveys for species identification, where echo identification was performed by using this trawl information and echotrace characteristics (Gutierrez et al., 2007). Each survey has a name code composed of eight characters: the symbol “Cr”, the year, and the starting and ending months (e.g. Cr160304 for a survey which started in March and ended in April 2016). After each survey, an an- chovy biomass estimate is obtained (for more details see IMARPE, 2016 and Simmonds et al., 2009)

2.2. Spatial model

Bayesian methods provide more realistic and accurate estimations of uncertainty because they allow the use of both observed data and model parameters as random variables (Benerjee et al., 2004). They also allow incorporation of the spatial component as a random-effect term in a natural way, thereby reducing its influence on estimates of effects of geographical variables (Gelfand et al., 2006). We used a Bayesian hierarchical model, which also permit to add non-Gaussian responses and to introduce sequentially the uncertainties associated with the entire sampling process as well as the spatial random effect to account for spatial autocorrelation. In addition to this, the Integrated Nested Laplace Approximation (INLA) methodology was employed, which provides a faster Bayesian inference than traditional Markov Chain Monte Carlo methods and more accurate numerical and de- terministic approximations to posterior marginal distributions (Rue et al., 2009). Additionally, the spatial component was achieved through a Stochastic Partial Differential Equations (SPDE) approach, which fits well into the INLA framework. This novel approach has been used in fisheries more frequently in the last years and is a suitable tool to deal with spatial ecology problems (Cosandey-godin et al., 2015; Muñoz et al., 2013; Paradinas et al., 2015; Quiroz et al., 2015).

To explain the model, the proxy of fish biomass (NASC) at an EDSU was assumed to have a non-Gaussian distribution. In this situation, Yi represents the NASC value plus an arbitrary value of one ( = +Y NASC 1i i ) in an EDSU i:

∼Y LogN μ σ( , )i i i 2 (1)

= = +log μ n β W( )i i i0 (2)

∼W N Q κ τ(0, ( , ) )i (3)

∼ ( )β N μ ρ,β0 0 (4) The logarithmic function links the linear predictor ni to the para-

meter =μ E Y( )i i . The linear predictor ni is important because it can capture the variability in the data caused by spatial autocorrelation. In the linear predictor (Eq. (2)), β0 is the intercept that represents the effects of the linear predictor and whose hyper-parameters are μ β0, the mean, and ρ, the precision parameter. On the other side, Wi represents the spatially structured random effect which is assumed to be a Gaus- sian Field (GF), accommodating the spatial random dependence with zero mean and a given Matérn covariance matrix Q. This matrix de- pends on the Euclidean distance between locations and the hy- perparameters κ and τ, which determinate the range of the effect

G. Moron et al. Deep-Sea Research Part II xxx (xxxx) xxx–xxx

2

( ≈φ τ8 / ) and the total variance ( ≈σ πκ τ1/(4 )w 2 2 2 ), respectively. The

matrix Q is computed internally by the SPDE approach and it represents the Gaussian Markov Random Fields (GMRF) approximation to the continuous Gaussian Field (GF) with Matérn covariance structure (Lindgren et al., 2011).

The spatial domain was defined by the construction of a Constrained Defined Delaunay triangulation (Chew, 1989), which, opposite to a regular grid, is a partition of the study area into triangles (Fig. S2). Initially, observations are treated as initial vertices for the triangula- tion, and extra vertices are added heuristically to minimize the number of triangles needed to cover the region subject. These extra vertices are used as prediction locations. The triangulation is denser in regions where there are more observations and consequently there is more in- formation, and therefore more detail is needed (Muñoz et al., 2013). One of the advantages of this approach is that it provides the posterior conditional distribution with more detail and saves computing time for the prediction of the response variable at any location (for more details see Muñoz et al., 2013).

As a result of the spatial model, the predictive mean posterior dis- tribution values of NASC (called as mNASC hereafter, which is in logarithmic scale) in the study area were calculated in each triangle vertex, and then they were linearly interpolated into a finer lattice of 1 × 1 nm grid, obtaining the final result of the model. This model was used for each survey, obtaining thus 41 anchovy spatial distribution maps. These distribution maps were the basis for the following analysis.

Each model was evaluated for goodness-of-fit computing the root mean squared estimation error (RMSEE), a common indicator used to assess the closeness between observed and estimated NASC values (Carson et al., 2017; Quiroz et al., 2015). The RMSEE is computed as follows:

∑= = −RMSEE n

d d log NASC mNASC 1

; ( ) obs

i i i i 2

where nobs is the number of NASC recorded in the survey. Clearly, smaller values are better. Since linearity is expected between the ob- served and predicted values, we also assessed the performance of the predicted values using Pearson's correlation coefficient (ρ), considering

>ρ 0.7 and p-values <0.01 as acceptable predictions (Rufener et al., 2017).

2.3. Identification of large-scale and mesoscale structures

We followed the large and mesoscale dimensions given in Bertrand et al. (2008a). The total area of extent was used as a metric of the large- scale structures (Fig. 1A). We referred to mesoscale structures as clus- ters of grids where the population is in high abundances (hereafter as ‘aggregation areas or structures’, Fig. 1B). There are many criteria to identify these clusters (e.g. Bartolino et al., 2011; Santora et al., 2011; Thakali et al., 2015), however, the rule followed in this study was simple, clustering grids with mNASC values greater than a common threshold among all surveys. This common threshold was the mean of the 95th percentiles of positive mNASC values of each survey. This simple approach has also been used in other studies to identify diversity hotspots (Williams et al., 1996).

2.4. Spatial indicators

We proposed 6 spatial indicators focused on large-scale and me- soscale structures. The following are the 3 large-scale indicators.

Presence area (PA, nm2): Indicates the total distribution area but does not take into account the level of abundance (Woillez et al., 2009).

Mean density ( +sA , m 2.nm−2): Indicates the density level (Gutiérrez

et al., 2007).

= ∑

∀ >+ =s mNASC n

mNASC, 0A i n

i i

1

Isotropy (I ): Indicates a measure of the elongation of the population in space (Woillez et al., 2009).

=I I I

min

max

where Imin and Imax are the minimal and maximal distance, respectively, of the population respecting to its center of gravity. A value near 1 means a spatial distribution composed by patches approximating a perfect circle and a value near 0 approximates an elongated ellipse.

There are also 3 mesoscale indicators. Aggregation area (, nm2): Indicates the total area occupied by ag-

gregation areas. Number of aggregation areas (N ). Center of aggregation (DC, nm): Indicates the mean distance to the

coast of aggregation areas weighted by their areas.

= ∑

∑

=

=

DC dc A

A i N

i i

i N

i

1

1

where dci is the distance to the coast of the hotspot area i and Ai is its area.

2.5. Temporal variations

We plotted the temporal series of each proposed spatial indicator. Then, a principal component analysis (PCA) was employed to sum- marize into few dimensions the variability of these indicators over surveys. After choosing an appropriate number of axes based on the variance explained, a breakpoints analysis (Bai, 1994) was performed on the time series of the scores of each axis. This allowed us to identify structural changes in the time series and thus possible changes in the aggregative behavior.

2.6. Spatial distribution scenarios

We clustered the surveys using the first axes of the previous PCA explaining more than 90% of the variance as variables. We obtained the

Fig. 1. Large- (A) and meso- (B) scale structures identified from the modeled spatial distribution. Large-scale structures are related to the total area of dis- tribution. Mesoscale structures are clusters of grids with values of mNASC greater than the common threshold.

G. Moron et al. Deep-Sea Research Part II xxx (xxxx) xxx–xxx

3

optimal number of clusters using the ‘NbClust’ package (Charrad et al., 2014) in the R software (R Development Core Team, 2011). This package obtains different optimal groups through different methodol- ogies based on Euclidean distance and k-means as the clustering method, and then suggests the more frequent number of optimal clus- ters. We considered to the optimal clusters obtained as the spatial scenarios since they shared a common spatial behavior. Then, we ob- tained an average map and a variability map from the spatial dis- tributions of each scenario. An average map is a summary of the area occupied by this species during a period (in this case a set of surveys), and it was calculated as the mean of values per grid. Similarly, a variability map shows the inter-survey variability of the anchovy spatial distribution and it was calculated as the standard deviation of values per grid (Bellier and Planque, 2007; Lelièvre et al., 2014; Saraux et al., 2014). Once obtained these two maps for each scenario, we classified areas as:

a) Recurrent - grids with high mean and low standard deviation; b) Occasional - grids with high mean and standard deviation; and c) Rare - grids with low mean and high standard deviation.

The threshold to consider a grid as ‘high’ or ‘low’ was the 85th percentile of the positive values of the average and variability map. Finally, in order to explore differences in biomass estimates for each spatial scenario, we performed an ANOVA test to look for differences among scenarios (anchovy biomass estimates of IMARPE can be found in: http://www.imarpe.pe/imarpe/archivos/informes/Biomasas_ metodo_hidroacustico_anchoveta_1985_2015_.pdf).

Analyses explained here were performed mostly using the R soft- ware (R Development Core Team, 2011), as well as ‘INLA’ package for the spatial model (Rue et al., 2014).

3. Results

3.1. Spatial distribution

In general, the model was able to reconstruct the spatial patterns of the Peruvian anchovy. Table S1 shows a numerical summary of the posterior distributions of the most important parameters for each model. RMSEE values ranged from 0.4 to 1.55 and Pearson coefficients between observed and predicted values were always greater than 0.8 and p-value ≤0.01 for all models, which confirms that predicted values follow the same pattern of observed values (Table S2).

The spatial distribution of this species showed high differences among surveys. Different types of distributions were obtained, from a constrained distribution area (e.g. Cr980809) to a wide-spread one (e.g. Cr000102) (Fig. S3). There was not a clear preference for a specific zone, being highly variable for each survey. Some surveys showed small areas with the highest abundances (e.g. Cr010304, Cr020203). The mean effective range for the structured spatial effect ranged from 10.2 nm (for Cr061112) to 25.2 nm (for Cr030203), with an average of 16.8 nm. Thus, we can interpret that the anchovy biomass for some sites depends upon its neighbor sites up to that distance. On the other side, the nominal variance ranged from 1.3 (for Cr980809) to 10 (for Cr020203), with an average value of 5.45. The mesoscale aggregation areas identified using our criterion ranged from 3.3 to 138 nm2 in average, which fit into the mesoscale classification employed.

3.2. Temporal variations

The series of presence area showed the lowest value for Cr980305 (7152 nm2), followed by Cr980809 (10,716 nm2) and Cr970910 (10,850 nm2). On the other hand, the highest value was determined for Cr130809 (43,025 nm2). After the 1997–1998 El Niño, a sudden in- crease in this indicator happened, with high values until 2005 and then lower values from 2006 to 2010 and from 2014 to 2016 (Fig. S4A). For

the case of +sA , low values were calculated for surveys in the winter and spring period; a continuous increase occurred from 1994 to 2001, fol- lowed by high values until 2013 (Fig. S4B). The last large-scale in- dicator, isotropy, did not show clear trends, attaining high values in 1994, 2000–2001 and 2012–2013, while the lowest one was de- termined for the summer of 2010. (Fig. S4C).

For mesoscale indicators, the lowest HA values were calculated in 1994–1998 and 2014–2016, while the highest ones for 1999–2004 (Fig. S4D). The number of aggregation areas (N) showed a continuous in- crease from 1994 to 2004, high values until 2013, and then a decrease (Fig. S4E). Finally, the DC indicator did not show a clear temporal trend, however, it is important to mention the low values during spring of 1997, summer 2008 and several periods after 2013 (Fig. S4F).

According to PCA results, the first two axes explained 81.8% of the total variance and we retained them to summarize the spatial indicators (Fig. 2). The first axis was dominated by N, HA and +sA indicators, therefore we will call it ‘the aggregation axis’, while the second one was represented by isotropy, the center of aggregation and the presence area, so that we will call it ‘the large-scale axis’.

The Bayesian Information Criterion (BIC) of the breakpoint analysis for the aggregation axis reported two optimal breakpoints (BIC = 162.2), dated at ~1999 and ~2013 (Fig. 3 and Table S3). The period 1994–1999 and 2013–2016 had negative scores while the period 1999–2013 had positive ones. This analysis was also performed con- sidering only scores obtained for surveys carried out during summer- autumn, thus avoiding possible seasonal effects, but we obtained the same results. For the case of the large-scale axis, no breakpoints were detected; however, the lowest values were calculated for Cr970809, Cr080204, Cr160304, and Cr160506, most of them influenced by El Niño events.

3.3. Spatial distribution scenarios

The first four axes of the PCA explained 94.5% of the variance and were used for the clustering analysis. We found four optimal number of

Fig. 2. PCA showing the relative importance of the six spatial indicators to explain the inter-survey variability. The two first axes explain 82% of the variance. The hierarchical cluster analysis of the surveys, according to spatial indicators, performed on the four first PCA dimensions shows four optimal clusters. Cluster 1: green, Cluster 2: sky-blue, Cluster 3: purple, and Cluster 4: red.

G. Moron et al. Deep-Sea Research Part II xxx (xxxx) xxx–xxx

4

clusters, which showed well-defined scenarios for the anchovy spatial distributions (Fig. 2, Fig. S5). In the first cluster (green group), we found three surveys characterized by high values for all spatial in- dicators, except Cr010304, which showed intermediate values for large- scale indicators (first quadrant). Surveys in the second cluster (sky-blue group), excepting Cr021011 and Cr081112, were carried out during summer or autumn, and they were located in the fourth quadrant. They were characterized by small presence areas, low values of isotropy and their aggregation areas placed close to the coast, but they had high fish density and larger and more mesoscale aggregation areas. The third cluster (purple group) was the ‘El Niño cluster’ since all of these surveys were carried out during moderate, strong or extraordinary El Niño events (ENFEN, 2012). These surveys were placed in the third quadrant, characterized by low values for all spatial indicators, opposite to the first cluster. The last cluster (red group) was mostly composed of sur- veys carried out during winter or spring, which were placed in the second quadrant of the PCA, therefore characterized by low values of density, smaller and fewer mesoscale structures, but large presence area. These clusters (or scenarios) were called 1, 2, 3 and 4 in the order just presented.

For Cluster 1, the average map showed that the main area where the anchovy was found was between 7 and 10°S. The variability map showed small areas with a high variability far from the coast and be- tween 7 and 9°S (Fig. 4). Recurrent areas were located in the northern zone and occupied 6184 nm2, while some rare and occasional areas were around them. In general, the preferred area for the anchovy during this scenario was placed north of 12°S.

In Cluster 2, there were six main areas with a high average abun- dance located along and near the coast north of 14°S. In the case of the variability map, the high values were not as concentrated as the average map and were placed until 90 nm from the coast. There were several recurrent areas which occupied 4621 nm2, but the main ones were located between 10° and 14°S near the coast. Moreover, there were occasional and small rare areas north of 9°S. Therefore, in this scenario, there was a slight preference for the southern part of the study area.

In Cluster 3, the average map only showed four small areas with high mean abundances. They were mostly placed between 12 and 14°S

and particularly near the coast. In the variability map, almost no high values were calculated, excepting some small areas. The combinations of these values resulted in two main types of areas: recurrent and oc- casional, both placed south of 12°S. The recurrent areas occupied 2489 nm2 and were mostly between 12 and 14°S and in the first 20 nm from the coast. The occasional areas were mostly between 14 and 16°S and near the coast, and some small rare areas were spread along the coast.

For Cluster 4, the average map showed only three main high abundances areas, while the variability map showed a high variability in the whole study area, spreading until ~110 nm from the coast. Occasional areas were very common due to the high variability, while recurrent areas occupied 4817 nm2 and they were not concentrated in a specific location. Rare areas were estimated far from the coast.

Finally, we found statistical differences between the mean biomass estimates for each scenario (ANOVA F = 10.4, p-value < 0.01, Fig. 5). The highest biomass values were estimated for Cluster 1, with an average biomass of 10.8 million tons, while the lowest ones for Cluster 3 (average biomass of 4 million tons). Furthermore, Cluster 2 had greater estimated biomass values (average biomass of 8.5 million tons) than Cluster 4 (average biomass of 6.7 million tons).

4. Discussion

The most important results were: (1) the detection of two structural changes in the time series of the mesoscale aggregative behavior in ~1999 and ~2013, (2) the identification of four spatial scenarios highly influenced by seasons and El Niño events, (3) differences in the locations of recurrent areas for each spatial scenario, and (4) differences in the biomass estimates between spatial scenarios. Unlike past studies, we modeled the whole spatial distribution of this species, allowing us to identify aggregations at specific spatial scales, the number of these aggregations, their areas, and their locations. This study expands the findings of previous investigations and, in general, the current knowl- edge of the spatial ecology of the Peruvian anchovy.

The non-Gaussian spatial model employed in this study displayed the capability to reproduce the spatial pattern of anchovy biomass during the 41 acoustic surveys at a fine spatial resolution, and the re- sults were similar to those based on geostatistical methodologies

Fig. 3. Time series of PCA scores for the first two axes. Two breakpoints are identified for the first axis (‘the aggregation axis’) in ~1999 and ~2013, identifying three periods (gray horizontal dashed lines).

Fig. 4. Maps of average, variability and type of area (red: recurrent, yellow: occasional, blue: rare) for each cluster or scenario.

G. Moron et al. Deep-Sea Research Part II xxx (xxxx) xxx–xxx

5

(Castillo et al., 2015). Quiroz et al. (2015) successfully estimated the biomass of this species using a Bayesian hurdle spatial model with the INLA methodology, accommodating zero and nonzero values as an in- tegrated process and using four geographic variables (distance to the coast, depth, latitude, and longitude) as covariates. Furthermore, other authors have also used this methodology to assess the spatial variability of fisheries resources in different ecosystems, being more common in marine ecology studies in recent years (Muñoz et al., 2013; Paradinas et al., 2015; Pennino et al., 2016). The model employed in this study is plausible to accomplish the proposed goals since predicted values follow the same spatial pattern as observed values, however, some modifications and covariates incorporation should be made if more precision is required (e.g. for biomass estimation).

4.1. Temporal variations

The first and second PCA axes were related to the aggregative and large-scale behavior of the Peruvian anchovy, respectively. The first breakpoint in the time series of the first axis was estimated in ~1999, when we found greater and more mesoscale structures. This year cor- responds to the beginning of the ‘full anchovy era’, a period when the anchovy had a complete dominance in terms of biomass and spatial occupancy of the northern HCS while the sardine (Sardinops sagax) disappeared (Gutierrez et al., 2007). Moreover, the 1994–1999 period is considered as a transition from a warm to a cold decadal regime in the Central and Eastern Pacific (Alheit and Niquen, 2004; Chavez et al., 2003). The Peruvian anchovy attained high abundances during 1999–2013, with estimates of up to 12.1 million metric tons in 2013. On average, the extension of the CCW was wider in this period in comparison to the warm regime (Swartzman et al., 2008) and the main prey of this species, euphausiids, increased its availability (Ayón et al., 2011; Espinoza and Bertrand, 2008). Feeding has been proposed to influence schooling behavior in small pelagic fish, forming more ag- gregations to feed areas with high abundances of prey (Bertrand et al., 2006; Brehmer et al., 2007). Therefore, the habitat expansion, in- creased prey availability and high abundances since 1999 might have triggered the formation of more mesoscale, and even finer, aggregation structures.

The second breakpoint identified in ~2013 coincides with the

beginning of two consecutive El Niño events in the southeastern Pacific: 2014 (moderate intensity) and 2015–2016 (strong intensity) (ENFEN, 2012, http://www.met.igp.gob.pe/elnino/lista_eventos.html). These events might be the main reason for the reduction of the number and dimensions of mesoscale aggregations during this period. The ag- gregative behavior during the last survey of 2016 (Cr160910, Fig. 3A), which was not affected by an El Niño event, supports this hypothesis since we observed a recovery in the aggregation indicators.

The variations in the aggregative behavior since 2013, a period highly influenced by consecutive El Niño events, have coincided with changes in the ecology of this species, especially diet, and might have caused a negative effect on its fishery. Regarding the ecology, the dominance of euphausiids in the diet was lower and other taxa in- creased its presence (e.g. copepods) (Espinoza, 2015). This was likely due to a more coastal habitat distribution during the last years, where euphausiids are less abundant (Ayón et al., 2011). Moreover, this might be heightened by changes of El Niño on the zooplankton community, which favors the presence of copepods during these events (Ayón et al., 2008). Regarding the fishery, the fishing efficiency had a decreasing tendency since 2013, and a slight recovery for the second semester of 2016 was observed (IMARPE, 2017), which agrees with the recovery of mesoscale structures. More studies about the consequences of the last El Niño events (2014–2016) on the HCS should be done to understand their all effects on the ecology of this species.

4.2. Spatial distribution scenarios

We found four spatial scenarios principally influenced by seasons and El Niño events. In a previous study, Joo et al. (2014) identified four ecosystem scenarios for the environment and the spatial behavior of the Peruvian anchovy, and our results are highly consistent with their findings. We refer to scenarios published in Joo et al. (2014) in Roman numerals hereafter.

4.2.1. Summer average scenario Scenario II is similar to Scenario 2 of this study since they share

spatial distributions occurring in summer seasons. These similarities allow us to deduce some environmental features. Scenario II was characterized by relatively high SST, high primary production, and shallow oxycline. Moreover, there is a reduction of the CCW extension and the number of fine-scale and sub-mesoscale physical structures are fewer during summer (Grados et al., 2016; Swartzman et al., 2008). Consequently, the anchovy preferred habitat is more constrained, and, for this reason, individuals tend to aggregate in areas closer to the coast, thus forming more and greater mesoscale structures as detected for this scenario.

The recurrent areas for Scenario 2 were located close to the coast and compactly formed. These areas had a slight preference for the southern part of the study area, possibly due to the intrusion of warmer waters from the equatorial and oceanic zone during summer in the northern zone, which might be the main reason of a higher presence of occasional areas in those locations (Checkley et al., 2008).

4.2.2. Winter and spring scenario Our Scenario 4 is similar to Scenario IV since both share the same

surveys, mostly carried out during winter or spring. This scenario is characterized by relatively low SST, low primary production and a deep oxycline (Joo et al., 2014). Furthermore, the CCW extension increases and there are more fine-scale and sub-mesoscale physical structures (Grados et al., 2016; Swartzman et al., 2008). Thus, during this sce- nario, the preferred habitat for the Peruvian anchovy is wider and the stock tends to occupy the whole area. For this reason, despite the ex- istence of more physical structures, the formation of mesoscale ag- gregation structures of anchovies is less, likely due to the tendency of individuals to avoid negative density dependence effects and exploit most of the available resources (Barra et al., 2015). This is a clear

Fig. 5. Differences in official estimated biomass values per cluster (ANOVA F=10.4, p-value < 0.01).

G. Moron et al. Deep-Sea Research Part II xxx (xxxx) xxx–xxx

6

example that more physical structures do not necessarily translate into more aggregation structures.

For Scenario 4, we detected almost no recurrent areas and the dominant area was the ‘occasional’ type. This high variability in the spatial distribution for winter and spring seasons might be influenced by the wide extension of the CCW (Swartzman et al., 2008), where it is less probable to find fidelity areas for the presence of anchovies (Castillo et al., 2015).

The estimated biomass values for Scenario 2 were slightly higher than the Scenario 4. Despite a wider extension of the favorable habitat during spring seasons, the aggregative behavior observed during summer seems to play a critical role in the current process of biomass estimation. Some authors have discussed possible effects of the wider dispersion during winter-spring surveys on the estimated biomass, which can lead to an underestimation due to a less anchovy availability for scientific vessels (Gutierrez et al., 2007). Therefore, it is possible that the real biomass during winter or spring could be higher than es- timated.

4.2.3. El Niño scenario The characteristics of Scenario 3 have not been identified in pre-

vious studies. Its identification here is possibly due to a longer time series and novel approach for the identification of aggregation struc- tures used in this research. This scenario is highly influenced by mod- erate, strong and extreme El Niño events, which negatively affect the spatial behavior at the two scales examined in this study. During these events the environment is characterized by high SST, lower primary productivity, a deep oxycline, and the extension of the CCW is ex- tremely reduced (Bertrand et al., 2011, 2004; Fiedler, 2002; Swartzman et al., 2008), therefore reducing the anchovy preferred habitat. It has been suggested that the number of sub-mesoscale physical structures is positively associated with the upwelling intensity (Grados et al., 2016; Nieto et al., 2012), therefore the number and dimensions of these structures would be more and greater during large extensions of the CCW and the opposite during its reduction. Thus, the number and di- mensions of physical structures might be fewer and smaller during El Niño events due to the drastic reduction of cold upwelled waters (Espinoza-Morriberon et al., 2017; Grados et al., 2016). These factors, combined with the lower presence of individuals due to a elevated natural mortality for this species during El Niño events (Pauly and Tsukayama, 1987), might induce a reduction in the number and di- mensions of aggregation structures.

The recurrent area identified was the most constrained of the four scenarios and was located between 12° and 14°S, near the coast. Some authors have hypothesized a high ability of this species to find refuges zones, or ‘loopholes’, during unfavorable events as El Niño (Bakun and Broad, 2003; Bertrand et al., 2004). For instance, it has been proposed that during El Niño 1997–1998, the distribution area of this species was limited to the area of the reduced CCW (Bertrand et al., 2004; Swartzman et al., 2008), where upwelled nutrient-rich waters were still observed (Bertrand et al., 2004). However, no area containing these loopholes has been localized until now. The recurrent area identified is consistent with a refuge zone and might be critical for the survival of the stock during these events; therefore, its protection might be ne- cessary.

The lowest biomass values have been estimated for this scenario, confirming the negative effects of El Niño on the stock (Pauly and Tsukayama, 1987). On the other hand, a downward bias in the biomass estimates has also been discussed, this due to the constrained and very coastal distribution of anchovies during these events, where scientific vessels can hardly access (Gutierrez, 2001). One fact that supports this hypothesis is the absence of adult individuals from August to September 1998 and their reappearance in catches in late 1998 and early 1999 (Bertrand et al., 2004). Further studies have to be carried out to un- derstand the main cause of these low abundances during these warm events: underestimation, a high natural mortality (Gutierrez, 2001), or

both.

4.2.4. Summer favorable scenario Scenario 1 is composed of three spatial distributions inferred from

surveys carried out during summer seasons and might be similar to the Scenario I since both have the Cr010304 as common survey. Scenario I was characterized by high SST, productivity, and a shallow oxycline (Joo et al., 2014). The main recurrent area was located in the northern zone, where one of the most important upwelling areas occurs (Gutiérrez et al., 2016). Moreover, it is possible that during these sur- veys the intrusion of warm waters in the study area has been scarce.

Furthermore, the highest biomass values were estimated during this scenario. This might be a consequence of the wide expansion of the area of presence and the existence of more and larger aggregation areas. These results are in agreement with the basin model, where the highest biomass estimates were associated with a large presence area and higher density (MacCall, 1990). This spatial behavior is known for this species and other small pelagic fish in the Mediterranean Sea (Barange et al., 2009; Barra et al., 2015; Saraux et al., 2014), which expand their distribution range with increasing biomass and limit their distribution to particular locations when biomass decreases (e.g. El Niño scenario in this study) (Checkley et al., 2008).

5. Conclusions

Using a novel approach to identify aggregation structures at two different spatial scales, we analyzed temporal and spatial variations of the aggregative behavior of the Peruvian anchovy. We identified four spatial scenarios, highly influenced by seasons and El Niño events. Each scenario differed in the aggregative behavior and locations of the re- current areas. For the El Niño scenario, the recurrent area was detected nearshore along the Peruvian coast (12–14°S). This area might be cri- tical for the survival of the stock during these unfavorable events and its protection might be required. We also detected two temporal changes in the formation of mesoscale aggregation structures. The first break- point was in 1999, influenced by a new decadal regime in the HCS. The second change was in 2013, influenced by consecutive El Niño events from 2014 to 2016. Finally, we detected differences in the biomass estimates for each scenario. The lowest biomass values were estimated for the El Niño scenario, characterized by a nearshore presence, and smaller and fewer mesoscale aggregation structures. On the other hand, the highest biomass estimates were calculated for the summer favorable scenario, characterized by a wider presence and more and larger me- soscale aggregation structures. These differences support the basin model proposed for this species in previous studies (Barange et al., 2009). Further investigation of the effects and causes of the variations of the anchovy aggregative behavior is needed, especially for the most recent years, in order to elucidate the biological and ecological impacts of these changes.

Acknowledgements

We would like to express our gratitude to Daniel Grados for valuable explanations about physical structures and to the two anonymous re- viewers for their comments to greatly improve the writing of this ar- ticle. We are grateful to IMARPE for the data for this work. This work was supported in part by the cooperation agreement among the Instituto del Mar del Peru (IMARPE) and the Institut de Recherche pour le Développement (IRD), France: International Join Laboratory - Dynamics of the Humboldt Current System (LMI-DISCOH) to one of the authors.

Appendix A. Supporting information

Supplementary data associated with this article can be found in the online version at doi:10.1016/j.dsr2.2018.11.009.

G. Moron et al. Deep-Sea Research Part II xxx (xxxx) xxx–xxx

7

References

Alheit, J., Niquen, M., 2004. Regime shifts in the Humboldt Current ecosystem. Prog. Oceanogr. 60, 201–222. https://doi.org/10.1016/j.pocean.2004.02.006.

Ayón, P., Criales-Hernandez, M.I., Schwamborn, R., Hirche, H.J., 2008. Zooplankton research off Peru: a review. Prog. Oceanogr. 79, 238–255. https://doi.org/10.1016/j. pocean.2008.10.020.

Ayón, P., Swartzman, G., Espinoza, P., Bertrand, A., 2011. Long-term changes in zoo- plankton size distribution in the Peruvian Humboldt Current System: conditions fa- vouring sardine or anchovy. Mar. Ecol. Prog. Ser. 422, 211–222. https://doi.org/10. 3354/meps08918.

Bai, J., 1994. Least squares estimation of a shift in linear processes. J. Time Ser. Anal. 15, 453–472. https://doi.org/10.1111/j.1467-9892.1994.tb00204.x.

Bakun, A., 1996. Patterns in the Oceans: Ocean Processes and Marine Population. California Sea Grant College System, National Oceanic and Atmospheric Administration, in cooperation with Centro de Investigaciones Biológicas del Noroeste, California, USA.

Bakun, A., Broad, K., 2003. Environmental ‘loopholes' and fish population dynamics: comparative pattern recognition with focus on El Niño effects in the Pacific. Fish. Oceanogr. 12, 458–473.

Barange, M., Coetzee, J., Takasuka, A., Hill, K., Gutierrez, M., Oozeki, Y., Lingen, C., van der, Agostini, V., 2009. Habitat expansion and contraction in anchovy and sardine populations. Prog. Oceanogr. 83, 251–260. https://doi.org/10.1016/j.pocean.2009. 07.027.

Barange, M., Coetzee, J.C., Twatwa, N.M., 2005. Strategies of space occupation by an- chovy and sardine in the southern Benguela: the role of stock size and intra-species competition. ICES J. Mar. Sci. 62, 645–654. https://doi.org/10.1016/j.icesjms.2004. 12.019.

Barra, M., Petitgas, P., Bonanno, A., Somarakis, S., Woillez, M., Machias, A., Mazzola, S., Basilone, G., Giannoulaki, M., 2015. Interannual changes in biomass affect the spatial aggregations of anchovy and sardine as evidenced by Geostatistical and spatial in- dicators. PLoS One 10. https://doi.org/10.1371/journal.pone.0135808.

Bartolino, V., Maiorano, L., Colloca, F., 2011. A frequency distribution approach to hotspot identification. Popul. Ecol. 53, 351–359. https://doi.org/10.1007/s10144- 010-0229-2.

Bellier, E., Planque, B., 2007. Historical fluctuations in spawning location of anchovy ( Engraulis encrasicolus) and sardine ( Sardina pilchardus) in the Bay of Biscay during 1967 – 73 and 2000 – 2004. Fish. Oceanogr. 16, 1–15. https://doi.org/10.1111/j. 1365-2419.2006.00410.x.

Benerjee, S., Carlin, B., Gelfand, A., 2004. Hierarchical Modeling and Analysis for Spatial Data. Chapman and Hall/CRC, Boca Raton, FL.

Bertrand, A., Barbieri, M.A., Gerlotto, F., Leiva, F., Córdova, J., 2006. Determinism and plasticity of fish schooling behaviour as exemplified by the South Pacific jack mackerel Trachurus murphyi. Mar. Ecol. Prog. Ser. 311, 145–156. https://doi.org/ 10.3354/meps311145.

Bertrand, A., Chaigneau, A., Peraltilla, S., Ledesma, J., Graco, M., Monetti, F., Chavez, F.P., 2011. Oxygen: a fundamental property regulating pelagic ecosystem structure in the coastal southeastern tropical pacific. PLoS One 6. https://doi.org/10.1371/ journal.pone.0029558.

Bertrand, A., Gerlotto, F., Bertrand, S., Gutiérrez, M., Alza, L., Chipollini, A., Díaz, E., Espinoza, P., Ledesma, J., Quesquén, R., Peraltilla, S., Chavez, F., 2008a. Schooling behaviour and environmental forcing in relation to anchoveta distribution: an ana- lysis across multiple spatial scales. Prog. Oceanogr. 79, 264–277. https://doi.org/10. 1016/j.pocean.2008.10.018.

Bertrand, A., Grados, D., Colas, F., Bertrand, S., Capet, X., Chaigneau, A., Vargas, G., Mousseigne, A., Fablet, R., 2014. Broad impacts of fine-scale dynamics on seascape structure from zooplankton to seabirds. Nat. Commun. 5, 1–9. https://doi.org/10. 1038/ncomms6239.

Bertrand, A., Segura, M., Gutiérrez, M., Vásquez, L., 2004. From small-scale habitat loopholes to decadal cycles: a habitat-based hypothesis explaining fluctuation in pelagic fish populations off Peru. Fish. Fish. 5, 296–316. https://doi.org/10.1111/j. 1467-2679.2004.00165.x.

Bertrand, S., Dewitte, B., Tam, J., Díaz, E., Bertrand, A., 2008b. Impacts of Kelvin wave forcing in the Peru Humboldt Current system: scenarios of spatial reorganizations from physics to fishers. Prog. Oceanogr. 79, 278–289. https://doi.org/10.1016/j. pocean.2008.10.017.

Brehmer, P., Gerlotto, F., Laurent, C., Cotel, P., Achury, A., Samb, B., 2007. Schooling behaviour of small pelagic fish: phenotypic expression of independent stimuli. Mar. Ecol. Prog. Ser. 334, 263–272. https://doi.org/10.3354/meps334263.

Carson, S., Shackell, N., Mills Flemming, J., 2017. Local overfishing may be avoided by examining parameters of a spatio-temporal model. PLoS One 12, 1–21. https://doi. org/10.1371/journal.pone.0184427.

Castillo, P.R., Madureira, L., Marangoni, J., Gerlotto, F., Guevara-carrasco, R., 2015. Variability in distribution and aggregation behavior of the Peruvian Anchovy (Engraulis ringens) analyzed using a Fifteen Year Long Series of Acoustic Surveys (2000–2014). RIO Acoust. 2015 (1), 1–9.

Castillo, R., Peraltilla, S., Aliaga, A., Flores, M., Ballon, M., Calderon, J., Gutierrez, M., 2009. Protocolo técnico para la evaluación acústica de las áreas de distribución y abundancia de recursos pelágicos en el mar peruano. Versión 2009. Callao.

Charrad, M., Ghazzali, N., Boiteau, V., Niknafs, A., 2014. NbClust: an R Package for Determining the. J. Stat. Softw. 61, 1–36.

Chavez, F.P., Bertrand, A., Guevara-Carrasco, R., Soler, P., Csirke, J., 2008. The northern Humboldt Current System: brief history, present status and a view towards the future. Prog. Oceanogr. 79, 95–105. https://doi.org/10.1016/j.pocean.2008.10.012.

Chavez, F.P., Ryan, J., Lluch-Cota, S.E., Niquen, C., M., 2003. From anchovies to sardines

and back: multidecadal change in the Pacific Ocean. Science 299, 217–221. https:// doi.org/10.1126/science.1075880.

Checkley, D.M., Alheit, J., Oozeki, Y., Roy, C., 2008. Climate Change and Small Pelagic Fish, 1st ed. Cambridge University Press, Cambridge.

Chew, L.P., 1989. Constrained Delaunay Triangulations. Algorithmica 4, 97–108. https:// doi.org/10.1007/BF01553881.

Cosandey-godin, A., Krainski, E.T., Worm, B., Flemming, J.M., 2015. Applying Bayesian spatiotemporal models to fisheries bycatch in the Canadian Arctic. Can. J. Fish. Aquat. Sci. 72, 1–12. https://doi.org/10.1139/cjfas-2014-0159.

Ekau, W., Auel, H., P̈ortner, H.O., Gilbert, D., 2010. Impacts of hypoxia on the structure and processes in pelagic communities (zooplankton, macro-invertebrates and fish). Biogeosciences 7, 1669–1699. https://doi.org/10.5194/bg-7-1669-2010.

ENFEN, 2012. Definición operacional de los eventos El Niño y La Niña y sus magnitudes en la costa del Perú. Lima.

Espinoza-Morriberon, D., Echevin, V., Colas, F., Tam, J., Ledesma, J., Vasquez, L., Graco, M., 2017. Impacts of El Niño events on the Peruvian upwelling system productivity. J. Geophys. Res. Ocean. 122, 1–22. https://doi.org/10.1002/2016JC012439.

Espinoza, P., 2015. Breve descripción de la dieta de la Anchoveta procedente de los cruceros de investigación realizados por el IMARPE entre verano 1996 y verano 2015. Callao, Peru.

Espinoza, P., Bertrand, A., 2008. Revisiting Peruvian anchovy (Engraulis ringens) tro- phodynamics provides a new vision of the Humboldt Current system. Prog. Oceanogr. 79, 215–227. https://doi.org/10.1016/j.pocean.2008.10.022.

Fiedler, P.C., 2002. Environmental change in the eastern tropical Pacific Ocean: review of ENSO and decadal variability. Mar. Ecol. Prog. Ser. 244, 265–283. https://doi.org/ 10.3354/meps244265.

Fréon, P., Arístegui, J., Bertrand, A., Crawford, R.J.M., Field, J.C., Gibbons, M.J., Tam, J., Hutchings, L., Masski, H., Mullon, C., Ramdani, M., Seret, B., Simier, M., 2009. Functional group biodiversity in Eastern Boundary Upwelling Ecosystems questions the wasp-waist trophic structure. Prog. Oceanogr. 83, 97–106. https://doi.org/10. 1016/j.pocean.2009.07.034.

Gelfand Jr., A.E., Silander, J.A., Wuz, S., Latimerx, A., Lewis, P.O., Rebelok, A.G., Holder, M., 2006. Explaining species distribution patterns through hierarchical modeling. Bayesian Anal. 1, 41–92.

Grados, D., Bertrand, A., Colas, F., Echevin, V., Chaigneau, A., Gutierrez, D., Vargas, G., Fablet, R., 2016. Spatial and seasonal patterns of fine-scale to mesoscale upper ocean dynamics in an Eastern Boundary Current System. Prog. Oceanogr. 142, 105–116. https://doi.org/10.1016/j.pocean.2016.02.002.

Gutierrez, D., Akester, M., Naranjo, L., 2016. Productivity and sustainable management of the Humboldt current large marine ecosystem under climate change. Environ. Dev. 17, 126–144. https://doi.org/10.1016/j.envdev.2015.11.004.

Gutiérrez, D., Akester, M., Naranjo, L., 2016. Productivity and sustainable management of the Humboldt current large marine ecosystem under climate change. Environ. Dev. 17, 126–144. https://doi.org/10.1016/j.envdev.2015.11.004.

Gutierrez, M., 2001. Efectos del evento El Niño 1997-98 sobre la distribución y abun- dancia de anchoveta (Engraulis ringens). El Niño en América Lat. impactos biológicos y Soc.

Gutierrez, M., Swartzman, G., Bertrand, A., Bertrand, S., 2007. Anchovy ( Engraulis ringens) and sardine ( Sardinops sagax) spatial dynamics and aggregation patterns in the Humboldt Current ecosystem, Peru, from 1983 – 2003. Fish. Oceanogr. 99, 1–14. https://doi.org/10.1111/j.1365-2419.2006.00422.x.

Gutiérrez, M., Swartzman, G., Bertrand, A., Bertrand, S., 2007. Anchovy (Engraulis ringens) and sardine (Sardinops sagax) spatial dynamics and aggregation patterns in the Humboldt Current ecosystem, Peru, from 1983–2003. Fish. Oceanogr. 16, 155–168. https://doi.org/10.1111/j.1365-2419.2006.00422.x.

Hixon, M.A., Anderson, T., Buch, K., Johnson, D., McLeod, J.B., Stallings, C., 2012. Density dependence and population regulation in marine fish: a large-scale, long- term field manipulation. Ecol. Monogr. 82, 467–489.

IMARPE, 2017. Situación del stock norte-centro de la Anchoveta peruana (Engraulis ringens) al 01 de noviembre de 2017 y perspectivas de explotación para la segunda temporada de pesca del 2017. Lima.

IMARPE, 2016. Elaboración de la Tabla de Decisión para la determinación del Límite Máximo de Captura Total Permisible para la pesquería del Stock Norte-Centro de la anchoveta peruana. Callao, Peru.

Joo, R., Bertrand, A., Bouchon, M., Chaigneau, A., Demarcq, H., Tam, J., Simier, M., Gutiérrez, D., Gutiérrez, M., Segura, M., Fablet, R., Bertrand, S., 2014. Ecosystem scenarios shape fishermen spatial behavior. The case of the Peruvian anchovy fishery in the Northern Humboldt Current System. Prog. Oceanogr. 128, 60–73. https://doi. org/10.1016/j.pocean.2014.08.009.

Lelièvre, S., Vaz, S., Martin, C.S., Loots, C., 2014. Delineating recurrent fi sh spawning habitats in the North Sea. J. Sea Res. 91, 1–14. https://doi.org/10.1016/j.seares. 2014.03.008.

Lindgren, F., Rue, H., Lindstrom, J., 2011. An explicit link between Gaussian fields and Gaussian Markov random fields. J. R. Stat. Soc. Ser. B 73, 423–498.

MacCall, A.D., 1990. Dynamic Geography of Marine Fish Populations, 1st ed. University of Washington Press, Seattle and London.

Muñoz, F., Pennino, M.G., Consea, D., López-Quílez, A., Bellido, J.M., 2013. Estimating and prediction of the spatial occurrence of fish species using Bayesian latente Gaussian models. Stoch. Environ. Res. Risk Assess. 27, 1171–1180. https://doi.org/ 10.1007/s00477-012-0652-3.

Nieto, K., Demarcq, H., McClatchie, S., 2012. Mesoscale frontal structures in the Canary Upwelling System: new front and filament detection algorithms applied to spatial and temporal patterns. Remote Sens. Environ. 123, 339–346. https://doi.org/10.1016/j. rse.2012.03.028.

Palomares, M.L., Muck, P., Mendo, J., Chuman, E., Gomez, O., Pauly, D., 1987. Growth of the Peruvian Anchoveta (Engraulis ringens), 1953 to 1982. In: Pauly, D., Tsukayama,

G. Moron et al. Deep-Sea Research Part II xxx (xxxx) xxx–xxx

8

I. (Eds.), The Peruvian Anchoveta and Its Upwelling Ecosystem: Three Decades of Change 15. International Center for Living Aquatic Resources Management ICLARM Studies and Reviews, Manila, pp. 117–144.

Paradinas, I., Conesa, D., Pennino, M., Muñoz, F., Fernandez, A., Lopez-Quilez, A., Bellido, J.M., 2015. Bayesian spatio-temporal approach to identifying fish nurseries by validating persistence areas. Mar. Ecol. Prog. Ser. 528, 245–255. https://doi.org/ 10.3354/meps11281.

Pauly, D., Tsukayama, L., 1987. The Peruvian anchovetta and its upwelling ecosystem: three decades of change, Instituto del Mar del Peru (IMARPE), Callao, Peru; Deutsche GesellIschaft fdr TechnLsche Zusammenarbeit (G72), GmbH, Eschbom, Federal Republic of Germany; and International Cente. ICLARM Studies and Reviews, Callao, Peru.

Pennino, M.G., Conesa, D., Lopez-Quilez, A., Muñoz, F., Fernandez, A., Bellido, J.M., 2016. Fishery-dependent and -independent data lead to consistent estimations of essential habitats. ICES J. Mar. Sci. 73, 2302–2310. https://doi.org/10.1093/ icesjms/fst048.

Perry, J.N., Liebhold, A.M., Rosenberg, M.S., Dungan, J., Miriti, M., Jakomulska, A., Citron-Pousty, S., 2002. Illustrations and guidelines for selecting statistical methods for quantifying spatial pattern in ecological data. Ecography (Cop.) 25, 578–600.

Petitgas, P., Alheit, J., Peck, M.A., Raab, K., Irigoien, X., Huret, M., Kooij, J., Van Der, Pohlmann, T., Wagner, C., Zarraonaindia, I., Dickey-collas, M., 2012. Anchovy po- pulation expansion in the North Sea. Mar. Ecol. Prog. Ser. 444, 1–13. https://doi.org/ 10.3354/meps09451.

Quiroz, Z.C., Prates, M.O., Rue, H., 2015. A Bayesian Approach to Estimate the Biomass of Anchovies Off the Coast of Peru. Biometrics 71, 208–217. https://doi.org/10.1111/ biom.12227.

R Development Core Team, 2011. R: a language and environment for statistical com- puting.

Rue, H., Martino, S., Chopin, N., 2009. Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. J. R. Stat. Soc. Ser. B 71, 319–392.

Rue, H., Martino, S., Simpson, D., Riebler, A., Teixeira-Krainski, E., 2014. INLA: Functions

which allow to perform full Bayesian analysis of latent Gaussian models using Integrated Nested Laplace Approximation.

Rufener, M.C., Kinas, P.G., Nóbrega, M.F., Lins Oliveira, J.E., 2017. Bayesian spatial predictive models for data-poor fisheries. Ecol. Modell. 348, 125–134. https://doi. org/10.1016/j.ecolmodel.2017.01.022.

Santora, J.A., Sydeman, W.J., Schroeder, I.D., Wells, B.K., Field, J.C., 2011. Mesoscale structure and oceanographic determinants of krill hotspots in the California Current: implications for trophic transfer and conservation. Prog. Oceanogr. 91, 397–409. https://doi.org/10.1016/j.pocean.2011.04.002.

Saraux, C., Fromentin, J.-M., Bigot, J., Bourdeix, J.-H., Morfin, M., Roos, D., Beveren, E. Van, Bez, N., 2014. Spatial structure and distribution of small pelagic fish in the Northwestern Mediterranean Sea. PLoS One 9, 1–12. https://doi.org/10.1371/ journal.pone.0111211.

Simmonds, E.J., Gutiérrez, M., Chipollini, A., Gerlotto, F., Woillez, M., Bertrand, A., 2009. Optimizing the design of acoustic surveys of Peruvian anchoveta. ICES J. Mar. Sci. 66, 1341–1348. https://doi.org/10.1093/icesjms/fsp118.

Simmonds, J., MacLennan, D., 2005. Fisheries Acoustics: Theory and Practice, 2nd ed. Blackwell Science Ltd, Great Britain.

Swartzman, G., Bertrand, A., Gutiérrez, M., Bertrand, S., Vasquez, L., 2008. The re- lationship of anchovy and sardine to water masses in the Peruvian Humboldt Current System from 1983 to 2005. Prog. Oceanogr. 79, 228–237. https://doi.org/10.1016/j. pocean.2008.10.021.

Thakali, L., Kwon, T.J., Fu, L., 2015. Identification of crash hotspots using kernel density estimation and kriging methods: a comparison. J. Mod. Transp. 23, 93–106. https:// doi.org/10.1007/s40534-015-0068-0.

Williams, P., Gibbons, D., Margules, C., Rebelo, A., Humphries, C., Pressey, R., 1996. A comparison of richness hotspots, rarity hotspots, and complementary areas for con- serving diversity of British birds. Conserv. Biol. 10, 155–174.

Woillez, M., Rivoirard, J., Petitgas, P., 2009. Notes on survey-based spatial indicators for monitoring fish populations. Aquat. Living Resour. 22, 155–164. https://doi.org/10. 1051/alr/2009017.

G. Moron et al. Deep-Sea Research Part II xxx (xxxx) xxx–xxx

9