Otolith chemical fingerprints of skipjack tuna (Katsuwonus pelamis) in the Indian Ocean: First insights into stock structure delineation

The chemical composition of otoliths (earbones) can provide valuable information about stock structure and connectivity patterns among marine fish. For that, chemical signatures must be sufficiently distinct to allow accurate classification of an unknown fish to their area of origin. Here we have examined the suitability of otolith microchemistry as a tool to better understand the spatial dynamics of skipjack tuna (Katsuwonus pelamis), a highly valuable commercial species for which uncertainties remain regarding its stock structure in the Indian Ocean. For this aim, we have compared the early life otolith chemical composition of young-of-the-year (<6 months) skipjack tuna captured from the three main nursery areas of the equatorial Indian Ocean (West, Central and East). Elemental (Li:Ca, Sr:Ca, Ba:Ca, Mg:Ca and Mn:Ca) and stable isotopic (δ13C, δ18O) signatures were used, from individuals captured in 2018 and 2019. Otolith Sr:Ca, Ba:Ca, Mg:Ca and δ18O significantly differed among fish from different nurseries, but, in general, the chemical signatures of the three nursery areas largely overlapped. Multivariate analyses of otolith chemical signatures revealed low geographic separation among Central and Eastern nurseries, achieving a maximum overall random forest cross validated classification success of 51%. Cohort effect on otolith trace element signatures was also detected, indicating that variations in chemical signatures associated with seasonal changes in oceanographic conditions must be well understood, particularly for species with several reproductive peaks throughout the year. Otolith microchemistry in conjunction with other techniques (e.g., genetics, particle tracking) should be further investigated to resolve skipjack stock structure, which will ultimately contribute to the sustainable management of this stock in the Indian Ocean.

Introduction Skipjack tuna (Katsuwonus pelamis) is a cosmopolitan species inhabiting tropical and subtropical waters of the Indian, Pacific and Atlantic Oceans [1]. This species is, by far, the most commonly caught tuna worldwide [2,3]. Currently, skipjack tuna stocks are considered to be in a healthy status in all oceans, both in terms of stock abundance and fishing mortality [3]. Five stocks of skipjack tuna are considered globally for management, two in the Atlantic Ocean, two in the Pacific Ocean and one in the Indian Ocean. Regional studies indicate a more complex stock structure than currently assumed for assessment and management purposes of skipjack tuna in the Indian Ocean [4,5]; however, lack of understanding of population dynamics and connectivity at oceanic scale do not allow separation of stocks [6].
Stock structure understanding is essential to determine a suitable spatial scale for management, as the way a stock will respond to management decisions cannot be accurately predicted if the boundaries that characterize a stock are not correctly defined [7]. Skipjack tuna spawns throughout the year in large warm water tropical areas, but as those areas are often poor feeding grounds, they move towards adjacent colder and more productive subtropical areas to feed [8][9][10]. However, not all skipjack tuna seem to perform these long-distance movements, but rather, some remain around the tropical spawning area [10,11]. As a result, the spatial dynamics of skipjack tuna and, hence, its implications on the stock structure of the species needs to be evaluated. Besides, skipjack tuna are fast growing, early maturing species, and have high reproductive potential, which makes this species more resilient to fishing pressure than many other tuna species [12][13][14]. However, fishing pressure may vary spatially due to factors such as fleet operational constrains, variations in population densities, and differences in the value of fish caught, among others [15]. Thus, excessive local fishing effort and catches, may lead to situations of local overfishing and depletion when skipjack movements between areas are limited [16]. The understanding of exchange rates and connectivity among recruits from different regions is therefore, essential to achieve a sustainable and effective management of this species [17].
There are several methods that can be used to study fish stock structure, which can provide information at different spatial and/or temporal scales [18]. Among them, analyses of otolith (i.e., calcified structures found in the inner ear of the fish) chemical composition is widely used to explore movements and habitat use of fish [19]. This method relies on the premise that ambient water chemistry and environmental conditions (but also other intrinsic factors such as physiology, diet or genetics) affect the elemental incorporation (at minor and trace quantities) into the concentric growth of the otolith during daily increment formation [20,21]. Since otoliths are both acellular and metabolically inert, material accreted during otolith formation is preserved as fish grows [22]. As such, the chemical composition of the material accreted during early life stages, serves as a natural marker to identify fish that have inhabited environments with distinct physicochemical characteristics [23]. Otolith microchemistry proved to be a powerful tool to discriminate among nursery fish groups for other tropical tuna species [24][25][26], and also to depict different movements and life history patters of skipjack tuna in the Pacific Ocean [27].
Here, we examined otoliths of young-of-the-year (YOY) skipjack tuna collected across the equatorial strip of the Indian Ocean, where high larvae concentrations and spawning activity of this species have been observed; off Seychelles, Somalia and Mozambique Channel in the western Indian Ocean, off Maldives in the central Indian Ocean, and off Sumatra in the eastern Indian Ocean [9,28,29]. Besides, those areas are also important fishery grounds of skipjack tuna in the Indian Ocean [6]. Our aim was to determine whether YOY skipjack tuna captured within these three regions could be spatially discriminated based on their early life otolith microchemistry composition. If so, the characterized nursery origin signature can then be used as a baseline sample to predict older skipjack tuna origin, by analyzing the same early life portion of the otolith [30,31]. This in turn, will contribute to improve our understanding of the connectivity and mixing rates and, ultimately, the stock structure of this species in the Indian Ocean.

