Otolith microchemistry: a useful tool for investigating stock structure of yellowfin tuna (Thunnus albacares) in the Indian Ocean

Stock identification is an essential component of modern fisheries management (Pita et al. 2016). Although several definitions of ‘stock’ have been used in the fisheries literature, the stock concept essentially describes demographically discrete units of fish assumed to be homogeneous for particular management purposes (Begg et al. 1999). A fundamental requirement for sustainable fisheries management is the match between management action units and biologically relevant processes (Reiss et al. 2009). Different stocks may possess specific genetic, physiological and behavioural traits that can influence life processes such as growth rates, fecundity, abundance and disease resistance (Stepien 1995), and thus different resilience to exploitation and environmental changes. As such, when stock structure is more complex than recognised, the sustainability of the resource and the profitability of the fishery can be jeopardised (Kerr et al. 2017; Taillebois et al. 2017).

Among the existing methods for fish stock structure assessment (e.g. genetic markers, tagging, life history parameters, parasites etc.), otolith microchemistry has been recognised as a useful stock discriminant (Kerr and Campana 2014; Tanner et al. 2016). These hard calcareous structures are found in the inner ear of vertebrates and are used for balance and hearing in teleost fish. Otoliths are acellular and metabolically inert (i.e. newly deposited material is neither reabsorbed nor reworked after deposition) and grow continually throughout the life of the fish (Campana and Neilson 1985). Ambient water chemistry and environmental conditions affect elemental incorporation into the otolith (Elsdon and Gillanders 2004; Walther and Thorrold 2006) and variations in the concentrations of selected elements and isotopes (also known as ‘chemical fingerprints’) can therefore be used as natural markers to discriminate among groups of fish for stock identification purposes (Kerr and Campana 2014). However, factors that affect the elemental incorporation into otoliths are not yet completely understood because various factors other than ambient water chemistry, including dietary sources, physiological processes and genetic variation, may also play important roles in elemental incorporation (Ranaldi and Gagnon 2008; Clarke et al. 2011; Sturrock et al. 2015; Izzo et al. 2018).

The potential of otolith chemistry to improve our understanding of population structure (e.g. Wells et al. 2015), natal origin (e.g. Rooker et al. 2003; Shiao et al. 2010; Wells et al. 2012) and movement patterns (e.g. Arai et al. 2005; Macdonald et al. 2013; Fraile et al. 2016) of different tuna species has been demonstrated. However, few studies have used otolith microchemistry to assess the stock structure of yellowfin tuna (Thunnus albacares (Bonnaterre, 1788)), and these studies are limited to the Pacific (Wells et al. 2012; Rooker et al. 2016) and Atlantic (Shuford et al. 2007; Kitchens et al. 2018) oceans.

Yellowfin tuna is a highly migratory species that inhabits the epipelagic zone of tropical and subtropical marine waters around the world (Froese and Pauly 1999). It is a significant component of the global fishing industry and economy (Galland et al. 2016) and constitutes the second largest tuna fishery worldwide (Food and Agriculture Organization of the United Nations 2016). The high levels of exploitation of yellowfin tuna render the species vulnerable to unsustainable harvest. In the Indian Ocean, in particular, increased harvest in recent years has led to overexploitation of the stock (International Seafood Sustainability Foundation 2017). The most recent stock assessment model by the Indian Ocean Tuna Commission (IOTC), which takes into account both stock abundance and fishing mortality, confirmed that the yellowfin tuna stock in the IOTC area is overfished (Indian Ocean Tuna Commission 2017a).

The current stock assessment model for yellowfin tuna in the Indian Ocean assumes a single stock (Indian Ocean Tuna Commission 2017a). This assumption is based on tag recoveries from the Regional Tuna Tagging Program of the Indian Ocean (RTTP-IO), which showed large and rapid movements of the species around the Indian Ocean (Fonteneau and Hallier 2015). However, uncertainties remain regarding the stock structure of yellowfin tuna in this region, because tag–recapture studies may not reveal important aspects of the movement dynamics of wide-ranging species (see Rooker et al. 2008a). Indeed, genetic studies using mitochondrial (mt)DNA suggest that a more complex stock structure may exist within yellowfin tuna populations. For example, genetically discrete groups of this species have been reported in the north Indian Ocean (Dammannagoda et al. 2008) and within Indian Ocean waters (Kunal et al. 2013).

