# Global patterns and inferences of tuna movements and trophodynamics from stable isotope analysis

J. M. Logan, H. Pethybridge, Anne Lorrain, C. J. Somes, V. Allain, N. Bodin, C. A. Choy, L. Duffy, N. Goni, B. Graham, et al.

► To cite this version:

J. M. Logan, H. Pethybridge, Anne Lorrain, C. J. Somes, V. Allain, et al. Global patterns and inferences of tuna movements and trophodynamics from stable isotope analysis. Deep Sea Research Part II: Topical Studies in Oceanography, 2020, 175, pp.104775. 10.1016/j.dsr2.2020.104775. hal-02987362

HAL Id: hal-02987362

https://hal.science/hal-02987362v1

Submitted on 5 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Global patterns and inferences of tuna movements and trophodynamics

Logan J.M. <fnref n="1" />,*, Pethybridge Heidi <fnref n="2" />, Lorrain Anne <fnref n="3" />, Somes C. <fnref n="4" />, Allain Valerie <fnref n="5" />, Bodin Nathalie <fnref n="6" />,
Choy C.A. <fnref n="7" />, Duffy L. <fnref n="8" />, Goñi N. <fnref n="9" />, Graham B. <fnref n="10" />, Langlais C. <fnref n="1" />, Ménard F. <fnref n="11" />, Olson R. <fnref n="8" />, Young J. <fnref n="2" />

## Abstract :

A global dataset of carbon stable isotope (δ13C) values from yellowfin, bigeye, and albacore tuna muscle tissue (n = 4275) was used to develop a novel tool to infer broad-scale movement and residency patterns of these highly mobile marine predators. This tool was coupled with environmental models and lipid content (C:N ratio) of tuna muscle tissues to examine ocean warming impacts on tuna ecology and bioenergetic condition across Longhurst provinces. Over a 16-year study period (2000–2015), latitudinal gradients in tuna δ13C values were consistent, with values decreasing with increasing latitude. Tuna δ13C values, reflecting modelled global phytoplankton δ13C landscapes (“isoscapes”), were largely related to spatial changes in oxygen concentrations at depth and temporal changes in sea surface temperature. Observed tuna isoscapes (δ13CLScorr), corrected for lipid content and the Suess effect (oceanic changes in CO2 over time), were subtracted from model-predicted baseline isoscapes (Δ13Ctuna-phyto) to infer spatial movement and residency patterns of the different tuna species. Stable isotope niche width was calculated for each Longhurst province using Δ13Ctuna-phyto and baseline-corrected nitrogen isotope (δ15Ntuna-phyto) values to further quantify isotopic variability as evidence of movements across isoscapes. A high degree of movement—defined as the deviation from expected range of Δ13Ctuna-phyto values—was evident in four Longhurst provinces: Guinea current coast, North Atlantic drift, Pacific equatorial divergence, and the North Pacific equatorial counter current. The highest level of population dispersal (variability in Δ13Ctuna-phyto values) was observed in Longhurst provinces within the western and central Pacific Oceans and in the Guinea current coast. While lipid content was low in yellowfin and bigeye, high and variable lipid stores in albacore muscle were consistent with seasonal movements between productive foraging and oligotrophic spawning habitats. Our ability to characterize tuna movement patterns without ambiguity remains challenged by uncertainty in trophic discrimination

<footnote n="1">¹ Massachusetts Division of Marine Fisheries, New Bedford, MA, 02744, USA</footnote>

<footnote n="2">² CSIRO Oceans & Atmosphere, GPO Box 158, Hobart, Tasmania, Australia</footnote>

<footnote n="3">³ IRD, Univ Brest, CNRS, Ifremer, LEMAR, F-29280, Plouzané, France</footnote>

<footnote n="4">⁴ GEOMAR Helmholtz Centre for Ocean Research Kiel, Kiel, Germany</footnote>

<footnote n="5">⁵ Pacific Community (SPC), Nouméa, New Caledonia</footnote>

<footnote n="6">⁶ Seychelles Fishing Authority (SFA), Fishing Port, Victoria, Mahé Island, Seychelles</footnote>

<footnote n="7">⁷ Scripps Institution of Oceanography, University of California, San Diego, Integrative Oceanography Division, La Jolla, CA, 92037, USA</footnote>

<footnote n="8">⁸ Inter-American Tropical Tuna Commission (IATTC), La Jolla, CA, 92037, USA</footnote>

<footnote n="9">⁹ AZTI-Tecnalia/Marine Research, Herrera Kaiaportualdea Z/g, 20110, Pasaia, Gipuzkoa, Spain</footnote>

<footnote n="10">¹⁰ National Institute of Water and Atmospheric Research, Ltd. (NIWA), Wellington, 6021, New Zealand</footnote>

<footnote n="11">¹¹ Aix Marseille Univ, Université de Toulon, CNRS, IRD, MIO, 13288, Marseille, France</footnote>

<footnote>* Corresponding author : J. M. Logan, email address : john.logan@mass.gov</footnote>
factors and ecological (e.g. diet variability) processes. However, this study illustrates that model-corrected δ13C values are a valuable, relatively cost-effective tool for identifying potential areas of mixing across management zones, particularly when electronic tagging studies are limited or absent. Stable isotope analyses of tuna tissues can therefore be an additional tool for guiding spatial stock assessments on top predator movement, dispersal patterns, and how they may be altered under a changing climate.

**Keywords** : Albacore tuna, Bigeye tuna, Yellowfin tuna, Carbon isotope analysis, Movement indicators, Global ocean, Fisheries management
59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

# 1. Introduction

Understanding how and why highly mobile species utilise and move between available habitats in time and space is vital to determining how successful stocks are managed (Evans et al., 2015). In recent years there has been considerable effort to empirically measure and forecast where and when highly mobile species are likely to occur. To date, this has largely been achieved through expensive tagging or bio-logging programs (Block et al., 2011) and habitat preference models that use environmental proxies (Eveson et al., 2015; Hobday et al., 2011; Lehodey et al., 2018; Sibert et al., 1999). Another important aspect in fisheries management and conservation biology is assessing population connectivity, including how much dispersal and mixing occurs among populations. Such information is currently acquired through extensive demographic datasets or using genetics, with the recent development of close-kin mark-recapture methods (Bravington et al., 2016). Many of these tools however still have foundational and logistical constraints (e.g. high cost, require large sample sizes, and limited number of experts and/or users). In contrast, a long-standing and widely utilized tool—analysis of spatial patterns of carbon stable isotope ($\delta^{13}C$) values in top predator tissues—is a potential complementary approach to assist in understanding both movement patterns and habitat use in marine predators (Bird et al., 2018a; Graham et al., 2010; Trueman et al., 2012; Trueman and St. John Glew, 2019).

Carbon stable isotope analysis has been used for over forty years to trace movements in estuarine and marine environments of a variety of consumers including crustaceans (Fry, 1981; Fry et al., 2003), cephalopods (Takai et al., 2000), seabirds (Jaeger et al., 2010), pinnipeds (Newsome et al., 2007), whales (Best and Schell, 1996; Schell et al., 1989), teleosts (Ménard et al., 2007), and elasmobranchs (Bird et al., 2018a). Early studies examined shrimp and whale movements based on tissue $\delta^{13}C$ values (Best and Schell, 1996; Fry, 1981; Schell et al., 1989). Shrimp movement detection was facilitated by the presence of primary producers with distinct $\delta^{13}C$ values at each foraging habitat while whale movements were inferred by spatial variability
in phytoplankton δ¹³C values. At broad oceanic scales, δ¹³C values of dissolved inorganic carbon (δ¹³C_DIC) are influenced by temperature (Laws et al., 1995), oceanic partial pressure of CO₂ (pCO₂), phytoplankton growth, and phytoplankton community structure (Fry and Wainright, 1991; Hinga et al., 1994; Popp et al., 1998). These variations in baseline δ¹³C values are driven by spatial differences in primary producers or abiotic conditions, which propagate up to top predators and provide a chemical tracer of movements across foraging grounds with different baseline values (isoscapes sensu West et al., 2010). Further studies have demonstrated success in linking isotopic variation to individual-level dietary differences and movement patterns (Araújo et al., 2007; Grecian et al., 2015; Shipley et al., 2019). Individual isotopic variability, quantifiable as isotopic niches based on scatter in bivariate isotope space (Jackson et al., 2011), can provide a relative measure of trophic and/or movement diversity within a group.

Before carbon isotope techniques can be applied to examine the feeding and behavioural ecology of consumers, spatiotemporal and biochemical factors need to be considered. Direct comparison of isotope values of consumers occupying different isoscapes requires some form of prior adjustment for spatial and temporal baseline variability. Over long time periods, plankton δ¹³C values are known to decline (typically by 0 – 0.08‰ per year) due to the Suess effect (ocean absorption of fossil fuel-derived CO₂) (Quay et al., 2017), which needs to be accounted for in temporal studies. Spatial baseline adjustments can be applied using site-specific modelled phytoplankton isotope values (e.g. Navarro et al., 2013; Pethybridge et al., 2018), normalization using a secondary consumer of a known trophic position (e.g. Olson et al., 2010; Post, 2002), and compound specific stable isotope analyses (e.g. Popp et al., 2007). Within fish muscle tissue, bulk δ¹³C values can also be influenced by lipid content, as lipids have a lower δ¹³C value than proteins (Abelson and Hoering, 1961; Park and Epstein, 1961). Lipid effects can be accounted for with chemical extractions (e.g. Pinnegar and Polunin, 1999) or a posteriori mathematical correction using carbon to nitrogen (C:N) ratios as lipid-content proxies (e.g. Logan et al., 2008).

Trophic influences on δ¹³C values also need to be accounted for in carbon isotope analyses. Biochemical processes alter δ¹³C values between a consumer and its prey (trophic discrimination factor; TDF), and this results in an increase by approximately 0.5–2.0‰ (Caut et al., 2009; DeNiro and Epstein, 1978; McCutchan et al., 2003; Post, 2002; Sweeting et al., 2007). Controlled feeding studies show that carbon TDFs for fish muscle are generally at the higher end of this range with estimates of 1.5‰ found across a variety of species (Sweeting et al., 2007) and 1.8‰ for *Thunnus* muscle specifically (Madigan et al., 2012). Cephalopods, an important prey group of tunas (Olson et al., 2016), are less studied but existing controlled studies showed limited (<0.5‰) increases in mantle δ¹³C values relative to diet (Hobson and Cherel, 2006). Estimates for coastal and open ocean zooplankton fall closer to 1‰ (France and Peters, 1997). Therefore, if the trophic position of yellowfin tuna (*Thunnus albacares*) in the equatorial Pacific Ocean is observed to be approximately 4.5‰ (Houssard et al., 2017; Lorrain et al., 2015; Popp et al., 2007), then the difference between the base of the food web and tuna δ¹³C values (Δ¹³C_tuna-phyto) should be approximately 2 to 8‰ given documented variability in TDFs (0.5-2.0‰ across taxa linking phytoplankton to tunas) and food chain length from tunas down through to lower trophic level species (Fig. 1). Variations from these values and/or a high degree of individual variability could be indicators of highly dynamic migration areas.

Most previous studies that have used stable isotopes to examine the movements of highly mobile marine animals have focused on regional scales and higher latitudes where regional gradients in δ¹³C values are more apparent (Cherel and Hobson, 2007; Cunjak et al., 2005; Hobson, 1999; Rubenstein and Hobson, 2004; Witteveen et al., 2009). A recent analysis of shark
muscle δ¹³C values examined latitudinal trends in relation to modelled phytoplankton δ¹³C values to establish global trends in shelf, oceanic, and deep sea trophic geography of sharks (Bird et al., 2018a). Our study applied a similar approach to tunas and focused primarily on the tropics and subtropics to examine animal movements, potential environmental drivers, and global-scale trends. Using a compiled dataset of δ¹³C values and C:N (a proxy for bioenergetic condition) values of yellowfin (*Thunnus albacares*), bigeye (*T. obesus*), and albacore (*T. alalunga*) tunas sampled over a 16-year period from the major ocean bodies, our study aims were to: (i) map spatial patterns at a global scale (i.e. develop carbon isoscapes), (ii) test the influence of biological and environmental factors on the temporal and spatial patterns, (iii) relate environmental, temporal, and spatial variables to tuna bioenergetic condition, and (iv) use modelled baseline and observed tuna δ¹³C values to identify likely transitional zones of tuna movements. Movement was further inferred using calculated isotopic niches (Jackson et al., 2011) based on baseline-corrected δ¹³C values (this study) and nitrogen stable isotope (δ¹⁵N) values (Pethybridge et al., 2018). Based on previous tagging and stable isotope studies that demonstrate limited movements of tunas (Evans et al., 2011; Graham et al., 2010; Houssard et al., 2017; Schaefer et al., 2007; Schaefer and Fuller, 2002), we hypothesized that the δ¹³C values of tropical tunas (i.e. yellowfin and bigeye) would track latitudinal baseline gradients. Given documented north-south migrations of the temperate albacore tuna (Nikolic et al., 2017; Williams et al., 2015), we hypothesized that δ¹³C values for this species would not track baseline values as extensive movements across isoscapes should create a disconnect between tuna δ¹³C values and the baseline values at the site of capture (Graham et al., 2010). Ultimately, identifying the biological and environmental mechanisms responsible for the large-scale patterns or isoscapes will improve our understanding of regional movements of tunas. Such data can inform stock assessment models and better assist fisheries managers and policy makers in managing ecologically dynamic and commercially valuable open-ocean predators.

158

## 2. Methods

160

### 2.1. Isotope data compilation

162

Overall, 4275 δ¹³C values were compiled from the muscle tissue of yellowfin (*Thunnus albacares*), bigeye (*T. obesus*), and albacore (*T. alalunga*) tunas from published and unpublished regional studies and assembled in a global database (Fig. 2; Table 1). Samples spanned a 16-year period from 2000 to 2015. Measurements of tuna size (fork length, cm) were taken, and dates, positions, and fishing gear type were provided for each sampling event. Data were classified by ocean basin (Pacific, Atlantic including Mediterranean Sea, Indian, and Antarctic Oceans) and Longhurst Biogeographical Province (Longhurst, 1998). Longhurst provinces are shown in Figure 2 and defined in Table S1.

171

### 2.2. δ¹³C lipid correction

173
Our dataset consists of both bulk and lipid-extracted samples. Generally, for aquatic samples, C:N ratios < 3.5 correspond to low lipid content and by association minor lipid bias (Hoffman and Sutton, 2010; Logan et al., 2008; Post et al., 2007). Given unlikely lipid bias and potential error created by variable C:N$_{lipid-free}$ among tuna species, no corrections were applied to samples with C:N values ≤ 3.5. For δ¹³C values with C:N values > 3.5, mathematical corrections were applied (Logan et al., 2008) to produce lipid corrected tuna δ¹³C values (δ¹³C$_{Lcorr\_tuna}$; Table 2). The correction equation was based on a dataset of bluefin tuna (T. thynnus) muscle with a maximum C:N of 6.5 so samples with C:N > 6.5 (n=47) were excluded from δ¹³C analyses. The equation used was derived from Fry (2002):

$$ \delta^{13}C_{Lcorr\_tuna} = P - \frac{P * F}{C: N} + \delta^{13}C_{Bulk} $$
where δ¹³C$_{Lcorr\_tuna}$ is the lipid corrected δ¹³C value, C:N is %C/%N, P is an estimate of the isotopic difference between protein and lipid (7.489), and F is an estimate of the C:N value of lipid-free tissue (3.097). Mathematical corrections were applied to 39.5% (yellowfin), 40.5% (bigeye), and 51.7% (albacore) of all samples used in δ¹³C analyses (n=4228). To assess the effect of lipid correction on δ¹³C values and determine if any potential lipid bias remained in the mathematically corrected values, bulk and lipid corrected δ¹³C values were plotted against C:N for the subset of samples that underwent mathematical lipid correction (3.5<C:N≤6.5) (Fig. S1). Equivalent comparisons were made for baseline corrected δ¹³C values (Δ¹³C$_{tuna-phyto}$) (Fig. S1).

## 2.3. Baseline δ¹³C values