Fish sampling
Skipjack tuna were collected from three distinct nursery areas in the Indian Ocean: West (10˚S-10˚N, 40˚E-60˚E), Central (0˚-10˚N, 65˚E-75˚E), and East (5˚S-10˚N, 85˚E-95˚E) (Fig 1). Samples were obtained directly by scientist or scientific observers on-board purse seine vessels or at port during two consecutive years (2018 and 2019), as part of a collaborative research project on Population Structure of Tuna, Billfish and Sharks of the Indian Ocean [32]. Fish size ranged between 24.5 and 35.0 cm fork length (FL) (Tables 1 and 2), and were estimated to be about 4-6 months old following the age-length relationship described by Eveson et al. [33]. Thus, we assumed that capture locations represent nursery areas where fish reside during the early juvenile stage. Note that fish were sampled at different years/time periods, which implies

Otolith preparation and analyses
Sagittal otoliths were extracted, cleaned of adhering organic tissue, rinsed with ultrapure water (Milli-Q), and stored dry in plastic vials. In the laboratory, one otolith of each specimen was embedded in two-part epoxy resin (Araldite 2020, Huntsman Advanced Materials, Switzerland). Blocks were polished using 3M 1 silicon carbide sandpaper (particle size = 220 μm) and a lapping wheel with a series of decreasing grain diameter (30, 15, 9, 3 and 1 μm) 3M 1 silicon carbide lapping discs, moistened with ultrapure water, to obtain a transverse section where the core was exposed. Sections were ultrasonically cleaned using ultrapure (Milli-Q) water for 10 minutes. Following sonication, otolith sections were left to air dry in loosely capped vials for 24 h before being glued in a sample plate using Crystalbond thermoplastic glue (Crystalbond 509; Buehler). When a single otolith was available only trace element analyses were conducted. When both pairs of otoliths were available, the second one was used for carbon and oxygen stable isotope analyses.
Trace element analyses. Otoliths (n = 128) were analysed for trace element chemistry using a high resolution inductively coupled plasma mass spectrometer (HR-ICP-MS, Element XR, Thermo Scientific, Bremen, Germany), coupled to a high repetition rate 1030 nm femtosecond laser (fs-LA) system (Alfamet, Neseya, Canejan, France) available at the Institut des Sci  direct daily microincrement counts, and thus to exclude the potential maternal effects the primordium (i.e., the initial structure of an otolith) may incorporate [34]. Laser ablation conditions were 200 Hz, a laser beam size of 15 μm and 30 μJ per pulse energy corresponding to a fluence of 14 J/cm2, until the depth limit of ablation (< 30 μm). During ablation the small laser beam was continuously and rapidly moved (500 μm s-1) at the surface of the sample due to a 2D galvanometric scanner in order to create a trajectory made of 6 concentric circles [35]. This resulted in the ablation of a crater 30 μm in diameter and 30 μm deep. The ablation cell was flushed with argon to transport laser-induced particles to the HR-ICP-MS. The fsLA-H-R-ICP-MS was tuned daily to reach optimal particle atomization conditions and minimal elemental fractionation. This was obtained for a U/Th signal ratio of 1 ± 0,05 using NIST 612. The mass spectrometer was used in the medium-resolution mode (R = 4000) to ensure a complete polyatomic interference removal for the isotopes of interest. Relative abundances of five isotopes ( 7 Li, 88 Sr, 138 Ba, 24 Mg and 55 Mn) were estimated, as well as 43 Ca, which was used as the internal standard. The concentration of 43 Ca in the otolith was assumed to be constant at 38.3% [36]. Data reduction including background subtraction, conversion to ppm and standardization to calcium (element:Ca μmol mol-1) was done using an in-lab developed software Stable isotopes analyses. For stable isotope analyses (n = 56), microsampling of otolith powder for carbon (δ 13 C) and oxygen (δ 18 O) stable isotope analysis was performed using a high-resolution computerised micromill (New Wave MicroMill System, NewWave Research G. C. Co., Ltd, Cambs, UK). The area of analysis on the smallest skipjack tuna otolith (24.5 cm FL) was used to create a standard template that was then applied to the remaining otoliths, to ensure that the same portion of the otolith was analyzed in every fish (S2 Fig). Therefore, the drill path covered the area of the otolith corresponding to the first~4 months of life (according to Eveson et al. [33] age-length relationship), with a larger time period of the otolith sampled for stable isotopes than trace elements due to differences in sample material requirements. Ten drill passes were run at a depth of 50 μm per pass over a preprogrammed drill path using a 300μm diameter carbide bit (Komet dental; Gebr. Basseler, Lemgo, Germany). Powdered material was then analysed for δ 13 C and δ 18 O on an automated carbonate preparation device (KIEL-III, Thermo-Fisher 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. The isotope ratio measurement was calibrated based on repeated measurements of NBS-19 and NBS-18 (International Atomic Energy Agency standards). Measurement precision was ± 0.10 ‰ for δ 18 O and ±0.08‰ for δ 13 C (1 sigma).

Statistical analyses
All statistical analyses were performed using open access R software [37]. Prior to all multivariate analyses, otolith chemistry data was standardized (i.e., for each element, the data was centered by subtracting the mean and scaled by dividing by the standard deviation) to give the same weight to all elements and stable isotopes.
Elements were first analyzed individually. Parametric assumptions were violated by Ba:Ca, Mg:Ca and Mn:Ca. Therefore, non-parametric tests were used for trace element analyses, to apply a consistent approach for all elements. To determine whether the otolith chemistry of YOY skipjack tuna varied spatially and/or temporally, a permutational multivariate analysis of variance (PERMANOVA) was used to test for differences in element:Ca ratios between nurseries and cohorts [38,39]. Nursery and cohort were fixed factors in the full factorial model. The resemblance matrix was based on euclidean distance dissimilarities and the number of unrestricted permutations was set to 999 random repeats using adonis {vegan}. Statistical significance was determined based on adjusted p-values after the Benjamini-Hochberg correction [40], and post hoc pairwise comparisons were applied to identify the source of differentiation between group means using lincon {WRS2}. Parametric assumptions of normality and homoscedasticity were met for δ 13 C and δ 18 O data, and therefore one-way ANOVA test was used for comparisons among nurseries, aov {stats}. When significant differences were detected, post hoc comparison were performed to the source of differences between means, using TukeyHSD {stats}.
Both trace elements and stable isotopes were them combined for multivariate analyses. The resultant dataset was limited by the number of individuals for which both types of data were available (n = 56). A multilevel pairwise comparison was performed using pairwise.adonis {pairwiseAdonis} to test for differences in the multielemental signature between nursery areas. The resemblance matrix was based on euclidean distance dissimilarities and the number of unrestricted permutations was set to 999 random repeats. Statistical significance was determined based on adjusted p-values after the Benjamini-Hochberg correction [40]. Then, multivariate data were reduced to two-dimensions and visualized by a canonical analysis of principal coordinates (CAP) using CAPdiscrim {BiodiversityR} function [41]. Random forest (RF) classification algorithm (number of trees = 500, mtry = 2) was implemented to test the ability of trace elements and stable isotope signatures to discriminate among fish belonging to different nursery areas [42]. Data was split into a training dataset (70%) and a testing dataset (30%), and this procedure was randomly repeated 1000 times. At each time, the rate of classification success (i.e., rate of correct predicted membership to nursery areas in which the fish were collected) was calculated, and mean values were extracted. Cohen´s Kappa (κ) statistic was also calculated, which is a method that accounts chance-corrected percentage of agreement between actual and predicted group memberships. Values of κ range from 0 to 1, where 0 indicates that the RF resulted in no improvement over chance, and 1 indicates perfect agreement [43]. Random forest was performed for each cohort, as well as for pooled data from both cohorts in the case of trace element data. For stable isotope data, and combined trace element and stable isotope data, pooled data from both cohorts was only used.

Individual elements
Significant differences were detected in the chemical signatures of skipjack tuna collected from different nurseries, but also between cohorts in the case of trace element data ( Table 3). The concentrations of Sr and Mn differed significantly between years, concentrations being higher in 2018 (Fig 2). Significant interactions between nursery and cohort were also detected for Ba, Mg and Mn, meaning that the pattern of variation was different among nurseries at each cohort (Table 3, Fig 2). Concentrations of Sr, Ba and Mg differed among nurseries in 2017, being lower for skipjack tuna collected from the West nursery (Fig 2). For fish belonging to the 2018 cohort no significant differences across regions were detected for any of the trace elements analyzed.
Otolith δ 18 O values also varied among nurseries (anova, p<0.001), with significantly higher values found in the West nursery with respect to the central and eastern nurseries (Fig 3). No significant differences were detected in otolith δ 13 C composition between nursery areas.

Multielemental signatures
When individual trace elements and carbon and oxygen stable isotopes were combined, significant regional differences were found among nursery areas (PERMANOVA, p = 0.005). Specifically, multielemental signatures from YOY skipjack tuna from the western nursery were differentiated from the central and eastern nurseries, as supported by posterior pairwise comparisons (Table 4, and Fig 4).
Classification accuracy was generally low regardless of the elemental combination used. Classification success using cohort pooled cohort trace element data was of 44% and κ = 0.16 (Table 5). Overall classification success of fish back to their nursery area was higher in 2017 (50%, κ = 0.33) than in 2018 (44%, κ = 0.15). Nursery specific signatures also varied considerably between cohorts. For YOY skipjack tuna from the West nursery classification accuracy decreased from 71% in 2017 to 49% in 2018, while for the Central nursery increased from 23% in 2017 to 40% in 2018. Cohort pooled classification accuracy of YOY skipjack tuna from East nursery remained similar (45% and 42% in 2017 and 2018, respectively).
Carbon and oxygen stable isotopes provided higher classification success than trace elements alone (51%, κ = 0.44, Table 5). The combination of trace elements and stable isotopes did not improve the overall classification accuracy with respect to the use of carbon and stable isotopes alone (51%, κ = 0.40). However, nursery specific classification did vary, YOY from the West nursery were better classified with stable isotopes alone, while for YOY from the Central nursery greater classification accuracy was obtained when stable isotope data was combined with trace element data (Table 5, Fig 4). Nevertheless, samples from the Central and East nurseries were often misclassified (<50%), regardless on the elemental combination used.

Discussion
Otolith chemical fingerprints are powerful discriminators of groups as long as differences exist, but of negligible value when differences cannot be detected, as the absence of difference does not necessarily imply that groups own a common origin [20]. The high overlap in the otolith chemical signatures of YOY from the central and eastern nurseries, may be at least partly explained by the homogeneity in the physicochemical properties of ambient water fish were exposed to during their early life. Depending on the period, the regions of Maldives (Central) and eastern Indian Ocean are quite homogeneous in terms of sea surface temperature (SST), salinity (SSS), and dissolved oxygen (DO) [44,45]. Another potential reason to the observed overlap of the chemical signatures, is that important transoceanic mixing may have occurred prior to capture. However, it is unlikely that transoceanic migrations occur before 4 to 6 months of life. To our knowledge, transoceanic movements of skipjack tuna within the first 6 months have not been documented, since tagging studies are generally not available for these size ranges [11]. Moreover, skipjack tuna predominantly shows diffusion type movements rather than long longitudinal migrations, as they move towards colder and richer subtropical areas adjacent to their spawning areas, where they can temporarily found higher productivity and prey availability than in warm equatorial waters [8,10]. Besides, the submesoscale activity present in the nursery areas can minimize long-distance dispersal, by retaining the larvae near their spawning area [46,47]. For all the above, we think that the observed results in the chemical signatures of YOY skipjack tuna are more likely attributed to the fact that biochemical properties of the water masses were relatively homogeneous during the early life of the fish, leading to a low discrimination capacity among fish captured in different nursery areas. Strong interannual effect was detected in the otolith trace element signatures of fish belonging to two different cohorts. While significant differences were detected among nurseries in Sr, Ba and Mg concentrations in 2017, differences were not significant in 2018. Barium and strontium are incorporated into the otolith via Ca substitution, and may be appropriate markers of environmental history [48][49][50]. Mg incorporation into the otolith by contrast, has been related with metabolic activity [51,52]. Due to sampling constrains, YOY skipjack tuna were hatched during different year periods (S1 Fig). North of 10˚S the Indian Ocean is characterized by two seasons with distinct wind regimes, the monsoon system, which drives the ocean circulation and climate in the tropics and northern subtropic region of the Indian Ocean [53,54]. The southwest monsoon, also known as summer monsoon, takes place from June to October with a maximum of intensity during the months of July, August, and September. During this season there is little difference in rainfall intensity within the equatorial region, while temperatures are colder west of 60˚E due to the coastal upwelling that takes place off Somalia and Oman Coasts [44,45,55]. The northeast Monsoon, or the winter monsoon, takes place between December and April with a maximum of intensity in the months of January, February, and March. Upwelling events take place off northwest Australian shelf and in open ocean at 5-10˚S region [53]. During the winter monsoon waters off Somalia are colder (26.5˚C), whereas SSTs are above 28˚C over the rest of the equatorial region [45]. Besides, during this season rainfall intensity is stronger in the west south (4.5˚S) of the equator off Madagascar and off Sumatra [55]. During the transition periods between the monsoons (i.e., May and November), a strong equatorial downwelling occurs. For skipjack tuna belonging to the 2017 cohort, those collected from the West nursery were estimated to be hatched during the summer monsoon, while those from Central and East nurseries were considered to be hatched during the winter monsoon and the autumn intermonsoon periods. All skipjack tuna belonging to the 2018 cohort were estimated to be hatched within the summer monsoon period and the two intermonsoons. As seasonal climatology strongly influences the biochemical variability throughout the northern Indian Ocean, it may be possible that observed regional differences in Sr, Ba and Mg in 2017 are due to the seasonal changes in water physicochemical properties that could have masked regional differences [56,57].
Otolith δ 18 O signatures were significantly higher from individuals of the West nursery than the Central and the East nurseries. The isotopic composition of ambient seawater is dependent on evaporation, and largely controlled by salinity [58]. As otolith δ 18 O is incorporated into otolith aragonite close to equilibrium with seawater, the isotopic composition of otolith reflects seawater δ 18 O values, being inversely correlated with ambient seawater temperature [59,60]. Otolith composition of YOY skipjack tuna from the Indian Ocean followed the expected trend for δ 18 O signatures, presenting higher values in the West (expected due to lower water temperatures, under the influence of the seasonal Somali upwelling) and decreasing towards the east (expected due to higher water temperatures, as confirmed by the overall sea surface temperature pattern in the Indian Ocean). Similar trends have been found in a preliminary study for yellowfin tuna from the Indian Ocean [61]. No differences were detected in otolith δ 13 C signatures between fish captured in the three different nursery areas. Otolith δ 13 C can be a proxy for metabolic rate in fish, and is influenced by the environment (dissolved inorganic carbon in water, DIC) and diet [62,63]. The observed absence of difference in δ 13 C signatures followed the expected trend for surface water carbon isotope composition, which is relatively homogenous within the equatorial Indian Ocean [64], and therefore providing little value for stock discrimination.
The overall low regional variability in chemical signatures from skipjack tuna resulted in a low classification of YOY individuals to the three nursery areas analyzed in the Indian Ocean. To date, there are no studies aiming to retrospectively determine skipjack tuna to their nursery origin based on otolith chemical signatures. However, studies done with other tuna species inhabiting tropical environments (e.g. yellowfin and bigeye) reported higher classification accuracies of fish back to their nursery area than those reported in this study [24][25][26]65]. The use of stable isotopes alone was sufficient for nursery discrimination of yellowfin tuna in the Pacific Ocean, whereas trace elements proved to be more efficient for discriminating among yellowfin nursery areas in the Atlantic and Pacific and Indian Oceans [24,26,65]. Here, the use of carbon and oxygen stable isotope values alone, proved to be effective to differentiate skipjack tuna from the western nursery, but the high overlap in the chemical signature of skipjack tuna from the Central and East nurseries resulted in an overall low classification. The use of trace element data or the combination of stable isotope and trace element data did not improve the overall classification success. It is also possible that differences in classification success between trace elemental ratios and stable isotopes is due to the different the time period of the otolith represented by each tracer (i.e.~day 13 to 15 for trace elements vs. 4-6 months for stable isotopes) [24]. One possible solution would be to analyze with LA-ICPMS the integrated signal of an ablation square that corresponds to the same otolith portion milled for stable isotopes analyses (i.e.~4-6 months of life). However, this data should be interpreted with caution as it has been proved that ontogeny strongly influences trace element uptake into the otolith [52,66].

