Ontogenetic and intraspecific variability in otolith shape of anchoveta (Engraulis ringens) used to identify demographic units in the Pacific Southeast off Chile

Additional keywords: elliptical Fourier analysis, generalised additive model otolith shape indices, stock.

A very complex issue in fisheries management is the uncertainty in the number, distribution and connectivity of demographic units stocks (Begg et al. 1999; Taylor and Dizon 1999; Luck et al. 2003; Kerr et al. 2014). In practice, it is not possible to estimate the productive surplus that can be harvested with acceptable precision and accuracy without knowing the effective distribution, abundance, degree of reproductive isolation, self-recruitment potential and dependence on immigration of a management unit (Secor 2013). To date, although control measures have been taken to achieve stock recovery in some overexploited fisheries (Hutchings et al. 2010; Murawski 2010; Petitgas et al. 2010), their poor success could be linked to a poor match between management units and biological stocks (i.e. demographic units).

Otolith morphology, and shape analysis in particular, have been widely used in fisheries science to determine demographic units (Campana and Casselman 1993; Begg and Brown 2000; Torres et al. 2000; Pothin et al. 2006; Canas et al. 2012). Otoliths are pairs of calcified structures, located in the inner ear of fishes, whose main functions are balance and hearing (Campana 1999). Although otolith shape is species specific (Hecht and Appelbaum 1982; Gaemers 1984), environmental variability (biotic and abiotic) may lead to intraspecific differences among geographical regions, probably mediated through differences in individual growth rates (Vignon 2012). The use of otolith shape analysis to differentiate demographic units is possible and adequate when enough and prolonged spatial isolation leads to detectable and consistent differences in otolith shape among demographic units (Neilson et al. 1985; Bird et al. 1986; Campana and Casselman 1993; Begg and Brown 2000; Turan 2000; Turan et al. 2006).

Otolith shape analysis includes shape indices (circularity, rectangularity, ellipticity, eccentricity, roundness), shape factors (Pothin et al. 2006) and otolith outline analyses, commonly based on elliptic Fourier and Wavelet transformations (Bird et al. 1986; Smith 1992; Campana and Casselman 1993; Begg and Brown 2000; Bergenius et al. 2005; Parisi-Baradad et al. 2005; Pothin et al. 2006; Turan et al. 2006). These approaches have benefited from increasing technological ability to produce reliable two-dimensional digital images in most fishery laboratories. Furthermore, otolith shape analysis is far less time consuming and less costly than other techniques used for stock discrimination, such as otolith microchemistry and otolith microstructure analysis (Chittaro et al. 2006). Consequently, otolith shape analysis can be used either complementary to or as an effective alternative for other techniques, particularly for long-term monitoring of the degree of separation or mixing of units previously established using multiple approaches.

Anchoveta (Engraulis ringens) is distributed from northern Peru (4°30′S) to southern Chile (42°30′S). Two main population units have been identified within this range: one off the northern to central coast of Peru and the other off southern Peru and northern Chile (Pauly and Tsukayama 1987). The northern unit supports important fishery activity throughout its distribution (i.e. from 16°00′S to 24°00′S), representing a substantial fraction of the global anchoveta fisheries (Bakun and Weeks 2008; Schreiber and Halliday 2013). Despite the tremendous ecological role anchoveta plays in the Humboldt Current system (Chavez et al. 2003; Espinoza and Bertrand 2008; Karstensen and Ulloa 2009) and its considerable contribution to regional economies from the 1950s to the present, there have been limited efforts to enhance our knowledge regarding the population structure of anchoveta (Galleguillos et al. 1996; Ferrada et al. 2002; Chávez et al. 2007; Valdivia et al. 2007; Rojas 2011; George-Nascimento and Moscoso 2013) and the degree of coherence between management and demographic units.

Until now, efforts to identify the number, distribution and connectivity of demographic units within Chilean waters have been limited. Genetic studies have failed to reject the null hypothesis of panmixia and have suggested the existence of a single evolutionary unit along the Chilean (Galleguillos et al. 1996; Ferrada et al. 2002) and Peruvian (Rojas 2011) coasts. Conversely, parasitological studies (Valdivia et al. 2007; George-Nascimento and Moscoso 2013) have indicated the possible existence of at least two separate demographic units, one in the north and another in southern Chile. Therefore, the aim of the present study was to investigate the spatial, temporal and ontogenetic variability in the shape of sagittal otoliths of anchoveta to contribute to efforts to reveal the population structure of this important fisheries resource.

