Combining microvolume isotope analysis and numerical simulation to reproduce fish migration history

Tracking the movement of migratory fish is of great importance for efficient conservation, although this has been technically difficult to achieve in small fish to which artificial tags cannot be attached. We show that migration history can be reproduced by combining high‐resolution otolith stable oxygen isotope ratio (δ18O) analysis and numerical simulation. High‐precision micromilling and microvolume carbonate analysing systems had the remarkable capability of extracting the otolith δ18O profiles with 10–30 days resolution. Furthermore, reasonable movements were reproduced by searching the routes consistent with the otolith δ18O profile, using an individual‐based model with random swimming behaviour. This method will be a valuable alternative to tagging and electronic loggers for revealing migration routes in early life stages, thereby providing crucial information to understand population structures and the environmental cause of recruitment variabilities, and to validate and improve fish movement models.


| INTRODUC TI ON
Marine fish stocks are an important food resource for humans (Crist, Mora, & Engelman, 2017) and marine predators (Crawford, Sabarros, Fairweather, Underhill, & Wolfaardt, 2008;Cury et al., 2011), although many have faced severe declines (Pauly et al., 2002). Understanding population structure is a basic and crucial requirement for efficient fisheries management (Cadrin, Karr, & Mariani, 2014) because, without it, fisheries can exert an unexpectedly high pressure on certain regional populations, leading to over-fishing of some stock components (Felix-Uraga et al., 2005;Hüssy et al., 2016). Therefore, revealing seasonal or annual changes in fish distribution at both individual and population levels has been a key topic in fisheries science (Block et al., 2005;Brennan et al., 2015;Tsukamoto & Nakai, 1998).
Several methods have been developed to track the movements of marine fish in open oceans. Electronic tagging has been valuable for reproducing the migration histories of large fish such as adult bluefin tuna based on recorded light levels and temperatures (Block et al., 2005). For early life stages in which swimming ability is negligible, numerical particle tracking experiments based on ocean models have allowed a better understanding of how fish disperse from spawning grounds (Allain, Petitgas, Lazure, & Grellier, 2007;Itoh et al., 2011). However, we still lack solid methods for small-sized, but actively swimming, juveniles that cannot be tagged, although survival rate during this stage is crucial for populations fluctuation (Anderson, 1988;Houde, 1987). Recently, numerous individual-based migration models coupled with fish bioenergetics have been developed to estimate migrations of this type of fish (Okunishi et al., 2012;Rice et al., 1993;Rose, Rutherford, McDermot, Forney, & Mills, 1999), although they still include assumptions that are very difficult to verify empirically, thereby presenting inherent uncertainty.
Otoliths are calcium carbonate crystals formed in the inner ear of fish, and their chemical composition is known to record ambient water chemistry (Campana, 1999;Sturrock, Trueman, Darnaude, & Hunter, 2012). In particular, the oxygen stable isotope ratio (δ 18 O), which reflects both temperature and δ 18 O of surrounding water (Kim, O'Neil, Hillaire-Marcel, & Mucci, 2007), has often been used as a temperature proxy to discuss fish migration (e.g., Carpenter, Erickson, & Holland, 2003;Shiao, Yui, Høie, Ninnemann, & Chang, 2009). It has mainly been applied to fish species with large otoliths because conventional isotope ratio mass spectrometry required large sample amounts. Recently, however, the analysing system MICAL3c has been developed. This system requires 1/100 of the sample quantity required for conventional mass spectrometry (Ishimura, Tsunogai, & Nakagawa, 2008;Nishida & Ishimura, 2017), thereby allowing an accurate analysis of narrow growth increments in small otoliths (Sakamoto et al., 2017). Despite this drastic technical improvement, interpreting otolith δ 18 O values to estimate fish migration history remains difficult. The variability in seawater δ 18 O, which is primarily related to salinity variation (Craig & Gordon, 1965), often confounds temperature estimation because it is essentially impossible to estimate the two parameters in a single equation. Moreover, because otolith δ 18 O can only distinguish water masses of different temperature or salinity, estimations of fish location become impractically rough when uniform water masses are distributed broadly (Torniainen et al., 2017), which is often the case in open oceans.
Here, we describe and demonstrate a new methodology to estimate the migration history of individual small fish, by using an interdisciplinary combination of isotopic analysis and numerical simulation.
The Japanese sardine was selected as the model for this demonstration because of its dynamic migration from the temperate Kuroshio region to the subarctic Oyashio region within 6 months after hatching (Kuroda, 1991), and economic and ecological importance in the western North Pacific (Ishida et al., 2009). To our knowledge, this is the first otolith microchemistry study to use numerical migration models to interpret the isotope profile. First, the distribution of the seawater δ 18 O and its relation with salinity in the habitat area of the Japanese sardine were investigated. Using micromilling and microvolume isotope analysis, accurate otolith δ 18 O profiles were obtained with high resolution. Furthermore, instead of estimating temperature histories from the otolith δ 18 O profiles, we searched for migration routes that should be followed to reproduce the δ 18 O profile, using an individual-based model in which a randomly swimming migration scheme was applied to a dataassimilated ocean model.  00-47.48°E, 138.48-172.38°N region, which roughly covers the whole distribution of larval and juvenile Japanese sardine. In all sampling points, the salinity profile was measured with a SBE 9 Plus CTD unit (Sea-Bird Electronics, Inc.), with precision better than 0.01 in all cruises. Seawater samples for δ 18 O analysis were taken from the sea surface using a roped plastic bucket on board, or at 5 m depth using CTD-attached Niskin bottles (General Oceanics, Inc.). Water samples were preserved in sealed glass vials to avoid evaporation. Samples obtained in 2012 and 2013 were measured at the Geo Science laboratory using the LWIA DLT-100 system (Los Gatos Research) with precision ±0.05‰.

| Distribution of seawater δ 18 O and its relationship with salinity
Samples from other cruises were analysed at the Atmosphere and Ocean Research Institute (AORI), University of Tokyo, using the Picarro L2120-i system with precision ±0.07‰. Samples analysed in AORI were membrane-filtered (pore size: 0.45 μm, Toyo Roshi Kaisha, Ltd.) before their introduction in the Picarro L2120-i system to reduce suspended particles and avoid blocking the sampling line.
The correlation between seawater δ 18 O and salinity was analysed by the least squares method using surface water δ 18 O and salinity measured at 10 m (or 5 m if available) depth by the CTD.

| Fish sampling and otolith micromilling and δ 18 O analyses
Six young-of-the-year Japanese sardine individuals were used. The position and number of daily rings in otoliths were examined along the axis of the postrostrum using an otolith measurement system (RATOC System Engineering Co. Ltd.). The daily age at capture was calculated by adding two to the total number of otolith daily rings (Takahashi, Nishida, Yatsu, & Watanabe, 2008), because the first daily ring is formed 2-3 days after hatching in Japanese sardine (Hayashi, Yamashita, Kawaguchi, & Ishii, 1989). At every 10 or 15 daily rings, a ring was tracked as far as possible from the axis and used as the boundary of the micromilling area. Date ranges corresponding to each milling area were back calculated from the capture date by subtracting the number of daily increments from the edge of the otolith. Marked pictures were input to the high-precision micromilling system GEOMILL 326 that comprises a micromill, a CMOS camera, a video monitor, and an image analyser controlled by a computer. This system allows accurate sampling at 1/1,000 mm scales (Sakai, 2009), and it has been used to accurately extract the edge of sardine otoliths (Oda, Tetsu, Sakai, & Ishimura, 2016;Sakamoto et al., 2017). Each area between marked daily rings was then sequentially extracted, and the resulting powders (0.3-11.4 μg) were collected into stainless microcups (YUKO-PARTS). After each milling, the otolith was cleaned with an air-duster to avoid cross-contamination between milling paths. The micromilling depth was set at 50 μm for the two areas nearest the core and 50-120 μm for the remaining areas.
The δ 18 O values were reported in δ-notation relative to the VPDB (Vienna Pee Dee Belemnite) scale, and given as a ‰ value. Analytical precision was better than ±0.10‰ throughout the entire analysis.
The acid fractionation factor of calcite was used to facilitate comparisons with isotopic values reported in previous studies (Amano, Shiao, Ishimura, Yokouchi, & Shirai, 2015). Although we could not find errors in the purifying and analysing processes, a sharp increase was observed in the δ 18 O profile (* marked in Figure 1); this value was excluded from the modelling analysis and considered as a failure of isotope analysis because we could not find any migration route that reproduced such a fluctuation.

| Estimation of migration history
We applied an individual-based model (IBM) with a randomswimming scheme to an ocean environment reproduced by the data-assimilated ocean model FRA-ROMS. FRA-ROMS is an ocean forecast and reanalysing system developed to realistically simulate mesoscale variations over the Kuroshio-Oyashio region (Kuroda et al., 2016). A random swimming scheme was chosen to find many, ideally all, possible routes that the fish can take to reach capture location, simply based on the distance that can be covered by advection and swimming. It would be interesting to add assumptions, such the fish preference to head the direction where they can yield higher growth rate, to further narrow route prediction. However, we avoided introducing such an assumption because the mechanism by which fish decide their heading direction is currently unknown and we wished to free our model from assumptions that are not empirically confirmed. For each of the six individual fish analysed, approximately 10 9 simulated individuals were released on the hatching date in the survey-based spawning grounds, and advected and diffused by a surface current field reproduced using FRA-ROMS. Individuals were set to swim at a speed of 3*SL/s after metamorphosis, randomly changing their swimming direction once a day. Otolith δ 18 O histories were calculated for the individuals that reached the capture point, based on ambient temperature and salinity histories. The possibility of each migration route was calculated as a product of egg abundance in the grid of origin and agreement between calculated otolith δ 18 O history and analysed profiles. Possibility-weighted mean position and temperature and water current velocity were calculated daily and considered representatives of estimated routes.

| FRA-ROMS configuration
The FRA-ROMS ocean forecast system was developed by the Japan where the Japanese sardine is distributed (Supporting Information Figure S1).

| IBM simulation design and route selection
The survey-based monthly egg abundance data (Oozeki, Takasuka, Kubota, & Barange, 2007), 1/4° grid scale for individuals from 2010 and 1/2° for individuals from 2014, were used to set the starting points of particle tracking. In the grid cells in which egg abundance was positive in the month that targeted individuals hatched, 120 × 120 particles were laid at regular intervals on the hatch date.
The number of egg-positive grid cells in the southern coast of Japan was 37 and 33 for 2010 and 2014, respectively. The horizontal movement of particles was regulated with advection, diffusion, and active swimming terms using the Euler scheme: where (x n , y n ) is the position of a particle at time step n, (u n , v n ) is the surface current velocity at (x n , y n ) interpolated linearly from the surface velocity reanalysed by FRA-ROMS, and (l x n , l y n ) is the random walk steps representing horizontal diffusion parameterized by the scheme of Smagorinsky (1963). The term (mig x n, mig y n ) corresponds to the active swimming velocity at time step n, which was updated daily as.
where C is a constant coefficient to estimate cruising speed from body length, BSL is a back calculated SL on that day, and θ is a random number between 0 and 1. Back calculated SL was based on the linear relationship between SL and otolith radius using a biological intercept method (Campana, 1990;Takahashi et al., 2008). Particles were tracked until the capture day and their position, ambient temperature, and salinity were output daily. This trial was repeated 2,000 times with different random number sequences. Hara (1987) reported that the cruising speed of the adult Japanese sardine is 1.2-4.1 body length (BL)/s. Although low swimming ability of sardine larvae would be negligible (Silva, Faria, Teodósio, & Garrido, 2014), we assumed that metamorphosed juveniles, compared to adults, swim equally or slightly faster in proportion to their body size because smaller sized fish have generally higher tail beat frequency (e.g., Bainbridge, 1958;Hunter & Zweifel, 1971). Therefore, C was set to 0 when BSL is under 30 mm, and changed to 3 when BSL exceeded 30 mm. Additional experiments were conducted with C for BSL > 30 mm equal to 2 or 4, aiming to test the sensitivity of the results to the value of C, using 1,000 replicates. As the scheme is essentially a random walk, the scale of dispersal is larger for smaller frequencies of direction changes.
To check the sensitivity of the results to the frequency of direction change, further experiments with C equal to 3 were conducted with frequency equal to 2 and 0.5 times/day, using 1,000 replicates.
The possibility of each route was calculated as follows. First, the particles with horizontal position out of the ±0.5° range from the captured location on the capture date were regarded as unrelated and excluded from further analysis. For particles that reached the capture location, egg abundance at the grid of origin was set as the initial weight. Calculation of the δ 18 O profile was conducted using averaged ambient temperature and salinity of the particle during the date corresponding to each milled otolith area. Salinity was converted into seawater δ 18 O using a regression between seawater δ 18 O and salinity in the habitat of the Japanese sardine (Equation 1, see Section 3). The overall possibility for each route was calculated as the product of initial weight and agreement for each date range. Routes with possibilities higher than 1/100 of the highest possibility were considered possible routes. For these possible routes, possibility-weighted mean positions, ambient water temperature and water current velocity were calculated daily to discuss the details of selected routes.  & Schmidt, 2006) or to the near shore Kuroshio-Oyashio system (Oba & Murayama, 2004 Figure 3b,d,f) took more offshore routes, being advected as far as ~165°E before heading north.