Conclusions
Overall, the differences in the chemical signature of YOY skipjack tuna from the three nursery areas in the Indian Ocean were not strong enough to describe a reference baseline for each nursery that allows the assignment of older individuals to their origin nursery. Temporal differences in the physicochemical characteristics of the Indian Ocean seem to strongly influence otolith trace element composition. Considering also that classification success for models based on otolith trace element data compared to those using only otolith stable isotopes was lower, the effort and cost to incorporate elemental ratios into Indian Ocean skipjack tuna otolith chemistry baseline delineation may not be worthwhile [25]. Until there is a proper understanding of the stock structure of skipjack tuna, the uncertainty in relation to the response of this stock to management decisions adopted considering a single stock will be maintained. Further research on skipjack population structure using otolith microchemistry should rely on younger (<4 months) individuals to reduce the possibility of movements between nursery areas and consider temporal stratification of sampling so that seasonal differences in oceanography can be disentangled from potential regional differences. Finally, an holistic approach may provide a more accurate overview to properly define management units, as single technical approaches may not be sufficient to delineate complex stock structures [67,68]. Recent studies combining different stock delineation techniques (e.g. otolith microchemistry, genetics, biophysical models. . .), proved to be more effective to fully understand the spatial ecology of highly migratory fish species, and, hence increased the resolving power of stock discrimination [69][70][71]. For instance, including passive drift trajectory simulations may also be helpful to understand potential patterns of skipjack larval dispersal in the Indian Ocean [72]. The combination of otolith chemistry data coupled with genetic analyses can increased details on connectivity patterns, as both techniques provide information on complementary timescales (individual for otolith chemistry, and evolutionary for genetics), unravelling otherwise hidden patterns [73]. Ultimately, a better understanding of the stock structure and spatial connectivity of skipjack tuna in the Indian Ocean will be essential to implement and enforce management strategies that ensure long-term sustainable fisheries of this important species.