The study area was divided in three sampling zones that corresponded to the three main fishing areas off Chile: Zone I, Arica to Antofagasta (18°24′–26°00′S); Zone II, Caldera to Coquimbo (26°01′–32°16′S); and Zone III, Valparaiso to Valdivia (32°17′–41°77′S; Fig. 1).

Fig. 1.  The anchoveta (Engraulis ringens) study area along the Chilean coast. Symbols show the fishing hauls in each year, and the dashed lines demarcate the study zones.

However, there are also oceanographic and topographic traits to justify this sampling design. First, there is a coastal transition zone (CTZ) between 34 and 39°S (Morales et al. 2010), which can act as natural boundary. Although there is no strict topographic boundary between the northern zones (I and II), there is a critical latitude described around Mejillones and Moreno bays (~22.8°S; Letelier et al. 2012). Furthermore, each of these zones has discrete catch and spawning areas, where catches and biological activity (e.g. spawning) are almost non-existent.

The anchoveta fishery in Zone I is concentrated from Arica (18°S) to Antofagasta (26°S; Böhm et al. 2018). The spatial distribution of spawning in this area is concentrated in two recurring foci, one from Arica (18°24′S) to Punta Patache (21°S), with another secondary focus located between Copaca (22°20′S) and Mejillones (23°S; Angulo et al. 2018). Historically, landings have been lower in Zone II (Caldera–Coquimbo) than Zone I, concentrated primarily in Caldera (26–28°S) and Coquimbo (29–30°S; Böhm et al. 2018). In these areas, spawning is distributed from Chañaral (16°10′S) to Tongoy (30°10′S; Reyes et al. 2018). In Zone III, anchoveta are distributed from 36 to 40°S (Aranis et al. 2018), and spawning occurs primarily between 35 and 40°S (Cubillos et al. 2017).

In all, 940 juveniles and adult anchoveta were collected between December 2015 and October 2016 in each of the three sampling zones (Table 1). Juvenile fish corresponded to specimens <12-cm total length (TL), sampled during the main recruitment period (December–February) from annual scientific surveys conducted during this period by the Chilean Instituto de Fomento Pesquero (IFOP). Samples of adult fish corresponded to individuals ≥12 cm TL, sampled during the main spawning period (August–October) from annual scientific surveys conducted during this period by IFOP. A complementary sample of adults was available from the 2013 annual reproductive survey and was included in the present study for an interannual comparison between 2016, a year with a strong El Niño–Southern Oscillation (ENSO), leading to an El Niño event, and a normal year (2013) without an El Niño event. TL to the nearest centimetre and total weight (g) were measured and sagittal otoliths were extracted under a stereomicroscope. Otoliths were cleaned in distilled water just after extraction, with any remaining tissue removed from the macula and vestibule using fine tweezers (Secor et al. 1992). Once cleaned, the otolith pairs were dried and kept dry in Eppendorf tubes.

Table 1.  Number (n) of sagittal otoliths of the anchoveta (Engraulis ringens) used in morphological analysis, classified by ontogenetic stages, area and locality of origin off Chile, for the years 2016 and 2013
The range of fish length is included for reference. TL, total length

Digital images of each left otolith were acquired using Image-Pro Plus software (Media Cybernetics, Bethesda, MD, USA) which also binarised the images and calculated otolith surface area, perimeter, Feret length (longest calliper length of the otolith) and Feret width (smallest calliper length of otolith). These morphometric variables were corrected to avoid any effects of fish size and the allometry that occurs during otolith growth (Lombarte and Lleonart 1993) using the method proposed by Lleonart et al. (2000), whereby the morphometric measures, used as dependent variables, were related to fish TL as an independent variable and then used to calculate five shape indices: form factor, roundness, circularity, rectangularity and ellipticity (Pothin et al. 2006). Finally, ~1000 Cartesian coordinates of each otolith outline were extracted using the ‘getaoipoints’ macro in Image-Pro Plus and were exported to a text file for further elliptic Fourier analysis.