Baseline δ¹³C predictions from a global model of ocean biogeochemistry and isotopes, including both δ¹³C and δ¹⁵N values (Schmittner and Somes, 2016), were used to assess and correct for spatial and temporal variability in isotopic baselines (Fig. S2). The model included a three-dimensional (1.8° x 3.6°, 19 vertical levels) ocean circulation model forced with prescribed (fixed) monthly climatological winds (Weaver et al., 2001). The two stable carbon isotopes, ¹²C and ¹³C, were included for dissolved inorganic carbon and organic carbon including phytoplankton, zooplankton, sinking detritus, and dissolved organic carbon. The main relevant processes incorporated in the model included air-sea gas exchange, physical transport, biological uptake, and remineralization of organic carbon. Model results were provided for each respective year of a hindcast from year 2000 of a pre-industrial simulation including increasing atmospheric CO₂ and decreasing δ¹³C-atmospheric CO₂ following Schmittner et al. (2013). All model results were averaged annually over 130 meters of surface ocean. To account for the Suess effect, model predicted phytoplankton δ¹³C values given for year 2000 were subtracted from the phytoplankton δ¹³C values that corresponded to the year of tuna sampling (δ¹³C$_{phyto}$), that is 2000-2015. The difference was used to correct for the Suess effect in tunas (δ¹³C$_{LScorr}$; Table 2) with an average decline of 0.05‰ per year over the 16-year period. Note that the model hindcast simulation also accounted for CO₂-driven warming over this period, which affects δ¹³C$_{phyto}$ values through changes in circulation, air-sea gas exchange, and uptake fractionation that is CO₂ and temperature-dependent. These processes contribute to a larger decrease in δ¹³C$_{phyto}$ than would be expected from the Suess effect alone. As additional biogeochemical and hydrological
corresponding diet data for tuna in the regions sampled in the present study (Duffy et al., 2017), this level of individual trophic variability is not expected and therefore, mixing across isoscapes was considered the most likely source. For these analyses, we assumed that model predictions are equally accurate in terms of both regional median and inter-quartile ranges and that any residual effect of C:N on δ¹³C values not accounted for with mathematical corrections can be ignored.

Movement was further inferred for each species and Longhurst region based on isotopic niches calculated using the Stable Isotope Bayesian Ellipse in R (SIBER; Jackson et al., 2011) and Stable Isotope Analysis in R (SIAR; Parnell et al., 2010) packages in R (R Core Team, 2019). Similar to our calculation table based on Δ¹³C<sub>tuna-phyto</sub> values, we assumed that increased intra-specific and intra-region isotopic variability and associated increased isotopic niche values were driven by movement across isoscapes. Niche calculations were made using δ¹³C<sub>Lcorr_tuna</sub> - δ¹³C<sub>Scorr_phyto</sub> (Δ¹³C<sub>tuna-phyto</sub>) and δ¹⁵N<sub>tuna</sub> - δ¹⁵N<sub>phyto</sub> (See Pethybridge et al., 2018 for nitrogen isotope baseline correction methods). Bayesian standard ellipse areas (SEA<sub>B</sub>), reported as ‰² values, were calculated alongside Bayesian credible intervals to provide a measure of isotopic variability of baseline-corrected δ¹³C and δ¹⁵N values of tuna muscle. The SEA contains approximately 40% of the data and provides a measure of the core niche area as it is effectively the equivalent to standard deviation for bivariate data (Batschelet, 1981; Jackson et al., 2011). In addition to SEA<sub>B</sub>, total area (TA) and SEA corrected for small sample sizes (SEA<sub>c</sub>) were also calculated as measures of isotopic niche width. SEA<sub>c</sub>, reported as ‰², provides a measure of the representative isotopic niche breadth of fish sampled from each region. TA (‰²) is the convex hull area contained within the smallest polygon containing all carbon and nitrogen Δ values and as such reflects total isotopic diversity (Layman et al., 2007). More detailed information for SEA<sub>c</sub> values was derived from calculations of eccentricity (E) and theta (θ°) (based on methods described by Reid et al., 2016; Shipley et al., 2019; Turner et al., 2010). Briefly, eccentricity is a measure of the relative contributions of δ¹³C and δ¹⁵N values to observed variance within the standard ellipse area. Values closer to 0 indicate contributions of both delta values and values approaching 1 reflect variability driven by a single isotope tracer. Theta, reported in degrees, provides a measure of the relative angle of isotopic dispersion relative to δ¹³C with values closer to 0° corresponding to δ¹³C as the primary driver and values closer to 90° instead suggesting δ¹⁵N as the primary influence. Negative theta values reflect opposing trends where δ¹³C values increase with decreasing δ¹⁵N values, which can provide evidence of baseline vs. trophic influences on intra-group isotopic variability as increases in trophic position would be expected to produce a positive relationship.

295

## 2.5. *Tuna bioenergetic condition estimates*

297

Tuna bioenergetic condition was compared across regions and in relation to environmental variables using muscle C:N as a proxy for lipid content and thus bioenergetic condition. Fish muscle C:N values increase linearly with lipid content due to the high percentage of carbon and the absence or minimal nitrogen content in lipids (Logan et al., 2015, 2008). A regression equation developed from Atlantic bluefin tuna muscle samples was used to convert C:N values to percent lipid content, % lipid = 3.9972 * C:N - 11.4663 (Logan et al., 2015). The C:N dataset used for the percent lipid estimates had a maximum C:N value of 7.3, which results in a percent lipid content estimate of 17.7%. Any samples in our C:N dataset with values > 7.3
are reported as >17.7%. Samples with C:N > 7.3 were converted to percent lipids only for mean estimates (Table 1). Chemical lipid extractions were applied to some of our samples (yellowfin: 36.9%, bigeye: 34.4%, albacore: 15.3%) prior to elemental composition analysis and therefore were not included in this condition analysis. Some muscle samples obtained from regions of the tuna other than caudal dorsal muscle (e.g. head, ventral) were also excluded from condition analyses since lipid content can vary across these different body regions in tunas (Balshaw et al., 2008). Following exclusion of lipid extracted and non-dorsal muscle samples, the remaining dataset available for C:N analyses contained 2205 samples.

314

## 2.6. Environmental data

316

