Shared fish stocks, namely fish resources with high migration ability that cross exclusive economic zone (EEZ) boundaries, are more prone to overexploitation than other stocks (e.g. McWhinnie, 2009; Kraska et al., 2015; Popova et al., 2019). Therefore, cooperative stock management among stakeholders is a necessary prerequisite for achieving sustainable fisheries of these species (e.g. Caddy, 1997; Miller and Munro, 2004). Appropriate understanding of migration routes and their individual and temporal variation is crucial for better cooperative management of shared fish stocks. Knowledge of a migration route is also essential for determining which countries share fishery resources from the same stock, as well as for assessing the effect of changes in the migration routes on the stock and fishing conditions in each country (Miller and Munro, 2004; Vatsov, 2016). Recently developed models (e.g. Carvalho et al., 2015; Sippel et al., 2015) can make use of migration data for stock assessment, thus demonstrating the benefits of understanding the migration routes of these fish.

In fishery science, one of the main approaches for studying fish migration is dynamic tracking using electronic technologies such as telemetry, radio, and computer networks (e.g. Priede et al., 1988; Thorstad et al., 2013). Although these methods have excellent spatial and temporal resolution, there are several intrinsic drawbacks, such as the high cost of bio-logging equipment (Hebblewhite and Haydon, 2010; Graves and Horodysky, 2015), the large effort required for capturing and tagging the fish, and body-size limitations for tagged fish: conventionally, the recommended tag weight is less than 3% of the fish body mass (e.g. Winter, 1983; Jepsen et al., 2002; Cooke et al., 2011). Stock assessment models are generally based on population-level data, but bio-logging methods only yield individual-level inferences, making them difficult to incorporate into overall stock condition assessments. In addition, equipment attached to a target fish can change its behavior, potentially resulting in biased migration data, although the effect is highly context-dependent and species-specific (Jepsen et al., 2015).

To overcome these drawbacks of bio-logging and make use of migration data for stock assessment of shared fish stocks, other approaches with lower cost and less effort are warranted. Recently, migration tracking techniques using isotope analysis of incremental-growth tissues such as otoliths, eye lenses, and vertebral bones (Matsubayashi et al., 2017; Sakamoto et al., 2019; Matsubayashi et al., 2020; Vecchio and Peebles, 2020) have proven promising for inferring the migration routes of individual fish. The isotopic composition of local environment and an animal’s tissue can be used as a natural tag to track its movements through isotopically distinct habitats (Graham et al., 2010), and we call this technique “iso-logging”. The major advantages of iso-logging over bio-logging are lower cost and effort, longer period of tracking (Vecchio and Peebles, 2020), and no influence on the behavior of target fish before death, although the temporal resolution and accuracy of estimated ranges are not equal to those from bio-logging.

In the case of marine fish, the ¹⁸O:¹⁶O ratio (δ ¹⁸O) in otoliths is the most frequently used isotope ratio for iso-logging (e.g. Trueman et al., 2012; Hsieh et al., 2019; Sakamoto et al., 2019). δ ¹⁸O can be retrospectively measured from otoliths and reflects the ambient temperature experienced by the fish. However, δ ¹⁸O is not effective for iso-logging of fish that 1) migrate in oceans with similar thermal gradients, 2) have a stable body temperature, like tuna (Thunnus spp.; Carey and Teal, 1969), or 3) often dive deep, crossing the thermocline. Skipjack tuna Katsuwonus pelamis comprises one of the most commercially important shared fish stocks inhabiting warm and temperate waters. Skipjack in western Pacific populations prefer sea-surface temperatures ranging from 20.5 to 26.0°C (Mugo et al., 2010) and can maintain its body temperature higher than that of ambient water (Barrett and Hester, 1964; Stevens et al., 1974), which hampers the use of otolith δ ¹⁸O for migration tracking. For this reason, it is necessary to test the applicability of other isotope ratios for tracking skipjack migration. Recent studies have suggested that stable carbon and nitrogen isotope ratios (δ ¹³C and δ ¹⁵N) can be excellent indicators for the migration histories of marine animals in some areas (Trueman et al., 2017; Matsubayashi et al., 2020), and these isotopes may be also applicable for tracking skipjack in the western Pacific (Coletto et al., 2021).