| Migration routes estimated using randomly swimming IBM
These routes seemed to be reasonable as, for most individuals analysed, the possibility-weighted mean routes from May to June, corresponding to 40-100 days post hatch, coincided with the area where late-larvae and juvenile sardines were caught during sampling surveys conducted in the corresponding period (Figure 4a,b).
Mean temperature and water current velocity were also calculated to examine key timings and features of the migration. Sardines were F I G U R E 3 High possibility routes for each individual according to random swimming model estimation. The thicker and darker coloured routes are the routes with higher possibilities. The possibility of each route was divided by the highest possibility for each individual. Capture locations are presented as green stars. The left upper panel within each picture shows calculated otolith δ 18 O history for each migration route (yellow-black) and the analysed otolith δ 18 O (green) initially transported along the Kuroshio until departing to the north in late May to mid-June (Figure 4c,d). They entered the Oyashio feeding ground dominated by subarctic water (surface water temperature <15°C) from the end of June to early July (Figure 4e,f), and reached the highest latitude (47-48°N) in late August (Figure 4c,d).
The substantial contribution of the northward water current for sardines reaching the feeding ground was suggested by the positive northward speed of this current in May and June (Figure 4g,h).
Selected routes were similar when swimming speed was changed, but moved slightly inshore with C = 2 (Supporting Information Figure S3). Similarly, selected routes were robust to changes in the frequency of swimming direction, but shifted inshore with high frequency of these changes (Supporting Information Figure S4).