In the present study we analysed stable isotopes and trace elements in otoliths of young-of-the-year (YOY) and age-1 yellowfin tuna, which are both considered immature (Zudaire et al. 2013). During the YOY period, it is assumed that movements are restricted to a small area around the natal site, whereas for age-1 individuals large-scale movements can be expected (Fonteneau and Hallier 2015). Fish were collected in two regions of the western Indian Ocean, the north-west (NW) Indian Ocean and Mozambique Channel, because major yellowfin tuna spawning grounds have been identified westward of 75°E (Zudaire et al. 2013). In particular, the equatorial area (0–10°S) and the Mozambique Channel have been identified as primary and secondary spawning areas (Reglero et al. 2014; Indian Ocean Tuna Commission 2017b). The western region accounted for the 80% of the Indian Ocean catches in 2017 (Indian Ocean Tuna Commission 2018), being the most important area for the Indian Ocean yellowfin tuna fishery. The aim of the study was to assess the feasibility of otolith microchemistry as a potential tool for investigating yellowfin tuna stock structure in the Indian Ocean. Specifically, our objectives were to: (1) explore otolith chemistry profiles of YOY yellowfin tuna; (2) seek the best methodology for fish origin analyses and test region-specific differences on chemical composition of YOY yellowfin tuna otoliths; (3) evaluate the usefulness of different classification methods in assigning YOY yellowfin tuna origin to their respective collection regions; and (4) investigate potential natal-origin differences in age-1 yellowfin tuna.

Samples were collected across two regions of the western Indian Ocean: (1) offshore western Madagascar (15°53′S–13°30′S, 43°33′E–47°05′E), hereafter ‘Mozambique Channel’; and (2) offshore Somalia and Seychelles regions (7°77′S–7°22′N, 48°57′E–58°37′E), hereafter ‘NW Indian Ocean’ (Fig. 1). These two regions belong to two distinct provinces according to the Longhurst (2007) biogeographical classification of the pelagic ecosystem: the Mozambique Channel belongs to the East African Coastal Province (EAFR) and the NW Indian Ocean is part of the Indian Monsoonal Gyre Province (MONS). These two regions are defined by different oceanographic characteristics that affect the biological production of the area (for further information, see Fig. S1–S5, available as Supplementary material to this paper). North of 10°S the Indian Ocean is characterised by two seasons with distinct wind regimes driving the ocean climate and circulation: the monsoon system (Schott and McCreary 2001). The south-west (summer) monsoon drives upwelling events that lead to high primary productivity off the coast of Somalia. By contrast, during the north-east (winter) monsoon, conditions are not favourable for upwelling and primary productivity is generally lower (Wiggert et al. 2006). The Mozambique Channel is characterised by the presence of anticyclonic eddies that generate an upward movement of nutrient-rich waters around their edge, and advection of nutrient-rich coastal waters when they run along the coast (Quartly and Srokosz 2004; Tew Kai and Marsac 2010).

Fig. 1.  Study area in the western Indian Ocean, where two Longhurst provinces, the Eastern Africa Coastal Province (EAFR) and the Indian Monsoon Gyres Province (MONS) are represented. Sampling locations of the different groups of yellowfin tuna (Thunnus albacares) used in this study are indicated as circles in the Mozambique Channel and triangles in the north-west Indian Ocean for young-of-the-year, and as crosses for age-1 fish of unknown origin.

Yellowfin tuna (n = 106) were collected in the western Indian Ocean by scientific observers on board Spanish commercial purse seiners between February 2009 and July 2010. Fish were measured (fork length (FL); to the nearest centimetre) and weighed (to the nearest 0.1 kg) on board. Sagittal otoliths were extracted, cleaned of adhering organic tissue, rinsed with ultrapure water and stored dry in plastic vials. The otolith collection available for this study comprised fish from different collection dates, age classes and cohorts. Following the growth curve described by Sardenne et al. (2015), fish were classified as YOY (<43 cm FL) or age-1 (43–75 cm FL). Birth year was then back-calculated from the catch date using the age–length relationship of Sardenne et al. (2015). In order to reduce potential variability due to temporal variation among samples, only fish born in the same monsoonal regime (i.e. summer monsoon) were selected for analysis. In addition, to minimise the possibility that YOY tuna had migrated outside their origin areas before sampling, only fish <40 cm were considered in the study. To date, no tagging experiments have considered fish of that size in movement range studies of yellowfin tuna (see Hallier and Fonteneau 2015). Nevertheless, our threshold size is even more conservative than those already considered in Atlantic Ocean yellowfin tuna stock structure analyses (Kitchens et al. 2018; Pecoraro et al. 2018). Thus, the final dataset consisted of 52 YOY and 15 age-1 fish (Table 1).

Table 1.  Summary of young-of-the-year (YOY) and age-1 yellowfin tuna (Thunnus albacares) used for otolith chemistry analyses from two regions of the western Indian Ocean: the Mozambique Channel and the north-west (NW) Indian Ocean
Size is fork length (FL). The cohort was calculated by back-calculation after age estimation following the growth curve described by Sardenne et al. (2015)

Whole otoliths were briefly soaked in nitric acid (1%) and cleaned with ultrapure (Milli-Q) water to eliminate any remaining biological material (Rooker et al. 2001). Left and right otoliths were embedded in two-part resin (EpoFix; Struers, Ballerup, Denmark) and sectioned using a Buehler low-speed Isomet saw (Buehler, Lake Bluff, IL, USA) fitted with a diamond-edged blade to obtain a transverse section that included the core. When both otoliths from the same fish were available, one otolith was randomly selected for stable C and O isotope analysis and the other was used for trace element analysis. When only one otolith was available, stable isotope analyses were conducted.