A suite of available environmental variables was selected to explore physical forcing on tuna δ¹³C values at the surface and at depth since these tunas are vertical foragers across epipelagic and mesopelagic zones (Olson et al., 2016). Surface variables included sea surface temperature (SST, °C) and chlorophyll-a (Chl-a, mg m⁻³). Monthly averaged NOAA Optimum Interpolation SST data were acquired from the NOAA National Centre for Environmental Information (http://www.esrl.noaa.gov/psd/) using AVHRR version 2 (Reynolds et al., 2007), available at https://ncics.org/data/obs4mips/. The spatial resolution of the SST data was a 0.25° grid, and this product was available for the years 2001–2014. For tuna sampling years outside this range, climatologies for the next closest year were used.

Monthly Chl-a concentrations were acquired from GlobColour data (http://globcolour.info) which has been developed, validated, and distributed by ACRI-ST, France and provides a continuous dataset of merged L3 Ocean Colour product from 1997 to 2016. Monthly means of net primary production (NPP, mg C m⁻²d⁻¹) were derived from a vertically generalized production model (VGPM, http://www.science.oregonstate.edu/ocean productivit y; Behrenfeld and Falkowski, 1997) available at a 9 km resolution which we interpolated onto a 1°x1° grid.

Variables at depth included the mixed layer depth (MLD, m), depth of the 12°C isotherms (Dᵢₐ₁₂, m), and oxygen concentrations (mL L⁻¹) at three depths (375 m O₂₃₇₅m, 100 m below the MLD O₂₁₀₀m, and 50 m below the MLD O₂₅₀m), which were taken from daily CARS climatology (Dunn and Ridgway, 2002; Ridgway et al., 2002) updated using the 2006-2015 Argo-only dataset. The MLD climatology defined the depth of the mixed layer based on an absolute temperature difference of 0.2°C from the surface (de Boyer Montégut et al., 2004). The variables at depth do not vary annually as they are based on a climatology.

For all environmental variables, spatial linear interpolation was used to collocate data to the tuna sample location. For each tuna location, environmental variables were also obtained at the matching tuna sampling date minus 3 and 6 months to account for time delays for prey tissues to be incorporated in tuna muscle (months; Madigan et al., 2012). If environmental data were not available for a particular year, data for the next closest year were used.

345

## 2.7. Data analyses

347

Boxplots were used to compare δ¹³Ccorr values, Δ¹³Ctuna-phyto, C:N ratios and isotopic niche areas between regions and tuna species. These exploratory analyses were performed only for Longhurst province-species datasets with n ≥ 5.

348

349

350
351

General trends in δ¹³Ccorr and Δ¹³Ctuna-phyto values in relation to latitude expressed as distance from the equator and C:N were assessed using linear and non-linear regression, respectively. Normality was assessed visually using QQ plots.

352

353

354

355

## 2.8. Generalized additive mixed models and isoscapes

356

357

Generalized additive mixed models (GAMMs) were used to test potential predictors of Δ¹³Ctuna-phyto, δ¹³Cscorr values, and lipid content. The predictors included spatial (latitude and longitude), morphometric (tuna fork length), and environmental (SST, Chl-a, NPP, MLD, Diso, and O₂) variables. Longhurst province was incorporated as a random effect to account for intra group correlations and overdispersal that frequently arises in spatial analyses (Wood, 2006). Before constructing GAMMs, matrix scatterplots and variance inflation factors (VIF) from the usdm package (Naimi, 2015) in R (R Core Team, 2019) were used to confirm that collinearity was not apparent among the explanatory variables. Individual GAMMs were fitted for each set of environmental variables at the matching sample date (D) and D minus 3, 6, and 9 months prior to tuna capture. For each tuna species, the best time-related environmental variable, determined by the percent deviance explained (DE%, or pseudo R-squared), was retained and used in subsequent models (Tables S3 and S4).

358

359

360

361

362

363

364

365

366

367

368

369

GAMMs were constructed in R using the gamm4 package (Wood, 2006) and the Gaussian family with an identity link function. Smoothing splines were used to examine non-linear fits for all marginal terms considered in the model while a scale-invariant tensor product was used to model the covariates for geographical location. Model selection was performed manually using a backward elimination approach, and the best model fit was selected using the Akaike information criterion for small sample sizes (AICc) and DE%. To determine the amount of variation explained by each predictor variable, we fitted a separate model for individual predictor variables. Subsequent models were then fitted, including each of the remaining smooth terms, one at a time, to assess additional variation explained. To test the performance of each model, standard diagnostic checks were performed on the residuals from the model, and constant variance was observed.

370

371

372

373

374

375

376

377

378

379

380

A GAMM based on two-dimensional thin plate regression splines fitted to the geographical location and interpolated at a ~2° resolution was used to visually explore broad-scale patterns in predicted Δ¹³Ctuna-phyto and lipid content. Due to the potential effect of year, isoscapes were produced for both the earlier (2000–2007) and later (2008–2015) periods of our analysis.

381

382

383

384

385

386

387

# 3. Results

388

389

## 3.1. Lipid corrections and tuna δ¹³C

390

391

Both bulk and δ¹³Ccorr_tuna datasets had negative correlations with C:N while Δ¹³Ctuna-phyto values only declined for the bulk δ¹³C dataset (Fig. S1). Bulk δ¹³C values were lower (~17.9±2.1‰) than lipid corrected (~16.3±1.6‰) values with a modelled maximum offset of 4.3‰ for high lipid samples (Fig. S1). Within the C:N range for which lipid corrections were applied

392

393

394
395 (3.5<C:N≤6.5), the modelled decrease in δ¹³C values with increasing C:N was 3.3 and 6.2‰ for
396 δ¹³C<sub>Lcorr_tuna</sub> and bulk δ¹³C values, respectively. For baseline-corrected δ¹³C values (Δ¹³C<sub>tuna-phyto</sub>),
397 the offset between tuna and baseline values decreased by 2.1‰ for bulk δ¹³C values while
398 modelled isotope spacing increased slightly (<1‰) for δ¹³C<sub>Lcorr_tuna</sub> values (Fig. S1).
399

400

3.2. *Global and ocean region observations*

401

402 Global mean δ¹³C<sub>Lscorr</sub> values were lowest in albacore (-17.2±1.0‰) and similar in
403 bigeye (-15.9±1.0‰) and yellowfin (-15.7±0.9‰) tunas (Table 1). The range of observed
404 δ¹³C<sub>Lscorr</sub> values was greatest in albacore (7.0‰) and smallest in bigeye (5.2‰). Across ocean
405 regions, highest mean δ¹³C<sub>Lscorr</sub> values were observed in Longhurst provinces within the eastern
406 Pacific Ocean and central tropical Pacific Ocean for yellowfin and bigeye and the western
407 Pacific Ocean for albacore (Figs 3-5; S3). The lowest mean δ¹³C<sub>Lscorr</sub> values were observed in
408 the AUSE Longhurst province of the western Pacific and NASW and GFST provinces in the
409 western North Atlantic Ocean for yellowfin and bigeye (Figs 3-4; S3). For albacore, the lowest
410 δ¹³C<sub>Lscorr</sub> values were in the eastern Atlantic in the NADR and CNRY provinces (Fig 5; S3).
411

411 Differences between δ¹³C values of tuna and phytoplankton (Δ¹³C<sub>tuna-phyto</sub>; Table 2), used
to report degrees of spatial movement and population dispersal, were globally similar among
412 species with means varying between -5.0 to -5.3‰ (Table 1). Median Δ¹³C<sub>tuna-phyto</sub> across
413 Longhurst provinces ranged from -4.6 to -8.3‰, -3.9 to -7.9‰, and -4.0 to -7.0‰ for yellowfin,
414 bigeye, and albacore, respectively (Figs. 3-5). Minimum differences between Δ¹³C<sub>tuna-phyto</sub> values
415 were within the anticipated -2 to -8‰ range, but maximum differences exceeded this range for
416 yellowfin and bigeye in several regions (Figs. 3,4). Elevated values (> -8‰) were observed in
417 the GUIN (Atlantic) and PEQD (Pacific) Longhurst provinces for both yellowfin and bigeye
418 translating to a “high degree of movement by some individuals” classification (Table 3). Similar
419 maximum values were observed for yellowfin in the PNEC province (Fig. 3) which were also
420 classified as “high degree of movement” (Table 3). For bigeye, the GFST (Atlantic) province
421 included high Δ¹³C<sub>tuna-phyto</sub> differences (maximum of -7.6‰) (Fig. 4). While none of the
422 maximum albacore values exceeded -8‰, values for bordering provinces in the Northeast
423 Atlantic (NADR and NECS) were -7.9 and -7.7‰, respectively, and therefore were classified as
424 “high degree of movement by some individuals” and “moderate movement”, respectively (Table
425 3). ARCH, TASM, and SANT, provinces in the western Pacific and Antarctic, had similarly high
426 (-7.4 to -7.7‰) maximum values (Fig. 5) and “moderate movement” classifications for albacore.
427

428 Regional variability in Δ¹³C<sub>tuna-phyto</sub> values, used to report levels of population dispersal,
differed between and within the Longhurst provinces with ranges among individual tuna of 0.0 to
429 -4.9‰ (yellowfin), 0.0 to -5.2 (bigeye), and -0.5 to -5.2‰ (albacore) (Figs. 3-5). The largest
430 ranges (> -3.6‰) among non-outlier values were found in the bordering WARM and PEQD
431 regions of the equatorial Pacific for both yellowfin and bigeye (Figs. 3,4), which were assigned a
432 classification of “very high dispersal” (Table 3). Yellowfin in the GUIN region also had “very
433 high dispersal” with a range of 4.7‰. AUSE (yellowfin), SPSG (both yellowfin and bigeye), and
434 ARCH (both yellowfin and bigeye) regions in the Pacific had similarly high ranges (Fig. 3) and
435 were assigned the “high dispersal” classification (Table 3). For albacore, the greatest ranges (> -
436 3.6‰) occurred in the bordering ARCH and SPSG Longhurst provinces (Fig. 4) and had a
437 dispersal classification of “very high” (Table 3). Relatively low regional variability was observed
438 in the Indian Ocean for all species (Figs. 3-5), and in albacore from the MEDI province (the
439
440 tropical yellowfin and bigeye tunas were not sampled from this region) (Fig. 5). These regions
441 all had population dispersal scores of “no”, “low”, or “moderate dispersal” (Table 3).

442 SEA<sub>B</sub>, SEA<sub>c</sub>, and TA values were highest for Longhurst provinces in the Pacific Ocean
443 for all three tuna species (Table S2). For albacore, SEA<sub>B</sub>, SEA<sub>c</sub>, and TA values from the SPSG
444 were elevated (e.g. SEA<sub>c</sub>>10) relative to all other sampled Longhurst regions (Fig. 5; Table S2).
445 For bigeye and yellowfin, highest SEA<sub>B</sub>, SEA<sub>c</sub>, and TA values were observed in the SPSG,
446 WARM, and PEQD regions (Table S2).

447 Eccentricity values fell close to 1 for most provinces (mean across Longhurst provinces ≥
448 0.8 for all three tuna species) while theta was closer to 90° (absolute mean across Longhurst
449 provinces > 70° for all three tuna species) indicating that most of the intra-specific variability
450 was driven by variation in δ<sup>15</sup>N values (Table S2). Theta values were negative for many
451 provinces (37.5%, 31.3%, and 44.4% of provinces for yellowfin, bigeye, and albacore,
452 respectively), which reflects increasing δ<sup>13</sup>C values with decreasing δ<sup>15</sup>N values in these regions
453 (Table S2).

454 C:N values were generally low for yellowfin and bigeye tunas but higher and more
455 variable for albacore tuna (Figs. 3-5). For yellowfin and bigeye, median values for all provinces
456 were < 4 (4.5% lipid content), but several outliers (with median C:N values ~ 5, 8.5% lipid
457 content) were observed in the western and central Pacific Ocean (WCPO) in the WARM
458 province (Figs. 3,4). In the WCPO, similar outliers with elevated C:N values were present in the
459 WARM (yellowfin) and NPTG (bigeye) provinces. For the Indian Ocean, yellowfin with
460 elevated C:N values (C:N~4-5) were present in the MONS province (Fig. 3). Based on
461 regressions of muscle C:N and lipid content, median values for yellowfin and bigeye correspond
462 to < 3% lipid content. For albacore, the north Atlantic Ocean had the highest C:N values with a
463 median value > 4.5 (6.5% lipid content) in the NADR and NECS and individual values as high
464 as 8 (>17.7% lipid content). While the regional median C:N was < 3.9 (4.1% lipid content) for
465 albacore in the ARCH and AUSE provinces of the Southwest Pacific Ocean, many outliers had
466 elevated C:N between 4 and 6 (4.5% to 12.5% lipid content) (Fig. 5).

467

3.3. *Model results*

469

470 The best model (i.e. lowest AIC<sub>c</sub> value) for δ<sup>13</sup>C<sub>LScorr</sub> values of all species combined
471 spatial, environmental, and year variables (DE% = 65.3% in yellowfin, 70.0% in bigeye, and
472 80.7% in albacore) (Table S3). According to AICc, year was one of the best standalone
473 predictive variables in all tuna species, with more negative δ<sup>13</sup>C<sub>LScorr</sub> values predicted in later
474 sampling years with intercepts reached at 2006 or 2007 (Fig. 6). Spatial variables were also
475 ranked highly for all species with an inverse relationship occurring between latitude and δ<sup>13</sup>C
476 <sub>LScorr</sub> values with higher values around the equator and lower values of δ<sup>13</sup>C <sub>LScorr</sub> at high latitudes
477 (Fig. 6; Fig. S3). There were no clear relationships with longitude.

478 Similarly, linear regressions of δ<sup>13</sup>C<sub>LScorr</sub> values with latitude expressed as distance from
479 the equator all had significant negative correlations (Fig. S4; All tuna: R<sup>2</sup>=0.43; F<sub>1,4226</sub>=3152;
480 p<0.001, yellowfin: R<sup>2</sup>=0.06; F<sub>1,1617</sub>=106.4; p<0.001, bigeye: R<sup>2</sup>=0.04; F<sub>1,932</sub>=36.5; p<0.001,
481 albacore: R<sup>2</sup>=0.52; F<sub>1,1781</sub>=1781; p<0.001). The intercept of the best fit lines for all tuna
482 combined was 5‰ above the intercept of modelled phytoplankton δ<sup>13</sup>C values while separation
483 for individual tuna species were 5.5‰ (yellowfin), 5.2‰ (bigeye), and 3.8‰ (albacore) (Fig.
484 S4). The slope of the best fit line for all tuna combined was -0.05 (‰ change in δ<sup>13</sup>C per degree

489

p<0.001, yellowfin: R²=0.06; F₁,₁₆₁₇=106.4; p<0.001, bigeye: R²=0.04; F₁,₉₃₂=36.5; p<0.001,

481 albacore: R²=0.52; F₁,₁₇₈₁=1781; p<0.001). The intercept of the best fit lines for all tuna

482 combined was 5‰ above the intercept of modelled phytoplankton δ¹³C values while separation

483 for individual tuna species were 5.5‰ (yellowfin), 5.2‰ (bigeye), and 3.8‰ (albacore) (Fig.

484 S4). The slope of the best fit line for all tuna combined was -0.05 (‰ change in δ¹³C per degree
485 of latitude) relative to -0.08 for associated modelled phytoplankton. Yellowfin and bigeye best fit
486 line slopes (-0.02) were similar to associated phytoplankton (-0.03 for both datasets) while the
487 slope for the albacore best fit line (-0.05) was less steep than the associated phytoplankton
488 baseline (-0.11). While all regressions showed a significant negative correlation between
489 modelled phytoplankton δ¹³C values and latitude (P<0.001 for all comparisons), R² values varied
490 across datasets. Lower values were observed associated with the sampling locations for
491 yellowfin (R²=0.10) and bigeye (R²=0.13) while modelled phytoplankton δ¹³C values associated
492 with the albacore dataset had an R² of 0.82.

493 For Δ¹³C_tuna-phyto, the full tuna dataset had a significant negative relationship with latitude
494 (R²=0.09; F₁,₄₂₃₂=443.5; p<0.001), but this relationship was mainly driven by albacore (R²=0.51;
495 F₁,₁₆₇₄=1778; p<0.001) (Fig. S4). For the tropical tunas, bigeye had a weak negative relationship
496 (R²=0.01; F₁,₉₃₄=9.5; p<0.01) while yellowfin Δ¹³C_tuna-phyto was not correlated with latitude
497 (R²=0.00; F₁,₁₆₂₀<0.01; p=0.98) (Fig. S4).

498 Smooth terms for most environmental variables were significant in all tuna species (Table
499 S3). Oxygen concentration at depth was consistently the best ranked environmental model
500 (according to AICc). For all species, predicted δ¹³C<sub>LScorr</sub> values declined after oxygen
501 concentrations reached 4 mL L⁻¹ (Fig. 6). The next best ranked environmental variable was SST
502 with δ¹³C<sub>corr</sub> values predicted to increase at higher temperatures with the predicted intercepts
503 reached at 25.8°C, 25.9°C, and 25.0°C for yellowfin, bigeye, and albacore, respectively (Fig. 6).
504 D<sub>Iso12</sub> consistently showed a parabolic curve in which δ¹³C<sub>LScorr</sub> values were highest at about 250,
505 300, and 550 m depth for yellowfin, bigeye, and albacore, respectively. As Chl-α and NPP
506 increased, δ¹³C<sub>LScorr</sub> values were predicted to increase in yellowfin and decrease in bigeye and
507 albacore. Trends with MLD were less clear.

508 Length was a less efficient predictor of δ¹³C<sub>LScorr</sub> values in all tuna species, globally.
509 Global patterns varied among species; decreasing in bigeye, increasing in albacore, and the trend
510 unclear in yellowfin (data not shown).

511 Model results of Δ¹³C<sub>tuna-phyto</sub> were largely similar to those for δ¹³C<sub>LScorr</sub> values with the
512 best models (i.e. those with the lowest AIC<sub>c</sub> values) including oxygen at depth and the spatial
513 and temporal variables (Fig. 7; Table S4). Predictive ability was slightly higher (DE% = 71.0%
514 in yellowfin, 71.5% in bigeye, and 84.0% in albacore). For all tuna species, the most significant
515 standalone variable was year followed by oxygen at depth. Annual changes in Δ¹³C<sub>tuna-phyto</sub> closely
516 matched corresponding δ¹³C<sub>LScorr</sub> patterns. Similar to δ¹³C<sub>LScorr</sub>, oxygen at depth was the best
517 environmental variable; values of Δ¹³C<sub>tuna-phyto</sub> declined with increasing oxygen at depth for all
518 three species with the caveat that values then increased at highest oxygen concentrations (> 4 mL
519 L⁻¹) for albacore. Relationships between Δ¹³C<sub>tuna-phyto</sub> and SST were the inverse of corresponding
520 δ¹³C<sub>LScorr</sub>; Δ¹³C<sub>tuna-phyto</sub> instead declined as SST increased. For latitude, Δ¹³C<sub>tuna-phyto</sub> increased in
521 northern latitudes relative to the equator for all three tuna species. In southern latitudes, albacore
522 Δ¹³C<sub>tuna-phyto</sub> showed a similar pattern while yellowfin and bigeye instead showed slight declining
523 trends with distance from the equator (Fig. 7).

524 The model predictive ability of C:N ratios was reasonably good for all tuna species (59.0
525 to 65.0 DE%) (Table S5). The best model (i.e. lowest AICc value) for yellowfin and bigeye was a
526 variable model that included year and the spatial variables while oxygen at depth was most
527 influential for albacore (Table S5).

528
## 3.4. Spatio-temporal patterns

530

531 All tuna species showed strong spatio-temporal patterns of δ¹³C<sub>LScorr</sub> values with distinct geographical regions identified in the GAMM predicted and interpolated isoscapes (Fig. S3).

532 Clear differences between the two selected time periods were observed with higher values predicted in the earlier (2000-2007) sampling years for each of the ocean basins. Largest temporal differences were observed in the Pacific Ocean and in the southwestern Indian Ocean.
533
534 In both sampling periods, predicted δ¹³C<sub>LScorr</sub> values were highest in the western central Pacific
535
536 Ocean and lowest in the North Atlantic Ocean and southwestern Pacific and Indian Oceans for
537
538 all tuna species.

539 Predicted Δ¹³C<sub>tuna-phyto</sub> varied spatially during both periods for the three tuna species (Fig. 8). During 2000-2007, the most negative Δ¹³C<sub>tuna-phyto</sub> values, suggesting higher degrees of spatial movement, for bigeye and yellowfin were in the eastern tropical Pacific while albacore values were most negative in the southern sampling extent of our western and central Pacific dataset.
540
541
542
543 During 2008-2015, bigeye and yellowfin data were only available for the Pacific and Indian Oceans, where most negative values occurred in the central Pacific. For albacore, most negative values were in the Northeast Atlantic region (Fig. 8).
544
545

546 Maps of model predicted C:N values for each tuna species sampled over the study period
547 showed few spatial gradients (Fig. 9). Yellowfin values were consistently low across ocean
548 regions and sampling periods. Bigeye C:N values were low and uniform throughout much of the
549 sampled area except for some elevated values in the central Pacific Ocean. For albacore, C:N
550 values were highest in the Atlantic Ocean (Fig. 9).
551

552

553

# 4. Discussion

554

555 An understanding of top predator movements and trophodynamics in the open ocean is
crucial to informing stock assessments of commercially important species, and in predicting and
556
determining shifts in a changing climate. Here we provide evidence of potential tuna movement
557
and residency patterns, illustrating how spatial isotope analyses can provide insight into the
558
habitat use of large migratory predators to directly inform conservation and management efforts.
559
Given stable isotope theory, we assumed greater regional residency where the difference
560
between corrected-tuna and model-derived phytoplankton δ¹³C values (Δ¹³C<sub>tuna-phyto</sub>, ‰) was
561
between -2 and -8‰. By correcting for known lipid and Suess effects and adjusting for baseline
562
variability our study demonstrates a novel and rigorous approach employing tuna δ¹³C values to
563
reflect spatiotemporal changes in movement patterns at a global scale. Predictive models of
564
Δ¹³C<sub>tuna-phyto</sub> values allowed us to determine which environmental variables were explaining
565
observed temporal, spatial, and ontogenetic patterns in tuna foraging habitats. Combining
566
baseline-corrected δ¹³C and δ¹⁵N data to calculate isotopic niche width allowed for further
567
assessment of individual isotopic variability as an indicator of population dispersal across
568
regional isoscapes. Assessing additional variability in percent lipid as derived from C:N values
569
further assisted in understanding tuna movements and trophodynamics.
570
571

571
## 4.1. Temporal, spatial, and environmental gradients

Tuna δ¹³C<sub>LScorr</sub> values declined across sampling years for all three tuna species, despite correcting for Suess and additional climate effects on phytoplankton. The decline was continuous across sampling years for albacore while yellowfin and bigeye declines were slightly delayed until 2003 (Fig. 6). This observed decline was thoroughly explored using time-series analysis and theoretical modelling (Lorrain et al., 2020). Analyses by Lorrain et al., (2020) showed that our observed temporal declines are most likely due to temporal shifts in phytoplankton community structure and/or physiology that differ across ocean basins (greater temporal decline in tuna muscle δ¹³C in the Pacific than Atlantic and Indian Oceans).

In general, phytoplankton and zooplankton baseline δ¹³C values increase with increasing SST toward the equator (Magozzi et al., 2017; McMahon et al., 2013). Here we demonstrate similar patterns for three tuna species indicating, in agreement with previous literature (e.g. Bird et al., 2018a), that marine carbon isoscapes can be directly linked with top predators. Overall adherence to these latitudinal trends also suggests relatively constrained movements within the timescale of muscle isotopic turnover (months; Madigan et al., 2012) as foraging across a broad latitudinal gradient was expected to blur baseline differences in phytoplankton and zooplankton δ¹³C values (exceptions discussed below). These baseline gradients would not be apparent in higher trophic level predators if e.g. the tuna sampled in the present study frequently traversed these latitudinal isoscapes within months prior to capture. A global analysis of shark muscle δ¹³C values found a similar trend for sharks foraging in continental shelf, but not pelagic, waters (Bird et al., 2018a). Similar patterns have been documented for yellowfin and bigeye in the Pacific based on latitudinal gradients in bulk δ¹⁵N values (Graham et al., 2010; Houssard et al., 2017; Olson et al., 2010). Popp et al. (2007) demonstrated that latitudinal patterns in yellowfin muscle bulk δ¹⁵N values were indicative of baseline effects based on amino acid stable isotope analysis (AA-SIA). Sharks inhabiting the continental shelf and the three tuna species included in our analysis—based solely on similarity in latitudinal trends in muscle δ¹³C and underlying baseline δ¹³C values—appeared to have constrained movements. Mean δ¹³C differences between phytoplankton and all three tuna species (~5‰) were similar to the 4.6‰ global average difference reported for sharks (Bird et al., 2018a) and consistent with a presumed average tuna TP of 4.5 to 5.0 based on collective expected TDFs in the trophic pathway to tunas (Caut et al., 2009; France and Peters, 1997; Hobson and Cherel, 2006; Madigan et al., 2012; Sweeting et al., 2007) (Fig. 1). Preservation of latitudinal baseline δ¹³C gradients at the global scale is not surprising for yellowfin or bigeye as these tropical tunas tend to have mainly east-west movements and more constrained north-south movements in most of their global distribution (Evans et al., 2011; Fonteneau and Hallier, 2015; Schaefer et al., 2007; Schaefer and Fuller, 2002, 2009).

Species undertaking extensive north-south migrations would be expected to have only a weak or no correlation with latitudinal baseline values as foraging across different isoscapes should cause tissue δ¹³C values to reflect a mix of baseline values rather than solely the baseline of their sampling location. Open ocean shark δ¹³C values, for example, did not closely match baseline latitudinal trends, which could be due to large scale movements and foraging across isoscapes and a consequent homogenization of latitudinal baseline values (Bird et al., 2018a). Albacore tuna are known or presumed to undertake long distance north-south migrations between foraging and spawning grounds (Nikolic et al., 2017; Williams et al., 2015) and therefore a similar homogenization and disassociation from baseline gradients was expected for
this species. While muscle tissue may not reach isotopic steady state at any given foraging ground for highly migratory species, some degree of resident feeding will reflect at least a fraction of the local δ¹³C baseline. Our assumption that muscle δ¹³C values of highly migratory species would blur isotopic baselines due to feeding and movement across isoscapes may have been overly simplistic with regards to species with directed movements between spawning and foraging habitats. Tunas have been described as "energy speculators" that undertake long migrations in search of productive foraging grounds (Brill, 1996). Upon reaching a productive foraging ground, tunas would begin to incorporate local δ¹³C baseline values in their tissues. For example, Atlantic bluefin tuna are known to undergo extensive north-south migrations throughout the North Atlantic (Mather et al., 1995), but their mean muscle stable isotope values still reflect local isoscapes (Logan, 2009; Logan et al., 2015). However, individual isotopic variability is high, likely due to past migrations from different isoscapes and different arrival times to the areas of capture (Estrada et al., 2005; Logan et al., 2015). While albacore δ¹³C values tracked baseline values across latitudes, further exploration of this relationship revealed that the rate of change for tuna δ¹³C values was far less than that observed across the same locations for modelled phytoplankton and so the isotopic spacing between tuna and local baselines varied negatively with latitude (Fig. S4). This latter pattern is consistent with albacore retaining some baseline isotope values from prior north-south migrations; fish foraging in higher latitudes with lower baseline δ¹³C values would have positively biased separation from local baselines while fish foraging near the equator with elevated phytoplankton baseline values would be negatively biased from previous foraging in high latitudes. Yellowfin and bigeye muscle δ¹³C relationships with phytoplankton were instead closely linked, which is consistent with spatially constrained movements.

The extent to which highly migratory predator muscle δ¹³C values will reflect isotope gradients is related to the duration of residency and feeding intensity within a given isoscape as well as tissue turnover rates. A fish would be expected to reflect local baselines if it remains resident in that region for a time sufficient for the sampled tissue to reach isotopic steady state. While migratory tunas may not reside and feed in any given area for a sufficient duration to reach isotopic steady state, muscle tissue from the relatively small (fork length ~ 80-90 cm) tunas sampled in our study should at least partly reflect local baselines given the elevated metabolic rates and regional endothermy attributes of these species (Dickson and Graham, 2004; Olson et al., 2016) and relationships between size, metabolic rate, and isotopic turnover rates (Thomas and Crowther, 2015; Vander Zanden et al., 2015). By contrast, pelagic sharks from the Bird et al. (2018a) analysis would presumably have slower isotopic turnover rates based on estimates from a representative shark species that exhibits ram ventilation (Logan and Lutcavage, 2010) and larger average size (178±61.0 cm TL) (Bird et al., 2018a, 2018b). These larger pelagic sharks with slower isotopic turnover rates would likely not reflect local baselines due to the mismatch between regional residency and isotopic turnover. For example, muscle δ¹⁵N values of shortfin mako sharks (*Isurus oxyrinchus*) of a similar size (178-232 cm FL) sampled in southern New England waters did not reflect known seasonal diet shifts and associated migrations between offshore and shelf waters (MacNeil et al., 2005). Shortfin makos are regional endotherms (Carey and Teal, 1969) with among the highest measured routine metabolic rates of any shark species (Sepulveda et al., 2007), so muscle tissue of other large pelagic sharks would also be unlikely to reflect local isoscapes from seasonal foraging grounds. These differences in size, metabolic rate, and associated isotopic incorporation rates may explain the disconnect from baseline trends for
663
664

pelagic shark δ¹³C values (Bird et al., 2018a) as opposed to the significant correlations observed in our tuna dataset.

665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698

Our models showed that spatiotemporal patterns in tuna δ¹³C<sub>LScorr</sub> and Δ¹³C<sub>tuna-phyto</sub> values were largely related to changes in oxygen concentrations at depth. The oxygen minimum layer was the most powerful environmental predictor of baseline-corrected tuna δ¹⁵N values for the same tuna sampled in this study (Pethybridge et al., 2018) of which tuna trophic position was shown to decline under reduced oxygen concentrations. Tuna δ¹³C<sub>LScorr</sub> and Δ¹³C<sub>tuna-phyto</sub> instead generally increased under reduced oxygen conditions (Figs. 6 and 8). Therefore, relationships between tuna isotope values and oxygen at depth could be influenced by vertical δ¹³C and δ¹⁵N isoscapes since baseline δ¹³C and δ¹⁵N values tend to decrease and increase with depth, respectively (Altabet et al., 1991; Miller et al., 2008; Richards et al., 2019). As the model-derived phytoplankton δ¹³C values were based on surface waters, model corrections would not account for vertical baseline shifts. If declining trends in tuna δ¹³C<sub>LScorr</sub> and Δ¹³C<sub>tuna-phyto</sub> values with increasing oxygen at depth are a result of vertical baseline differences, this relationship would suggest that tuna vertical feeding is influenced by oxygen concentration, which in turn would provide support for foraging habitat alteration with ocean warming. These relationships would be expected for yellowfin as this species is more vertically constrained and limited by oxygen concentrations (>3.5 ml O₂ l⁻¹ required; Brill et al., 2005). Bigeye instead routinely feed in the mesopelagic zone (Lam et al., 2014; Schaefer and Fuller, 2002), can tolerate low oxygen concentrations (< 1.5 ml O₂ l⁻¹; Brill et al., 2005), and thus would not be expected to vary foraging habitat across a higher oxygen concentration gradient. The decline in Δ¹³C<sub>tuna-phyto</sub> with increasing dissolved oxygen concentrations at depth follows a steeper gradient for bigeye than yellowfin, which is consistent with bigeye only being impacted by the lowest oxygen concentrations (i.e. < 1.5 ml O₂ l⁻¹) and yellowfin vertical distribution and foraging being influenced across a broader range of low dissolved oxygen values (i.e. < 3.5 ml O₂ l⁻¹). Albacore tolerance of low dissolved oxygen is less than bigeye and similar to yellowfin (Arrizabalaga et al., 2015; Sund et al., 1981). While albacore δ¹³C<sub>LScorr</sub> values decline with increasing oxygen at depth, the more parabolic Δ¹³C<sub>tuna-phyto</sub> relationship with oxygen concentration may reflect combined influences of vertical δ¹³C baselines, trophic position, and body condition across different vertical foraging habitats (Table S5; Pethybridge et al., 2018). The apparent tracking of vertical δ¹³C baseline gradients in tuna muscle with oxygen concentration at depth suggests feeding at depth follows the relative physiological tolerances of these tuna species. Consequently, foraging habitat could be impacted by habitat compression resulting from the predicted expansion of the oxygen minimum layer with ocean warming (Oschlies et al., 2008; Stramma et al., 2012).

699
700
701
702
703
704
705
706
707
708

## 4.2. *Spatial movement and population dispersal patterns*

While variability and uncertainty in trophic discrimination factors (TDFs) for tunas and underlying lower trophic level taxa limit our ability to interpret modelled differences between tuna and phytoplankton δ¹³C values, some general patterns emerged. High Δ¹³C<sub>tuna-phyto</sub> values near or beyond the extremes of the range presumed based on trophic enrichment (-2 to -8‰) and high individual variability (≥3‰) provided evidence of potential movements and dispersal across isoscapes (Table 3). While trophic diversity cannot be completely ruled out as an alternative explanation for these patterns, Δ¹³C<sub>tuna-phyto</sub> values near or outside of our assumed -2 to -8‰ would require TDFs across trophic levels to be consistently at the higher or lower range of TDFs
709 expected for these species. The higher tuna-phytoplankton δ¹³C offsets observed in several

710 regions could be the result of higher TDFs across trophic linkages in these food webs, but

711 presuming tuna TPs of ~4.5, these greater offsets would require all trophic linkages to be

712 enriched in ¹³C relative to diet by almost 2‰. Several controlled feeding studies have shown a

713 positive relationship between fish muscle TDFs and temperature (Bloomfield et al., 2011;

714 Martino et al., 2019), which could bias spatial trends in Δ¹³C<sub>tuna-phyto</sub> values. However, Δ¹³C<sub>tuna-</sub>

715 phyto values in our dataset were lower in the warmer SST regions near the equator relative to

716 higher latitudes, suggesting factors other than TDF variability are driving these patterns.

717 Similarly, variability in Δ¹³C<sub>tuna-phyto</sub> values consistent with our “high dispersal” and “very high

718 dispersal” classifications (Table 3) would require ≥ two trophic levels of dietary segregation

719 among individuals in a given foraging location over the timescale of muscle isotopic turnover

720 (months), a degree of individual variability not consistent with past diet studies (Duffy et al.,

721 2017). Tuna size was not included in our best GAMMs (i.e. those with lowest AIC<sub>c</sub> values),

722 which, given that tuna diet and TP often vary ontogenetically (Olson et al., 2016), suggests

723 observed δ¹³C variability was explained more by spatial baseline than trophic differences.

724 Moderate to high degrees of movement were observed for all three tuna species in several

725 regions of the Atlantic and Pacific oceans (Table 3). The elevated values relative to planktonic

726 baselines are not all due simply to higher tuna TP in these regions, as corresponding TP

727 estimates based on model-corrected δ¹⁵N were actually below global mean estimates for most of

728 the oceanic regions containing these Longhurst provinces (Pethybridge et al., 2018). Even with

729 some degree of individual variability in TP and inherent variability in TDFs, the wide ranges of

730 Δ¹³C<sub>tuna-phyto</sub> values observed in our dataset for certain regions exceeded expected ranges (3–4‰)

731 and likely relate to movements across carbon isoscapes. Large isotopic ellipse areas and

732 variability in baseline-corrected δ¹³C values further reflect movements across isoscapes as

733 observed ranges would require trophic separation of greater than two trophic levels across

734 individuals based on a presumed tuna δ¹³C TDF of 1.8‰ (Madigan et al., 2012). This shows that

735 our classification of both movement and dispersal based on baseline-corrected δ¹³C<sub>LScorr</sub> values is

736 a useful tool, particularly when combined with results from complementary global analyses of

737 tuna stomach contents and trophic position and verified with existing tagging data.

738 Moderate to high movement and dispersal were consistently reported among the tuna

739 species within the WCPO. The WCPO is a large and dynamic area that includes both

740 oligotrophic and nutrient-rich waters (Brown et al., 2014; Le Borgne et al., 2011) therefore both

741 variable baseline-corrected δ¹³C and bioenergetic condition are consistent with movements and

742 feeding across these differing water masses. High variability in the diet diversity and TP of these

743 tuna species has been reported in this area (Duffy et al., 2017; Pethybridge et al., 2018),

744 consistent with mixing across isoscapes. In addition, bigeye in the WCPO also showed high

745 individual C:N variability (Fig. 4). Based on previous tagging studies within this region, bigeye

746 undertake long distance movements (Fonteneau and Hallier, 2015; Schaefer et al., 2015), while

747 yellowfin display higher residency, particularly in the WARM province (Rooker et al., 2016), a

748 result similar to our study where yellowfin were assigned a classification of “low movement” in

749 the same area (Table 3). Following model baseline corrections, δ¹³C values remained most

750 variable for tunas in the WCPO, particularly the SPSG province for which albacore were

751 classified as having “very high dispersal.” Phytoplankton carbon isoscapes are highly variable

752 across this region of the Pacific (Magozzi et al., 2017) while particulate organic matter (POM)

753 δ¹⁵N isoscapes vary by as much as 10‰ (Houssard et al., 2017). Based on similarities in baseline

754 POM and tuna δ¹⁵N isoscapes from the WCPO, Houssard et al. (2017) concluded that yellowfin
and bigeye had constrained movements within this ocean region. While latitudinal baseline isoscapes generally matched tuna $\delta^{15}$N in Houssard et al. (2017) and $\delta^{13}$C values in our study, individual variability within Longhurst provinces is sufficiently high for both $\delta^{13}$C and $\delta^{15}$N values to suggest a degree of mixing and foraging within the WCPO.

Another pattern consistent among tuna species was moderate degrees of movement reported in the GFST province in the Northwest Atlantic Ocean based on tuna $\Delta^{13}C_{tuna-phyto}$ median values of -4 to -6‰. The GFST province borders southerly waters with zooplankton-based isoscapes with higher $\delta^{13}$C and lower $\delta^{15}$N values (McMahon et al., 2013). Electronically tagged bigeye have been shown to travel between the northern waters represented by our isotope dataset and more southerly waters extending as far south as the Brazilian shelf (Lam et al., 2014). Recent movements from these southerly isoscapes would produce the elevated $\Delta^{13}C_{tuna-phyto}$ and negative theta values observed in our dataset as well as lower TP estimates reported in Pethybridge et al. (2018).

Only albacore were sampled from the NADR province in the Northeast Atlantic Ocean, of which the reported “*high degrees of movement*” agree with existing tagging data showing that albacore travelled to southerly and offshore waters mainly southwest off the Iberian Peninsula (Cosgrove et al., 2014). Zooplankton carbon isoscapes for this region are elevated relative to NECS and NADR waters while corresponding $\delta^{15}$N isoscapes are fairly homogenous (McMahon et al., 2013). Similar to the elevated $\Delta^{13}C_{tuna-phyto}$ observed in GFST tunas, the elevated median $\Delta^{13}C_{tuna-phyto}$ observed for NADR albacore could be due to movements from these offshore waters with higher baseline $\delta^{13}$C values. Despite the elevated $\Delta^{13}C_{tuna-phyto}$ values, ranges of these values in addition to SEA values were more constrained suggesting fairly homogenous movements for these NADR albacore. Albacore diet diversity generally decreases with increasing latitude and so more constrained isotopic niches for this species in the Northeast Atlantic could be partly due to a corresponding lack of individual variability in diet in this region (Duffy et al., 2017).

Yellowfin and bigeye sampled from the GUIN province in the eastern Atlantic had the highest estimated degrees of movement observed in this study, based on $\Delta^{13}C_{tuna-phyto}$ values. The global model used for $\delta^{13}$C baseline corrections could potentially underestimate phytoplankton $\delta^{13}$C near this coastal upwelling zone, which would result in artificially high estimates of $\Delta^{13}C_{tuna-phyto}$ values. Bigeye in the GUIN showed little $\Delta^{13}C_{tuna-phyto}$ variability and tagging data have shown fairly constrained movements for this region (Carruthers et al., 2010), which further suggests the observed difference may not be indicative of movements across isoscapes. In contrast, yellowfin in the GUIN province had the greatest variability in $\Delta^{13}C_{tuna-phyto}$ values (4.7‰) which corresponds to “*very high dispersal*” and high SEA$_c$. While absolute baseline-corrected $\delta^{13}$C values may be inaccurate due to modelling deficiencies, the high variability in estimates across individuals does suggest some influence of movement across isoscapes for yellowfin. Diet diversity is low for yellowfin in this and other coastal upwelling regions (Duffy et al., 2017), therefore wide isotopic niche space is not likely due to trophic diversity. Many conventionally-tagged yellowfin in the Atlantic travelled to the eastern central Atlantic Ocean (ECAO) region from distant waters in the western Atlantic and Gulf of Mexico while others displayed more localized movements within the ECAO (Carruthers et al., 2010). Yellowfin from the ECAO also showed a high degree of variability in estimated TP using baseline-corrected $\delta^{15}$N (Pethybridge et al., 2018), which is consistent with influences of movements from different isoscapes combined with more constrained movements for other individuals within this region.

Another consistent pattern among species was the classification of “*limited movement*” in the EAFR province in the coastal waters in the southwestern Indian Ocean. In fact, limited to low
degrees of movement were found for all tunas sampled in southwestern Indian Ocean provinces (ARAB, MONS, EAFR, ISSG). Limited tagging data exist for these provinces, but diet diversity and TP have been shown to be low for the same tuna species sampled in the present study (Duffy et al., 2017; Pethybridge et al., 2018). Levels of population dispersal and SEA values were also relatively low within the southwestern Indian Ocean, suggesting residency.

806

### 4.3. Condition

808

Total lipid content is used as a proxy for bioenergetic condition and can reflect a requirement of energy for life history events such as reproduction or migration (Marshall et al., 1999). For albacore, the higher and more variable lipid content likely reflects its migration range over a larger spread of latitudes than yellowfin and bigeye (Childers et al., 2011; Williams et al., 2015). Albacore encounter a more diverse mix of both productive temperate waters with abundant, energy-rich prey as well as more oligotrophic waters in the tropics where prey is less abundant and lower energetic value (Duffy et al., 2017). As yellowfin and bigeye movements are more constrained latitudinally, they also forage in less productive waters. Albacore in the eastern North Atlantic undertake seasonal migrations to foraging grounds in the Bay of Biscay with muscle lipid stores increasing across the summer to fall foraging period (Goñi and Arrizabalaga, 2010). Muscle lipid stores also increase with albacore size in the Bay of Biscay (Goñi et al., 2011), with a peak at 80 cm FL followed by a declining trend (Goñi and Arrizabalaga, 2010). High but variable lipid stores in albacore sampled from the eastern North Atlantic Ocean (NAO) likely reflect this documented gradient as individuals were sampled across different sizes, months, and locations. Previous analyses of muscle lipid content in Pacific albacore found elevated values (13–19%) in fish sampled in higher latitudes (>40°N), a region not included in our analysis but consistent with the greater lipid stores observed for our most northerly samples in the Atlantic Ocean (Rasmussen et al., 2006). While the median value was low, many albacore sampled from the Southwest Pacific Ocean (SWPO) had elevated C:N values. Muscle lipid stores for albacore in this region vary with SST and like eastern Atlantic albacore (Goñi et al., 2011) with length (Pethybridge et al., 2015a, 2015b). Lower lipid values in many SWPO individuals could relate to recent energy expenditures associated with long distance movements for reproduction. More detailed analyses of this broader SWPO region revealed spatial C:N gradients with highest values in the Tasman Sea and Central offshore bioregions as well as in inshore waters off the south-central coast of Australia, an area of higher productivity (Pethybridge et al., 2015b). Albacore in lower latitudes off New Zealand feed on energy-rich prey like krill and have elevated muscle lipid stores while albacore feeding in more tropical waters have a more diverse but less energy-rich diet (Williams et al., 2015) and correspondingly low muscle lipid stores (Fig. 9). As documented by Dhurmeea et al. (2018), albacore in the central and southwestern Indian Ocean had low muscle lipid content, which could relate to a capital-income breeder strategy in this region in which fish rely more on local foraging than stored energy from previous feeding migrations.

809

810

811

812

813

814

815

816

817

818

819

820

821

822

823

824

825

826

827

828

829

830

831

832

833

834

835

836

837

838

839

840

841

842

843

844

845

846

Yellowfin and bigeye tunas in the Atlantic, Indian, and Pacific Oceans all have been reported to have low lipid content in white muscle tissue (<5%) (Aubourg et al., 1996; Medina et al., 1995; Peng et al., 2013; Sardenne et al., 2018, 2016; Zudaire et al., 2014). Elevated lipid stores were instead concentrated mainly in liver tissue in these species (Sardenne et al., 2018, 2016; Zudaire et al., 2014). Some outliers with elevated C:N values may have been biased by sampling location. While some samples were identified as being taken from non-dorsal muscle
locations and consequently excluded from C:N analyses, sample location was not verified for all remaining samples. Some of the higher lipid content samples may have been derived from regions with higher lipid content than dorsal muscle in tuna (Balshaw et al., 2008). For albacore, the best model (i.e. lowest AIC<sub>c</sub> value) included only dissolved oxygen (D.O.) at 375 m depth, and the deviance explained was high (65%), which could relate to differences in vertical distribution and associated dietary differences. Goñi and Arrizabalaga (2010) noted ontogenetic increases in vertical range in the water column and associated dietary differences as a possible explanation for observed pronounced seasonal changes in older (age 3 and 4 years) relative to younger albacore. Pethybridge et al. (2018) observed an increase in albacore TP in relation to D.O. at 375 m while C:N instead followed a declining trend, suggesting a possible reduction in foraging efficiency or prey quality for albacore feeding at higher TPs and/or at greater depths.

858

### 4.4. Caveats and Next Steps

860

Regional estimates of tuna residency based on differences between tuna muscle and phytoplankton δ¹³C are limited by several sources of uncertainty and variability related to trophic and physiological processes. As δ¹³C TDFs vary across taxa and environmental conditions (McCutchan et al., 2003), potential variability in tuna-phytoplankton Δ¹³C values is relatively high. Even fish held under controlled environmental conditions and fed a common diet have some degree of inherent variability in tissue δ¹³C values (Barnes et al., 2008) that will contribute to uncertainty in relating individual isotopic variability to movement or other ecological processes. Thus, our resolution of residency evaluation is fairly coarse. The uncertainty considered in our analysis for δ¹³C spacing between local phytoplankton and tuna (~2-8‰) exceeds spatial isoscapes across our sampling regions. In addition to variability associated with TDFs for tunas and underlying food web components, TP of tunas, as generalist predators, will also vary across regions and size classes (Houssard et al., 2017; Olson et al., 2016), further contributing to variability in tuna-phytoplankton δ¹³C spacing unrelated to movement.

862

863

864

865

866

867

868

869

870

871

872

873

874

875

876

877

878

879

880

881

882

883

884

885

886

887

Lipid biases on δ¹³C values were largely addressed using mathematical corrections rather than chemical extractions for our dataset, which can result in error (Abrantes et al., 2012; Sardenne et al., 2015). While our lipid-corrected δ¹³C values were negatively related to C:N, this relationship did not persist for values that were both lipid and baseline corrected (Fig. S1) suggesting the former relationship was due to underlying negative relationships between tuna condition and baseline δ¹³C gradients rather than lipid bias. Indeed, albacore, the species with the highest muscle lipid content in our dataset (Table 1), had highest C:N values in regions with low baseline δ¹³C values (e.g. Northeast Atlantic) (Magozzi et al., 2017) and low C:N values in regions with elevated baseline δ¹³C values (e.g. equatorial Pacific, Magozzi et al., 2017). The accuracy of lipid correction models largely relates to the similarity between the type of samples used to generate equation parameters and the dataset to which the equation is applied (Logan et al., 2008). Our samples were corrected using parameters that were tissue and genus specific, which were shown to give reliable predictions for yellowfin and bigeye tuna white muscle (Sardenne et al., 2015) and thus likely did not bias our δ¹³C analyses.

888

889

Residency based on isoscapes is also relative to the proximity and magnitude of bordering isoscape gradients. Evidence of movement for all three tuna species in the WCPO
based on high $\Delta^{13}$C$_{tuna-phyto}$ variability and wide isotopic niches is enhanced by the more distinct isoscapes in this region. Foraging and movements across isoscapes in this region are more easily detected and distinguished from more constrained variability potentially resulting from trophic differences, relative to other ocean regions with less pronounced isoscape gradients.

Consequently, a carbon isoscape approach to estimating residency and movements of tunas and other highly migratory species may only be suitable for assessing movements across broader isoscapes much as isotope-based food web analyses are best applied in distinguishing among isotopically distinct prey (Fry, 2006). Potential applicability of isotope-based evaluations of movement would also be limited to cases where management zones occupy multiple isoscapes.

Relatedly, the resolution of movement distance estimates will be limited to the spatial scale of bordering isoscapes. Relatively constrained movements could still impart high intraspecific isotopic variability if fish are foraging across the boundaries of bordering isoscapes.

Additionally, temporal variability in baseline, lower trophic level, and tuna $\delta^{13}$C values relative to the timing of sampling can also influence residency estimates. Timing of sampling will influence the ability to detect any movements across isoscapes with declining detections for samples collected following greater durations of residency in a given region. Movements of prey species across isoscapes rather than movements of the tunas themselves could further impart observed individual variability in tuna $\delta^{13}$C values. Seasonal variability in baseline $\delta^{13}$C values driven by changes in primary producer populations, blooms, or abiotic conditions (Goering et al., 1990; Ostrom et al., 1997) could also impart variability in tuna $\delta^{13}$C values in the absence of movements across isoscapes.

Finally, various model assumptions and an unbalanced sampling design also may have influenced our regional movement assessments. Our assessment of movement and dispersion assumed that the model is equally accurate across all Longhurst provinces, which is likely not the case. While a modelled carbon isotope baseline correction approach allows for evaluations of isotope patterns across isoscapes, modelled estimates are prone to biases, which potentially created unique values for our datasets in the GUIN province. Sample sizes also varied by species and Longhurst province, which would affect our dispersal classifications as range would be expected to increase with sample size. Dispersal estimates would be negatively biased for regions with more limited sample sizes. The model derived correction for the Suess effect also represents higher values than those reported in the literature, as it incorporates other climate effects on carbon isotope fractionation in the water column. More observations are needed to better validate the spatial and temporal changes of $\delta^{13}$C$_{DIC}$ and $\delta^{13}$C$_{POC}$ over the global ocean to have better confidence in baseline variation corrections on tuna $\delta^{13}$C values.

A variety of additional data as well as a more hypothesis-driven sample collection approach would further refine the general observations derived from the current study which took a meta-data analysis approach. A balanced sampling design across Longhurst provinces would aid in direct comparisons of movement and dispersal classifications across regions. Compound-specific stable isotope analysis (CSIA) is an emerging tool for distinguishing trophic and movement-based influences on bulk stable isotope values (McMahon and McCarthy, 2016; McMahon and Newsome, 2019) and smaller samples of AA- $\delta^{13}$C and $\delta^{15}$N data could be paired with larger bulk datasets like our tuna muscle dataset to further resolve trophic vs. movement influences on bulk patterns (e.g. Houssard et al., 2017; Popp et al., 2007). Inclusion of data from more rapid turnover tissues (e.g. liver) would also help define duration of residency through
comparisons with the slower turnover white muscle tissue used in our study and local baseline values (MacNeil et al., 2005; Sardenne et al., 2018). Comparing isotopic offsets between tissue types among individuals in a given region would provide another means of assessing homogeneity of movements and trophic status (Bean and Logan, 2019; Fry et al., 2003). Rapid turnover tissues would also more likely reflect local baselines and trophic information, consequently reducing potential biases of local baseline corrections for fish that recently migrated prior to sampling. Combining isotope data with other techniques (e.g. electronic tagging) or tracers (e.g. organochlorines, Dickhut et al., 2009; or mercury, Houssard et al., 2019) would also strengthen results. The current dataset could be better refined in future work using existing isotope or tagging data to develop a sampling plan more directly designed to quantify dynamics of residency, movements, and bioenergetic condition.

945

946

## 5. Conclusions

948

Despite limitations associated with our dataset, our analysis was able to document several important spatial-temporal and environmental relationships. Our GAMM analyses showed large-scale latitudinal gradients in tuna muscle δ¹³C consistent with previously documented spatial variability at the base of these food webs. Preservation of these gradients at higher trophic levels is consistent with a general degree of residency within the carbon isoscape regions occupied by the tuna species in the present study, but a more detailed analysis is required to fully assess residency and movement across isoscapes. For example, while albacore muscle δ¹³C values preserved underlying latitudinal baseline patterns, the slope of δ¹³C change with distance from the equator differed from the underlying baseline suggesting some degree of mixing across isoscapes. High individual variability in δ¹³C spacing between phytoplankton and tuna muscle in several regions also suggested exceptions to the general residency pattern as observed variability (e.g. 5‰ for YFT in the WARM province) was likely due to mixing across isoscapes rather than any purely trophic differences. The WCPO showed particularly high variability in baseline-corrected δ¹³C and wide isotopic niche space beyond what would be expected from trophic and inherent variability, reflective of movements across isoscapes in this region. Carbon, and through inclusion in isotopic ellipses, nitrogen stable isotopes of tropical and temperate tunas provided insight into both global and regional patterns of movements as well as biological, temporal, and environmental associations. Continued application of these approaches will contribute to our understanding of tuna movement and ecology as climate change continues to modify habitats (Erauskin-Extramiana et al., 2019) and underlying food webs (Le Borgne et al., 2011) for these mobile top predators.

949

950

951

952

953

954

955

956

957

958

959

960

961

962

963

964

965

966

967

968

969

970

971

## Acknowledgments

973

974

975

976

977

978

We thank the observer programs and the many observers who collected the samples in each of the ocean basins. We also thank the two anonymous reviewers that provided insightful comments on earlier drafts of this manuscript. This is a contribution to IMBeR-CLIOTOP Task Teams 2016-2018/01.
979

980

981

982

983

984

985

986

987

988

989

990

991

992

993

994

995

996

997

998

999

1000

1001

1002

1003

1004

1005

1006

1007

1008

1009

1010

1011

1012

1013

1014

1015

1016

1017

1018

1019

1020

1021

1022

1023

1024

References

Abelson, P.H., Hoering, T.C., 1961. Carbon isotope fractionation in formation of amino acids by photosynthetic organisms. Proc. Natl. Acad. Sci. U. S. A. 47, 623–632.

Abrantes, K.G., Semmens, J.M., Lyle, J.M., Nichols, P.D., 2012. Normalisation models for accounting for fat content in stable isotope measurements in salmonid muscle tissue. Mar. Biol. 159, 57–64.

Altabet, M.A., Deuser, W.G., Honjo, S., Steinen, C., 1991. Seasonal and depth-related changes in the source of sinking particles in the North Atlantic. Nature 354, 136–139.

Araújo, M.S., Bolnick, D.I., Machado, G., Giaretta, A.A., dos Reis, S.F., 2007. Using δ¹³C stable isotopes to quantify individual-level diet variation. Oecologia 152, 643–654.

Arrizabalaga, H., Dufour, F., Kell, L., Merino, G., Ibaibarriaga, L., Chust, G., Irigoien, X., Santiago, J., Murua, H., Fraile, I., Chifflet, M., Goikoetxea, N., Sagarminaga, Y., Aumont, O., Bopp, L., Herrera, M., Fromentin, J.M., Bonhomeau, S., 2015. Global habitat preferences of commercially valuable tuna. Deep Sea Res. Part II Top. Stud. Oceanogr. 113, 102–112. https://doi.org/10.1016/j.dsr2.2014.07.001.

Aubourg, S.P., Medina, I., Pérez-Martín, R., 1996. Polyunsaturated fatty acids in tuna phospholipids: distribution in the sn-2 location and changes during cooking. J. Agric. Food Chem. 44, 585–589.

Balshaw, S., Edwards, J.W., Ross, K.E., Daughtry, B.J., 2008. Mercury distribution in the muscular tissue of farmed southern bluefin tuna (*Thunnus maccooyi*) is inversely related to the lipid content of tissues. Food Chem. 111, 616–621.

Barnes, C., Jennings, S., Polunin, N.V.C., Lancaster, J.E., 2008. The importance of quantifying inherent variability when interpreting stable isotope field data. Oecologia 155, 227–235.

Batschelet, E., 1981. Circular Statistics in Biology. Academic Press, London.

Bean, S.B., Logan, J.M., 2019. Stable isotope analyses of cold-stunned Kemp's ridley (*Lepidochelys kempii*) sea turtles at the northern extent of their coastal range. Mar. Biol. 166, 64. https://doi.org/10.1007/s00227-019-3516-2.

Behrenfeld, M.J., Falkowski, P.G., 1997. Photosynthetic rates derived from satellite-based chlorophyll concentration. Limnol. Oceanogr. 42, 1–20.

Best, P.B., Schell, D.M., 1996. Stable isotopes in southern right whale (*Eubalaena australis*) baleen as indicators of seasonal movements, feeding and growth. Mar. Biol. 124, 483–494.

Bird, C.S., Veríssimo, A., Magozzi, S., Abrantes, K.G., Aguilar, A., Al-Reasi, H., Barnett, A., Bethea, D.M., Biais, G., Borrell, A., Bouchoucha, M., Boyle, M., Brooks, E.J., Brunnschweiler, J., Bustamante, P., Carlisle, A., Catarino, D., Caut, S., Cherel, Y., Chouvelon, T., Churchill, D., Ciancio, J., Claes, J., Colaço, A., Courtney, D.L., Cresson, P., Daly, R., de Necker, L., Endo, T., Figueiredo, I., Frisch, A.J., Hansen, J.H., Heithaus, M., Hussey, N.E., Itembu, J., Juanes, F., Kinney, M.J., Kiszka, J.J., Klarian, S.A., Kopp, D., Leaf, R., Li, Y., Lorrain, A., Madigan, D.J., Maljković, A., Malpica-Cruz, L., Matich, P., Meekan, M.G., Ménard, F., Menezes, G.M., Munroe, S.E.M., Newman, M.C., Papastamatiou, Y.P., Pethybridge, H., Plumlee, J.D., Polo-Silva, C., Quaeck-Davies, K., Raoult, V., Reum, J., Torres-Rojas, Y.E., Shiffman, D.S., Shipley, O.N., Speed, C.W., Staudinger, M.D., Teffer, A.K., Tilley, A., Valls, M., Vaudo, J.J., Wai, T.-C., Wells, R.J.D., Wyatt, A.S.J., Yool, A., Trueman, C.N., 2018a. A global perspective on the
1025 trophic geography of sharks. Nat. Ecol. Evol. 2, 299–305.

1026 https://doi.org/10.1038/s41559-017-0432-z.

1027 Bird, C.S., Veríssimo, A., Magozzi, S., Abrantes, K.G., Aguilar, A., Al-Reasi, H., Barnett, A.,

1028 Bethea, D.M., Biais, G., Borrell, A., Bouchoucha, M., Boyle, M., Brooks, E.J.,

1029 Brunnschweiler, J., Bustamante, P., Carlisle, A., Catarino, D., Caut, S., Cherel, Y.,

1030 Chouvelon, T., Churchill, D., Ciancio, J., Claes, J., Colaço, A., Courtney, D.L., Cresson,

1031 P., Daly, R., de Necker, L., Endo, T., Figueiredo, I., Frisch, A.J., Hansen, J.H., Heithaus,

1032 M., Hussey, N.E., Iitembu, J., Juanes, F., Kinney, M.J., Kiszka, J.J., Klarian, S.A., Kopp,

1033 D., Leaf, R., Li, Y., Lorrain, A., Madigan, D.J., Maljković, A., Malpica-Cruz, L., Matich,

1034 P., Meekan, M.G., Ménard, F., Menezes, G.M., Munroe, S.E.M., Newman, M.C.,

1035 Papastamatiou, Y.P., Pethybridge, H., Plumlee, J.D., Polo-Silva, C., Quaeck-Davies, K.,

1036 Raoult, V., Reum, J., Torres-Rojas, Y.E., Shiffman, D.S., Shipley, O.N., Speed, C.W.,

1037 Staudinger, M.D., Teffer, A.K., Tilley, A., Valls, M., Vaudo, J.J., Wai, T.-C., Wells,

1038 R.J.D., Wyatt, A.S.J., Yool, A., Trueman, C.N., 2018b. A global perspective on the

1039 trophic geography of sharks. Dryad, Dataset, https://doi.org/10.5061/dryad.d1f0d.

1040 Block, B.A., Jonsen, I.D., Jorgensen, S.J., Winship, A.J., Shaffer, S.A., Bograd, S.J., Hazen,

1041 E.L., Foley, D.G., Breed, G.A., Harrison, A.L., Ganong, J.E., Swithenbank, A., Castleton,

1042 M., Dewar, H., Mate, B.R., Shillinger, G.L., Schaefer, K.M., Benson, S.R., Weise, M.J.,

1043 Henry, R.W., Costa, D.P., 2011. Tracking apex marine predator movements in a dynamic

1044 ocean. Nature 475, 86–90.

1045 Bloomfield, A.L., Eldson, T.S., Walther, B.D., Gier, E.J., Gillanders, B.M., 2011. Temperature

1046 and diet affect carbon and nitrogen isotopes of fish muscle: can amino acid nitrogen

1047 isotopes explain effects? J. Exp. Mar. Biol. Ecol. 399, 48–59.

1048 Bravington, M.V., Grewe, P.M., Davies, C.R., 2016. Absolute abundance of southern bluefin

1049 tuna estimated by close-kin mark-recapture. Nat. Commun. 7, 13162.

1050 Brill, R.W., 1996. Selective advantages conferred by the high performance physiology of tunas,

1051 billfishes, and dolphin fish. Comp. Biochem. Physiol. 113, 3–15.

1052 Brill, R.W., Bigelow, K.A., Musyl, M.K., Fritsches, K.A., Warrant, E.J., 2005. Bigeye tuna

1053 (*Thunnus obesus*) behavior and physiology and their relevance to stock assessments and

1054 fishery biology. Col Vol Sci Pap ICCAT 57, 142–161.

1055 Brown, J.N., Langlais, C., Maes, C., 2014. Zonal structure and variability of the Western Pacific

1056 dynamic warm pool edge in CMIP5. Clim. Dyn. 42, 3061–3076.

1057 Carey, F.G., Teal, J.M., 1969. Mako and porbeagle - warm-bodied sharks. Comp. Biochem.

1058 Physiol. 28, 199–204.

1059 Carruthers, T.R., McAllister, M.K., Ahrens, R.N.M., 2010. Simulating spatial dynamics to

1060 evaluate methods of deriving abundance indices for tropical tunas. Can. J. Fish. Aquat.

1061 Sci. 67, 1409–1427.

1062 Caut, S., Angulo, E., Courchamp, F., 2009. Variation in discrimination factors (Δ¹⁵N and Δ¹³C):

1063 the effect of diet isotopic values and applications for diet reconstruction. J. Appl. Ecol.

1064 46, 443–453.

1065 Cherel, Y., Hobson, K.A., 2007. Geographical variation in carbon stable isotope signatures of

1066 marine predators: a tool to investigate their foraging areas in the Southern Ocean. Mar.

1067 Ecol. Prog. Ser. 329, 281–287.

1068 Childers, J., Snyder, S., Kohin, S., 2011. Migration and behavior of juvenile North Pacific

1069 albacore (*Thunnus alalunga*). Fish. Oceanogr. 20, 157–173.
1070 Choy, C.A., Popp, B.N., Hannides, C.C.S., Drazen, J.C., 2015. Trophic structure and food resources of epipelagic and mesopelagic fishes in the North Pacific Subtropical Gyre ecosystem inferred from nitrogen isotopic compositions. Limnol. Oceanogr. 60, 1156–1171.

1074 Cosgrove, R., Arregui, I., Arrizabalaga, H., Goñi, N., Sheridan, M., 2014. New insights to behaviour of North Atlantic albacore tuna (*Thunnus alalunga*) observed with pop-up satellite archival tags. Fish. Res. 150, 89–99.

1077 Cunjak, R.A., Roussel, J.M., Gray, M.A., Dietrich, J.P., Cartwright, D.F., Munkittrick, K.R., Jardine, T.D., 2005. Using stable isotope analysis with telemetry or mark-recapture data to identify fish movement and foraging. Oecologia 144, 1–11.

1080 de Boyer Montégut, C., Madec, G., Fischer, A.S., Lazar, A., Iudicone, D., 2004. Mixed layer depth over the global ocean: An examination of profile data and a profile-based climatology. J. Geophys. Res. Oceans 109, C12003.

1083 DeNiro, M.J., Epstein, S., 1978. Influence of diet on the distribution of carbon isotopes in animals. Geochim. Cosmochim. Acta 42, 495–506.

1085 Dhurmeea, Z., Pethybridge, H., Appadoo, C., Bodin, N., 2018. Lipid and fatty acid dynamics in mature female albacore tuna (*Thunnus alalunga*) in the western Indian Ocean. PLoS ONE 13, e0194558. https://doi.org/10.1371/journal.pone.0194558.

1088 Dickhut, R.M., Deshpande, A.D., Cincinelli, A., Cochran, M.A., Corsolini, S., Brill, R.W., Secor, D.H., Graves, J.E., 2009. Atlantic bluefin tuna (*Thunnus thynnus*) population dynamics delineated by organochlorine tracers. Environ. Sci. Technol. 43, 8522–8527.

1091 Dickson, K., Graham, J., 2004. Evolution and consequences of endothermy in fishes. Physiol. Biochem. Zool. 77, 998–1018.

1093 Duffy, L.M., Kuhnert, P.M., Pethybridge, H.R., Young, J.W., Olson, R.J., Logan, J.M., Goñi, N., Romanov, E., Allain, V., Staudinger, M.D., Abecassis, M., Choy, C.A., Hobday, A.J., Simier, M., Galván-Magana, F., Potier, M., Ménard, F., 2017. Global trophic ecology of yellowfin, bigeye, and albacore tunas: Understanding predation on micronekton communities at ocean-basin scales. Deep Sea Res. II 140, 55–73.

1098 Dunn, J., Ridgway, K., 2002. Mapping ocean properties in regions of complex topography. Deep Sea Res. Part Oceanogr. Res. Pap. 49, 591–604.

1100 Erauskin-Extramiana, M., Arrizabalaga, H., Hobday, A.J., Cabré, A., Ibaibarriaga, L., Arregui, I., Murua, H., Chust, G., 2019. Large-scale distribution of tuna species in a warming ocean. Glob. Change Biol. 25, 2043–2060.

1103 Estrada, J.A., Lutcavage, M., Thorrold, S.R., 2005. Diet and trophic position of Atlantic bluefin tuna (*Thunnus thynnus*) inferred from stable carbon and nitrogen isotope analysis. Mar. Biol. 147, 37–45.

1106 Evans, K., Patterson, T., Pedersen, M., 2011. Movement patterns of yellowfin tuna in the Coral Sea region: defining connectivity with stocks in the western Pacific Ocean region. Final Report 2008/804 to the Australian Fisheries Management Authority, Canberra.

1109 Evans, K., Young, J.W., Nicol, S., Kolody, D., Allain, V., Bell, J., Brown, J.N., Ganachaud, A., Hobday, A.J., Hunt, B., Innes, J., 2015. Optimising fisheries management in relation to tuna catches in the western central Pacific Ocean: A review of research priorities and opportunities. Mar. Policy 59, 94–104.

1113 Eveson, J.P., Hobday, A.J., Hartog, J.R., Spillman, C.M., Rough, K.M., 2015. Seasonal forecasting of tuna habitat in the Great Australian Bight. Fish. Res. 170, 39–49.
1115 Fonteneau, A., Hallier, J.P., 2015. Fifty years of dart tag recoveries for tropical tuna: A global
1116 comparison of results for the western Pacific, eastern Pacific, Atlantic, and Indian
1117 Oceans. Fish. Res. 163, 7–22.

1118 France, R.L., Peters, R.H., 1997. Ecosystem differences in the trophic enrichment of ¹³C in
1119 aquatic food webs. Can. J. Fish. Aquat. Sci. 54, 1255–1258.

1120 Fry, B., 2006. Stable isotope ecology. Springer, New York, NY, USA.

1121 Fry, B., 2002. Stable isotope indicators of habitat use by Mississippi River fish. J. North Am.
1122 Benthol. Soc. 21, 676–685.

1123 Fry, B., 1981. Natural stable carbon isotope tag traces Texas shrimp migrations. Fish. Bull. 79,
1124 337–345.

1125 Fry, B., Baltz, D.M., Benfield, M.C., Fleeger, J.W., Gace, A., Haas, H.L., Quiñones-Rivera, Z.J.,
1126 2003. Stable isotope indicators of movement and residency for brown shrimp
1127 (*Farfantepenaeus aztecus*) in coastal Louisiana marshscapes. Estuaries 26, 82–97.

1128 Fry, B., Wainright, S.C., 1991. Diatom sources of ¹³C-rich carbon in marine food webs. Mar.
1129 Ecol. Prog. Ser. 76, 149–157.

1130 Goering, J., Alexander, V., Haubenstock, N., 1990. Seasonal variability of stable carbon and
1131 nitrogen isotope ratios of organisms in a North Pacific Bay. Estuar. Coast. Shelf Sci. 30,
1132 239–260. https://doi.org/10.1016/0272-7714(90)90050-2.

1133 Goñi, N., Arrizabalaga, H., 2010. Seasonal and interannual variability of fat content of juvenile
1134 albacore (*Thunnus alalunga*) and bluefin (*Thunnus thynnus*) tuna during their feeding
1135 migration to the Bay of Biscay. Prog. Oceanogr. 86, 115–123.

1136 Goñi, N., Logan, J.M., Arrizabalaga, H., Jarry, M., Lutcavage, M., 2011. Diet of albacore tuna
1137 (*Thunnus alalunga*) in the Bay of Biscay and Mediterranean Sea using stomach content
1138 and stable isotope analyses. Mar. Biol. 158, 1057–1073.

1139 Graham, B.S., Koch, P., Newsome, S., McMahon, K.W., Aurioles, D., 2010. Using isoscapes to
1140 trace the movements and foraging behavior of top predators in oceanic ecosystems, in:
1141 West, J.B., Bowen, G.J., Dawson, T.E., Tu, K.P. (Eds.), Isoscapes: Understanding
1142 Movement, Pattern and Processes on Earth through Isotope Mapping. Springer, New
1143 York, NY, USA.

1144 Grecian, W.J., McGill, R.A.R., Phillips, R.A., Ryan, P.G., Furness, R.W., 2015. Quantifying
1145 variation in δ¹³C and δ¹⁵N isotopes within and between feathers and individuals: Is one
1146 sample enough? Mar. Biol. 162, 733–741.

1147 Hinga, K.R., Arthur, M.A., Pilson, M.E.Q., Whitaker, D., 1994. Carbon isotope fractionation by
1148 marine phytoplankton in culture: the effects of CO₂ concentration, pH, temperature, and
1149 species. Glob. Biogeochem. Cycles 8, 91–102.

1150 Hobday, A.J., Hartog, J.R., Spillman, C.M., Alves, O., 2011. Seasonal forecasting of tuna habitat
1151 for dynamic spatial management. Can. J. Fish. Aquat. Sci. 68, 898–911.

1152 Hobson, K.A., 1999. Tracing origins and migration of wildlife using stable isotopes: a review.
1153 Oecologia 120, 314–326.

1154 Hobson, K.A., Cherel, Y., 2006. Isotopic reconstruction of marine food webs using cephalopod
1155 beaks: new insight from captively raised *Sepia officinalis*. Can. J. Zool. 84, 766–770.

1156 Hoffman, J.C., Sutton, T.T., 2010. Lipid correction for carbon stable isotope analysis of deep-sea
1157 fishes. Deep Sea Res. Part Oceanogr. Res. Pap. 57, 956–964.

1158 Houssard, P., Lorrain, A., Tremblay-Boyer, L., Allain, V., Graham, B.S., Menkes, C.E.,
1159 Pethybridge, H., Couturier, L.I.E., Point, D., Leroy, B., Receveur, A., Hunt, B.P.V.,
1160 Vourey, E., Bonnet, S., Rodier, M., Raimbault, P., Feunteun, E., Kuhnert, P., Munaron,
1161 J.M., Lebreton, B., Otake, T., Letourneur, Y., 2017. Trophic position increases with thermocline depth in yellowfin and bigeye tuna across the Western and Central Pacific Ocean. Prog. Oceanogr. 154, 49–63.

1162 Houssard, P., Point, D., Tremblay-Boyer, L., Allain, V., Pethybridge, H., Mashou, J., Ferriss, B.E., Baya, P.A., Lagane, C., Menkes, C.E., Letourneur, Y., Lorrain, A., 2019. A model of mercury distribution in tuna from the Western and Central Pacific Ocean: influence of physiology, ecology and environmental factors. Environ. Sci. Technol. 53, 1422–1431.

1163 Jackson, A.L., Inger, R., Parnell, A.C., Bearhop, S., 2011. Comparing isotopic niche widths among and within communities: SIBER - Stable Isotope Bayesian Ellipses in R. J. Anim. Ecol. 80, 595–602.

1171 Jaeger, A., Lecomte, V.J., Weimerskirch, H., Richard, P., Cherel, Y., 2010. Seabird satellite tracking validates the use of latitudinal isoscapes to depict predators' foraging areas in the Southern Ocean. Rapid Commun. Mass Spectrom. 24, 3456–3460.

1174 Lam, C.H., Galuardi, B., Lutcavage, M.E., 2014. Movements and oceanographic associations of bigeye tuna (*Thunnus obesus*) in the Northwestern Atlantic. Can. J. Fish. Aquat. Sci. 71, 1529–1543.

1177 Laws, E., Popp, B.N., Bidigare, R.R., Kennicutt, M.C., Macko, S.A., 1995. Dependence of phytoplankton carbon isotopic composition on growth rate and [CO₂]aq: Theoretical considerations and experimental results. Geochim. Cosmochim. Acta 59, 1131–1138.

1180 Layman, C.A., Arrington, D.A., Montaña, C.G., Post, D.M., 2007. Can stable isotope ratios provide for community-wide measures of trophic structure? Ecology 88, 42–48.

1182 Le Borgne, R., Allain, V., Griffiths, S.P., Matear, R.J., McKinnon, A.D., Richardson, A.J., Young, J.W., 2011. Chapter 4, Vulnerability of open ocean food webs in the tropical Pacific to climate change, in: Bell, J.D., Hobday, A.J. (Eds.), Vulnerability of Tropical Pacific Fisheries and Aquaculture to Climate Change. Secretariat of the Pacific Community, New Caledonia.

1187 Lehodey, P., Senina, I., Wibawa, T.A., Titaud, O., Calmettes, B., Conchon, A., Tranchant, B., Gaspar, P., 2018. Operational modelling of bigeye tuna (*Thunnus obesus*) spatial dynamics in the Indonesian region. Mar. Pollut. Bull. 131, 19–32.

1190 Logan, J., 2009. Tracking diet and movement of Atlantic bluefin tuna (*Thunnus thynnus*) using carbon and nitrogen stable isotopes. University of New Hampshire, Ph.D. Thesis. University of New Hampshire, Durham, NH, USA.

1193 Logan, J.M., Golet, W., Lutcavage, M.E., 2015. Diet and condition of Atlantic bluefin tuna (*Thunnus thynnus*) in the Gulf of Maine, 2004-2008. Environ. Biol. Fishes 98, 1411–1430.

1196 Logan, J.M., Jardine, T.D., Miller, T.J., Bunn, S.E., Cunjak, R.A., Lutcavage, M.E., 2008. Lipid corrections in carbon and nitrogen stable isotope analyses: comparison of chemical extraction and modelling methods. J. Anim. Ecol. 77, 838–846.

1199 Logan, J.M., Lutcavage, M.E., 2010. Stable isotope dynamics in elasmobranch fishes. Hydrobiologia 644, 231–244.

1201 Longhurst, A.R., 1998. Ecological Geography of the Sea. Academic Press, San Diego, CA, USA.

1202 Lorrain, A., Graham, B.S., Popp, B.N., Allain, V., Olson, R.J., Hunt, B.P.V., Potier, M., Fry, B., Galván-Magaña, F., Menkes, C.E.R., Kaehler, S., Ménard, F., 2015. Nitrogen isotopic baselines and implications for estimating foraging habitat and trophic position of yellowfin tuna in the Indian and Pacific Oceans. Deep-Sea Res. II 113, 188–198.
1206 Lorrain, A., Pethybridge, H., Cassar, N., Receveur, A., Allain, V., Bodin, N., Bopp, L., Choy, C.A., Duffy, L., Fry, B., Goñi, N., Graham, B.S., Hobday, A.J., Logan, J.M., Ménard, F., Menkes, C., Olson, R.J., Point, D., Revill, A.T., Somes, C.J., Young, J.W., 2020. Trends in tuna carbon isotopes reflect global changes in pelagic phytoplankton communities. Glob. Change Biol. 26(2), 458-470.

1211 MacNeil, M.A., Skomal, G.B., Fisk, A.T., 2005. Stable isotopes from multiple tissues reveal diet switching in sharks. Mar. Ecol. Prog. Ser. 302, 199-206.

1213 Madigan, D.J., Litvin, S.Y., Popp, B.N., Carlisle, A.B., Farwell, C.J., Block, B.A., 2012. Tissue turnover rates and isotopic trophic discrimination factors in the endothermic teleost, Pacific bluefin tuna (*Thunnus orientalis*). PLoS ONE 7(11): e49220.
doi:10.1371/journal.pone.0049220.

1217 Magozzi, S., Yool, A., Vander Zanden, H.B., Wunder, M.B., Trueman, C.N., 2017. Using ocean models to predict spatial and temporal variation in marine carbon isoscapes. Ecosphere 8, e01763.

1220 Marshall, C.T., Yaragina, N.A., Lambert, Y., Kjesbu, O., 1999. Total lipid energy as a proxy for total egg production by fish stocks. Nature 402, 288-290.

1222 Martino, J.C., Doubleday, Z.A., Gillanders, B.M., 2019. Metabolic effects on carbon isotope biomarkers in fish. Ecol. Indic. 97, 10-16.

1224 Mather, F.J., Mason Jr., J.M., Jones, A.C., 1995. Historical document: life history and fisheries of Atlantic bluefin tuna. NOAA Technical Memorandum NMFS-SEFSC-370. Miami, Fla. 165 pp.

1227 McCutchan, J.H., Lewis Jr, W.M., Kendall, C., McGrath, C.C., 2003. Variation in trophic shift for stable isotope ratios of carbon, nitrogen, and sulfur. Oikos 102, 378-390.

1229 McMahon, K.W., Hamady, L.L., Thorrold, S.R., 2013. A review of ecogeochemistry approaches to estimating movements of marine animals. Limnol. Oceanogr. 58, 697-714.

1231 McMahon, K.W., McCarthy, M.D., 2016. Embracing variability in amino acid δ¹⁵N fractionation: mechanisms, implications, and applications for trophic ecology. Ecosphere 7, e01511.

1234 McMahon, K.W., Newsome, S.D., 2019. Amino acid isotope analysis: a new frontier in studies of animal migration and foraging ecology, in: Hobson, K.A., Wassenaar, L.I. (Eds.), Tracking Animal Migrations with Stable Isotopes, Second Edition. Academic Press, San Diego, CA, USA pp. 173-190.

1238 Medina, I., Aubourg, S.P., Pérez Martín, R., 1995. Composition of phospholipids of white muscle of six tuna species. Lipids 30, 1127-1135.

1240 Ménard, F., Lorrain, A., Potier, M., Marsac, F., 2007. Isotopic evidence of distinct feeding ecologies and movement patterns in two migratory predators (yellowfin tuna and swordfish) of the western Indian Ocean. Mar. Biol. 153, 141-152.

1243 Miller, T.W., Brodeur, R.D., Rau, G.H., 2008. Carbon stable isotopes reveal relative contribution of shelf-slope production to the northern California Current pelagic community. Limnol. Oceanogr. 53, 1493-1503.

1246 Naimi, B., 2015. usdm: Uncertainty analysis for species distribution models (R package version 1.1-12). R Documentation http://www.r documentation. org/packages/usdm.

1248 Navarro, J., Coll, M., Somes, C.J., Olson, R.J., 2013. Trophic niche of squids: Insights from isotopic data in marine systems worldwide. Deep Sea Res. II 95, 93-102.

1250 Newsome, S.D., Etnier, M.A., Gifford-Gonzalez, D., Phillips, D.L., van Tuinen, M., Hadley, E.A., Costa, D.P., Kennett, D.J., Guilderson, T.P., Koch, P.L., 2007. The shifting baseline
1252 of northern fur seal ecology in the northeast Pacific Ocean. Proc. Natl. Acad. Sci. 104, 9709-9714.

1254 Nikolic, N., Morandeau, G., Hoarau, L., West, W., Arrizabalaga, H., Hoyle, S., Nicol, S.J., Bourjea, J., Puech, A., Farley, J.H., Williams, A.J., Fonteneau, A., 2017. Review of albacore tuna, *Thunnus alalunga*, biology, fisheries and management. Rev. Fish Biol. Fish. 27, 775-810.

1258 Olson, R.J., Popp, B.N., Graham, B., López-Ibarra, G.A., Galván-Magaña, F., Lennert-Cody, C.E., Bocanegra-Castillo, N., Wallsgrove, N.J., Gier, E., Alatorre-Ramírez, V., Ballance, L.T., Fry, B., 2010. Food-web inferences of stable isotope spatial patterns in copepods and yellowfin tuna in the pelagic eastern Pacific Ocean. Prog. Oceanogr. 86, 124-138.

1262 Olson, R.J., Young, J.W., Ménard, F., Potier, M., Allain, V., Goñi, N., Logan, J.M., Galván-Magaña, F., 2016. Bioenergetics, trophic ecology, and niche separation of tunas. Adv. Mar. Biol. 74, 199-344.

1265 Oschlies, A., Schulz, K., Riebesell, U., Schmittner, A., 2008. Simulated 21st century's increase in oceanic subboxia by CO₂-enhanced biotic carbon export. Glob. Biogeochem. Cycles 22, GB4008.

1268 Ostrom, N.E., Macko, S.A., Deibel, D., Thompson, R.J., 1997. Seasonal variation in the stable carbon and nitrogen isotope biogeochemistry of a coastal cold ocean environment. Geochim. Cosmochim. Acta 61, 2929-2942. https://doi.org/10.1016/S0016-7037(97)00131-2.

1272 Park, R., Epstein, S., 1961. Metabolic fractionation of C¹³ & C¹² in plants. Plant Physiol. 36, 133-138.

1274 Parnell, A.C., Inger, R., Bearhop, S., Jackson, A.L., 2010. Source partitioning using stable isotopes: coping with too much variation. PLoS ONE 5(3): e9672. doi:10.1371/journal.pone.0009672.

1277 Peng, S., Chen, C., Shi, Z., Wang, L., 2013. Amino acid and fatty acid composition of the muscle tissue of yellowfin tuna (*Thunnus albacares*) and bigeye tuna (*Thunnus obesus*). J. Food Nutr. Res. 1, 42-45.

1280 Pethybridge, H., Choy, C.A., Logan, J.M., Allain, V., Lorrain, A., Bodin, N., Somes, C.J., Young, J., Ménard, F., Langlais, C., Duffy, L., Hobday, A.J., Kuhnert, P., Fry, B., Menkes, C., Olson, R.J., 2018. A global meta-analysis of marine predator nitrogen stable isotopes: Relationships between trophic structure and environmental conditions. Glob. Ecol. Biogeogr. 27, 1043-1055.

1285 Pethybridge, H.R., Parrish, C.C., Morrongiello, J., Young, J.W., Farley, J.H., Gunasekera, R.M., Nichols, P.D., 2015a. Spatial patterns and temperature predictions of tuna fatty acids: tracing essential nutrients and changes in primary producers. PLoS ONE 10, e0131598. doi:10.1371/journal.pone.0131598.

1289 Pethybridge, H.R., Young, J.W., Kuhnert, P.M., Farley, J.H., 2015b. Using stable isotopes of albacore tuna and predictive models to characterize bioregions and examine ecological change in the SW Pacific Ocean. Prog. Oceanogr. 134, 293-303.

1292 Pinnegar, J.K., Polunin, N.V.C., 1999. Differential fractionation of δ¹³C and δ¹⁵N among fish tissues: implications for the study of trophic interactions. Funct. Ecol. 13, 225-231.

1294 Popp, B., Graham, B., Olson, R., Hannides, C., Lott, M., López-Ibarra, G., Galván-Magaña, F., Fry, B., 2007. Insight into the trophic ecology of yellowfin tuna, *Thunnus albacares*, from compound-specific nitrogen isotope analysis of proteinaceous amino acids, in:
1297 Dawson, T., Seigwolf, R. (Eds.), Isotopes as Tracers of Ecological Change. Elsevier

1298 Academic Press, San Diego, CA, USA pp. 173–190.

1299 Popp, B., Laws, E., Bidigare, R., Dore, J., Hanson, K., Wakeham, S., 1998. Effect of phytoplankton cell geometry on carbon isotopic fractionation. Geochim. Cosmochim.

1300 Acta 62, 69–77.

1301 Post, D.M., 2002. Using stable isotopes to estimate trophic position: models, methods, and assumptions. Ecology 83, 703–718.

1302 Post, D.M., Layman, C.A., Arrington, D.A., Takimoto, G., Quattrochi, J., Montaña, C.G., 2007.

1303 Getting to the fat of the matter: models, methods and assumptions for dealing with lipids in stable isotope analyses. Oecologia 152, 179–189.

1304 Quay, P., Sonnerup, R., Munro, D., Sweeney, C., 2017. Anthropogenic CO₂ accumulation and uptake rates in the Pacific Ocean based on changes in the ¹³C/¹²C of dissolved inorganic carbon. Glob. Biogeochem. Cycles 31, 59–80.

1305 R Core Team, 2019. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL http://www.R-project.org/.

1306 Rasmussen, R.S., Morrissey, M.T., Carroll, S., 2006. Effect of seasonality, location, and size on lipid content in North Pacific troll-caught albacore tuna (*Thunnus alalunga*). J. Aquat. Food Prod. Technol. 15, 73–86.

1307 Reid, W.D.K., Sweeting, C.J., Wigham, B.D., McGill, R.A.R., Polunin, N.V.C., 2016. Isotopic niche variability in macroconsumers of the East Scotia Ridge (Southern Ocean) hydrothermal vents: What more can we learn from an ellipse? Mar. Ecol. Prog. Ser. 542, 13–24.

1308 Reynolds, R.W., Smith, T.M., Liu, C., Chelton, D.B., Casey, K.S., Schlax, M.G., 2007. Daily high resolution blended analyses for sea surface temperature. J. Clim. 20, 5473–5496.

1309 Richards, T.M., Gipson, E.E., Cook, A., Sutton, T.T., Wells, R.J.D., 2019. Trophic ecology of meso- and bathypelagic predatory fishes in the Gulf of Mexico. ICES J. Mar. Sci. 76, 662–672.

1310 Ridgway, K., Dunn, J., Wilkin, J., 2002. Ocean interpolation by four-dimensional weighted least squares-Application to the waters around Australasia. J. Atmospheric Ocean. Technol. 19, 1357–1375.

1311 Rubenstein, D.R., Hobson, K.A., 2004. From birds to butterflies: animal movement patterns and stable isotopes. Trends Ecol. Evol. 19, 256–263.

1312 Sardenne, F., Chassot, E., Fouché, E., Ménard, F., Lucas, V., Bodin, N., 2016. Are condition factors powerful proxies of energy content in wild tropical tunas? Ecol. Indic. 71, 467–476.

1313 Sardenne, F., Dunn, J., Wilkin, J., 2002. Ocean interpolation by four-dimensional weighted least squares-Application to the waters around Australasia. J. Atmospheric Ocean. Technol. 19, 1357–1375.

1314 Oceanogr. 25, 277–291.

1315 Rooker, J.R., Wells, R.J.D., Itano, D.G., Thorrold, S.R., Lee, J.M., 2016. Natal origin and population connectivity of bigeye and yellowfin tuna in the Pacific Ocean. Fish.

1316 Oceanogr. 25, 277–291.

1317 Rubenstein, D.R., Hobson, K.A., 2004. From birds to butterflies: animal movement patterns and stable isotopes. Trends Ecol. Evol. 19, 256–263.

1318 Sardenne, F., Chassot, E., Fouché, E., Ménard, F., Lucas, V., Bodin, N., 2016. Are condition factors powerful proxies of energy content in wild tropical tunas? Ecol. Indic. 71, 467–476.

1319 Sardenne, F., Diaha, N.C., Amandé, M.J., Zudaire, I., Couturier, L.I.E., Metral, L., Le Grand, F., Bodin, N., 2018. Habitat and size influence on the trophic niche of co-occurring tropical tunas in the eastern Atlantic Ocean. Can. J. Fish. Aquat. Sci. 76, 69–80.

1320 Sardenne, F., Ménard, F., Degroote, M., Fouché, E., Guillou, G., Lebreton, B., Hollanda, S., Bodin, N., 2015. Methods of lipid-normalization for multi-tissue stable isotope analyses in tropical tuna. Rapid Commun. Mass Spectrom. 29, 1253–1267.
1341 Schaefer, K., Fuller, D., 2002. Movements, behavior, and habitat selection of bigeye tuna (*Thunnus obesus*) in the eastern equatorial Pacific, ascertained through archival tags.
1342
1343 Fish. Bull. 100, 765–788.

1344 Schaefer, K., Fuller, D., Block, B., 2007. Movements, behavior, and habitat utilization of
1345 yellowfin tuna (*Thunnus albacares*) in the northeastern Pacific Ocean, ascertained
1346 through archival tag data. Mar. Biol. 152, 503–525.

1347 Schaefer, K., Fuller, D., Hampton, J., Caillot, S., Leroy, B., Itano, D., 2015. Movements, 1348 dispersion, and mixing of bigeye tuna (*Thunnus obesus*) tagged and released in the 1349 equatorial Central Pacific Ocean, with conventional and archival tags. Fish. Res. 161, 1350 336–355.

1351 Schaefer, K.M., Fuller, D.W., 2009. Horizontal movements of bigeye tuna (*Thunnus obesus*) in 1352 the eastern Pacific Ocean, as determined from conventional and archival tagging 1353 experiments initiated during 2000–2005. Inter-Am. Trop. Tuna Comm. Bull. 24, 191– 1354 247.

1355 Schell, D.M., Saupe, S.M., Haubenstock, N., 1989. Bowhead whale (*Balaena mysticetus*) growth 1356 and feeding as estimated by δ¹³C techniques. Mar. Biol. 103, 433–443.

1357 Schmittner, A., Gruber, N., Mix, A.C., Key, R.M., Tagliabue, A., Westberry, T.K., 2013. 1358 Biology and air-sea gas exchange controls on the distribution of carbon isotope ratios 1359 (δ¹³C) in the ocean. Biogeosciences 10, 5793–5816.

1360 Schmittner, A., Somes, C.J. 2016. Complementary constraints from carbon (¹³C) and nitrogen 1361 (¹⁵N) isotopes on the glacial ocean's soft-tissue biological pump. Paleooceanogr. 1362 Paleoclimatology 31, 669-693.

1363 Sepulveda, C.A., Graham, J.B., Bernal, D., 2007. Aerobic metabolic rates of swimming juvenile 1364 mako sharks, *Isurus oxyrinchus*. Mar. Biol. 152, 1087–1094.

1365 Shipley, O.N., Gallagher, A.J., Shiffman, D.S., Kaufman, L., Hammerschlag, N., 2019. Diverse 1366 resource-use strategies in a large-bodied marine predator guild: evidence from differential 1367 use of resource subsidies and intraspecific isotopic variation. Mar. Ecol. Prog. Ser. 623, 1368 71–83.

1369 Sibert, J.R., Hampton, J., Fournier, D.A., Bills, P.J., 1999. An advection-diffusion-reaction 1370 model for the estimation of fish movement parameters from tagging data, with 1371 application to skipjack tuna (*Katsuwonus pelamis*). Can. J. Fish. Aquat. Sci. 56, 925–938.

1372 Stramma, L., Prince, E.D., Schmidtko, S., Luo, J., Hoolihan, J.P., Visbeck, M., Wallace, D.W.R., 1373 Brandt, P., Körtzinger, A., 2012. Expansion of oxygen minimum zones may reduce 1374 available habitat for tropical pelagic fishes. Nat. Clim. Change 2, 33–37.

1375 Sund, P.N., Blackburn, M., Williams, F., 1981. Tunas and their environment in the Pacific 1376 Ocean: a review. Oceanogr. Mar. Biol. Annu. Rev. 19, 443–512.

1377 Sweeting, C.J., Barry, J.T., Polunin, N.V.C., Jennings, S., 2007. Effects of body size and 1378 environment on diet-tissue δ¹³C fractionation in fishes. J. Exp. Mar. Biol. Ecol. 352, 165– 1379 176.

1380 Takai, N., Onaka, S., Ikeda, Y., Yatsu, A., Kidokoro, H., Sakamoto, W., 2000. Geographical 1381 variations in carbon and nitrogen stable isotope ratios in squid. J. Mar. Biol. Assoc. U. K. 1382 80, 675–684.

1383 Thomas, S.M., Crowther, T.W., 2015. Predicting rates of isotopic turnover across the animal 1384 kingdom: a synthesis of existing data. J. Anim. Ecol. 84, 861–870.

1385 Trueman, C.N., MacKenzie, K.M., Palmer, M.R., 2012. Identifying migrations in marine fishes 1386 through stable-isotope analysis. J. Fish Biol. 81, 826–847.
1387 Trueman, C.N., St. John Glew, K., 2019. Isotopic tracking of marine animal movement, in:

1388 Hobson, K.A., Wassenaar, L.I. (Eds.), Tracking Animal Migrations with Stable Isotopes,

1389 Second Edition. Academic Press, San Diego, CA, USA.

1390 Turner, T.F., Collyer, M.L., Krabbenhoft, T.J., 2010. A general hypothesis-testing framework for

1391 stable isotope ratios in ecological studies. Ecology 91, 2227–2233.

1392 Vander Zanden, M.J., Clayton, M.K., Moody, E.K., Solomon, C.T., Weidel, B.C., 2015. Stable

1393 isotope turnover and half-life in animal tissues: a literature synthesis. PLoS ONE 10,

1394 e0116182. doi:10.1371/journal.pone.0116182.

1395 Weaver, A.J., Eby, M., Wiebe, E.C., Bitz, C.M., Duffy, P.B., Ewen, T.L., Fanning, A.F.,

1396 Holland, M.M., MacFadyen, A., Matthews, H.D., Meissner, K.J., Saenko, O., Schmittner,

1397 A., Wang, H., Yoshimori, M., 2001. The UVic earth system climate model: Model

1398 description, climatology, and applications to past, present and future climates.

1399 Atmosphere-Ocean 39, 361–428.

1400 West, J.B., Bowen, G.J., Dawson, T.E., Tu, K.P., 2010. Isoscapes: Understanding movement,

1401 pattern, and processing on Earth through isotope mapping. Springer, New York, NY,

1402 USA.

1403 Williams, A.J., Allain, V., Nicol, S.J., Evans, K.J., Hoyle, S.D., Dupoux, C., Vourey, E., Dubosc,

1404 J., 2015. Vertical behavior and diet of albacore tuna (*Thunnus alalunga*) vary with

1405 latitude in the South Pacific Ocean. Deep Sea Res. Part II Top. Stud. Oceanogr. 113,

1406 154–169.

1407 Witteveen, B.H., Worthy, G.A.J., Roth, J.D., 2009. Tracing migratory movements of breeding

1408 North Pacific humpback whales using stable isotope analysis. Mar. Ecol. Prog. Ser. 393,

1409 173–183.

1410 Wood, S., 2006. Generalized additive models: An introduction with R. CRC Press, Boca Raton,

1411 FL, USA.

1412 Zudaire, I., Murua, H., Grande, M., Pernet, F., Bodin, N., 2014. Accumulation and mobilization

1413 of lipids in relation to reproduction of yellowfin tuna (*Thunnus albacares*) in the Western

1414 Indian Ocean. Fish. Res. 160, 50–59.

1415
**Figure 1.** Diagram of representative food chain with tuna occupying a trophic position of 5, showing upper and lower bounds of estimated $\Delta^{13}$C between tuna and phytoplankton. The upper estimate (blue) assumes all trophic transfers impart a 2‰ increase in $\delta^{13}$C to the consumer, which translates to a $\Delta^{13}C_{tuna-phyto}$ value of -8‰. The lower estimate (red) assumes all trophic transfers impart a 0.5‰ increase in $\delta^{13}$C to the consumer, which translates to a $\Delta^{13}C_{tuna-phyto}$ value of -2‰. The 0.5‰ and 2‰ trophic discrimination factor (TDF) estimates are based on the lower and upper ends, respectively, of reported TDFs for representative taxa in pelagic food webs that include tuna. Values in black are best estimates based on controlled feeding experiments on taxa representative of each trophic step:

ªFrance and Peters 1997; ᵇSweeting et al. 2007; ᶜCaut et al. 2009; ᵈHobson and Cherel 2006;

ᶜMadigan et al. 2012). Phytoplankton and tuna images: Tracey Saxby, IAN Image Library (http://ian.umecs.edu/imagelibrary/). Zooplankton image: Jane Hawkey, IAN Image Library.

Anchovy image: Dieter Tracey, IAN Image Library. Squid image: Jane Thomas, IAN Image Library.

**Figure 2.** Map of Longhurst provinces and overlapping sampling locations for yellowfin, bigeye, and albacore tuna. See Supplementary Table S1 for code names and related sampling details.

**Figure 3.** Boxplots of yellowfin tuna (A) lipid and Suess corrected carbon isotopes values of tuna ($\delta^{13}C_{LScorr}$); (B) difference between lipid corrected carbon isotope values of tuna and model-derived Suess corrected phytoplankton carbon isotope values ($\Delta^{13}C_{tuna-phyto}$); (C) C:N ratios; and (D) standard ellipse area (SEA) as an indicator of isotopic niche area, for each Longhurst province (Table S1) in this study. Biome: C – Coastal, T – Trade wind; W – Westerly. For boxplots (A-C), the dotted line corresponds to the global median value, the black line is the region median, the shaded box includes the interquartile range (IQR), whiskers represent 1.5 IQR, and individual points are outliers outside of the IQR. For (B) the dark and light grey band highlight $\Delta^{13}C_{tuna-phyto}$ values that fall between -4 to -6‰ and -2 to -8‰, respectively (that define limited and moderate degrees of tuna movement; see Table 2). SEA density plots (D) represent the credibility intervals of the Bayesian SEA with shaded boxes encompassing the 50th, 75th, and 95th credibility intervals for each Longhurst province. Black circles and red x-marks represent the mode SEA and SEA corrected for small sample sizes ($SEA_c$).

**Figure 4.** Boxplots of bigeye tuna (A) lipid and Suess corrected carbon isotopes values of tuna ($\delta^{13}C_{LScorr}$); (B) difference between lipid corrected carbon isotope values of tuna and model-derived Suess corrected phytoplankton carbon isotope values ($\Delta^{13}C_{tuna-phyto}$); (C) C:N ratios; and (D) standard ellipse area (SEA) as an indicator of isotopic niche area, for each Longhurst province (Table S1) in this study. Biome: C – Coastal, T – Trade wind; W – Westerly. For boxplots (A-C), the dotted line corresponds to the global median value, the black line is the region median, the shaded box includes the interquartile range (IQR), whiskers represent 1.5 IQR, and individual points are outliers outside of the IQR. For (B) the dark and light grey band highlight $\Delta^{13}C_{tuna-phyto}$ values that fall between -4 to -6‰ and -2 to -8‰, respectively (that define limited and moderate degrees of tuna movement; see Table 2). SEA density plots (D) represent the credibility intervals of the Bayesian SEA with shaded boxes encompassing the 50th, 75th, and 95th credibility intervals for each Longhurst province. Black circles and red x-marks represent the mode SEA and SEA corrected for small sample sizes ($SEA_c$).

**Figure 5.** Boxplots of albacore tuna (A) lipid and Suess corrected carbon isotopes values of tuna ($\delta^{13}C_{LScorr}$); (B) difference between lipid corrected carbon isotope values of tuna and model-derived Suess corrected phytoplankton carbon isotope values ($\Delta^{13}C_{tuna-phyto}$); (C) C:N ratios; and (D) standard ellipse area (SEA) as an indicator of isotopic niche area, for each Longhurst province (Table S1) in this study. Biome: C – Coastal, T – Trade wind; W – Westerly. For boxplots (A-C), the dotted line corresponds to the global median value, the black line is the region median, the shaded box includes the interquartile range (IQR), whiskers represent 1.5 IQR, and individual points are outliers outside of the IQR. For (B) the dark and light grey band highlight $\Delta^{13}C_{tuna-phyto}$ values that fall between -4 to -
6‰ and -2 to -8‰, respectively (that define limited and moderate degrees of tuna movement; see Table 2). SEA density plots (D) represent the credibility intervals of the Bayesian SEA with shaded boxes encompassing the 50th, 75th, and 95th credibility intervals for each Longhurst province. Black circles and red x-marks represent the mode SEA and SEA corrected for small sample sizes (SEA_c).

**Figure 6.** Relationships between predicted lipid- and Suess-corrected (δ¹³Cₐₗₐ) values of tunas and observations of standalone variables: year, latitude, oxygen at depth, and SST for (A) yellowfin, (B) bigeye, and (C) albacore tunas. The dashed horizontal red lines represent the intercept (or zero line) in each plot. The solid blue areas bracketing the response curves show the confidence limits of the model and are twice the standard error.

**Figure 7.** Relationships between difference between lipid corrected carbon isotope values of tuna and model-derived Suess corrected phytoplankton carbon isotope values (Δ¹³Cₐₗₐ) and observations of standalone variables: year, latitude, oxygen at depth, and SST for (A) yellowfin, (B) bigeye and (C) albacore tunas. The dashed horizontal red lines represent the intercept (or zero line) in each plot. The solid blue areas bracketing the response curves show the confidence limits of the model and are twice the standard error.

**Figure 8.** GAMM produced isoscapes of the difference between lipid corrected carbon isotope values of tuna and model-derived Suess corrected phytoplankton carbon isotope values (Δ¹³Cₐₗₐ) in (A) yellowfin, (B) bigeye, and (C) albacore tunas sampled between 2000 – 2007 and 2008 – 2015.

**Figure 9.** GAMM produced isoscapes of C:N ratios representing condition of (A) bigeye, (B) yellowfin, and (C) albacore tunas sampled between 2000 – 2015.

**Figure S1.** (A) Tuna muscle δ¹³C and (B) baseline-corrected δ¹³C (Δ¹³Cₐₗₐ) values (‰) in relation to bulk C:N values pre (green squares) and post (blue circles) mathematical correction for lipid content.

**Figure S2.** Model-data comparison of δ¹³C of particulate organic carbon (POC) in the surface layer (0-50 meters) of the ocean: (A) Observations of δ¹³C-POC in the data compilation from (Goericke and Fry, 1994), which were collected during 1970-1990; (B) Model δ¹³C-POC (Schmittner and Somes, 2016) averaged over the same time period (1970-1990) in a hindcast simulation forced by observed increasing atmospheric CO₂ and decreasing δ¹³CO₂ (i.e. Suess effect) as described in Schmittner et al. (2013); and (C) δ¹³C-POC versus latitude of observations in the Atlantic (red squares), Pacific (green circles), and Indian (blue diamonds), along with zonally-averaged model simulations from 1970-1990 when the observations were sampled (solid line) as well as from 2000-2009 (dotted line), which shows lower values due to the Suess effect.

**Figure S3.** GAMM produced isoscapes of lipid and Suess corrected δ¹³C values (δ¹³Cₐₗₐ) in (A) yellowfin, (B) bigeye, and (C) albacore tunas sampled between 2000–2007 and 2008–2015.

**Figure S4.** Tuna muscle (lipid and Suess-corrected) and modelled phytoplankton (Suess corrected) (A) δ¹³C values and (B) δ¹³C separation (Δ¹³Cₐₗₐ) in relation to distance from the equator.
<table><caption>**Table 1.** Global observations of tuna number (n) per time period, length, δ¹³C values of tuna corrected for lipid content and the Suess effect (δ¹³C<sub>corr_tuna</sub>), and differences between δ¹³C<sub>Lcorr_tuna</sub> and δ¹³C<sub>Scorr_phyto</sub> values (Δ¹³C<sub>tuna-phyto</sub>), C:N ratio and lipid content. Mean ± standard deviations are indicated. Sample sizes in parentheses refer to C:N ratio and lipid content data. A breakdown of sample sizes by species and region, and associated references can be found in Table S1.</caption><thead><tr><th></th><th>Unit</th><th>Yellowfin</th><th>Bigeye</th><th>Albacore</th></tr></thead><tbody><tr><td>2000-2007, n</td><td>-</td><td>1118 (790)</td><td>561 (455)</td><td>343 (247)</td></tr><tr><td>2008-2015, n</td><td>-</td><td>503 (32)</td><td>378 (30)</td><td>1372 (651)</td></tr><tr><td>Length</td><td>cm</td><td>94.3±37.9</td><td>84.8±34.2</td><td>82.4±16.8</td></tr><tr><td>δ¹³C<sub>LScorr_tuna</sub></td><td>‰</td><td>-15.7±0.9</td><td>-15.9±1.0</td><td>-17.2±1.0</td></tr><tr><td>δ¹³C<sub>LScorr_tuna</sub> min/max</td><td>‰</td><td>-19.3/-12.9</td><td>-18.4/-13.2</td><td>-20.0/-13.2</td></tr><tr><td>Δ¹³C<sub>tuna-phyto</sub></td><td>‰</td><td>-5.3±1.1</td><td>-5.0±1.1</td><td>-5.2±1.2</td></tr><tr><td>C:N ratio</td><td>-</td><td>3.6±0.3</td><td>3.5±0.3</td><td>4.1±1.1</td></tr><tr><td>Lipid content</td><td>%</td><td>2.8±1.0</td><td>2.6±1.2</td><td>4.8±4.3</td></tr></tbody></table>
<table><caption>**Table 2.** Summary of variables used in the global tuna analysis, based on δ¹³C values.</caption><thead><tr><td>Variable</td><td>Definition</td><td>Application</td></tr></thead><tbody><tr><td>δ¹³C<sub>Lcorr_tuna</sub></td><td>Tuna muscle δ¹³C values corrected for lipids by chemical extraction and/or mathematical correction or deemed to not require correction (i.e., C:N&lt;3.5)</td><td>Lipid corrected values to be normalized with baseline phytoplankton δ¹³C values for global tuna comparisons</td></tr><tr><td>δ¹³C<sub>Scorr phyto</sub></td><td>Model-derived baseline phytoplankton δ¹³C values corrected for the Suess effect</td><td>Baseline correction for global tuna comparisons</td></tr><tr><td>Δ¹³C<sub>tuna-phyto</sub></td><td>δ¹³C<sub>Scorr phyto</sub> - δ¹³C<sub>Lcorr tuna</sub></td><td>Movement and residency classifications and isotopic niche calculations</td></tr><tr><td>δ¹³C<sub>LScorr</sub></td><td>Tuna muscle δ¹³C values corrected for lipids and the Suess effect</td><td>Global tuna isoscapes</td></tr></tbody></table>
<table><caption>**Table 3.** Classification of spatial movement (or occupancy) and population dispersal in yellowfin, bigeye, and albacore tuna sampled in the different Longhurst provinces (Table S1) based on trends in Δ¹³C_tuna-phyto values presented in this study (Figures 2-3 B). Spatial movement is scaled from 0 to 4 with 0: limited movement (median and range falls within -4 to -6‰); 1: moderate movements by some individuals (median within -4 to -6‰, with range within -2 to -4 or -6 to -8‰); 2: moderate movement (median within either -2 to -4 or -6 to -8‰ with range between -2 and -8‰); 3: high degree of movement by some individuals (median within either -2 to -4 or -6 to -8‰, but range <-2 or >-8‰); and 4: high degree of movement (median is outside <-2 or >-8‰). Population dispersal was scaled from 0-4 with 0: no dispersal (range < 1‰); 1: low dispersal (range 1-2‰); 2: moderate dispersal (range 2-3‰); 3: high dispersal (range 3-4‰); and 4 very high dispersal (range >4‰). Blank spaces indicate areas where sample size was ≤ 5. Most limited (0) and highest (4) movement scores are highlighted in red and bold, respectively. Range corresponds to 1.5 of the interquartile range represented by the boxplot whiskers.</caption><thead><tr><th rowspan="2">Ocean</th><th rowspan="2">Biome</th><th rowspan="2">Longhurst Province</th><th colspan="3">Spatial movement</th><th colspan="3">Population dispersal</th></tr><tr><th>Yellowfin</th><th>Bigeye</th><th>Albacore</th><th>Yellowfin</th><th>Bigeye</th><th>Albacore</th></tr></thead><tbody><tr><td>Atlantic</td><td>W</td><td>MEDI</td><td></td><td></td><td>1</td><td></td><td></td><td>1</td></tr><tr><td>Atlantic</td><td>C</td><td>NECS</td><td></td><td></td><td>2</td><td></td><td></td><td>1</td></tr><tr><td>Atlantic</td><td>C</td><td>CNRY</td><td></td><td></td><td>1</td><td></td><td></td><td>1</td></tr><tr><td>Atlantic</td><td>C</td><td>GUIN</td><td>4</td><td>3</td><td></td><td>4</td><td>1</td><td></td></tr><tr><td>Atlantic</td><td>TW</td><td>ETRA</td><td>1</td><td>1</td><td></td><td>2</td><td>2</td><td></td></tr><tr><td>Atlantic</td><td>TW</td><td>WTRA</td><td></td><td>0</td><td></td><td></td><td>1</td><td></td></tr><tr><td>Atlantic</td><td>TW</td><td>NATR</td><td></td><td>0</td><td></td><td></td><td>1</td><td></td></tr><tr><td>Atlantic</td><td>W</td><td>NADR</td><td></td><td></td><td>2</td><td></td><td></td><td></td></tr><tr><td>Atlantic</td><td>W</td><td>GFST</td><td>2</td><td>2</td><td>1</td><td>1</td><td>2</td><td>0</td></tr><tr><td>Atlantic</td><td>W</td><td>NASW</td><td>0</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>Pacific</td><td>TW</td><td>SPSG</td><td>1</td><td>1</td><td>1</td><td>3</td><td>3</td><td>4</td></tr><tr><td>Pacific</td><td>TW</td><td>PEQD</td><td>3</td><td>3</td><td></td><td>4</td><td>4</td><td></td></tr><tr><td>Pacific</td><td>TW</td><td>PNEC</td><td>3</td><td>2</td><td></td><td>2</td><td>1</td><td></td></tr><tr><td>Pacific</td><td>C</td><td>CCAL</td><td>2</td><td></td><td></td><td>0</td><td></td><td></td></tr><tr><td>Pacific</td><td>TW</td><td>NPTG</td><td>2</td><td>1</td><td></td><td>2</td><td>3</td><td></td></tr><tr><td>Pacific</td><td>TW</td><td>WARM</td><td>1</td><td>1</td><td>1</td><td>4</td><td>4</td><td>1</td></tr><tr><td>Pacific</td><td>TW</td><td>ARCH</td><td>1</td><td>1</td><td>1</td><td>3</td><td>3</td><td>4</td></tr><tr><td>Pacific</td><td>C</td><td>SUND</td><td>1</td><td></td><td></td><td>1</td><td></td><td></td></tr><tr><td>Pacific</td><td>W</td><td>TASM</td><td></td><td></td><td>2</td><td></td><td></td><td>1</td></tr><tr><td>Pacific</td><td>C</td><td>AUSE</td><td>1</td><td>1</td><td>0</td><td>3</td><td>2</td><td>2</td></tr><tr><td>Indian</td><td>TW</td><td>MONS</td><td>1</td><td>1</td><td>1</td><td>2</td><td>1</td><td>0</td></tr><tr><td>Indian</td><td>W</td><td>ARAB</td><td>1</td><td>0</td><td></td><td>2</td><td>0</td><td></td></tr><tr><td>Indian</td><td>TW</td><td>ISSG</td><td></td><td></td><td>1</td><td></td><td></td><td>1</td></tr><tr><td>Indian</td><td>C</td><td>EAFR</td><td>0</td><td>0</td><td>0</td><td>1</td><td>1</td><td>0</td></tr><tr><td>Antarctic</td><td>W</td><td>SSTC</td><td></td><td></td><td>2</td><td></td><td></td><td>1</td></tr><tr><td>Antarctic</td><td>W</td><td>SANT</td><td></td><td></td><td>2</td><td></td><td></td><td>2</td></tr></tbody></table>
<figure><img src="image_2.png" /></figure>

$\Delta^{13}C_{\text{Tuna-Phyto}} = -8\permil$

$\Delta^{13}C_{\text{Tuna-Phyto}} = -2\permil$

<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>
<figure><img src="image_15.png" /></figure>
<figure><img src="image_16.png" /></figure>
**Conflict of Interest**

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.