Isotopic Tracers Suggest Limited Trans-Oceanic Movements and Regional Residency in North Pacific Blue Sharks (Prionace glauca)

Blue sharks (Prionace glauca) are globally distributed, large-bodied pelagic sharks that make extensive migrations throughout their range. In the North Pacific, mark-recapture studies have shown trans-Pacific migrations, but knowledge gaps in migration frequency hinder understanding of regional connectivity and assessments of regional demography for stock assessments. Here, we use oceanographic gradients of stable isotope ratios (i.e., regional isoscapes) to determine exchange rates of blue sharks between the East and West North Pacific Ocean (EPO and WPO). We generated regional δ13C and δ15N distributions for blue sharks from published values in the North Pacific (n = 180; both sexes, juveniles and adults combined). Discriminant analysis suggested low trans-Pacific exchange, categorizing all western (100%) and most eastern (95.3%) blue sharks as resident to their sampling region, with isotopic niche overlap of WPO and EPO highly distinct (0.01–5.6% overlap). Limited trans-Pacific movements suggest that other mechanisms maintain genetic mixing of the North Pacific blue shark population. Potential finer scale movement structure was indicated by isotopic differences in sub-regions of the eastern and western Pacific, though application of mixing models are currently limited by aberrantly low blue shark δ13C values across studies. Our results suggest that blue shark population dynamics may be effectively assessed on a regional basis (i.e., WPO and EPO). We recommend further studies to provide size- and sex-specific movement patterns based on empirical isotopic values with large sample sizes from targeted regions. Strategically applied stable isotope approaches can continue to elucidate migration dynamics of mobile marine predators, complementing traditional approaches to fisheries biology and ecology.