Transverse sections were polished with a series of grinding and polishing films moistened with ultrapure water, and then with a microcloth and 0.3-mm alumina powder to ensure a smooth surface. Sections were glued in a sample plate using Crystalbond thermoplastic glue (Crystalbond 509; Buehler). For YOY, 41 otoliths were used for trace element chemistry and 52 for stable isotope analyses (Table 1). For age-1 fish, trace element composition of 15 yellowfin tuna was analysed (Table 1).

Laser ablation–inductively coupled plasma–mass spectrometry (LA-ICP-MS) was used to determine the trace element chemistry of the otolith. Analyses were performed using a 7700X quadrupole inductively coupled plasma mass spectrometer (ICP-MS, Agilent Technologies, Santa Clara, CA, USA) coupled to a laser ablation system with a Laurin Technic 2 volume cell (RESOlution LA system, Resonetics, Nashua, NH, USA) available at the School of Earth Sciences, The University of Melbourne (Melbourne, Vic., Australia). The laser was operated using a spot size of 26 μm in diameter with laser energy at 2.4 J cm^–2 and a repetition rate of 10 Hz. Otoliths were placed in the sealed chamber and ablated continuously across the growth axis from the primordium to proximal edge at a scan speed of 6 µm s^–1 (Fig. 2a). Ablation occurred in an atmosphere of pure He to minimise condensation of ablated material. The ablated material was transported to the ICP-MS mixed with argon carrier gas. A preablation step was implemented to minimise potential surface contamination. National Institute of Standards and Technology (NIST) 612 glass standard with known chemical composition was used for calibration. Standards were analysed twice at the beginning and end of each session, and once after one-third and two-thirds of the duration of the run. In addition, measurement precision was determined based on the analyses of NIST610 and the calcium carbonate pressed pellet standard MACS-3 (US Geological Survey, Reston, VA, USA; Wolf and Wilson 2007; Limburg et al. 2011). The variation in ablation yields between standards and unknown samples was corrected with the use of calcium as an internal standard (Craig et al. 2000). Relative abundances of 10 isotopes (⁷Li, ²⁵Mg, ⁵⁵Mn, ⁵⁷Fe, ⁵⁹Co, ⁶⁰Ni, ⁶³Cu, ⁶⁶Zn, ⁸⁸Sr and ¹³⁸Ba) were measured along the whole transect using the LA-ICP-MS system, and elemental concentrations (ppm) were then calculated. Data reduction and processing were performed with the software package Iolite (ver. 2.2, School of Earth Sciences, The University of Melbourne), which operates within IGOR Pro (ver. 6.2, WaveMetrics, Lake Oswego, OR, USA; Paton et al. 2011; Paul et al. 2012). The limit of detection (LOD) was calculated based on the ablation yield equivalent to three times the standard deviation of the blank signal. Averaged LOD values (ppm) of each element across all samples were as follows: Li, 0.049; Mg, 0.423; Mn, 0.079; Fe 0.920; Co, 0.020; Ni, 0.030; Cu, 0.023; Zn, 0.120; Sr, 0.0145; and Ba, 0.004. Five elements (Ba, Mg, Mn, Fe and Sr) were consistently present at concentrations above LODs within the whole transects. However, low measurement precision and accuracy was reported for Fe in MACS-3 and NIST610 (Table S1, available as Supplementary material to this paper), due to time-dependant increases in apparent concentrations. This is likely due to decreasing oxide–hydroxide interferences (e.g. ⁴⁰Ca-¹⁶O-H) with time affecting the low Fe/Ca NIST612 calibration standard proportionally more than the higher Fe/Ca NIST610 and MACS-3 standards. Thus, Fe was discarded for analyses. In addition, an ablation transect was performed on a test sample composed only of resin and thermoplastic glue. This analysis showed that resin was enriched in Zn, whereas high levels of Cu were found on the thermoplastic glue. Therefore, unusually high concentrations of Cu and Zn throughout the transect were used as indicators of potential sample contamination. Aragonitic calcium carbonate stoichiometry (400 000 mg Ca g^–1 otolith) was used for calcium concentration calculation, and hereafter the concentrations of the rest of the elements are expressed as element : Ca ratios (Ludsin et al. 2006).

Fig. 2.  Transverse section of a young-of-the-year yellowfin tuna (Thunnus albacares) sagittal otolith showing (a) the transect of the laser ablation spots performed along the growth axis, (b) the position of the data slices represented as the core and near-core and (c) pattern of the portion milled for stable isotope analysis.

In addition to whole transects, two data slices 78 µm in length were isolated from the primordium (hereafter ‘core’) and 200 µm beyond the primordium (hereafter ‘near-core’), by averaging all trace element concentration estimates across each data slice (Fig. 2b). Daily increment formations were counted along the first 550 μm of the transect in a subset of YOY otoliths (n = 8) under a light microscope (Fig. S6, available as Supplementary material to this paper). It was determined that the two 78-µm slices represent the otolith material accreted during the first ~10 days and first ~2–3 weeks after hatching respectively, when the fish are still in the larval phase (Kaji et al. 1999).

