NWP area is referred to as FAO Major Fishing Area 61, which is compassed by the north of 20°N and the west of 175°W across the mainland coast of countries, including China, North Korea, South Korea, Japan, and Russia. Note that SKJ catch in the southwestern corner of FAO Area 61, i.e., west of 115°E and south of 20°N, is excluded in the study as it is near zero and the study area is limited to the south of 45°N of NWP area, considering the distribution of SKJ (Figure 1).

Monthly aggregated SKJ catch data using fishery gear types (pole and line, and purse seine) were provided from WCPFC (Western and Central Pacific Fisheries Commission, http://www.wcpfcnt/public-domain, last accessed in November 2020). The primary result of retrieved data for the NWP area showed that pole and line fishery was the main operating gear with its catch nearly six times to the purse seine fishery during 1972–2019. Finally, this study uses the pole and line data with a resolution of 5°, which comprises year, month, catch, effort (operating days), latitude, and longitude. CPUE is assumed to be a common proxy for fish relative abundance. In this study, the CPUE of SKJ is defined as SKJ catch in metric tons per operating day.

Five climate indices associated largely with the NWP and SKJ migration (Figure 1), and SST data of the NWP with a resolution of 5°are employed in this study (Table 1). All data were monthly time series with a period of 1972–2019.

Convergent cross mapping (CCM) is an approach that can distinguish causality from mirage correlation in time series from dynamical (i.e., nonlinear) systems, which is nonparametric and rooted on state space reconstruction (Sugihara et al., 2012). The fundamental CCM principle is that the cause time series left its footprints on the effect time series at a state space level. Before applying CCM analysis, CPUEda is created to represent the anomaly fluctuations of SKJ relative abundance by detrending CPUE first and subsequently monthly anomaly of the detrended CPUE. The relationship between CPUEda and climate indices were preliminarily analyzed using the Spearman rank correlation, which provides a measure of a monotonic relationship between two continuous random variables and is useful with non-normal data, the 95% significance level is set on the correlation analysis. Results showed no prominent relationship between climate indices and SKJ existed except for a weak and negative correlation between NPGO and SKJ (Table S1). Here, time series are transformed into a dynamic system of low-dimension and nonlinear to represent its space state before performing CCM by the following steps:

Generalized additive models (GAM) and threshold generalized additive models (TGAM) were used to determine the relationship type (stationary or non-stationary) between SKJ and the causal climate index (AMO, from the above results). Specifically, a “stationary” relationship is better fitted by a single function throughout the period of the time series and is typically formulated using a GAM (Ciannelli et al., 2004):

where Y is the response variable (the square root of CPUE that conforms to normal distribution), X₁ is the predictor (AMO), X₂ is the month variable (as a categorical variable), and s, α, and ε are smooth function (with k ≤ 3 to avoid overfitting), intercept, and error terms, respectively.

Different functions better fits a “non-stationary” relationship for different periods, and the responses to drivers have an abrupt change over a threshold year (Litzow et al., 2018), which is formulated using a TGAM (with specific to two time periods) (Puerta et al., 2019):

where y is the threshold year that separates two periods with varying responses to drivers, set as between the 0.1 lower and the 0.9 upper quantiles of the time series at a month scale, the designations of X₁, X₂, and Y are the same as them in the GAM. The nonstationary optimization model is selected by minimizing the model’s generalized cross-validation score (GCV) (Casini et al., 2009). The superior model was selected and further conducted based on the minimum Akaike information criterion (AIC). The analyses in this section were conducted using “mgcv” package in R.Analyses flow is shown in Figure 2.

The fluctuations of SKJ catch present earlier increase and later decrease in trend during 1972–2019: An increasing trend was determined from 1972 to 1984 with the highest catch (151,649-mt) occurring in 1984; a relative stable segment during 1985–2005 is characterized by a slight decline; at last, an obvious declining trajectory appeared since 2006, where the lowest catch (25,267-mt) occurred in 2014. (Figure 3A). The CPUE of SKJ fluctuates with a rising trend over the study period at the decadal scale. It remains a line before the 1980s, but increases sharply in the mid-1990s, reaching the first peak (5.61-mt/day) in 1993 and gently fluctuating during 1994–2014, then the maximum peak (5.84-mt/day) occurred in 2015. (Figure 3A).

The high-catch areas distribute along the western boundary of the NWP, which overlaps the path of the Kuroshio Current fitly. Among these, the biggest patchy area of high catch is near the mixed area between Kuroshio Current front and Oyashio Current (Figure 3B). By contrast, the spatial distribution of SKJ CPUE is characterized by many meridional strips, where two intensified trips of CPUE are located at the most northern area (along 35°–40° N) and the southeast corner (150°E across on the 20°N), respectively (Figure 3C).

An obvious convergence of the AMO causing on CPUEda and a declined trend of the opposite predicting skill (CPUEda causing on the AMO) are found with the library size larger (Figure 4), which shows the interaction between the AMO and SKJ is unidirectional, and indicates that the abnormal fluctuation of SKJ CPUE in the NWP depended on the long-term variation of the AMO, despite the insignificant correlation (r = −0.004, Table S1). The ECCM results further attest that the strongest causation occurs at a lag of 1 year at an annual scale. (Figure 5A). More specifically, the maximum lag on the month scale is identified at 15 months (Figure 5B).

For other climate indices, the predicting skill of the ONI causing on CPUEda shows a weak convergence at the library size (L = 400) with a declined trend after the library size reached 400 (Figure S6). Considering the maximum predicting skill is less than 0.1, the ONI effect on SKJ is almost neglected. The predicting skills between PDOfd and NPGOfd cross-with CPUEda presented with negative values make no sense (Figure S6), which shows no causal interaction between SKJ with PDO and NPGO at the statistical level.

For the comparison between stationary and nonstationary models, the TGAM model obviously causes lower AIC (Figure 6), showing better model performances than the GAM model. Consequently, the relationship between SKJ and the AMO during the whole period is nonstationary. Following the variations in the GCV of TGAM model, the threshold years that distinguished eras for fitting AMO-CPUE relationship separately was around 1991/92 (Figure 7). Two eras were found with various relationships between CPUE and the AMO (Table 2), where the AMO is significantly related to CPUE in the era1 (1972–1991), while this relationship was insignificant in the era2 (1992–2019, Figure 8).

As mentioned earlier, three groups migrated from distinct areas constitute most SKJ in the NWP (Figure 1). Because of the least contribution of the northeast Pacific SKJ catch to its catch in Pacific during 1972–2019 (Figure S1), and the scared long-scale SKJ movement between the northeast Pacific and NWP through the conventional tagging (Moore et al., 2020), it is assumed that the SKJ from the northeast Pacific (i.e., group 1) hardly contributes to the SKJ in the NWP. Besides, no causal linkage between NPGO and SKJ in this paper supports the assumption from another perspective (Figure S6).

In contrast, the second group from the western tropical Pacific is considered as an important component of the SKJ in the NWP. SKJ catch in the western tropical Pacific has ranked the first since 1950 (Figure S1) and was significantly related to ENSO events (Lehodey et al., 1997). During the ENSO episodes, the SKJ purse seine fishery in the western tropical Pacific showed profound displacement between the warm and cold pools with catch fluctuations. Tag-release programs and otolith studies showed that there are one or two various migratory routes between the NWP and western tropical Pacific (Arai et al., 2005; Kiyofuji et al., 2019). However, the ONI makes a weak linkage with SKJ in this study. Accounting for a weaker effect on the high-latitude extratropical area (like the NWP) of ENSO than on low-latitude tropical areas (Alexander et al., 2004), ENSO records may be faded in the NWP, which may have little effect on SKJ. Besides, recent studies found that the fluctuation of SST in the western tropical Pacific is largely in phase with AMO on the decadal scales through atmospheric teleconnection and the interaction of air-ocean dynamics (Sun et al., 2017; Sun et al., 2021). When the Atlantic is in a positive phase with positive SST anomalies, so is the western tropical Pacific, which fitly overlaps and affects the migratory environment of the second group to reach the NWP area (Figure 11). Therefore, the production of the second migratory group may largely depends on the spawning parents linked with ENSO in the western tropical Pacific, while the survival rate of the second group migrated to the NWP is most likely to rest on the AMO.

More notably, the subsurface layer of the NWP was evidenced to correlate with the AMO without lag (Wu et al., 2020b). Since the 21st century, studies have indicated that the AMO dominates the multidecadal variability of the Atlantic Ocean and exhibits significant footprints on transoceanic basins, especially on the Pacific Ocean, which reflects in modulating the variability of Pacific prominent modes (such as ENSO, PDO et al.) and inducing the synchronous response of physical environment (such as the anomalous warm of western Pacific SST, the increased upper ocean heat content and the intensified tropical cyclones during positive AMO phases) (Enfield et al., 2001; d’Orgeville & Peltier, 2007; Levine et al., 2017; Sun et al., 2017; Wang et al., 2017; Gong et al., 2020; Sun et al., 2020; Wu et al., 2020b; Zhang & Delworth, 2007). Moreover, a novel discovery is that the AMO played a decisive role in SKJ in the NWP during 1972–2019, which is consistent with the notion that the AMO could explain fish catch variability in the NWP by dominating the subtropical mode water (Figure 11) (Wu et al., 2020b).Indeed, an obvious and significant causal affection of the AMO on SST in the NWP is identified in this paper using CCM and ECCM techniques (Figures 9A,B), due to the negative optimal SST lag caused (predicted) by the AMO and the positive optimal AMO lag caused by SST (Figure 9B), which indicates that SST should be the essential oceanographic medium between the AMO and SKJ, and which is further verified using CCM test between SST and CPUEda who showed an obvious convergent trait of the prediction skill from CPUEda to SST (Figure 10). Accordingly, the SKJ inhabiting in the NWP, especially for the residents (i.e., the third group) and the immigrants from the second group, are possibly forced by the AMO through SST (Figure 11).

Previous studies have shown that environmental changes caused by climate variability could affect tunas in recruitment, spawning, distribution, etc. For instance, SST anomaly concurred with ENSO events in the Indian Ocean could reflect in the tuna fishery (e.g., yellowfin tuna and bigeye tuna) in the Indian Ocean (Syamsuddin et al., 2013; Báez et al., 2020). Although PDO is the proxy of SST anomaly signals in the north Pacific, this study shows that PDO makes no causality with SKJ in the NWP, which is consensus with the point that the NWP ecosystem could not be completely explained by PDO, of which many species fluctuations could not be associated with PDO yet on a decadal scale (Tian et al., 2014; Ma et al., 2020; Wu et al., 2020b; Ma et al., 2021). Conversely, in this study, SKJ in the NWP over 1972–2019 is caused by the AMO through SST. The plausible explanation for this is that SST variation first affects the habitat suitability of SKJ and second triggers the SKJ relative abundance anomalies. Similar phenomena was confirmed in Atlantic bluefin tuna in the north Atlantic and SKJ in the tropical Pacific (Lehodey et al., 1997; Faillettaz et al., 2019). However, the fishery ground covered the whole study area without obvious changes all the way in this paper. Therefore, the possible explanation is that the SST variability involves in the growth and development of SKJ during its lifetime. SKJ could adopt various breeding strategies depending on the external SST; it could spawn all over the year in the tropical ocean, but spawn seasonally in the extratropical water (Ashida, 2020). Moreover, the early-life ontogenetic SKJ development (juvenile) is positive and significantly related to SST (Ashida et al., 2018) too. Combined with the maturity age (1.1–1.6 year) of SKJ in the NWP (Ashida, 2020), and the optimal causal relationship between the AMO and SKJ at the lag of 15 months (Figure 5B), this paper implies the AMO effects on the recruitment process of one generation of SKJ cohort through SST. Specifically, the second group leading a bold and long-distance migration for feeding ground in the NWP (Moore et al., 2020) would experience SST variability caused by the AMO in their migratory processes, which directly influenced their survival and growth rates. On the other hand, the fecundity of the local spawning cohort and the subsequent growth and development rate of SKJ would depend on SST controlled by the AMO. Consequently, the early life of SKJ in the NWP is affected by the AMO through SST.

TGAM results show that SKJ CPUE correlates linearly with AMO at the former era (1972–1991), while the correlation disappears at the latter era (1992–2019). Although the most famous and climatic regime shift documented in the North Pacific Ocean was in 1976–1977, recent major alterations in some components of the North Pacific ecosystem were identified a further shift in 1989 in the northeast Pacific and a nonlinear and non-stationary with threshold years in the 1990s in the NWP, respectively (Hare & Mantua, 2000; Litzow et al., 2018; Ma et al., 2021). This time-dependent nonstationary responses between AMO and SKJ occurred around the early 1990s as well, which may indicate a complex ecological effect of the ecosystem. As the opportunistic feeder, SKJ could consume various preys, such as squid, crustaceans, and some mackerel or perciform fish during migration (Aoki et al., 2017). Such region-dependent feeding could probably moderate the down-top biotic effect on SKJ through the food chain, while the inverse control from predators (like sharks and billfishes in the eastern tropical Pacific Ocean, (Hunsicker et al., 2012; Chang et al., 2022) is unclear in the NWP.

Studies to date have explored the teleconnection mechanism between the AMO and Pacific ocean through the atmosphere bridge and ocean-atmosphere interaction, of which the direct oceanographic variations in the western pacific surface triggered by the AMO is SST (Sun et al., 2020; Wu et al., 2020b), In this sense, this paper explored the causation between local SST and SKJ, which is consistent with many studies that highlighted the importance of SST on tuna abundance (Zainuddin, 2011; Tangke et al., 2020; Salazar et al., 2021). The distinct relationships between the AMO and SKJ during two eras in this study may respond to the various stages in the thermal niche response curve of SKJ with warmer SST. SST before 1991 may be within the ascending stages of the niche response curve, which is slowly more suitable for the SKJ growth. While SST since 1991 could match the gentle slope stages of the curve (relative decrease but still in positive). Apart from SST, other biotic factors (such as Mixed layer depths and Chlorophyll-a) could also represent the primary habitat environment in defining SKJ suitable habitat (Mugo et al., 2020; Salazar et al., 2021). Although no obvious shift of the SKJ fishery ground is found during different eras in current study, future work is needed to explore the multi-effect from environment factors on the SKJ distribution and abundance. Additionally, the temporal evolution of AMO has shown a trend of increasing volatility in recent decades, indicating that the AMO is non-stationary (Beyraghdar Kashkooli & Modarres, 2020). Chen et al. (2019)reported that modulating effect from the AMO on the relationship between wintertime North Pacific oscillations and ENSO is different during various AMO episodes, which may further show the nonstationary relationship between the AMO and Pacific. Therefore, the nonstationary trait between AMO and SKJ may be attributed to a comprehensive process coupled with the AMO variability, the SKJ population dynamics affected by biotic and abiotic factors.

Compare with the spearman correlation results (Table S1), CCM method shows a stronger ability for determining a distinct causation between SKJ and climate variability (AMO) in this study, and provides the key evidence of the causation direction from AMO via SST effecting on SKJ. This is consistent with the notion that the CCM is superior to traditional correlation techniques on determining the caused signal with big datasets (Chang et al., 2017). With the high development of observed tools and computer science, big ocean data are birthed and how to mine the potential information under these common datasets becomes a hot issue. The CCM method is of great benefit to prompt this analysis. Moreover, the extension of CCM methods involved with multiple factors, such as Multivariate CCM and Multiview CCM (Ye & Sugihara, 2016; Hu et al., 2021) are also recommended understanding the comprehensive mechanism between organisms with climate variability and local physical factors in the ecosystem.

This study proposes that the trans-basin interaction between the AMO and Pacific Ocean has affected the NWP SKJ fishery at 15 months in advance. Similar phenomena occurred in the Indian ocean as well. Such as, the tuna fishery in the Indian ocean could suffer an obvious decline were out of phase with PDO (Báez et al., 2020; Wu et al., 2020a). Tunas or other large pelagic fish (such as swordfishes, sharks) who migrated at the meso-large spatial scale experiences large environmental variation during its whole life history, large-scale climate indices that combine many physical variables would likely serve as an appropriate proxy for explaining and predicting the long-term abundance variation of the species. With the deeper understanding of the teleconnection pattern among oceans and atmosphere under the global climate change, the trans-basin effect of climate patterns on the prediction of the large-scale migratory species should be highlighted as the same as the basin climatic oscillations in the future.