To use δ ¹³C and δ ¹⁵N for iso-logging of marine fish, at least three types of information are required: 1) the isoscapes across the habitat, 2) any ontogenetic shifts of isotope ratios of target fish, and 3) isotopic offsets between isoscapes and target fish. Unfortunately, it is generally difficult to obtain all of this information for δ ¹³C and δ ¹⁵N. The generation of isoscapes requires thorough field sampling of proxy organisms with low swimming ability (Trueman et al., 2017; Matsubayashi et al., 2020). Understanding isotopic offsets between isoscape and target animal requires feeding experiments over several months (e.g. Matsubayashi et al., 2019; Canseco et al., 2021). Furthermore, in the case of highly migratory piscivorous species, understanding the ontogenetic shifts of δ ¹³C and δ ¹⁵N is particularly difficult because isoscapes and ontogenetic dietary shifts have mutual interactions. To overcome these difficulties in iso-logging, we conceived of using a geostatistical model to simultaneously estimate isoscapes and ontogenetic shifts of isotope ratios using the fish targeted for migration tracking as a proxy organism for isoscapes. With this approach, isotopic offsets between isoscapes and target fish can be ignored. Although such modeling was technically difficult up until a few years ago, recent advances in geostatistical analysis (e.g. Rue et al., 2009; Thorson, 2019; Anderson et al., 2022) allow us to account for these factors in a single model with a user-friendly interface.

Here, we test the utility of geostatistical analysis for generating isoscapes using δ ¹³C and δ ¹⁵N of skipjack in the western Pacific, considering ontogenetic shifts in isotope ratios of this species. In this region, skipjack is mainly distributed from 20°S to 40°N and from 110°E to 180°E (Bartoo, 1987; Kiyofuji and Ochi, 2016), and some individuals exhibit dynamic northward migration from the subtropics into waters off the Pacific coast of Japan (Kiyofuji et al., 2019a). The water temperature range preferred by skipjack is from 18.8°C to 28.1°C and the 18.0 isotherm is the lowest thermal limit of skipjack tuna distribution (Kiyofuji et al., 2019a). We sampled the muscle tissue of skipjack of various body sizes and determined the δ ¹³C and δ ¹⁵N values of bulk tissue (δ ¹³C_Bulk and δ ¹⁵N_Bulk). We then applied a geostatistical model to account for variations in the δ ¹³C and δ ¹⁵N of skipjack muscle with fork length, several environmental variables, and a spatial random field. Furthermore, we performed compound-specific nitrogen isotope analysis of amino acids (δ ¹⁵N_AAs) for five samples with different δ ¹⁵N_Bulk values to evaluate whether the differences in bulk values resulted from a baseline shift in isotope ratios or from trophic differences (see Matsubayashi et al., 2020). On the basis of these data, we discuss the applicability of δ ¹³C_Bulk and δ ¹⁵N_Bulk for assessing the migration of skipjack in the western Pacific.

For generation of isoscapes, previous studies have used species with lower mobility like jellyfish and copepods as proxy organisms (Trueman et al., 2017; Matsubayashi et al., 2020). The advantage of using other species with low swimming ability to create isoscape is that location-specific isotope ratios can be detected with high representativeness. On the other hand, the ontogenetic shift of isotope ratios in target fish and the offset of isotope ratios between target fish and isoscapes should be estimated separately in this method. This study used a target fish as a proxy organism for isoscapes, and successfully estimated both isoscapes ( Figure 4 ) and the ontogenetic shifts of isotope ratios ( Figure 3 ) of skipjack without compound-specific isotope analysis, through the use of geostatistical models. Furthermore, we can apply iso-logging of skipjack without considering the isotopic offset between isoscapes and target organism, because the isoscapes were generated based on the isotope ratios of the skipjack. Thus, our approach can be a useful alternative to the previous method, although our method has one drawback in that the uncertainty of isoscapes will increase with the number of fish captured that had migrated across regions with different isoscapes in the several months before they were caught. We addressed this issue by increasing the number of adult fish at each sampling point and period ( Figure 1 ), and this may have contributed somewhat to decreasing the uncertainty of isoscapes given the high prediction ability of our models. Another solution would be to use juvenile fish only for generating isoscapes and then to construct another model to account for the ontogenetic shift of isotope ratios. However, this is not applicable when the distribution of juvenile and adult fish does not completely overlap, as with skipjack in the western Pacific where juvenile fish are scarce in the northern area (Lehodey et al., 2013). To minimize the uncertainty associated with a high migration ability of proxy organisms, it is important to use tissues with a high turnover rate. As of now, the fish tissue with the shortest stable isotope turnover time is epidermal mucus (Winter and Britton, 2021), but this is difficult to sample from juvenile fish smaller than 100 mm. Therefore, it is important to further explore turnover rates in other tissues to discover better ways to sample proxy organisms for generation of isoscapes.