| D ISCUSS I ON
Here, we reproduced the migration history of the small pelagic Japanese sardine using the combination of high-resolution otolith is initiated controls the period of larvae advection by the Kuroshio, thereby determining the transported distance. Furthermore, when fish have low swimming ability but need to move north, they require strong northward cross-jet currents. In the Kuroshio Extension (KE), a cross-frontal northward current is frequently observed between the trough and the crest of the meander (Sainz-Trápaga & Sugimoto, 2000), as is the case of the Gulf Stream (Bower & Rossby, 1989).
Warm streamers are also spread intermittently from the KE (Sainz-Trápaga & Sugimoto, 2000) and quasi-stationary jets extend northeastward in the transitional area between the KE and the Oyashio fronts (Isoguchi, Kawamura, & Oka, 2006), probably providing conditions for longitudinal selection. Although these processes seem unique for the KE, similar processes may occur in other regions as many fish larvae depend on regional flows, such as the coastal jets in the Benguela Current (Hutchings et al., 2002) or the eddies in the California Current (Logerwell, Lavaniegos, & Smith, 2001), to effectively reach food-rich feeding grounds.
The method used here revealed new aspects of feeding migration of the Japanese sardine. Our estimations suggested that it is from the middle to end of May that sardines start to move northward departing from KE, and it is from the end of June to early July that they enter the Oyashio feeding ground north of 42°N. Although the spawning grounds and transportation during the larval stage have been well documented through a number of surveys (e.g., Kuroda, 1991), these key timings in the feeding migration are the first to be clarified because this type of information can only be obtained by continuously tracking fish movement. The estimation provided here also showed that sardines utilise the northward current during their northward migration, supporting the hypothesis of Isoguchi et al. (2006), who pointed out that warm tongues and meanders from KE might be the direct routes for the northward migration of pelagic fishes based on their analyses of surface flow field. Considering all individuals analysed in this study are the survivors that reached the feeding ground, riding the northward current near KE might be an important condition for the sardine to recruit and thus the strength of such currents can significantly affect the recruitment rate.
Between the 2 years analysed, the selected routes in 2014 were more offshore than those in 2010. This shift may be related to the seven-times larger population size in 2014 (Yukami et al., 2017), as sardines tend to expand their habitat area offshore when biomass is increased (Barange et al., 2009). All these detailed descriptions of migration patterns would be useful for narrowing down the season and areas that are important for recruitment variability, studying how fish decide their swimming direction, and testing and improving migration models, which will be the scopes of future studies.
It should be noted that the accuracy of route estimation using such processes largely relies on the reproducibility of the ocean model. The model we used was developed to realistically simulate the environment by assimilating numerous satellite and in situ observation data, although the number of observations is still limited, especially for salinity, as this is not measured by satellites. Further accurate estimation requires both expansion of the monitoring network and improvement of simulation models. It is also noteworthy that model schemes require some assumptions on the swimming depth. Because age-0 Japanese sardine vertically migrate mainly within the surface mixed layer (e.g., Takahashi  assumed they experience an environment similar to that on the water surface. When applying this methodology to species distributed in wider depth range, some adjustments, such as averaging the temperature and salinity within the water column, or introducing vertical migrations into IBM, may be needed. In addition, due to spatially and temporally dense observations of Japanese sardine spawning, these data were used as starting points for IBM estimation, but this is not always the case. In the absence of such survey datasets, modelled spawning grounds can be used (e.g., Zwolinski, Emmett, & Demer, 2011).
Overall, the method presented here would be a great alternative to electronic tags understanding movements of fish during early life history or for fish that are otherwise too small to accommodate electronic tags. It would also be useful not only to reveal population structure but also to validate and improve movement models that have been lacking detailed reference data to discuss their accuracy.
Moreover, this method has the unique and strong advantage of examining the environmental history of successfully recruited individuals. This is crucial information to understand the environmental conditions necessary for fish to survive, providing clues on how environmental variabilities are driving fish population fluctuations. Because otolith δ 18 O is generally recognized as temperature dependent (Høie, Otterlei, & Folkvord, 2004;Kim et al., 2007;Kitagawa et al., 2013;Sakamoto et al., 2017;Storm-Suke, Dempson, Reist, & Power, 2007), the method presented here will notably improve the knowledge on the survival and migration ecology of early life stages of numerous fish species.

ACK N OWLED G EM ENTS
The authors gratefully acknowledge the support from Grant-in-Aid for JSPS Fellows 17J00556, JSPS KAKENHI grants number JP15H04541, 16H02944, and 18H04921, and from grants for marine fisheries stock assessment and evaluation for Japanese waters from the Fisheries Agency and Fisheries Research and Education Agency of Japan. We thank S. Itoh for technical advices on numerical simulation and insightful comments to improve the manuscript.
We thank Mr. Zhang who provided the seawater δ 18 O data in 2012 and 2013. There are no conflicts of interest to declare.

AUTH O R S' CO NTR I B UTI O N S
T.S. and K.K. designed and conceived the study. T.S., T.H., and T.I. performed isotope analyses. T.S. and K.K. wrote the program and performed the numerical simulations. T.Setou provided ocean model data and validated its accuracy. Y.K., C.W., and A.K. provided fish samples. All authors contributed to manuscript editing. T.S., K.K., and K.S. wrote the manuscript.