Microsampling of otolith powder for carbon (δ¹³C) and oxygen (δ¹⁸O) stable isotope analysis was performed using a high-resolution computerised micromill (New Wave MicroMill System, New Wave Research G. C. Co., Ltd, Cambs, UK) only in YOY yellowfin tuna. The portion of the smallest fish (29 cm FL) in our sample was used to create a standard template that was then used with the remaining otoliths (Fig. 2c) to ensure that the same portion of the otolith was analysed in each individual. This is approximately related to the first 5–6 months of life according to the age–length relationship described by Sardenne et al. (2015).

Fourteen drill passes were run at a depth of 50 µm per pass over a preprogrammed drill path using a 200-µm diameter carbide bit (Komet dental; Gebr. Basseler, Lemgo, Germany). Powdered material was then analysed for δ¹³C and δ¹⁸O on an automated carbonate preparation device (KIEL-III, ThermoFisher Scientific, Waltham, MA, USA) coupled to a gas-ratio mass spectrometer (Finnigan MAT 252, ThermoFisher Scientific) at the Environmental Isotope Laboratory of the University of Arizona. All isotope values were reported according to standards of the International Atomic Energy Agency in Vienna (see https://www-pub.iaea.org/MTCD/publications/PDF/te_825_prn.pdf, accessed May 2012). Values of δ¹³C and δ¹⁸O represent ratios of ¹³C/¹²C and ¹⁸O/¹⁶O in the sample relative to the Vienna Pee Dee Belemnite (VPDB) scale.

All statistical analyses were performed using open access R software (R Foundation for Statistical Computing, Vienna, Austria, see http://www.R-project.org/, accessed September 2014). Prior to all multivariate analyses, trace element data were scaled to give the same weight to all elements. An α level of 0.05 was used for all two-tailed statistical tests.

Whole transects of YOY yellowfin tuna were transformed into equally spaced observations, generating a regular time series of elemental concentration from 0 to ~550 µm (i.e. approximately the first month after hatch). Data transformation and ordination were performed using zoo (ver. 1.8-6, see http://zoo.R-Forge.R-project.org/, accessed 6 November 2019; Zeileis and Grothendieck 2005) and vegan (J. Oksanen, F. Blanchet, M. Friendly, R. Kindt, P. Legendre, D. McGlinn, P. R. Minchin, R. B. O’Hara, G. L. Simpson, P. Solymos, M. Henry, H. Stevens, E. Szoecs, and H. Wagner, see https://cran.r-project.org/package=vegan, accessed 13 March 2018) packages. To examine leading patterns and temporal variations for the elements, regular time series data were simplified by averaging concentrations every 25 µm. Estimates of the s.d. were also calculated. Mean estimates of concentrations and profiles of element : Ca ratios were obtained for each region. In addition, a principal component analysis (PCA) was applied to each element in the time series data. To determine whether differences existed between regions, a permutational multivariate analysis of variance (PERMANOVA) test was performed on the scores of the first two axes.

YOY yellowfin tuna were used to determine whether the fish collected in Mozambique Channel and NW Indian Ocean origin could be differentiated based on otolith stable isotope and trace element concentrations measured in the otolith core and near-core. Normality and homoscedasticity were tested using Shapiro–Wilks and Fligner–Killeen tests respectively. Stable isotope data did meet these assumptions, so Student’s t-test was used for inter-region comparisons. Trace element data violated parametric assumptions; consequently, a robust statistical method was used to investigate these data, specifically a bootstrapped Yuen t-test (trimmed level for the mean = 0.2, number of bootstrap samples = 599). Robust statistics are not unduly affected by outliers and violations of parametric assumptions, reducing the probability of making a Type I error and improving test power when data are not normal or homogeneous (Daszykowski et al. 2007; Erceg-Hurn and Mirosevich 2008; Wilcox 2012). Accordingly, a Hotelling T² test was used to test the significance of differences in δ¹³C and δ¹⁸O values among regions and PERMANOVA was used to test the significance of differences in element : Ca ratios. Finally, core and near-core trace element data were examined using non-metric multidimensional scaling (nMDS) to obtain a reduced two-dimension graphical representation of the relationships between samples (Hout et al. 2013). The dissimilarity matrix construction was based on Euclidean distances. The otolith portion showing greatest variation among regions was used as the origin indicator in subsequent analyses.

Discrimination capacity between fish collected in the two regions was assessed using YOY otolith near-core microchemistry data using different classification methods. Four different methodologies were used, two classical methods (linear discriminant analysis (LDA) and quadratic discriminant analysis (QDA)) and two data-mining methods (random forest (RF) and artificial neural network (ANN)). LDA and QDA were completed using lda and qda functions from the MASS library; RF was implemented using the randomForest function from the randomForest library, using default settings (Liaw and Wiener 2002) and ANN was performed using the nnet function from the nnet library, using default settings and 30 hidden neurons at the intermediate layer (Venables and Ripley 2002). Detailed descriptions of each technique and their assumptions can be found in Jones et al. (2017) and Mercier et al. (2011). In all cases, data were randomly split into a training dataset (75%) and a testing dataset (25%) to perform a cross-validation procedure and measure prediction quality. To avoid sampling effects, the training subset was resampled 10 times and mean prediction accuracy (i.e. percentage of correct assignment of the fish to their known region) was obtained. In addition, the kappa index was calculated as a measure of classification success. Kappa values are always ≤1, with a value of 1 implying perfect agreement. The kappa index is unbiased when the number of samples differs between groups, and can present more reliable results because it considers the possibility of the agreement occurring by chance (Fielding and Bell 1997; Jones et al. 2017). Following Mercier et al. (2011), the elemental combination maximising habitat discrimination was selected by performing a stepwise forward variable selection procedure using the Wilk’s lambda criterion for LDA and QDA. In the case of RF, the best elemental combination was obtained by stepwise procedure using the Gini index (Cutler et al. 2007), and in the case of ANN following Olden et al. (2004).

RF clustering analysis of the near-core trace element data was performed in age-1 fish to determine the number of larval sources contributing to the sample. This approach proved to be useful for natal origin-related investigations when there is no previous baseline of potential larval sources (Gibb et al. 2017; Régnier et al. 2017; Wright et al. 2018). Moreover, the RF clustering approach does not require continuous multivariate data with normal distribution (Shi and Horvath 2006), an assumption that is often violated when working with otolith microchemistry data. The RF classification process was applied following Shi and Horvath (2006). First, a synthetic dataset was created by random sampling from the product of empirical marginal distributions of the elements. Then, unsupervised RF was used to separate the synthetic data (used as a class factor) from the observed data, and a dissimilarity matrix was created. Second, the dissimilarity matrix was used as input for partitioning around the medoid (PAM) clustering. Finally, the appropriate number of clusters was determined using the gap statistic (Tibshirani et al. 2001), with the fviz_nbclust function from the factoextra library (bootstrap number = 100). Canonical variable coefficients of YOY and age-1 yellowfin tuna trace element composition were plotted to visualise relationships (Hamer et al. 2005).

Patterns of element : Ca variation in YOY along the analysed transect were very similar between the two regions, except for Mn : Ca (Fig. 3, 4). Ba : Ca concentrations declined with increasing time after hatch, with YOY from the NW Indian Ocean having higher mean Ba values than those from the Mozambique Channel along most of the transect. As of the third week after hatch, the decreasing trend smoothed and became more similar between regions (Fig. 3a). Mg : Ca concentrations of both regions declined almost linearly within the first month of life. However, YOY from the Mozambique Channel had higher mean concentrations than those from the NW Indian Ocean along the whole transect (Fig. 3b). Mean Mn : Ca concentrations showed different trends between regions. For the Mozambique Channel, they were fairly constant during the first month after hatch. However, for NW Indian Ocean YOY, Mn : Ca increased sharply in the first ~2.5 weeks after hatch, with a progressive decrease thereafter (Fig. 3c). Sr : Ca showed little variation both among regions and along the transect, although slightly decreasing trends of concentrations were observed (Fig. 3d).

Fig. 3.  Mean (±s.d.) element : Ca concentrations across the first month of life of young-of-the-year yellowfin tuna (Thunnus albacares) collected in the Mozambique Channel (circles) and north-west Indian Ocean (triangles). Shaded areas represent the signatures that will correspond to the core and near-core portions of the otolith.

Fig. 4.  Principal component (PC) analysis correlation biplot on element : Ca time series data of the first month of life in young-of-the-year yellowfin tuna (Thunnus albacares) collected in the Mozambique Channel (circles) and north-west Indian Ocean (triangles). Ellipses represent the 95% confidence intervals for each region (dotted ellipse, Mozambique Channel; solid ellipse, north-west Indian Ocean). P-values are for differences in PC1 and PC2 scores according to permutational multivariate analysis of variance (PERMANOVA).

The PCA showed that all element : Ca time series were related among regions with the exception of Mn : Ca (Fig. 4). For this element, the scores of the first two axes of the PCA varied significantly between fish from the two regions (PERMANOVA, F_1,39 = 14.96, P = 0.01). The first two axes explained 39% of the total variation.

Significant differences in multi-elemental trace signatures (i.e. Ba : Ca, Mg : Ca, Mn : Ca and Sr : Ca) were found between the Mozambique Channel and NW Indian Ocean in the near-core portion of the otolith (PERMANOVA, F_1,39 = 4.09, P = 0.010) but not in the core (PERMANOVA, F_1,39 = 2.32, P = 0.050). Similarly, differences in trace element composition were more marked in the two-dimension nMDS plot of the otolith near-core than in the core (Fig. 5). Indeed, only Mg : Ca ratios were significantly different in the otolith core (Yuen test, t = 2.22, d.f. = 11.36, P = 0.046), with YOY from the Mozambique Channel having higher values than those from the NW Indian Ocean (Fig. 6b). YOY from the Mozambique Channel also had higher Mg concentrations than those from the NW Indian Ocean in the near-core portion (Fig. 6b; Yuen test, t = 4.05, d.f. = 12.06, P = 0.001). Ba and Mn also varied significantly in the otolith near-core among regions (Yuen test, t = 2.11 (d.f. = 24.89; P = 0.045) and t = 2.36 (d.f. = 16.60; P = 0.030) respectively), but in this case YOY from the NW Indian Ocean had higher values of both elements (Fig. 6a, c). There were no significant variations in Sr concentrations between regions in near-core signatures (Fig. 6d; t = 1.18, d.f. = 16.08, P = 0.256).

Fig. 5.  Non-metric multidimensional scaling (nMDS) of young-of-the-year yellowfin tuna (Thunnus albacares) otolith (a) core and (b) near-core trace elements (Ba : Ca, Mg : Ca, Mn : Ca and Sr : Ca). Each point represents an individual fish belonging to the Mozambique Channel (circles) and north-west Indian Ocean (triangles) regions.

Fig. 6.  Mean (±s.d.) element : Ca ratios in otolith cores and near-cores of young-of-the-year yellowfin tuna (Thunnus albacares) collected at two different regions, Mozambique Channel (left bars) and NW Indian Ocean (right bars). Asterisks indicate the significance of differences among regions according to bootstrapped Yuen t-tests (*, P < 0.05; **, P < 0.01).

Stable isotope composition of YOY yellowfin tuna did not vary between fish collected in the Mozambique Channel and NW Indian Ocean (Fig. 7; Hotelling T² test, T_2,41 = 0.76, P = 0.474). Mean (±s.d.) otolith δ¹³C values were –10.27 ± 0.43 and –10.13 ± 0.35‰ for YOY from the Mozambique Channel and NW Indian Ocean respectively (Student’s t-test, t_26,84 = –1.15, P = 0.259). The mean otolith δ¹⁸O value for fish collected in the Mozambique Channel was –2.32 ± 0.14‰, whereas for fish collected in the NW Indian Ocean it was –2.29 ± 0.23‰ (Student’s t-test, t_41,72 = –0.53, P = 0.600).

Fig. 7.  Scatter plot of otolith δ¹³C and δ¹⁸O of young-of-the-year yellowfin tuna (Thunnus albacares) captured in the Mozambique Channel (circles) and north-west Indian Ocean (triangles) regions. Marginal distribution curves of δ¹³C and δ¹⁸O are also presented for the Mozambique Channel (dotted lines) and north-west Indian Ocean (solid lines).

Cross-validated classification success based on YOY yellowfin tuna near-core trace element data ranged from 80 to 91% depending on the statistical method selected. However, results were more conservative according to the kappa index (κ = 0.53–0.72). LDA had the lowest differentiation capacity (Table 2). QDA and ANN showed similar results, with ANN presenting slightly higher accuracy and kappa value than QDA (Table 2). RF outperformed the other methods and produced the highest classification success (91%; κ = 0.72). Regardless of the method used for prediction, the best elemental combination always included Mg : Ca and Mn : Ca.

Table 2.  Comparison of classification methods for assigning young-of-the-year (YOY) yellowfin tuna (Thunnus albacares) collected in the Mozambique Channel and north-west (NW) Indian Ocean to their origin location based on near-core otolith microchemistry
Best elemental combination and performance, as a measure of maximal accuracy and the kappa index, are shown for each method. LDA, linear discriminant analysis; QDA, quadratic discriminant analysis; ANN, artificial neural network; RF, random forest

RF clustering of age-1 fish from unknown origin suggested a common larval source among the analysed samples. The canonical variable plot displaying spatial differences of the multi-elemental data showed that age-1 yellowfin near-core elemental composition was more similar to YOY of the NW Indian Ocean than to YOY of the Mozambique Channel, although there was considerable overlap among the regions (Fig. 8).

Fig. 8.  Canonical variable plots comparing trace element tags (Ba : Ca, Mg : Ca, Mn : Ca and Sr : Ca) of age-1 yellowfin tuna (Thunnus albacares) from unknown origin (crosses) with those of young-of-the-year (YOY) yellowfin tuna captured in the Mozambique Channel (circles) and north-west Indian Ocean (triangles) regions. Confidence ellipses cover the 95% of each group (dashed ellipse, age-1 tuna; dotted ellipse, Mozambique Channel YOY; solid ellipse, north-west Indian Ocean YOY).

Improving understanding of the stock structure of exploited species is essential for defining appropriate management units and optimising the use of fishery resources. In this context, this study used, for the first time, an otolith microchemistry approach for yellowfin tuna in the Indian Ocean and revealed the potential of this technique to provide valuable information about its stock structure. The main limitation of the study was that the YOY yellowfin tuna collected in the two regions were sampled from different cohorts. Therefore, observed differences between regions could result from the interannual variability in the oceanic system rather than variability among regions. However, this hypothesis seems unlikely, because interannual variability in environmental conditions was negligible for the observed periods (Table S2, available as Supplementary material to this paper).

Almost all elements followed a similar pattern of variation within the first month after hatch in both regions. Element profile analyses of natural otolith tags can be used for environmental and life history reconstructions (Elsdon et al. 2008). For example, Ba and Sr incorporation into otoliths has been shown to be correlated with ambient concentrations (Webb et al. 2012; Izzo et al. 2018). These alkaline metals are incorporated into the mineral component of the otolith (Izzo et al. 2016; Loewen et al. 2016), and thus appear to be reliable ‘geographical markers’ (Thomas et al. 2017). Although both elements have been related to salinity (e.g. Elsdon and Gillanders 2005; Macdonald and Crook 2010), the fact that Ba : Ca varied along the transect whereas Sr : Ca did not may suggest that observed Ba variation did not reflect changes in salinity. Indeed, no noticeable seasonal variation of sea surface salinity was detected in the western Indian Ocean (Fig. S1). Elevated Ba : Ca concentrations in the primordium relative to adjacent parts of the otolith have been previously described in the literature (Ruttenberg et al. 2005). However, this does not exclude other external influences of Ba incorporation. For example, because dissolved Ba in seawater possesses a nutrient-like profile (i.e. very low concentration in surface waters and increased values with deeper water), an increase in the Ba : Ca concentration may reflect periods of residence near upwelling zones (Patterson et al. 2004; Kingsford et al. 2009; Wang et al. 2009). During the summer monsoon (the period when YOY were born), elevated phytoplankton blooms have been described to extend off the Somali coast (Fig. S2; Wiggert et al. 2006), which could also be the reason for increased Ba concentration. Equally, classic eddy upwelling occurs throughout the Mozambique Channel, enhancing chlorophyll-a concentrations in the centre of the cyclones (Lamont et al. 2014). Magnesium concentrations decreased linearly with time after hatch. This pattern matches with previously described Mg : Ca profiles in the literature, with higher levels at early life (near otolith primordium) followed by a transition to lower Mg levels (Limburg et al. 2018). The uptake rate of Mg into otoliths has been related to somatic growth and otolith precipitation (Bath Martin and Thorrold 2005; Limburg and Casini 2018). The incorporation of this element into otoliths seems to be strongly regulated by physiological processes between the otolith and water (Hamer and Jenkins 2007). Therefore, the concentration of Mg in fish otoliths has been questioned as a reliable environmental indicator (Woodcock et al. 2012; Thomas et al. 2017). Mn : Ca profiles proved to be significantly different between the NW Indian Ocean and Mozambique Channel YOY yellowfin tuna. Manganese is a transition metal and, as such, its incorporation into the otolith can be physiologically regulated, which may caution against its use as an environmental marker (Thomas et al. 2017). Manganese concentrations have also been found to be high near the otolith primordium, possibly due to maternal transfer (Ruttenberg et al. 2005). However, no higher Mn concentrations in the first weeks of life were observed for YOY from the Mozambique Channel relative to the rest of the transect. Although sensitive to growth effects, environmental influences seem also to affect otolith Mn concentrations (Limburg and Casini 2018). In particular, otolith Mn : Ca concentrations have been shown to be a reliable proxy for fish exposure to hypoxia (Limburg et al. 2011, 2015). During the summer monsoon, dissolved oxygen concentration is lower in the NW Indian Ocean than Mozambique Channel (Fig. S3). Thermal vents, sediments, rivers and atmospheric dust are also known sources of manganese in seawater (van Hulten et al. 2017). The influence of aeolian deposition of Mn is greater in the NW Indian Ocean because of its proximity to the arid Arabian subcontinent (van Hulten et al. 2017). Therefore, changes in wind patterns may be another plausible cause of observed differences in Mn : Ca profiles of the two regions. Life history otolith chemistry profiles would be valuable for resolving stock structure of the species (e.g. Fowler et al. 2017), but there is still a considerable need for further empirical studies regarding elemental incorporation processes into otoliths for a comprehensive interpretation of the results.

The otolith near-core portion showed greater variability in trace elements among the two regions than did the core. Higher levels of Zn and Cu were found in core portions compared with the near-core (Fig. S7, available as Supplementary material to this paper). It is possible that some thermoplastic glue or resin residues had entered the pits and crevices occasioned in the primordium of some samples due to polishing (Di Franco et al. 2014). It is also possible that the chemical signature of this zone, which is formed during egg development, could have been affected by maternal inheritance (e.g. Campbell et al. 2002; Ruttenberg et al. 2005; Thorrold et al. 2006). Therefore, this signature can be affected by the environment that is experienced by the maturing female (Volk et al. 2000), hindering its use as an indicator of fish origin. Given these circumstances, trace element signatures from the primordium may not be appropriate for determining YOY yellowfin tuna origin, and the use of otolith material from a few weeks after hatch could represent a better alternative for early life history comparisons (see also Macdonald et al. 2008).

Region-specific chemical signatures were detected in the otolith near-core composition, with three elements (Ba, Mg and Mn) showing regional differences in otolith near-core chemistry of YOY yellowfin tuna. As stated above, in the NW Indian Ocean, upwelling events occurring during the south-west (summer) monsoon generate an upward movement of nutrient-rich waters from below the thermocline to shallower depths, which possibly implies the observed higher Ba : Ca concentrations in otolith near-core compared with YOY from the Mozambique Channel. YOY from the NW Indian Ocean also had higher concentrations of Mn : Ca in the otolith near-core than YOY from the Mozambique Channel, possibly due to their proximity to an arid region or an oxygen minimum zone (McCreary et al. 2013). Differences in Mg : Ca ratios in otoliths of fish from the Mozambique Channel and NW Indian Ocean may also reflect physiological differences rather than environmental differences (Woodcock et al. 2012). Nonetheless, each of the above elements differed significantly between groups, and they may be useful as tracers of nursery origin regardless of the cause of variation between regions.

Owing to the lack of difference in otolith isotope composition between the two regions in the first 5–6 months of life, δ¹³C and δ¹⁸O seem to be of negligible value for investigation of YOY yellowfin tuna origin in the western Indian Ocean. This is likely explained by the lack of variation in oxygen isotopic composition in the seawater of this region (LeGrande and Schmidt 2006).

In this study, RF was the best method to predict natal origin, with highest accuracy and kappa index. When using otolith microchemistry for identification of fish origin, the choice of the optimal classification method is important in order to avoid misinterpreting the results. Traditional approaches (LDA and QDA) have been commonly used as classification tools in other tuna otolith chemistry studies (e.g. Rooker et al. 2008b, 2016; Wells et al. 2012; Macdonald et al. 2013). However, in recent years, machine learning techniques have emerged as promising classification tools in otolith-related studies (Zhang et al. 2016; Tournois et al. 2017; Bouchoucha et al. 2018). The accuracy of each technique will depend on the nature of the data analysed, and the results of this study agree with those of recent studies encouraging the use of machine learning methods when otolith chemical data are not multivariate normal or exhibit skewed distributions (Mercier et al. 2011; Jones et al. 2017). The reasonably easy use of RF, together with a lack of distributional assumption requirements and the robustness of RF against overfitting (Breiman 2001), make this technique a useful approach for population structure analyses of yellowfin tuna in the Indian Ocean.

Given the high migratory behaviour and dispersal capacity of yellowfin tuna (Fonteneau and Hallier 2015), obtaining a sufficiently representative baseline at such a large scale would have been logistically unfeasible in the present study. Alternatively, using unsupervised learning methods to cluster larval elemental signatures is well suited for studies where no representative baseline samples are available (Gibb et al. 2017). In the present study, the RF clustering approach used on the near-core otolith portion inferred the presence of a cluster for age-1 yellowfin tuna of unknown origin. This implies either that the samples share a common origin, with a common larval source contributing to the samples, or a lack of difference in elemental composition between different source locations. Although the approach does not directly provide information on the location of this potential common larval source, a canonical covariates plot revealed higher similarity with NW Indian Ocean YOY signatures than with Mozambique Channel YOY. The sizes of the age-1 individuals in this study (i.e. 52–64 cm FL) were of a similar range as those from the RTTP-IO (i.e. 56–75 cm FL), which reported large movements and a low average rate of homing fidelity in yellowfin tuna of the Indian Ocean (Fonteneau and Hallier 2015). Therefore, the presence of different larval sources would have not been surprising in the age-1 sample composition. However, it has also been shown that distances travelled by yellowfin tuna captured around Seychelles Islands (NW Indian Ocean) were lower than for tuna from other areas (Hallier and Million 2012). This could possibly explain the origin and retention of age-1 samples in the region. This would highlight the importance of local production of yellowfin tuna recruits in the NW Indian Ocean fishery. A high degree of local production and retention of yellowfin tuna has also been reported in other oceans (Rooker et al. 2016; Wells et al. 2012).

This study supports the utility of otolith chemical analyses to provide valuable information about the stock structure of yellowfin tuna in the western Indian Ocean and suggests a more complex stock structure in the area than previously assumed. However, future work should delve into analysis at a whole-ocean scale, preferably sampling distinct regions at the same period. Moreover, further development and application of this approach is now needed, in particular by quantifying spatial and temporal variation in oceanic water chemistry at appropriate scales so that this information can help the interpretation of otolith chemistry data. Otolith chemistry analysis is one of a suite of techniques that can be used to study stock structure in marine fishes. Integrating the information gained from otolith chemistry analysis with complementary techniques, such as genetics, tagging and morphometrics (e.g. Abaunza et al. 2008; Taillebois et al. 2017), will provide a powerful means to define appropriate management units for the sustainable exploitation and conservation objectives of yellowfin tuna in the Indian Ocean.

The authors declare that they have no conflicts of interest.

This work was supported by a grant from the Department of Agriculture, Fishing and Alimentation from the Basque Government to Iraide Artetxe-Arrate for her doctoral studies.