Together with the results of the geostatistical models, our δ ¹⁵N_AAs data also show that the variation in isotope ratios of skipjack in the western Pacific was generally derived from the baseline shift of the isoscape and not from differences in prey items. Adult skipjack in the western Pacific are mainly piscivorous and occasionally consume cephalopods, with low selectivity for their prey (Iizuka et al., 1989). Therefore, regional shifts in prey abundance can cause regional differences in δ ¹⁵N_Bulk. If the observed fluctuations of δ ¹³C_Bulk and δ ¹⁵N_Bulk in skipjack muscle resulted from differences in trophic position, the regional trend of δ ¹⁵N_Bulk should be similar to that of TP_Glu–Phe ( Figure 5 ). However, our δ ¹⁵N_AAs data show that skipjack in the western Pacific have small individual variation in TP_Glu–Phe (maximum difference: 0.4). Assuming that δ ¹⁵N_Bulk increases 3.4‰ with trophic level (Minagawa and Wada, 1984), differences in TP_Glu–Phe can only explain 1.4‰ of the observed variation of δ ¹⁵N_Bulk out of 10.2‰ ( Table 2 , Figure 5 ). Given that the δ ¹⁵N_Phe, which generally has a small trophic discrimination factor (about 0.4%; Chikaraishi et al., 2009), was more correlated with δ ¹⁵N_Bulk of skipjack muscle than TP_Glu–Phe, individual differences in δ ¹⁵N_Bulk of adult skipjack largely reflect the variation of δ ¹⁵N at the base of the food web.

TP_Glu–Phe of skipjack estimated by commonly used TDF_Glu-Phe (7.6 ‰) were somewhat lower (2.7–3.1) than the trophic position estimated by stomach-content analysis (TP_SCA = 3.8; Bradley et al., 2015). Although the optimal TDF_Glu–Phe to use for estimating the trophic position of skipjack remains controversial (e.g., Bradley et al., 2015; McMahon and McCarthy, 2016; Ohkouchi et al., 2017; Coletto et al., 2022), differences in TP_Glu–Phe of skipjack do not affect intraspecific variation in TDF_Glu–Phe, and thus would not affect our conclusion that regional variation in TDF_Glu–Phe are unlikely to explain δ ¹⁵N_Bulk differences among regions.

Regional differences in δ ¹⁵N of skipjack in the western Pacific were largely consistent with the δ ¹⁵N distribution of nitrate in this region (Rafter et al., 2019). Low δ ¹⁵N around 10°N to 35°N ( Figure 4 ) is related to extensive N₂ fixation in this region (Gruber, 2016), because N₂ fixers typically exhibit lower δ ¹⁵N values than non-N₂ fixers (McClelland et al., 2003). On the other hand, our isoscape suggests that the low-δ ¹⁵N area extends as far as 40°N, which is not consistent with known distributions of N₂ fixers or low nitrate δ ¹⁵N (Gruber, 2016; Rafter et al., 2019). Presumably, this inconsistency is attributable to the relatively longer turnover time of skipjack muscle. The area around 35–40°N is the northern limit of skipjack migration and they can stay there at most four months a year (from June to September; Mugo et al., 2011; Kiyofuji et al., 2019a). Given that the turnover time of adult skipjack white muscle is about four months, isotope ratios of the muscle of skipjack captured around 35–40°N are likely to reflect the low δ ¹⁵N values south of 35°N, and thus our isoscape shows low δ ¹⁵N up to 40°N. On the other hand, the high δ ¹⁵N south of 10–15°N can be explained by high nitrate utilization and limited N₂ fixation. In the surface layer of the tropical Pacific Ocean, a strong halocline above the thermocline generates a barrier layer, which impedes vertical mixing and nutrient supply from subsurface to surface waters (Sprintall & Tomczak, 1992). The limited nutrient influx results in high nitrate utilization (i.e. low nitrate concentration)—the ratio of nitrate assimilation by phytoplankton to nitrate supply in the euphotic layer—and nitrate δ ¹⁵N increases with nitrate depletion (Yoshikawa et al., 2018). Furthermore, N₂ fixation is limited in this region (Gruber, 2016; Tanita et al., 2021) mainly due to the depleted dissolved iron concentration (Tanita et al., 2021). Such biogeochemical factors result in a clear latitudinal trend in δ ¹⁵N, rendering it useful for segregating the habitat use of skipjack in the western Pacific.