INTRODUCTION
Blue sharks (Prionace glauca) are large-bodied, highly migratory sharks with a global distribution extending throughout temperate and subtropical waters (Nakano and Stevens, 2008;Coelho et al., 2018). Blue shark populations have declined broadly, with high longline bycatch and mortality rates due to extensive overlap with commercial fisheries across much of their global range (Queiroz et al., 2016(Queiroz et al., , 2019. This is true for the North Atlantic (50-79% decline over 30 years; ICCAT, 2015) and .8% decline since the early 19th century; Ferretti et al., 2008) populations as a result of both targeted fisheries (i.e., for fins, meat, squalene) and bycatch (Clarke et al., 2006a,b;Cardeñosa et al., 2020), though the North Pacific population has recently been assessed as not overfished (ISC, 2017). While studies have challenged model-based inferences of shark population declines (Burgess et al., 2005), blue sharks are the major bycatch species in high-seas fisheries regionally (McKinnell and Seki, 1998;Francis et al., 2001) and perhaps globally (Clarke et al., 2006b;Campana et al., 2009). Limited genetic structure has been observed across populations sampled from disparate oceanic regions Taguchi et al., 2015;Veríssimo et al., 2017;Bailleul et al., 2018), though regional populations are managed separately. Mechanisms for high genetic homogeneity across global blue shark populations remain mostly speculative, mainly because studies are lacking in adequate sample sizes and sufficient demographic coverage to allow for robust conclusions to be made (Veríssimo et al., 2017). This limitation spans not only genetic information, but also a robust understanding of other aspects of the species' life history, particularly movements and migrations in relation to proposed mating and parturition grounds.
In the North Pacific Ocean, data suggest that mature blue sharks migrate to a latitudinal band spanning ∼20-30 • N for mating during the early summer months, with pupping typically occurring the following summer after a ∼12-month gestation. Gravid females are believed to move further north to parturition grounds located in sub-arctic waters between 35 and 45 • N (Nakano, 1994). Males and females typically segregate spatially prior to mating events, but some overlap between immature and mature individuals of opposite sexes can occur (Maxwell et al., 2019). Broadscale movements in the Pacific have been described with conventional and electronic tagging studies (Musyl et al., 2011;Maxwell et al., 2019), and although some blue sharks exhibit long-distance migrations over thousands of kilometers (Maxwell et al., 2019), most tracking data shows predominately latitudinal movements. However, recent conventional markrecapture information has shown that some individuals migrate across the northern Pacific Ocean from the western Pacific Ocean (WPO) into the eastern Pacific Ocean (EPO) and vice versa (Sippel et al., 2011), a behavior that also has been observed in the North Atlantic (Howey et al., 2017). The extent of connectivity between North Pacific sub-populations (e.g., WPO and EPO), however, remains unquantified. For example, it remains unclear whether genetic mixing and homogeneity  is maintained by trans-Pacific migrations of juveniles or adults, or by some other mechanism such as pupping and recruitment dynamics. Thus, the proportion of blue sharks performing trans-Pacific migrations warrants investigation, as quantifying movement connectivity can clarify regional sourcesink dynamics, inform spatial scales of management, and help explain the mechanisms that facilitate genetic homogeneity.
While electronic and conventional tagging approaches have provided useful information on blue shark movement and migration dynamics, studies are limited by high cost (for electronic tagging) and often require high sample sizes and protracted study duration to yield necessary ecological information . Consequently, complementary approaches are required for rapid assessment of blue shark migration and movement connectivity. Intrinsic chemical tracers measured in animal tissues, such as stable isotope (SI) ratios, are useful for reconstructing prior animal migrations (Graham et al., 2010;Madigan et al., 2021). Stable isotope analysis (SIA)-based movement studies utilize the distinct isotopic composition of prey baselines (i.e., regional isoscapes) across oceanic sub-regions, driven by local oceanographic and biogeochemical regimes (McMahon et al., 2013;Brault et al., 2018;Espinasse et al., 2020). In the North Pacific Ocean, pelagic prey fields in the EPO and WPO are isotopically distinct, particularly for nitrogen isotope ratios (δ 15 N) (Matsubayashi et al., 2020). In the EPO, upwelling of nutrient-rich water in the California Current promotes larger nitrate metabolizing primary producers such as diatoms, which creates a 15 N-enriched isotopic composition of regional prey (Altabet et al., 1999;Montoya, 2007;Madigan et al., 2012aMadigan et al., , 2017. Comparatively, WPO waters are nutrient-poor, which promotes dominance of nitrogen-fixing picophytoplankton at the base of the food web and lower (i.e., 15 N depleted) δ 15 N composition of regional prey pools (Takai et al., 2007;Fujinami et al., 2018;Ohshimo et al., 2019). Knowledge of regional isotopic baseline variation can be combined with measured predator stable isotope ratios and tissue-specific incorporation rates to identify individuals that have recently migrated from one system into another (Madigan et al., 2021;Shipley et al., 2021). Tissues of recent migrants will reflect an isotopic mix of prey baselines from the prior and current regions, provided that sampling has occurred prior to the consumer reaching isotopic steady-state with prey from their current region (Heady and Moore, 2013;Moore et al., 2016;Madigan et al., 2021). Residents are then defined as those individuals that are at isotopic steady-state with prey baselines in their current region (Madigan et al., 2014). This technique has been applied to characterize trans-Pacific migrations of Pacific bluefin tuna (Thunnus orientalis), using machine learning algorithms to define migrants vs. residents, and thus predicting the extent of mixing between WPO and EPO populations (Madigan et al., 2014). More recently, SIA approaches have been used to determine movement transitions across systems spanning marine, brackish, and freshwater habitats (Moore et al., 2016;Shipley et al., 2021), and have proven to be a robust and insightful approach for clarifying aspects of animal migration.
In the current study, we use regional δ 13 C and δ 15 N values from sampled sharks and prey to quantify potential trans-Pacific exchange rates of blue sharks between the WPO and EPO. We use multiple analytical methods with isotopic data to examine the extent of regional residency and foraging connectivity between oceanic sub-regions. This approach provides a tracerbased assessment of habitat use that is complementary to traditional tagging approaches and provides a framework than can be adopted across other study taxa and ecosystems. Inferred movements from isotopic signatures aid in constraining the extent of blue shark movements in the North Pacific basin, clarify the migratory mechanisms that may drive a mixed genetic stock, and can inform appropriate multi-national or regional management strategies for North Pacific blue sharks.

Data Compilation
Blue Shark SI Data A literature search was performed to obtain all SI data to-date for blue sharks in the North Pacific. Based on conventional tagging data and available SI data, we categorized North Pacific data into the WPO and EPO. From all studies, mean, reported error (SD or SE), and minimum and maximum values (when reported) of blue shark δ 13 C and δ 15 N values were tabulated. The literature search showed studies in four discrete regions of the EPO, so EPO blue sharks were further categorized into these four EPO sub-regions: Northern California Current (NCC), Southern California Bight (SCB), southern Baja (SBaja), and the eastern Tropical Pacific (ETP) (see section "Results"). We only used studies that accounted for lipid and urea effects on δ 13 C and δ 15 N through either chemical extraction or arithmetic correction and DI rinsing. While we recognize that different treatments can affect isotopic values, lipid content in most shark species is low (Hussey et al., 2012) and correction for urea in available data was not feasible. As a result, values were used as reported in published studies. Blue shark δ 13 C and δ 15 N values were estimated to represent their past foraging behavior for ∼0.5-1.5 years before sampling, based on published turnover rate estimates (Madigan et al., 2012b) and blue shark body size ranges across published studies (Thomas and Crowther, 2015;Vander Zanden et al., 2015).

Estimating Population-Wide Blue Shark SI Values
We used an iterative bootstrapping approach to resample each published δ 13 C and 15 N data distribution (n = 5 studies) to generate estimates of population-wide blue shark SI values for each ocean region. Blue shark δ 13 C and 15 N estimates were bootstrapped (1000×) by randomly sampling from mean (±SD) values published for each study region. This resulted in 1 × 10 3 estimates for the WPO (n = 1 study) and 4 × 10 3 estimates for the EPO (n = 4 studies), with separate, region-specific resampled populations for sub-regions of the EPO.

Prey SI Data
We searched the literature for published studies reporting SI values of epi-and mesopelagic prey in the WPO and EPO, based on known regional blue shark diets of epi-and mesopelagic forage fish, cephalopods, and crustaceans (Preti et al., 2012;Fujinami et al., 2018). Across relevant studies, prey-specific δ 13 C and δ 15 N means (±SD) were tabulated for subsequent analyses. As with blue sharks above, prey values were also obtained and regional means (±SD) calculated for the four defined sub-regions of the EPO (NCC, SCB, SBaja, and ETP). As with blue shark data, we only used studies that accounted for lipid effects on δ 13 C through chemical extraction or arithmetic correction, and values were used as they were reported in published studies.
Prey δ 13 C and δ 15 N values were used to generate regional mean diet δ 13 C and δ 15 N values in the EPO and WPO for subsequent analyses (see sections "Discriminant Analysis" and Isotopic Mixing Models"). For prey data in each region, we accounted for associated error of each prey mean δ 13 C and δ 15 N values by bootstrapping 1000 values from reported means (±SD) for each prey item. We then randomly selected from these prey distributions (1000×) to generate a mean blue shark diet value for the WPO (n = 145 prey species) and EPO (n = 75 prey species). Each prey item was weighed equally based on demonstrably broad and opportunistic blue shark diets. This resulted in 1 × 10 3 estimated mean (±SD) diet δ 15 N values for both the WPO and EPO. We also used the above approach to generate mean δ 13 C and δ 15 N diet estimates in the four EPO sub-regions. Combining data from studies that span different time periods accepts potential temporal variation in isotopic baselines and consumers, but given the observed distinction between isotope values between the WPO and EPO (see section "Results"), we deemed it unlikely that this would significantly impact overall results.

Data Analyses
Three analytical approaches were applied to blue shark and prey SI data to characterize blue shark movements (explained in detail below). Isotopic niche overlap of blue shark values was used to obtain general, quantitative metrics of likely blue shark exchange between sub-regions. We then used discriminant analysis to explicitly categorize individual shark isotope values as indicative of prior use of WPO or EPO waters. Finally, mixing models were used to estimate blue shark use of region-specific prey based on the isotopic composition of sharks and prey in each region.

Regional Variability and Overlap of Blue Shark SI Values
For the WPO and four EPO sub-regions (NCC, SCB, SBaja, ETP), integrated δ 13 C and δ 15 N niche areas (Bayesian ellipses: SEA B ) were determined, generating ellipses for each sub-region that incorporate 40% of the available data (Jackson et al., 2011). Isotopic overlap between each sub-region was then inferred using a Bayesian approach implemented in the R package "nicheROVER." Overlap estimates were generated from 1000 posterior draws based on 95% probabilistic niche regions (Swanson et al., 2015).

Diet-Dependent Diet-Tissue Discrimination Factors (DTDFs)
Discriminant analysis and isotopic mixing models required application of blue shark δ 15 N and δ 13 C diet-tissue discrimination factors (DTDFs). Because trophic discrimination factors can vary due to a suite of environmental and physiological processes (Hussey et al., 2012;Shipley and Matich, 2020) and have been shown to co-vary with diet δ 15 N and δ 13 C values (Caut et al., 2009;Hussey et al., 2014), we calculated diet-specific DTDFs for blue sharks within each region. For each calculated prey mean, we calculated mean DTDFs from algorithms reported in two studies: (1) Caut et al. (2009): (1) and (2) Hussey et al. (2014): where 15 N and 13 C represent diet-derived DTDFs for δ 15 N and δ 13 C, and δ 15 N diet and δ 13 C diet represent mean diet δ 15 N and δ 13 C values, respectively. For each regional diet δ 15 N mean, we calculated two DTDF values following Eqs 1, 2 and used the mean of both DTDF values (i.e., estimated DTDF values from equations) for applied diet-based blue shark 15 N and 13 C values for each sub-region.

Discriminant Analysis
Only δ 15 N values were used in discriminant analysis due to non-differentiation of WPO and EPO δ 13 C values (see sections "Results" and "Discussion"), following methods in Madigan et al. (2014). We generated training data for discriminant analysis using regional prey means and regional prey-based DTDFs. Specifically, regional DTDFs were added to regional prey means to generate 1 × 10 3 estimated blue shark δ 15 N values for the WPO, EPO, and EPO sub-regions. We then applied these training data to discriminant analysis of WPO and EPO blue shark δ 15 N data to classify individual sharks as recent migrants or long-term (≥1 year, based on ectotherm isotopic turnover rates; Thomas and Crowther, 2015;Vander Zanden et al., 2015) residents to the region in which they were sampled. Discriminant analysis reported an error value for the classification of unknown data, which estimates the percentage of individuals that were classified incorrectly (Klecka, 1980).

Isotopic Mixing Models
We applied isotopic mixing models to assess their efficacy for describing foraging across sub-regions in both the EPO and WPO. Bayesian isotope mixing models (R package "simmr"; Parnell, 2020) were used to estimate regional prey inputs to EPO and WPO blue sharks, providing estimates of regional connectivity. For the WPO, mixing models were run for a single blue shark population (due to n = 1 study in the WPO) and regional prey were based on reported prey values for four WPO sub-regions (Eastern Japan, Kuroshio-Oyashio, Sea of Japan, and offshore Taiwan) in Madigan et al. (2015), with additional prey data from Ohshimo et al. (2019). For the EPO, mixing models were run for regional blue shark populations and the pooled EPO population; regional prey endmembers were the NCC, SCB, SBaja, ETP, and the WPO. Endmember values and DTDFs were based on compilations of regional prey fields ( Table 1).
We estimated the accuracy of regional endmembers (n = 4 for WPO, n = 5 for EPO) by simulating 10,000 prey mixing polygons (Smith et al., 2013) and quantifying the probability of each individual consumer falling outside of the 95% prey mixing space. A relatively high proportion of blue sharks had >95% probability of falling outside the mixing space due to low δ 13 C values (see section "Results"), as has been observed in other studies (Rabehagasoa et al., 2012;Li et al., 2014;Kiszka et al., 2015). We consequently evaluated the two-isotope mixing model results to assess the effects of potential bias toward low δ 13 C diet inputs. We also performed a single isotope (δ 15 N TABLE 1 | δ 13 C and δ 15 N values (mean ± 1 SD) of published prey items by North Pacific Ocean sub-region, and calculated diet-dependent diet-tissue discrimination factors (DTDFs; 13 C and 15 N) estimated for blue sharks (Prionace glauca).

Region in North Pacific
Sub-regions a n b Years of collection Frontiers in Marine Science | www.frontiersin.org only) mixing model for both the WPO and EPO to compare to the two-isotope model results. The probable contributions of regional prey endmembers to blue shark diet were inferred from 10,000 model iterations, with a burn-in period of 1000 and a thinning interval of 100. Model convergence was evaluated based on inspection of Gelman-Rubin diagnostics, where values for each parameter should equal ∼1.0 (Phillips et al., 2014;Parnell, 2020).

Data Compilation
Blue Shark SI Data We obtained δ 13 C and δ 15 N values for blue shark muscle from one WPO study (Fujinami et al., 2018) (total n = 120 individuals) and four EPO studies (Miller et al., 2010;Madigan et al., 2012a;Li et al., 2014;Hernández-Aguilar et al., 2015) (total n = 60) from different regions of the EPO (NCC, SCB, SBaja, and ETP) (Figure 1). All studies included male and female sharks across juvenile and adult size ranges, though reporting of size and sex metadata did not allow for this information to be associated with individual δ 13 C and δ 15 N values. We assumed that published blue shark δ 13 C and δ 15 N values were generally normally distributed, based on reporting of mean ± SD or SE. Blue shark δ 13 C values showed high overlap between the WPO and EPO, with the majority of data in both regions falling between −20.0 and −17.0 (Figure 2); this was also observed in WPO and EPO prey that were used to calculate isoscape "baselines" for both regions (see Figure 2 and section below). In contrast, blue shark δ 15 N values were distinct between the WPO and EPO (Figure 2).
As with WPO and EPO blue shark data, prey δ 13 C values showed high overlap while prey δ 15 N values were discrete between the WPO and EPO (Table 1 and Figure 2). For this reason, we used only δ 15 N values for discriminant analysis, consistent with previous studies that used δ 15 N to quantify exchange rates of WPO and EPO Pacific bluefin tuna (Madigan et al., 2017;Tawa et al., 2017). Overall WPO prey δ 13 C  . Prey data (diamonds; species mean ± SD) were selected based on general prey groups (e.g., mesopelagic squids, forage fish) reported in published diet studies, and include epi-and mesopelagic forage fish, squids, and crustaceans. Lines represent linear fits to eastern (dashed line) and western (solid line) prey data.

Regional Variability and Overlap of Blue Shark SI Values
Across the five sub-regions, isotopic variability measured in terms of niche space (which was uniform across SEA and SEA B estimates) was greatest for individuals captured in the ETP (1.8 2 ), and lowest in the NCC (0.6 2 ). SCB, SBaja, and WPO individuals exhibited intermediate variability (1.1 2 -1.2 2 , Table 2 and Figure 4).
Isotopic niche overlap was extremely low between all possible combinations of EPO sub-regions and WPO blue sharks (<6%, Table 2 and Figure 4). Overlap between sub-regions of the EPO was substantial, but highly variable, with overlap estimates ranging from 12 to 94% (Table 2 and Figure 4). Blue sharks sampled in the NCC and SCB generally overlapped significantly with the isotopic niches of those from the ETP (>70%), but there was less overlap with SBaja (<48%). SBaja blue sharks overlapped minimally with NCC and SCB sharks (<39%), but overlapped highly with ETP sharks (94% , Table 2 and Figures 4, 5). Overlap of ETP sharks was relatively high with SCB and SBaja sharks (61 and 75%, respectively), but lower with the NCC (30% ,  Table 2 and Figures 4, 5). Overall, bilateral isotopic niche overlap (see Table 2) suggested that the highest level of WPO↔EPO connectivity was in the NCC (12%, 28%; Table 2), and the lowest level of WPO↔EPO connectivity in SBaja (<1%; Table 2 and Figure 5). Within the EPO, the highest δ 15 N-inferred connectivity of blue sharks was between the NCC and SCB (72%, 94%; Table 2) and SBaja and the ETP (94%, 75%; Table 2), and the lowest connectivity between SBaja and the NCC (12%, 28%; Table 2 and Figure 5).

Diet-Dependent DTDFs
Across all regions (i.e., EPO and WPO values, including all subregions), calculated blue shark 13 C values ranged from 0.8 to 1.2 and for 15 N from 1.7 to 3.9 . Due to substantial variation in WPO diet δ 15 N values, overall WPO DTDFs generated in bootstrapped estimates ranged from 1.8 to 4.8 (3.1 ± 0.5 ), and EPO DTDFs from 1.4 to 1.9 (1.5 ± 0.1 ) (see Table 1 for all regional DTDFs). Sub-region diet-dependent DTDFs varied substantially within both the EPO and WPO, based on differences in sub-region prey baselines (Table 1).

Discriminant Analysis
Discriminant analysis classified all WPO blue sharks as residents to the WPO (0% EPO migrants). In the EPO, 95.3% of blue sharks were categorized as residents to the EPO (∼5% WPO migrants), which varied by EPO sub-region: NCC (11% migrants), ETP (5% migrants), SCB (2% migrants), and SBaja (0% migrants), with <1% classification error across all regional discriminant analyses. The value representing the cutoff point in discriminant analyses (threshold δ 15 N value between WPO-and EPO-classified sharks) was similar across all analyses conducted, at ∼14.0 . Most δ 15 N values of WPO migrants in the EPO (δ 15 N < 14.0 ) were in the tails of bootstrapped population δ 15 N values, though one empirical EPO value (δ 15 N = 13.8 ; minimum reported in ETP) was classified as a WPO migrant.
Training data for discriminant analysis (estimated blue shark δ 15 N values from WPO and EPO prey, calculated from regionspecific diet estimates and DTDFs), were (WPO) 13.1 ± 1.1 and (EPO) 15.8 ± 1.1 . WPO and EPO training data were highly discrete between the two regions (13% overlap; Figure 3) and highly coherent with bootstrapped blue shark δ 15 N distributions (Figure 3).

Isotopic Mixing Models
Substantial blue shark SI data (9-32%) fell outside of the prey mixing space due to low shark δ 13 C values ( Supplementary  Figure 1 and Supplementary Table 2), leading to unreliable mixing model results. There were large discrepancies between dual (δ 13 C and δ 15 N) and single (δ 15 N) isotope mixing models in both the EPO and WPO (Supplementary Figures 2, 3 and Supplementary Tables 2, 3). In the EPO, the dual isotope model suggested high SCB inputs to all regions (66-100%), likely biased by the SCB having the lowest prey δ 13 C values  In the WPO, a similarly high proportion of blue shark SI data (55%) fell outside of prey mixing space due to low shark δ 13 C values ( Supplementary Figure 1 and Supplementary Table 2). The dual isotope model suggested blue sharks foraged most on prey from Kuroshio-Oyashio (65%), though model results were likely biased by that region's relatively low prey δ 13 C values (Supplementary Figure 3). When δ 13 C values were excluded from mixing models, the single-isotope δ 15 N model suggested greater contributions to blue shark diet from all sub-regions, with the southern region (offshore Taiwan) contributing the most (56%; Supplementary Figure 3 and Supplementary Table 3). The 95% credible intervals for all mixing models results are shown in Supplementary Tables 2, 3.

DISCUSSION
Region-specific stable isotope values of blue sharks and prey allowed for inferences of prior movement patterns and regional connectivity in the North Pacific Ocean, quantifying dynamics that have been observed (e.g., trans-Pacific migrations), but for which frequency and exchange rates are unknown. Our results suggest minimal trans-Pacific movements and indicate potential finer-scale movement and residency dynamics within sub-regions of the EPO and WPO, though key caveats were evident in estimates of finer scale movements. High coherence of prey-estimated blue shark δ 15 N values with empirical SI values demonstrated the effectiveness and predictive value of North Pacific isoscapes. Overall, our results demonstrate the efficacy of our multi-analytical stable isotope approach to identify movements of a highly migratory pelagic species in the North Pacific, and the potential for future analyses with consideration of sample treatment and quantification of crucial species-specific isotopic parameters.
Trans-Pacific migrations in North Pacific blue sharks were first demonstrated by conventional tagging efforts in the eastern Pacific (Sippel et al., 2011), and we used SI data to assess the frequency of EPO↔WPO movements. In the WPO, SI datasets revealed no evidence of migrations from the EPO, compared to conventional tagging studies that reported 4 of 205 (∼2%) EPOtagged (NCC, SCB, and SBaja) blue sharks migrating to the WPO (Sippel et al., 2011). In the EPO, discriminant analysis indicated trans-Pacific migration from the WPO in low proportions (∼5%) of our EPO bootstrapped population estimates, including at least one empirical migrant value (δ 15 N = 13.8 ; ETP). However, conventional tagging found no trans-Pacific migration to the EPO from the WPO (n = 207) (Sippel et al., 2011), and more recent satellite tagging in the EPO (n = 47) and WPO (n = 21) showed movements only within those respective ocean regions (Maxwell et al., 2019;Fujinami et al., 2021). While tagging studies have thus provided quantitative metrics of trans-Pacific migrations, tag-inferred movements have limitations, including simplistic movement information and potential non-reporting of tag recovery (conventional tags), limited sample size (electronic tags), and short-term tracks that do not capture long-distance movements (both tag types) (Siskey et al., 2019). Tagging data is also prospective, capturing future rather than prior movements that are potentially biased by tagging location. In contrast, isotopic measurements coupled with isotopic turnover rates and DTDFs can provide retrospective, quantifiable timeframes of prior movements. While SI-inferred movements are limited by isotopic turnover rates of analyzed tissue, prior movements can be characterized due to mobile predators integrating prey isotopic signatures during movements through isotopically distinct regions (Graham et al., 2010;MacKenzie et al., 2012;Carlisle et al., 2015;Madigan, 2015;Trueman and Glew, 2019). As such, isotope-and tag-inferred migration patterns are effective complementary techniques (Carlisle et al., 2012(Carlisle et al., , 2015Madigan et al., 2015Madigan et al., , 2018Shipley et al., 2021). Here, comparing isotopic estimates to past conventional tagging suggests that blue sharks make both eastward and westward trans-Pacific migrations, though the number of sharks that make these migrations appears to be low (<5%).
The low trans-Pacific exchange inferred here improve understanding of blue shark population dynamics in the North Pacific. Genetic analyses have shown global panmixia across regional blue shark populations, with minimal evidence of regional population structure Veríssimo et al., 2017;Bailleul et al., 2018). This lack of observed population structure requires some mechanism of regional population mixing, and while maintenance of population genetic homogeneity does not necessarily require high mixing (Bremer et al., 2005;Waples and Gaggiotti, 2006), the minimal trans-Pacific exchange we observed here is unlikely to be the primary driver. Current understanding of blue shark life history provides several alternative scenarios for mixing of EPO and WPO populations. If WPO and EPO sharks mate (central-southern waters; ∼20-30 • N) and pup (northern waters; ∼35-45 • N) (Nakano, 1994;Nakano and Stevens, 2008) in their respective ocean basins, EPO-and WPO-origin young-of-the-year (YOYs) could subsequently recruit to either the WPO or EPO. In addition, the Central Pacific Ocean (CPO) (i.e., waters around Hawaii) could serve as a mixing region for WPO and EPO sharks, as some exchange of blue sharks in this region has FIGURE 5 | Summarized estimates of blue shark (Prionace glauca) migratory exchange between eastern North Pacific Ocean regions, inferred from niche overlap of regional δ 13 C and δ 15 N values. Migratory exchange was inferred quantitatively by calculating overlap between 95% probabilistic isotopic niches of each sub-region, using a Bayesian approach ("nicheROVER"; Swanson et al., 2015). Arrows are scaled to the degree of isotopic niche overlap between regions; note that exchange rates are relative and approximate.
been observed previously (Sippel et al., 2011). Consequently, dispersive YOY recruitment and/or partial, temporary mixing of adults in the CPO (if for mating or parturition) could drive genetic mixing while maintaining regional WPO and EPO isotopic signatures that indicate regional residency of juveniles and adults. Alternatively, genetic analysis has demonstrated the possibility of "genetic time-lag" effects in blue sharks (Bailleul et al., 2018), with the possibility of discrete sub-populations despite genetic homogeneity; thus, WPO and EPO separation observed here could indicate discrete, minimally mixed subpopulations. As such, the degree to which YOY recruitment and adult mixing contribute to North Pacific blue shark population dynamics warrants further study.
Comparison of blue shark and prey values between the WPO and EPO demonstrated that δ 15 N, and not δ 13 C, serves as a regional diagnostic tracer in the North Pacific Ocean as has been previously observed in Pacific bluefin tuna (Madigan et al., 2017). We observed high overlap of WPO and EPO prey δ 13 C values, but almost no overlap of prey δ 15 N values (Figure 2). Mechanisms for this have been demonstrated, with photosynthetic pathways in pelagic primary producers (C 3 photosynthesis) varying more with latitude (i.e., due to variable productivity, temperature, and seawater pCO 2 regimes), rather than longitude (Bowen, 2010;Magozzi et al., 2017;Brault et al., 2018;Ohshimo et al., 2019). Coherence of prey-based estimated δ 15 N blue shark values with shark-derived values (Figure 3) further supports both the robustness of the δ 15 N isoscape approach and the marked separation between δ 15 N values of WPO and EPO sharks. Similar patterns of distinct EPO and WPO isotopic signatures have been found in Pacific bluefin tuna, a large-bodied pelagic teleost in the WPO and EPO that makes seasonal migrations through the same WPO and EPO regions (Boustany et al., 2010;Madigan et al., 2017;Tawa et al., 2017). In Pacific bluefin tuna, δ 15 Nbased estimates of trans-Pacific migration were demonstrably effective based on coherence with other chemical tracers, while non-differentiation of tuna δ 13 C made it an ineffective migration tracer (Madigan et al., 2017). This differs from more typical isotopic applications to ecology, which generally use δ 15 N to estimate trophic dynamics and δ 13 C to trace energy source and foraging location(s), due to lower trophic fractionation of δ 13 C than of δ 15 N (Post, 2002). However, more recent studies have demonstrated that in some systems and for certain predators, baseline δ 15 N values across ecoregions can result in predators acquiring regional δ 15 N signatures that outweigh trophic effects (Graham et al., 2010;Hobson et al., 2010;Madigan et al., 2017;Shipley et al., 2021). Isoscapes will be basin-and ecosystemdependent. In the North Pacific, the δ 15 N gradient is caused by upwelling-driven enrichment of 15 N in the EPO and low δ 15 N values due to oligotrophic N-fixation in the WPO (Takai et al., 2007;Madigan et al., 2012aMadigan et al., , 2017Fujinami et al., 2018;Ohshimo et al., 2019). The finer scale structure of δ 13 C and δ 15 N values in the WPO and EPO shown here demonstrate the utility of SIA in these ecosystems to track inter-and intra-basin predator movements.
We used niche overlap of blue shark δ 13 C and δ 15 N values to quantify finer scale connectivity dynamics within EPO and WPO sub-regions. In the EPO, isotopic niche overlap suggested variable mixing between the NCC, SCB, SBaja, and ETP (Figure 5). High exchange was suggested between the NCC and SCB; this is supported by conventional (Sippel et al., 2011) and electronic tag studies, which also found sexual segregation between these two regions (Maxwell et al., 2019). Lower overlap between SBaja and the NCC/SCB suggests that sharks may be more resident to this region. Regional proximity likely plays a role in exchange dynamics, as in general, more proximate regions showed greater migratory exchange (Figure 5). Collectively, estimated connectivity in the ETP and SBaja supports the premise that coastal and pelagic waters off SBaja may serve as both an overwintering ground for juveniles and a potential reproductive function for adults (Vögler et al., 2012), with isotopic overlap suggesting migration from SBaja to the ETP (Figure 5). Importantly, movement dynamics related to size and sex structure have been observed across these regions (Nakano and Stevens, 2008;Vögler et al., 2012;Maxwell et al., 2019), which we could not evaluate here due to a lack of size-and sex-specific information to match with blue shark isotopic values. Improved insight into size-and sex-specific movement dynamics could be accomplished with robust sampling of blue sharks for empirical isotopic measurements, coupled with size and sex metadata across the defined study regions. This could easily be achieved through sampling of fisheries [by]catch, potentially incorporating tissues of different turnover rates (i.e., fast [plasma, liver]; slow [muscle]; Thomas and Crowther, 2015;Vander Zanden et al., 2015).
While mixing models provided exploratory and potentially informative estimates of regional connectivity, results of finer scale sub-regional movements appeared biased and thus unreliable. Although the δ 15 N values of sharks and DTDFcorrected prey highly overlapped, shark δ 13 C values were low relative to prey in both the EPO and WPO (Supplementary  Figures 2, 3). Furthermore, prey mixing polygons (Smith et al., 2013) showed up to 32% and 55% of blue shark values falling outside of the simulated prey mixing space, due to these low shark δ 13 C values (Supplementary Figure 1 and Supplementary  Tables 2, 3). This likely biased mixing model estimates toward the lowest δ 13 C prey input(s); this confounding factor was observed in the two-isotope mixing models in both the WPO (high input of Kuroshio-Oyashio) and EPO (high input of SCB) (Supplementary Figures 2, 3). A similar effect of low blue shark δ 13 C values was observed in a previously published estimate of diet based on isotope mixing models in the SBaja sub-region, in which low blue shark δ 13 C values relative to prey resulted in the two lowest δ 13 C prey items (pelagic octopus Argonauta spp. and pelagic red crab P. planipes) dominating diet estimates (Hernández-Aguilar et al., 2015). Consequently, mixing model results here demonstrate that biased results will likely occur when predator SI values are substantially offset from prey inputs. While this is unsurprising and has previously been identified in the literature, studies continue to adopt this approach and report findings without identifying potential bias. It is imperative that investigators use appropriate diagnostic tools (i.e., quantitative assessment of prey vs. predator data following Smith et al., 2013) to ensure results are accurate or that potential bias is reported. For blue sharks, low δ 13 C values that are irreconcilable with local prey is consistent across studies and requires clarification.
There are three possible explanations for the consistently low δ 13 C values observed in blue sharks: (i) feeding in a region, within timeframes, or on prey not represented in regional studies; (ii) estimated 13 C values for blue sharks are too high, and fractionation between prey and blue shark muscle is lower than our estimated DTDF values; or (iii) shark sample preparation for SIA resulted in artificially low δ 13 C values across studies. Missing or isotopically misrepresented prey sources (scenario [i]) is a common issue across studies (Smith et al., 2013); here, sampling too close to coastlines could result in unrealistically high prey δ 13 C and/or δ 15 N values for a pelagic shark. However, most studies used in these analyses explicitly sampled in offshore epi-and mesopelagic zones. In addition, prey data are available for two offshore regions that were not included here, the TZCF and CPO, and prey δ 13 C values from those regions (−18.2 ± 0.8 and −17.7 ± 0.8 , respectively) (Gould et al., 1997;Choy et al., 2015) are also not low enough to explain the observed shark δ 13 C values. Temporal variability in isotopic values, while possibly influential, is also an unlikely explanation here, as studies spanned multiple seasons and years and different EPO and WPO values are driven by coarsely consistent oceanographic conditions (Madigan et al., 2017). It is possible that our DTDF estimates are imprecise (scenario [ii]), as they were calculated as diet-based DTDFs following Caut et al. (2009) andHussey et al. (2014), rather than empirically derived from laboratory experiments. Differential amino acid composition and subsequent 13 C fractionation in blue sharks could drive atypical DTDFs . Laboratoryderived 13 C DTDFs available for other elasmobranch species (Hussey et al., 2010;Logan and Lutcavage, 2010;Kim et al., 2012) are rarely <1.0 (the mean DTDF applied here), and a DTDF of ≤0 would be necessary for blue shark values to be highly coherent with prey. Negative 13 C values for sharks are rare in available studies, though one negative value (−0.5 ) has been estimated for one prey type (of 8 total diet items in natural diet) in the catshark Scyliorhinus canicula, so this DTDF in wild blue sharks cannot be ruled out (Caut et al., 2013). Similar observations in other ocean basins of low blue shark δ 13 C values relative to prey and/or other sharks (Rabehagasoa et al., 2012;Li et al., 2014;Kiszka et al., 2015) support the possibility that a unique aspect of blue shark muscle composition or physiology could drive atypically low 13 C, but this likely can only be validated with captive studies. Finally, it is possible that sample preparation partially contributed to low δ 13 C values in at least some studies (scenario [iii]). Studies included here performed lipid and/or urea extraction, but not always both (Madigan et al., 2012a;Li et al., 2014;Hernández-Aguilar et al., 2015;Fujinami et al., 2018), or mathematically corrected for lipids (Miller, 2006). Many of these studies preceded thorough published analyses demonstrating the importance of both lipid and urea extraction in elasmobranch tissues, which prevents artificially low δ 13 C and/or δ 15 N values (Carlisle et al., 2016;Li et al., 2016;Arostegui et al., 2019); recent work has also shown that water rinses for urea extraction may actually increase δ 13 C values (Carlisle et al., 2016). Our analyses, which draw upon published studies, underscore the need for standardized, consistent sample treatment in ongoing and future isotopic studies (Wolf et al., 2009;Shipley and Matich, 2020), while also noting the value of large archival datasets to address questions at an ocean basin scale. While laboratoryderived DTDFs will be difficult to obtain directly for blue sharks in captivity, improved understanding of 13 C dynamics in this species will be necessary to refine the accuracy of mixing model results, an issue that may also be applicable across other elasmobranch species.
While the data used here came from known high-use regions for blue sharks in the North Pacific, there are other relevant regions that could not be included in our analyses. In particular, the North Pacific Transition Zone (NPTZ) (also referred to as the Transition Zone Chlorophyll Front) has been shown as a region of high blue shark abundance (Pearcy, 1991;Polovina et al., 2001;Kubodera et al., 2007;Vögler et al., 2012). Recently, the NPTZ has been demonstrated as a migratory corridor for pregnant females that were satellite tagged in the WPO (Fujinami et al., 2021), has been suggested as a nursery ground for young sharks (Nakano and Stevens, 2008), and may be a migratory corridor for other blue shark life stages/sexes as well as tunas, swordfish, and turtles (Block et al., 2011). While blue shark isotope data were not available from the NPTZ, prey δ 13 C and δ 15 N appear to be similar to the WPO (Gould et al., 1997). As such, blue sharks with long-term residency in the NPTZ may be isotopically indistinguishable from WPO residents. It is currently unknown whether blue sharks in the NPTZ are residential to the frontal region for adequate timeframes (i.e., months to >1 yr) to acquire the NPTZ isotopic signal or whether they use the region temporarily during inter-region migrations; the WPO migrants we observed here in the NCC could reasonably be entering the California Current via the NPTZ. Isotopic sampling of NPTZ sharks, likely accessible from bycatch in offshore fisheries, could ascertain migration dynamics in NPTZ sharks. Similarly, blue sharks are common in the Central Pacific (i.e., waters around Hawaii; CPO), where conventional tagging suggests mixing with both the EPO and CPO (Sippel et al., 2011). With no isotopic characterizations of blue sharks in the CPO, we could not include this region here; however, isotopic analysis of CPO sharks could reveal the extent to which these sharks are residential to the region, as local prey seem to be distinctive from other regions, particularly prey δ 15 N values (Choy et al., 2015). Collecting samples across the spatial range of North Pacific fisheries is tractable, especially since a small biopsy sample from a subsequently released, live shark is sufficient for SIA.
This study was limited to the use of bootstrap-estimated blue shark δ 13 C and δ 15 N values rather than using direct empirical measurements. While our approach allowed for populationwide estimates based on empirically derived data distributions, some regional blue shark studies had relatively low sample sizes (e.g., n = 9 in SCB, n = 10 in NCC), and in this context, the "tails" of bootstrapped distributions likely lead to unrealistically high and low SI values (Madigan et al., 2017). Bootstrapped estimates can also lead to de-coupled shark δ 13 C and δ 15 N values (see shark data ellipses in Figures 2, 4 and Supplementary  Figure 2), when empirical δ 13 C and δ 15 N values are often positively correlated. This could affect mixing model estimates of prey/regional contributions to blue shark diet. The lack of size and sex data also limited our ability to undertake more detailed reconstruction of blue shark residency and movement dynamics across life history. Finally, isotopic identification of migrants is dictated by tissue turnover to steady-state conditions (here, likely ∼0.5 to 1.5 year before sampling), precluding identification of long-distance migrations that occurred prior to these timeframes. As such, results here should be taken in the context of these limitations and be viewed as a preliminary framework for regionally focused, empirical investigation.

CONCLUSION
Our results, drawing upon published δ 13 C and δ 15 N data for blue sharks and prey sampled at multiple locations in the EPO and WPO, provide a new and replicable means to assess blue shark residency and migration dynamics in the North Pacific. The analyzed data provide strong evidence for limited direct migrations between the WPO and EPO and reiterate the utility of δ 15 N isoscapes for the reconstruction of migratory predator movements in the North Pacific Ocean. Limited trans-Pacific migrations suggest that other mechanisms maintain genetic homogeneity of the North Pacific blue shark population, including YOY movements and/or partial mixing of adults in the Central Pacific. Regional structure in δ 13 C and δ 15 N data have promise for further quantification of finer-scale blue shark movements, increasing the resolutions of movement patterns suggested here, but consideration of isotopic parameters (e.g., accurate species-specific DTDFs), appropriate sample preparation of shark tissues, and length/sex metadata of sampled sharks are necessary. With emerging research showing varying residency and trans-regional movements in migratory predators, isoscapes can employ high sample sizes across a breadth of animal life stages, regions, and timeframes to reconstruct habitat use of highly mobile marine animals. Through these isotopic approaches, population-level estimates of movement dynamics are feasible on scales that may not be readily available from conventional tagging or telemetry studies.

DATA AVAILABILITY STATEMENT
The data analyzed in this study is subject to the following licenses/restrictions: Data came from published manuscripts, which are only publicly available through purchase of publications. Requests to access these datasets should be directed to daniel.madigan@stonybrook.edu.

AUTHOR CONTRIBUTIONS
DJM, ONS, NEH, and ABC conceived the study. DJM, ONS, and ABC performed the analyses. DJM and ONS wrote the manuscript, with input from all authors. All authors interpreted results and refined analyses.

FUNDING
This project was supported by WWF Canada and an NSERC Discovery Grant (#04922-2017) to NEH.

ACKNOWLEDGMENTS
Information and comments regarding conventional tagging data were generously provided by T. Sippel, and K. James provided valuable comments on the manuscript.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021. 653606/full#supplementary-material Supplementary Figure 1 | Comparisons of blue shark (Prionace glauca) δ 13 C and δ 15 N values to regional prey. Prey fields (filled gray forms) were generated by simulating 10,000 polygons using prey δ 13 C and δ 15 N values, following Smith et al. (2013). All prey values are adjusted by the addition of calculated diet-based DTDFs to allow quantification of overlap with blue shark δ 13 C and δ 15 N values (black circles). Proportion of blue shark values falling outside prey polygons, mostly due to low shark δ 13 C values, were 55% in the WPO (upper left panel) and 9-32% in regions of the EPO (other panels). These data highlight the importance of quantitatively assessing prey/predator isotope dynamics to ensure accurate interpretation of mixing models results and/or to determine (and report) the level of potential bias.
Supplementary Figure 2 | Isotopic overlap of regional blue shark (Prionace glauca) data with regional prey, and exploratory mixing model estimates of regional prey contributions, in sub-regions of the eastern Pacific Ocean. (A) Bootstrapped blue shark δ 13 C and δ 15 N values (small circles, colored by EPO sampling sub-region) and regional prey means (large circles; error bars ± SD), from the western (WPO) and eastern (EPO) Pacific Ocean. Mean prey δ 13 C and δ 15 N values are adjusted by the addition of calculated diet-dependent diet-tissue discrimination factors (DTDFs) (Caut et al., 2009;Hussey et al., 2014). After prey mean adjustment for DTDF, most blue shark δ 13 C values were left-shifted (lower δ 13 C) relative to prey δ 13 C values. (B) Estimated regional prey inputs to EPO blue shark diet from Bayesian mixing models. Left panel shows results from the dual isotope model (δ 13 C and δ 15 N), which were biased toward the regional prey with lowest δ 13 C values, and right panel the single isotope (δ 15 N) model.
Supplementary Figure 3 | Isotopic overlap of blue shark (Prionace glauca) data with regional prey, and exploratory mixing model estimates of regional prey contributions, in the western Pacific Ocean. (A) Bootstrapped blue shark δ 15 N values (small gray circles) and regional prey means (large circles, colored by WPO sub-region) from the western Pacific Ocean (WPO). Mean prey δ 13 C and δ 15 N values are adjusted by the addition of calculated diet-dependent diet-tissue discrimination factors (DTDFs) (Caut et al., 2009;Hussey et al., 2014). After prey mean adjustment for DTDF, most blue shark δ 13 C values were left-shifted (lower δ 13 C) from expected prey-based values. (B) Estimated regional prey inputs to WPO blue shark diet from Bayesian mixing models. Left panel shows results from two isotope model ( 13 C and δ 15 N), which were biased toward the regional prey with lowest δ 13 C values, and right panel shows a single isotope (δ 15 N) model. Table 1 | Isotopic niche metrics generated from carbon and nitrogen stable isotope values of blue sharks (Prionace glauca). Standard ellipse (SEA) areas ( 2 ) are derived from Northern California Current, Southern California Bight, Southern Baja, Eastern Tropical Pacific, and West Pacific Ocean (WPO). SEA estimates represent maximum likelihood (SEA) and Bayesian (SEA B ; [75% Cis]) derived estimates based on 40% of the data.

Supplementary
Supplementary Table 2 | Reliance of EPO blue shark populations on regional prey groups as inferred from Bayesian isotope mixing models. Results are median estimates (95% credible intervals [CIs]) derived from the posterior distributions of dual (δ 13 C and δ 15 N) and single (δ 15 N) isotope models. For dual isotope models, the percentage of individuals with >95% probability of falling outside of the simulated prey mixing space is shown.
Supplementary Table 3 | Reliance of WPO blue shark populations on regional prey groups as inferred from Bayesian isotope mixing models. Results are median estimates (95% credible intervals [CIs]) derived from the posterior distributions of dual (δ 13 C and δ 15 N) and single (δ 15 N) isotope models. For dual isotope models the percentage of individuals that had a >95% probability of falling outside of the simulated prey mixing space is indicated.