Descriptive analyses of shape indices were organised into six length-interval groups linked to ontogenetic stages: two pre-recruit intervals of 4.0–5.5 and 6.0–7.5 cm TL; two recruit intervals of 8.0–9.5 and 10.0–11.0 cm TL; and two adults intervals, one from 2016 (13.5–14.0 cm TL) and another from 2013 (14.0–14.5 cm TL).

The relationship between shape indices and both sampling zone and fish length (as a covariate) was modelled and analysed through a generalised additive model (GAM) using the R package ‘mgcv’ (see https://cran.r-project.org/web/packages/mgcv; Wood 2006). Both Gaussian and gamma distribution models were compared using the Akaike information criterion (AIC; Akaike 1974) to select the most informative representation of the distribution data. Explanatory variable effects were then evaluated on the selected model through deviance analysis (Venables and Ripley 2002). The general equation of the evaluated models was as follows:

where φ is the link function that links the mean to the model, SI corresponds to the shape indices, β₀ is the intercept of model, TL is total fish length, s is the smoothing function and Z (zone) is a dummy variable that represents the origin of the sample with j = 1,2 and i is the ith fish.

Elliptic Fourier analysis is a method that allows a closed curve to be described, characterised by equidistant (x,y) coordinates such as the otolith contour, through the infinite summation of ellipses with different amplitudes and angles. However, the objective of the analysis is to describe the otolith outline using the minimum number of ellipses. Each ellipse of the Fourier analysis is called an elliptical Fourier descriptor (EFD) or simply a ‘harmonic’, because of its functional form (sines and cosines). Because the shape described by the EFDs is sensitive to the orientation, size and starting point of the ordered pairs (x,y), a normalisation process was performed following Kuhl and Giardina (1982). Although normalisation reduced bias in EFDs due to correlation with other morphological variables, such as fish length, in order to avoid any undesired effect of fish size on otolith outlines, the analyses were conducted according to fish size classes. Hence, separated elliptical Fourier analyses were performed for fish ranging in size from 6.5 to 9.5 cm TL in the case of juveniles and from 13.5 to 14.0 cm in the case of adults collected in 2016; for adults collected in 2013, EFDs were obtained for adult fish ranging in size from 14.0 to 14.5 cm TL. A total of 200 EFDs was initially calculated for each otolith according to the algorithm implemented by Claude (2008) in the R package (R Foundation for Statistical Computing, Vienna, Austria). After a power analysis was performed, only 20 harmonics were sufficient to explain 95% of the variance in otolith outlines. The first EFDs, which corresponded to the ellipse, used to normalise the data were discarded for further analyses.

Stepwise canonical discriminant analysis (CDA) was used to determine whether otoliths collected in different sampling zones could be distinguished based on the 20 selected EFDs. CDA is a standard method where data belonging to known groups (sampling zones) are used to find linear combinations of descriptors that maximise Wilks’ lambda (λ; Ramsay and Silverman 2005; Pothin et al. 2006). Wilks’ λ is the ratio between the intragroup variance and total variance, and provides an objective means of calculating the chance-corrected percentage of agreement between real and predicted group membership. Values of Wilks’ λ range from 0 to 1: the closer λ is to 0, the better the discriminating power of the CDA (Lord et al. 2012).

There were differences in area, perimeter, Ferret length and Ferret width among length groups (ontogenetic state). Overall, shape indices showed a latitudinal gradient for all ontogenetic stages (i.e. the ellipticity and circularity were higher in Zone I and lower in Zones II and III). Conversely, rectangularity, form factor and roundness were higher in Zone III than in the other zones, at least for pre-recruits and recruits. In adults, mean shape indices showed closer values between Zones I and II than between each of these zones and Zone III (Table 2).

Table 2.  Mean (±s.d.) values for basic morphometric variables and shape indices of sagittal otoliths of pre-recruits, recruits and adults of the anchoveta (Engraulis ringens) by sampling zone
Zone I, Arica–Antofagasta; Zone II, Caldera–Coquimbo; Zone III, Valparaíso–Valdivia. TL, total length

Roundness and form factor decreased as fish length increased, whereas the opposite was found for circularity and ellipticity, which tended to increase with fish length. GAM analysis of the relationship between shape indices and fish length (Year 2016) across sampling zones showed an overall trend towards a continuous change in otolith shape with fish length, at least until 14 cm TL. Within this general pattern, more pronounced changes with size were noticeable between 4 and 8 cm TL, and there was slight evidence of stabilisation between 8 and 10 cm TL (Fig. 2).

Fig. 2.  Left-hand panels show generalised additive models for form factor, roundness, ellipticity, circularity and rectangularity as a nonlinear function of the total length of the anchoveta (Engraulis ringens) for Zone I (Arica–Antofagasta), Zone II (Caldera–Coquimbo) and Zone III (Valparaíso–Valdivia) off the coast of Chile. In the left-hand panels, the solid black line shows the tendency of partial residual of dependent variable (smooth term) against the fish length as predictor variable, and the grey shaded areas indicate the 95% confidence interval. Panels on the right show the standardised estimated effect, partial residual of dependent variable, for each level of factor zone. In the right-hand panels, the solid line shows the mean value of standardised dependent variable to each level factor zone and the dashed line shows the 95% confidence interval.

GAM deviance analysis showed sampling zone and fish length explained over 88% of the deviance in all shape indices except rectangularity, for which these parameters only explained 62.8% of deviance (Table 3). Considering fish length covariance, there were significant differences for all indices between sampling Zones I and III, as well as between Zones II and III. However, significant differences between Zones I and II were limited to form factor and ellipticity.

Table 3.  Deviance analysis from generalised additive models for each shape index by sampling zone
Zone I, Arica–Antofagasta; Zone II, Caldera–Coquimbo; Zone III, Valparaíso–Valdivia

A two-function discriminant model fit using otolith EFDs of juvenile fish (6.5–9.5 cm TL) explained 79 and 21% of variance, with canonical correlations of 0.88 and 0.69 respectively. Only the first discriminant function presented a significant value for Wilks’ λ of 0.117 (P < 0.0001). The lower Wilks’ λ of the second function (0.521; P < 0.2), suggested overlap between zones regarding this function. These canonical discriminant functions produced an accurate self-classification of 89% of samples, withe accurate classification in 95% of cases for Zone III. The total percentage of correct classifications of cross-validated grouped cases in their original group was 74%, and >76% in Zones I and III (Table 4; Fig. 3a).

Table 4.  Classification matrix of sagittal otoliths of juvenile anchoveta collected in 2016
Bold numbers indicate the numbers and percentage of fish correctly classification to the origin zone. Zone I, Arica–Antofagasta; Zone II, Caldera–Coquimbo; Zone III, Valparaíso–Valdivia

Fig. 3.  Scatterplots of the first and second discriminant function scores for Fourier coefficients of sagittal otolith outlines of anchoveta (Engraulis ringens) for (a) juveniles in 2016, (b) adults in 2016 and (c) adults in 2013. Zone I, Arica–Antofagasta; Zone II, Caldera–Coquimbo; Zone III, Valparaíso–Valdivia.

The two-function discriminant model fit using otolith EFDs of 2016 adults (13.5–14.0 cm TL) explained 67 and 33% of variance, with canonical correlations of 0.78 and 0.57 respectively. The first discriminant function was significant, with Wilks’ λ of 0.222 (P < 0.000). The Wilks’ λ of the second function was higher, at 0.568, at the edge of significance (P < 0.053). Classification reached an overall accuracy of 85%, with all sampling zones exhibiting crude accuracies >81%. The total percentage of correct classifications of cross-validated grouped cases in their original group was 64%, and 68% in Zone III (Table 5; Fig. 3b).

Table 5.  Classification matrix of sagittal otoliths of adult anchoveta collected in 2016
Bold numbers indicate the numbers and percentage of fish correctly classification to the origin zone. Zone I, Arica–Antofagasta; Zone II, Caldera–Coquimbo; Zone III, Valparaíso–Valdivia

Discriminant functions for 2013 adults explained 63 and 37% of the variance, with canonical correlations of 0.83 and 0.76 respectively. Wilks’ λ for the first discriminant function (0.130) was highly significant (P < 0.000), whereas the second function showed a much lower Wilks’ λ value (0.425; P = 0.072), indicating a greater overlap between sampling zones. The accuracy for 2013 adults reached 90% of all samples, with all sampling zones reaching accuracies >87%. The total percentage of correct classifications of cross-validated grouped cases in their original group was 66%, and 70% in Zone III (Table 6; Fig. 3c).

Table 6.  Classification matrix of sagittal otoliths of adult anchoveta collected in 2013
Bold numbers indicate the numbers and percentage of fish correctly classification to the origin zone. Zone I, Arica–Antofagasta; Zone II, Caldera–Coquimbo; Zone III, Valparaíso–Valdivia

Fourier reconstructed outlines (from the first 20 EFDs) for the six different ontogenetic stages suggested important differences among sampling zones (Fig. 4).

Fig. 4.  Mean otolith shapes from elliptical Fourier coefficients for four length ranges, representing three ontogenetic stages of anchoveta (Engraulis ringens). Shape indices were organised into six length interval groups linked to ontogenetic stages: two pre-recruit intervals of 4.0–5.5- and 6.0–7.5-cm total length (TL); two recruit intervals of 8.0–9.5 and 10.0–11.0 cm TL; and two adults intervals, one from 2016 (13.5–14.0 cm TL) and another from 2013 (14.0–14.5 cm TL).

For a long time, the anchovy fishery in Chile was managed using a purely administrative criterion without knowing the population structure of the species. In this context, the present study has made a substantive contribution, identifying at least two completely independent demographic units (I and III) of anchoveta (E. ringens) off the Chilean coast through of otolith shape analysis, where both basic shape indices and EFDs showed consistent results.

These results match with the stocks or management units used today, which support stock assessment of this species based on separated stocks. However, the lower discriminating power between Zones I and II suggests some level of mixing between these zones, which should be elucidated in further studies in order to specify the exact limits between populations to improve the management of this resource. Some of the characteristic features associated with these new findings are discussed below.

Otolith shape indices (form factor, circularity, rectangularity, roundness, ellipticity) showed high variability and significant differences among sampling zones. Although the magnitude of these morphological differences tended to increase with age and size, we found evidence indicating the effects of fish length on shape indices decreased with size, leading to relative stability in morphological differences among sampling zones after fish reached sexual maturity (~14 cm TL). Such a condition has been reported for the clupeoid Strangomera bentincki by Curin-Osorio et al. (2012), as well as for other species (e.g. Tuset et al. 2003). Consequently, evaluation of otolith shape stability is important to determine when the otolith outline is not affected by fish ontogeny, but rather by genetic and environmental variability. Furthermore, it would be very interesting to evaluate whether the pattern of otolith shape stability found in E. ringens in the present study is also applicable to other engraulid species. For example, Zengin et al. (2015) and Jemaa et al. (2015) identified population units of Engraulis encrasicolus in the Black and Mediterranean seas respectively, although in both these studies only immature fish (11–12 cm TL) were used, when the shape or otolith outline has not completely yet. In this sense, it is reasonable to infer that results derived from otolith shape analysis of adult fish in the present study are reliable.

A distinctive finding of the present study was the significant discriminatory capacity of EFDs to distinguish anchoveta from the three localities, even when the analysis was conducted in two different years (i.e. 2013 v. 2016). These results support previous findings where a high discriminatory capacity of EFDs has been reported in several species of teleost fishes (i.e. Popper et al. 2005; Brophy et al. 2016; Afanasyev et al. 2017; Duncan et al. 2018). Moreover, similar original classification accuracies (>81%) were found for E. encrasicolus when this technique was used to differentiate four populations in the Mediterranean Sea (Jemaa et al. 2015). These are promising results for fisheries resource monitoring purposes, because otolith shape is a cost-effective technique that would continuously be updated and enhanced because of advances in image-based platforms.

A further finding of the present study was that Fourier reconstructed otolith outlines showed differences between zones that appear to be larger for sampling Zone III (Valparaíso–Valdivia), where at the same fish length otoliths were smaller than those sampled in Zones I and II, at least at the pre-recruit stage (<8 cm TL). The smaller and rounder otoliths of Zone III could be associated with growth differences among juveniles, because a recent study reported much slower growth of pre-recruits in southern (zone III) than northern zones (zone I) (Cerna and Plaza 2015).

The marked differences found in otolith morphology between the two extreme localities in the present study are strong evidence of an environmental spatial heterogeneity of anchoveta along the Chilean coast. This evidence matches well with previous studies, which have proposed that environmental variability can produce important differences in otolith shape (Lombarte and Lleonart 1993; Cardinale et al. 2004; Vignon 2012). Indeed, juveniles and adult fish inhabiting the extreme northward and southward zones analysed for anchoveta in this study must cope with important environmental differences in biological productivity (e.g. chlorophyll-a), upwelling intensity, temperature and both physical and topographic changes. An important feature here is the CTZ defined between 34 and 39°S (Morales et al. 2010). Although the northern zones (Zones I and II) do not have a strict topographic boundary between them, there is a critical latitude (~22.8°S) described around Mejillones and Moreno bays that has higher concentrations of chlorophyll throughout the year and stronger water retention dynamics related to local topography and particular Ekman transport characteristics of adjacent waters (Letelier et al. 2012). Hence, it is reasonable to hypothesise that this critical latitudinal zone could enhance retention of early stages of anchoveta and limit the flow of early juveniles between Zones I and II.

Phenotypic differences found among sampling zones were large and consistent between years and through ontogenetic stages, which is indicative of a prolonged and demographically relevant separation among units (Reiss et al. 2009). However, these differences do not provide direct evidence of genetic (i.e. nearly complete) isolation between sampling areas. Nonetheless, although classical genetic studies performed for this area and species have failed to reject the conventional hypothesis of panmixia (Galleguillos et al. 1996; Ferrada et al. 2002; Rojas 2011), new evidence from Ferrada et al. (2018) suggests progressive differentiation with distance and may change the current view regarding the population genetics of anchoveta along the Humboldt Current. In this sense, in the present study, otolith-based morphometric differences were larger for sampling Zone III (Valparaiso–Valdivia), than for samples from Zones I and II. The high differences between Zone III and both Zones I and II suggest that genetic or environmental differences are greater between these areas, but lower between Zones I and II, which is consistent with the geographic distribution of the species (Castillo et al. 1998; Leiva et al. 2016), as well as with previous parasitological studies (Valdivia et al. 2007; George-Nascimento et al. 2011).

The similar classification accuracies found in 2013 and 2016 provide the first evidence that otolith morphological differences among zones were not sensitive to the El Niño event, which dominated the oceanographic conditions during 2014–16 in the Southern Hemisphere (Vera and Osman 2018). This finding suggests the absence of a generalised movement and mixing of schools in a southward direction, which could be expected as a result of movement of warm water masses from north to south. In fact, latitudinal stability in fish distribution under El Niño events has been observed in recruiting acoustic cruises (Castillo et al. 1998; Leiva et al. 2016), where a change in bathymetric distribution seems to be the dominant response to warm water conditions due to the El Niño event, where anchoveta move to deeper waters and are less vulnerable to purse seine fishing gear.

The existence of clear otolith-based morphological differences in both juvenile and adult fish among the two extreme localities supports the hypothesis that individuals have remained segregated throughout the entire (or at least a fraction of) their life cycle, supporting the increasing scientific evidence over the past two decades of homing in marine fishes (e.g. Thorrold et al. 2001; Rooker et al. 2008; Sólmundsson et al. 2015). Hence, further studies focused on extending both the spatial and temporal scale, as well as incorporating additional ecological markers (e.g. otolith microchemistry and microstructure) would provide additional insights to reveal homing and the level of mixing, if any, in this important fishery resource in the Humboldt Current system.

The authors declare that they have no conflicts of interest.

The sampling for this study was performed as part of the hydroacoustic recruitment survey of anchoveta and the Spawning Monitoring Project, funded by the Chilean Ministry of Economy, Promotion and Tourism (Asesoría Integral en Pesca y Acuicultura, ASIPA, grants 2015 and 2016). Morphological analyses were funded by the Chilean Fund for Fisheries and Aquaculture Research (Fondo de Investigación Pesquera y Acuícola, FIPA, grant 2015–22).