Interpretation of spatial variation in δ ¹³C is much more difficult than that of δ ¹⁵N. The spatial variations in δ ¹³C of skipjack muscle should be ultimately determined by that of phytoplankton (δ ¹³C_phy), which is strongly influenced by the δ ¹³C of dissolved inorganic carbon (DIC) and the isotopic fractionation associated with carbon uptake (Tuerena et al., 2019). In the case of the western Pacific, there are no clear latitudinal differences in δ ¹³C of DIC (Quay et al., 2003), and therefore the spatial variation in the δ ¹³C isoscape would be driven by the uptake fractionation by phytoplankton. Until the early 1990s, the δ ¹³C of phytoplankton was thought to be determined by the DIC concentration in the ambient water (Sackett et al., 1965; Rau et al., 1991). Our geostatistical model also predicted an inverse relationship between DIC and δ ¹³C ( Supplementary Figure 2 ), and the belt-shaped low-δ ¹³C zone around the equator ( Figure 4 ) may be mainly due to the high spCO₂ in this region ( Supplementary Figure 3 ). However, there is a large uncertainty associated with the equatorial low-δ ¹³C zone and further sampling in this region would provide important clarification ( Figure 4 ). On the other hand, the clear δ ¹³C gradient around 20°N is not likely explained by spatial variation in spCO₂. Recent studies have shown that local differences in phytoplankton species composition can be important controls on δ ¹³C_phy (e.g. Francois et al., 1993; Bidigare et al., 1997; Popp et al., 1998; Tuerena et al., 2019) as well as the DIC concentration, specifically in regions where the DIC is not high and is less variable, as in the western Pacific ( Supplementary Figure 3 ). Thus, we presume that local differences in environmental variables altered the phytoplankton species composition, and ultimately generated spatial variation in δ ¹³C of skipjack muscle in this region. It is worth noting that anthropogenic CO₂ emission can be another factor influencing the spatial variation in δ ¹³C_phy (Tuerena et al., 2019).

In this study, we demonstrated that δ ¹³C and δ ¹⁵N measurements of skipjack in the western Pacific can help to identify their latitudinal migration. Because the dorsal muscle, which is a major tissue used for isotope analysis of fishes, only represents isotope ratios at a specific point in time, the utility of isoscapes will be maximized when they are used along with the isotopic analysis of incremental growth tissues such as eye lenses and vertebral bone (e.g. Quaeck-Davies et al., 2018; Matsubayashi et al., 2019). Iso-logging using such incremental growth tissue would allow us to understand how skipjack in this region utilize different habitats though their lifetime without using electronic tags. For example, if time series isotopic data for eye lens are available, the migration pattern can be tracked by the following steps: i) calculate the fork length when each layer is formed from the equation relating lens diameter to fork length, ii) using the equation relating fork length to isotopic ratios ( Figure 3 ), correct the isotope ratios of eye lens so that the fork length is at zero; iii) subtract Δδ ¹⁵N_{muscle - eye lens} from the value obtained in ii), iv) corrected lens isotope ratios can then be compared to the isoscape ( Figure 4 ). We also suggest that the isoscapes generated in this study can be applied to other fish species that migrate in the western Pacific Ocean, such as albacore Thunnus alalunga (Nikolic et al., 2017), Pacific bluefin tuna Thunnus orientalis (Fujioka et al., 2018), and Japanese eel Anguilla japonica (Aoyama et al., 2014). For these other target species, however, it will be necessary to apply an appropriate correction to the trophic discrimination factor for differences between skipjack muscle and specific tissues of these other species. By using this method, we can better understand the immigration/emigration rates and connectivity in each region, which is crucial for cooperative management of shared fish stocks like skipjack. Further studies to generate precise isoscapes of δ ¹³C and δ ¹⁵N based on ocean models hold promise for improving the spatiotemporal resolution of isotope tracking and, together with the results of this study, will facilitate its use.

The original contributions presented in the study are included in the article/ Supplementary Material . Further inquiries can be directed to the corresponding author.

JM and HM conceived this study. HM and YT provided skipjack samples and contextual information. HM, YC, NO, KK, NFI, and NOO performed bulk and compound-specific stable isotope analyses. JM performed statistical analyses and wrote the first draft of the paper. All authors contributed to the article and approved the submitted version.

This study was partially supported by the Foundation for the Japan Society for the Promotion of Science (JSPS) KAKENHI Grants 19H04247, 21H03579, 21H04784, and 22H05028 to JM, and 20H00208 to NO, and the Research and Assessment Program for Fisheries Resources of the Fisheries Agency of Japan.

This study was conducted with the support of the Joint Research Grant for the Environmental Isotope Study of the Research Institute for Humanity and Nature. We thank H. Ashida, S. Ohashi, F. Tanaka, and H. Kiyofuji at the Fisheries Resources Institute, Japan Fisheries Research and Education Agency, for their help with sample collection. We thank C. Yoshikawa at the Japan Agency for Marine-Earth Science and Technology for her helpful comments on a draft of this manuscript. We also thank the reviewers for their prompt and careful review of our draft and for their helpful comments.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.