The fishery data of O. bartramii and D. gigas were collected from the Chinese Squid-jigging Technology Group of National Data Centers of Chinese Distant-Water Fisheries, Shanghai Ocean University. Fishing activities were operated from 2003 to 2020. The data comprised the time, location (latitude and longitude), and catch (unit: ton) of each squid fishing fleets’ operation. The stock/group studied in the analysis was chosen based on the fishery history and biological information described as follows.

According to catch data from 2003 to 2020, the Chinese main fishing ground of O. bartramii in the Northwest Pacific Ocean was in the western water of 170°E, 35°N-45°N (Area 1, Figure 1. The boundary is based on the historical fishing location for O. bartramii. The latitude and longitude points for drawing the boundary of Area 1 are presented in Supplementary Table 1). During 2003∼2018, Area 1 accounted for 82.35–99.95% of the overall catch (Figure 2A). These proportions only had a decrease in 2019 (79.35%) and 2020 (63.36%). Previous studies revealed that the O. bartramii had two cohorts: the western winter-spring cohort and the eastern autumn cohort (Yatsu et al., 1997; Nagasawa et al., 1998). The stock of the winter-spring cohort is mainly located at the western water of 170°E, while the autumn cohort is primarily distributed in the eastern area of 170°E (Chen and Chiu, 2003; Wang and Chen, 2005; Yu et al., 2016a). Through the biological analysis based on hard tissue (beak and statolith), significant differences in the population structure had been observed (such as sex ratio, growth, and age, Jin, 2015), as well as trophic position and migration patterns (Fang et al., 2016a,b) for the individuals captured between the eastern and western waters of 170°E. Based on the fishery history and the biological differences, the fishery data of O. bartramii within Area 1 were used, and all catches were simply considered as the stock of the winter-spring cohort in our analysis.

Before 2012, the main fishing grounds of D. gigas included the waters off PEEZ, off the Exclusive Economic Zone of Costa Rica or Chile (Chen et al., 2008; Liu et al., 2013a). After 2012, some squid-jigging fleets transferred to the equatorial waters (EW, 5°N-5°S, 95°-120°W), which has become a new fishing ground for D. gigas. Catch proportions in different regions exhibited high fluctuations in the 2000s due to the exploration of the new fishing grounds. Nevertheless, from 2011 to 2019, catch in the waters off PEEZ (Area 2, Figure 1, the boundary is based on the historical fishing location for D. gigas. The latitude and longitude points for drawing the boundary of Area 2 are listed in Supplementary Table 1), accounted for 61.02–88.60% of the overall catch (Figure 2B). D. gigas in different geographic regions have various size structures, age compositions (Liu et al., 2013a,2015b; Gong et al., 2017), and migration patterns (Liu et al., 2015a). Considering the complicated population classification, spawning ground locations, migration routes, and the fishery history of D. gigas, we only chose the fishery data of D. gigas within Area 2 for analysis in this study.

Cooper (2006) proposed that quantifying the fishery stock abundance could be determined from both the fishery-dependent data and the fishery-independent survey. Catch, effort or CPUE data derived from the fishery-dependent data, such as landing records, logbooks, vessel monitoring systems, and onboard observers, are commonly regarded as a crucial signal of fishery stock abundance dynamics (Cooper, 2006; Tian et al., 2009b; Pauly et al., 2013). A fishery-independent survey utilizing standardized sampling gears, including trawls, seines, hydroacoustics, and video, is expensive, particularly for the high sea fishery stocks (Cooper, 2006; Dennis et al., 2015). For example, for O. bartramii, only Japanese scientists conducted the driftnet surveys for the stock of the eastern autumn cohort before 2021 (Ichii et al., 2011; Igarashi et al., 2017). Currently, most studies of environment- and climate-induced oceanic squid stock abundance dynamics based on historical Chinese Squid-jigging fleets operation data consider the nominal CPUE as the squid abundance index (for example, Chen et al., 2007; Cao et al., 2009; Yu et al., 2017a,2021a,2021b). We followed their choice. And the nominal CPUE was defined as:

where CPUE denotes the yearly nominal CPUE (ton/days) at Area 1 or Area 2; Catch represents the total yearly catch for all squid-jigging fleets at Area 1 or Area 2; Effort indicates total yearly fishing days at Area 1 or Area 2. The nominal CPUE could be regarded as a credible index of squid abundance for three reasons: (1) According to the collected squid-jigging fleets’ parameters and historical reports (Wang and Chen, 2005), Chinese squid-jigging fleets for O. bartramii and D. gigas had almost similar fishing powers, lamps, and engines; all fleets operated one time in a night during fishing season. (2) Study areas were restricted at the main fishing grounds for these two squids (Area 1 and Area 2) and without bycatch and (3) previous studies have done the CPUE standardization works to investigate if the catchability is constant for the O. bartramii and D. gigas abundance based on the same fishery data source in our analysis and suggested that the standardized CPUE and nominal CPUE were similar and had the same trends (Tian et al., 2009a; Xu et al., 2017). Thus, although catchability can be impacted by fishers’ behaviors, fishing gear type, ocean environment, management and other factors (Campbell, 2004), it is reasonable to infer that that the nominal CPUE is proportional to the squid abundance in this study.

Monthly SST data in the Pacific Ocean from 2002 to 2020 were available at the Oceanwatch database from NOAA.¹ The data were provided by the NOAA Coral Reef Watch (CRW) program with 5 km × 5 km spatial resolution.

ONI and PDOI were applied as the representatives for the climate events El Niño and La Niña and PDO. The ONI indicated the 3-month running mean of SST anomalies (NOAA ERSST.V5) in the Niño 3.4 region (5°N∼5°S, 120°∼170°W). The PDOI is the first principal component of monthly SST anomalies in the North Pacific Ocean (poleward of 20°N) based on the empirical orthogonal function. These data covered the period 2002–2020 and had the time resolution of a month, available in the NOAA Physical Sciences Laboratory.² Since the study aims to explore the abundance fluctuations of squids at annual and interannual scales, monthly ONI and PDOI were calculated as annual mean in the analysis (Supplementary Figure 1).

The annual and interannual abundance trends of O. bartramii and D. gigas and their PC_i were quantified by the Hodrick—Prescott filter analysis. Given a time series (TS), the method can be considered an additive decomposition: TS = ts + cs, where ts and cs denote the trend series and the cyclical series, respectively. It is a data-smoothing technique by minimizing an objective function that penalizes the sum of squares for cs with the sum of squares of second-order differences for ts in a time series (Hodrick and Prescott, 1997; Adiabouah, 2020). The objective function can be written as (Hodrick and Prescott, 1997):

where T represents the time size, and λ refers to the smooth parameter and is set to 6. Through this method, it can be assumed that the annual fluctuations are short-term variations that belong to the cs and can be removed. Thus, the long-term trend (ts) defined as an interannual trend can be better observed.

We hypothesized that impacts of habitat temperature variabilities separately correspond to abundance fluctuations at different time scales (annual and interannual, respectively). And considering that the associations might consist of the non-linear effects, the GAM (Zhao et al., 2014; Wood, 2017) was used to detect the relationships between the ts of abundance (TA) and its PC_i (TPC_i) and between the cs of abundance (CA) and its PC_i (CPC_i) so as to validate this assumption. The model can be written as follows (Zhao et al., 2014):

where n denotes the number of PC_i we determined; α refers to the model intercept; s () indicates the smoothing function of TPC_i and CPC_i in which k represents the effective degrees of freedom determining the degree of non-linearity. k was restricted to a maximum of 3 to avoid over-fitting. ε suggests the random error with the normal distribution N (0, σ²). The Restricted Maximum Likelihood Estimation (RMLE) method was applied to fit the model. Moreover, the double penalty approach was used for model fitting (Marra and Wood, 2011), in which a second penalty on each smoothing function was added and only affected the functions in the null space, that is, one smoothing function might be flat (means no effect).

Regarding model validation and choice of the key TPC_i or CPC_i that influence the abundance fluctuations, a cross-validation test would be conducted (the 18-year data were divided into two parts: training data for 16 years and validation data for 2 years). Training data were selected to build the model. Then, the model was used to predict the validation data. Whenever possible 16 training and 2 validation data set combinations within 2003∼2020 were used at the cross-validation test. Therefore, 153 sub-models were established. After the construction of all sub-models, model validation was evaluated by the performance of the predictability to the actual TA or CA and its changing trends. Specifically, statistical measures and regression analysis were applied in the model validation step.

(1) Statistical measures

The predictability to the actual TA or CA was verified by the mean absolute error (MAE) and root mean squared error (RMSE) between the validating and actual TA or CA. The predictability to the changing trend of TA or CA was validated by the mean absolute percentage change error (MAPCE) and root mean squared percentage change error (RMSPCE) according to the actual TA or CA values in the last year. The MAPCE and RMSPCE can be written as:

where X_m,i and X_i denote the model predicted and actual TA or CA in year i, respectively, and X_i–1 represents the actual TA or CA in the last year i-1. Particularly, if i = 2003, the TA or CA in 2003 and 2004 (X_m,i–X_i–1 = X_m,2003–X₂₀₀₄ and X_i–X_i–1 = X₂₀₀₃–X₂₀₀₄) was used to calculate the MAPCE and RMSPCE.

(2) Regression analysis

Regression analysis is a common approach to evaluate the model goodness of fit. In the analysis, a simple regression model is used as follows:

where A and M denote the actual and model-predicted TA, CA, and their percentage change; α and β are the regression parameters indicating their slope and intercept; if α and β are not significantly different from 1 and 0, this implies that the model has an unbiased predictive skill (Chang et al., 2010; Zhao et al., 2014; Li et al., 2017).

After the model validation, key TPC_i and CPC_i were determined by the proportion of one TPC_i or CPC_i identified as significant (significant level: 0.05) in the sub-models for the 153 runs.

Key TPC_i or CPC_i was investigated by its corrections with ONI and PDOI. Previous studies demonstrated that El Niño and La Niña and PDO could affect the habitat and abundance fluctuation of O. bartramii and D. gigas in 0∼1-year effect (Chen et al., 2007; Yu et al., 2021a,b). Therefore, ONI, PDOI, and the values in 1-year lag (ONI_–1 and PDOI_–1) were included in the analysis. The best linear combination of climate indices (ONI, ONI_–1 PDOI, and PDOI_–1) was identified based on the step regression method with the minimum Akaike information criterion (AIC) to explain the variations of key TPC_i or CPC_i (Venables and Ripley, 2002).

All analyses in this study were implemented by the statistical software “R,” version 4.1.0 with package “mgcv” (version 1.8-35, Wood, 2017) and “mFilter” (version.1-5, Balcilar, 2019). The important abbreviations and their meaning for the statement of Results and Discussion are provided in Table A1.

The overall abundance, TA, and CA of O. bartramii are illustrated in Figure 3A. Within 2003∼2020, the overall abundance of O. bartramii changed in 1.25∼4.11 t/day. Its TA indicated that the abundance of O. bartramii had a downward interannual trend from 2003 to 2013 and tended to increase after 2014. Its CA suggested that its abundance at an annual scale had about 3∼4-year cycles (defined as a period during relative lowest or highest value) such as 2012–2014 (3 years), 2006–2009 (4 years). The amplitude (differences between the lowest and highest value in one cycle) for each cycle was not similar. For example, the CA changed in –0.032∼0.1 t/day during the cycle of 2003∼2005, but in –0.89∼1.24 t/day within the next cycle (2006∼2009).

Within 2003∼2020, the overall abundance of D. gigas changed in 2.99∼7.10 t/day (Figure 3B). TA reflected that its interannual abundance trend kept decreasing from 2003 to 2018. But it remained steady during 2007∼2014 and started to decrease clearly after 2014. Moreover, the CA of D. gigas had a 3-year cycle within 2003∼2007 and a 5∼year cycle within 2007∼2020. Similar to O. bartramii, the amplitude of D. gigas for the 3∼year cycle was not similar. Nonetheless, the amplitude was similar (±1t/day) within the 5∼year cycle (2007∼2020).

Briefly, results showed different but not synchronous trends of abundance for O. bartramii and D. gigas. Concerning O. bartramii, the interannual abundance first decreased (2003∼2013) and then increased (after 2014). Regarding D. gigas, the interannual abundance kept decreasing within 2003∼2020. Their annual trends presented large fluctuations, with the changing amplitude greater than the increase or decrease speed of the interannual trend. There were no clear synchronous variations for the overall abundance, TA, and CA between O. bartramii and D. gigas (The correlation coefficients for the overall abundance, TA, and CA were 0.09, 0.33, and –0.07, respectively).

The scree plot of eigenvalue and cumulative variance proportion of the SST at the spawning and feeding ground for O. bartramii indicated that six PC_i could be retained (Figure 4A). The first six PC_i totally explained 94.68% of the total variation. The results of D. gigas suggested that the first five PC_i need to be retained, accounting for 92.30% of the total variation (Figure 4B). Details of the PCA results, such as the retained loadings to each PC_i, are presented in Supplementary Tables 2, 3.

Figure 5 exhibits the overall, annual, and interannual trends of the first six PC_i (TPC_i and CPC_i) for O. bartramii from the Hodrick—Prescott filter analysis. Each TPC_i had its independent characteristics within 2003∼2020. TPC₁ possessed one and a half cycles with its increase within 2003∼2009 and 2016∼2020 and decrease within 2010∼2016 (Figure 5A). TPC₂ first decreased and then increased after 2012 (Figure 5B). TPC₃, TPC₄, and TPC₅ first increased, decreased, and increased, and then decreased within 2003∼2020, respectively (Figures 5C,E). TPC₆ changed little compared with the amplitude of other TPC_i (Figure 5F). There was no particular pattern of change for each CPC_i. Cycles appeared irregularly and could continue about 3∼5 years such as 2007∼2019 for CPC₂ (a cycle in 3 years, Figure 5B) and 2007∼2012 for CPC₄ (a cycle in 5 years, Figure 5D). Moreover, the amplitude in different cycles was not common. For example, a cycle of CPC₂ changed from –1.72 to 2.38 within 2015∼2018, while the range narrowed to –0.32∼0.43 for the cycle within 2018∼2020 (Figure 5B).

Overall, annual and interannual trends of the first five PC_i for D. gigas are exhibited in Figure 6. Similar to the results of O. bartramii, each TPC_i had its trend with independent shape through time. TPC₁ had an entire cycle within 2011∼2020 while slightly decreasing before 2011 (Figure 6A). TPC₂ had a cycle of increasing in 2003∼2007 and decreasing after 2011 (Figure 6B). Within 2003∼2020, TPC₃, TPC₄, and TPC₅ generally increased, decreased, and slightly decreased, respectively (Figures 6C–E). There was no particular pattern of change for each CPC_i. Irregular cycles and different amplitudes could also be observed on the cycles for each CPC_i (Figure 6).

Briefly, each TPC_i and CPC_i for these two squids had its independent characteristics. CPC_i had larger fluctuations than those of TPC_i by the irregular cycles and different amplitudes along the years.

Sub-models of GAM for TA to O. bartramii and D. gigas could perform acceptable predictions for the validation data (Figures 7, 8). MAE and RMSE for O. bartramii were 0.031 t/day and 0.046 t/day, respectively. Regarding D. gigas, MAE and RMSE were 0.11t/day and 0.14 t/day, respectively. MAPCE and RMAPCE for the TA of O. bartramii were 1.39 and 2.04%, respectively. MAPCE and RMAPCE for the TA of D. gigas were 2.48 and 3.04%, respectively. The errors mainly occurred when TA_i-TA_i–1 was near 0% (–1∼3%, Figures 7C, 8C). On the contrary, sub-models of GAM for CA to O. bartramii and D. gigas presented poor predictions for the validation data (Figures 9, 10). MAE and RMSE for O. bartramii were 0.49 t/day and 0.67 t/day, respectively. For D. gigas, MAE and RMSE were 0.89 t/day and 1.05 t/day, respectively. MAPCE and RMAPCE for the CA of O. bartramii were 477.24 and 1054.93%, respectively, and 439.05 and 262.48% for the CA of D. gigas.

Regression analysis (Table 1) showed the model-predicted TA of O. bartramii positively changed with actual TA. The slope α was closer 1 (0.989), but α and intercept β were significantly different from 1 and 0 (p < 0.05). However, the coefficients α and β for the changes in the TA percentage did not present significant difference with 1 and 0 (p > 0.05). Regarding D. gigas, the model-predicted TA positively changed with the actual TA and the coefficients α and β were not significantly different from 1 and 0 (p > 0.05). But coefficients α and β for the TA percentage change had significant difference with 1 and 0 (p < 0.05). For the CA models, the null hypothesis (H₀: α = 1, β = 0) was believed to be rejected because P-values for the coefficients α and β were less than 0.05 (except for the β in the regression of the CA for O. bartramii).

In summary, GAM performed well on the relationship between the abundance trend and its habitat temperature variabilities at the interannual scale while presenting poor results at the annual scale for both O. bartramii and D. gigas. In terms of the MAE, MSE, MAPCE, and RMAPCE values, GAM exhibited a better result for O. bartramii in the relationship of TA∼TPC_i than D. gigas. But from the results of regression analysis, GAM gave better predictions for the TA changing trend but had the slightly weaker predictive ability for the TA value of O. bartramii than D. gigas.

Given the poor results of CA∼CPC_i for both two squids, only key TPC_i would be determined. According to the proportion of each TPC_i identified as significant in the 153 sub-models of GAM (Table 2), TPC₁, TPC₄, and TPC₅ were considered the primary TPC_i affecting the abundance of O. bartramii with all the sub-model identifying them as significant (p < 0.05). TPC₃ was the second essential TPC_i with the proportion confirming as significant (p < 0.05) in sub-models at 92.81%. TPC₂ and TPC₆ were regarded as the least essential TPC_i with only 32.68% and 18.95% of sub-models identifying it as significant. For D. gigas, TPC₅ was taken as the most primary TPC_i with all the sub-models identifying it as significant (p < 0.05, Table 2), followed by TPC₁ (99.35%) and TPC₄ (85.62%). The least important TPC_i was TPC₂ and TPC₃ with the proportion as significance in sub-models as 54.90 and 29.41%, respectively. Therefore, TPC₁, TPC₃, TPC₄, and TPC₅ for O. bartramii and TPC₁, TPC₄, and TPC₅ for D. gigas were believed to be the key TPC_i affecting the interannual abundance trend of O. bartramii and D. gigas. Thus, they were used for further analysis.

Additionally, GAM was operated with all-year data (2003∼2020) to explore the impact of each TPC_i on TA (Supplementary Table 2). The impact of TPC₁, TPC₃, TPC₄, and TPC₅ on TA for O. bartramii decreased, increased, decreased, and increased, respectively. The impact of TPC₁, TPC₄, and TPC₅ on TA for D. gigas decreased, increased, stopped changing, and increased, respectively (Supplementary Table 2 and Supplementary Figure 2). Details of smooth function plots are provided in Supplementary Figure 2.

According to step regression analysis with minimum AIC value (Table 3), TPC₁ for O. bartramii could be explained by PDOI. TPC₁ and PDOI had a negative relationship with the decrease that occurred when PDO changed from the warm phase to the cold phase, and vice versa (Figure 6A and Supplementary Figure 1). TPC₅ could be explained by PDOI and PDOI_–1. Nonetheless, the regression result was not convincing because no significant relationship was found from the ANOVA test (p > 0.05), lower R² (0.19), and no significant parameters. No combination of parameters was demonstrated to be better compared to the null model for TPC₃ and TPC₄ (TPC₃/TPC₄ ∼ Intercept). Regression results for D. gigas suggested that TPC₁ could be explained by PDOI. TPC₁ and PDOI had a positive relationship with the value generally along the trend of PDOI (Figure 6A and Supplementary Figure 1). TPC₅ could be explained by PDOI and PDOI_–1, although the regression was not significant (P = 0.07). Furthermore, TPC₄ was related to ONI and ONI_–1, and both the model and its parameters were significant (p < 0.05).

It was assumed in this study that the squids’ habitat temperature annual and interannual variability was separately related to the abundance fluctuation at the corresponding scale. The hypothesis was tested by the GAM model and cross-validation test. Results demonstrated that the model only performed well on modeling the TA (Figures 7–10). This indicated that the abundance fluctuations of O. bartramii and D. gigas could be investigated by their habitat temperature variabilities at the interannual scale instead of the annual scale. Potential mechanisms might be induced by the complicated relationship between ocean environment and oceanic squids’ unique life histories. The oceanic squids are short-lived species with a life span of 1 year for O. bartramii and 1∼2 years for D. gigas (O’Dor, 1998; Ichii et al., 2009; Arkhipkin et al., 2015b). During their short life histories, the abundance fluctuations are highly sensitive to the large- and small-scale environmental changes (Anderson and Rodhouse, 2001). The abundance-environment relationship below the annual scale has been widely investigated in previous studies. For example, Ichii et al. (2011) proposed the stock abundance of the autumn cohort of O. bartramii coincided with the Chl-a concentration in their spawning ground and feeding ground. At the feeding ground of the O. bartramii, Nishikawa et al. (2014) discovered that the higher CPUE usually occurred at the location with higher SST or lower SST anomaly. Yu et al. (2016c) observed that food availability on the spawning ground, SST, and Chl-a concentration on the fishing ground together influenced the stock abundance of the winter-spring cohort of O. bartramii. With respect to D. gigas, Waluda et al. (2006) and Waluda and Rodhouse (2006) reported that the D. gigas abundance off Peru waters was associated with the mesoscale SST variability controlled by El Niño and La Niña. In the Gulf of California, the concentration of D. gigas was controlled by SST and Chl-a (Robinson et al., 2013). Yu et al. (2019) indicated that the seasonal CPUE was affected by the SST, net primary production (NPP), and sea surface height anomaly (SSHA). In short, many results of previous studies confirmed that the intra- and annual abundance fluctuations of O. bartramii and D. gigas were sensitive to environmental changes, including, not only temperature, but also other factors such as physical oceanic environmental parameters (such as Chl-a, NPP, and SSHA). It might be the reason for a poor relationship found on CA∼CPC since we only used habitat temperature as the explanatory factor (Figures 9, 10).

On the contrary, the results suggested that the interannual abundance fluctuations of O. bartramii and D. gigas could be investigated by their habitat temperature variabilities. The CA is observed to have large fluctuations influenced by many environmental factors. In fact, for squid species with 1-year life span, these environmental factors and recruitment dynamics are only related to their abundance in the same year (in other words, abundance is not associated with its spawning biomass or the abundance in previous years, O’Dor, 1998; Ichii et al., 2009; Arkhipkin et al., 2015a). Yu et al. (2016c) verified that the oceanic squids could respond quickly to regional and large-scale environmental conditions. Therefore, TA may reflect only the effects of long-term environmental and climate variabilities on the squids’ abundance fluctuations when the short-term CA is removed.

Interestingly, why is the habitat temperature essential to the effects of long-term environmental and climate variabilities? The annual and interannual abundance fluctuations of these two squids had been widely reported by Yatsu et al. (2000), Igarashi et al. (2017), Alabia et al. (2020) and Yu et al. (2021a,b) for O. bartramii and Taipe et al. (2001) and Yu et al. (2021a,b) for D. gigas. All studies considered SST and revealed powerful relationships wherein the habitat temperature factor is the most imperative one in squids’ abundance dynamics. From the biological perspective, Sakai et al. (2004) discovered that higher temperatures within optimum ranges could improve squids’ recruitment by reducing the hatching time and benefiting embryonic development. The abundance of the stock of the winter-spring cohort of O. bartramii could be largely controlled by the proportion of favorable-SST areas at the spawning ground (Cao et al., 2009). Additionally, squid growth responds to temperature mostly in nature (Jackson and Domeier, 2003; Forsythe, 2004). Temperature is also essential for the maturity and reproduction of O. bartramii and D. gigas (Bower and Ichii, 2005; Liu et al., 2013a). In short, temperature can be the most important one affecting the entire life history of O. bartramii and D. gigas (Yatsu et al., 2000; Ichii et al., 2009; Yu et al., 2021a). From another perspective, SST, as commonly used data in the fishery oceanography, is also adopted to describe the meso-scale and large-scale oceanic processes, such as fronts (Redelsperger et al., 2019), eddys (Jia et al., 2019), El Niño and La Niña (Hayashi et al., 2020), PDO (Mantua et al., 1997), and ocean warming (Cheng et al., 2019). Thus, temperature and temperature-related oceanic processes related to their abundance might be consequently reflected on their interannual abundance fluctuations, considering the temperature plays a key role in squids’ life history.

As demonstrated by our cross-validation results, the TA and its changing trends of validation data could not be all predicted well (Table 1 and Figures 7C,D, 8C,D). Based on our fishery data, an additional hindcast for the TA of O. bartramii in 1995∼2002 was made by the GAM model for the data over 2003∼2020 (Supplementary Figure 3). The results indicated the overall shape of TA between the hindcast and actual value was similar; however, an error occurred on the TA-changing trend in 1998∼1999. Thus, a conclusion cannot be easily extended to the interannual scale that the temperature plays a critical role in the interannual abundance fluctuations, owing to some potential mechanisms. Previous studies also indicated that other factors (partly) independent from temperature were also important. For example, according to the commercial driftnet fishery monitoring data, Yatsu et al. (2000) observed the interannual stock abundance of the autumn cohort of O. bartramii and discovered relatively low-level abundance in 1978∼1992 and an increasing trend after 1992. They revealed that the interannual variability might be caused by fishing mortality and SST. Nishikawa et al. (2014, 2015), Igarashi et al. (2017), and Yu et al. (2016a) presented that the wind-induced Ekman transport and mixed layer depth could be associated with the abundance of O. bartramii. CPUE of the autumn cohort of O. bartramii was considered to be related to the location of the transition zone chlorophyll front (TZCF, Ichii et al., 2011). Medellín-Ortiz et al. (2016) demonstrated that the population size of D. gigas along Baja California’s west coast was significantly affected by the primary productivity, catches of Pacific sardine (Sardinops caeruleus), and North Pacific Gyre Oscillation (NPGO). To summarize, temperature can be connected with other oceanic processes, while these factors would not be completely included in the information of habitat temperature variabilities, resulting in the errors of the model predictive performance. Thus, other potential factors should be considered in the future to explore the differences in their importance. This can deepen our understanding of the interaction between the interannual abundance fluctuations and environmental variabilities for O. bartramii and D. gigas.

Our results have three findings related to PDO regarding these two squids. (1) Based on the GAM results (Supplementary Table 2 and Supplementary Figure 2), the impacts of TPC₁ to TA of both O. bartramii and D. gigas were similar (decrease); (2) TPC₁ of O. bartramii was negatively related to PDOI, while TPC₁ of D. gigas had a positive correlation with PDOI; (3) the PDO-related TPC₁ between O. bartramii and D. gigas had a significant negative relationship (r = –0.76, p < 0.01). Thus, the mechanism of the effect of PDO on the O. bartramii and D. gigas is that PDO could induce opposite conditions of habitat temperature variabilities for these two squid species while acting similarly to the interannual abundance fluctuations for O. bartramii and D. gigas through habitat temperature variabilities. This is consistent with the conclusions of Yu et al. (2021b) who investigated the habitat variations of O. bartramii and D. gigas in 1950∼2015 and their relationship with PDO. They revealed that, during the cold PDO phase, SST anomaly increased but led to a decrease for the area with suitable SST for O. bartramii, while the decrease of SST anomaly led to the formation of areas with more suitable SST for D. gigas, and vice versa. What is more, SST and other environmental factors were usually included in the habitat suitability index model (Yu et al., 2017b,2021b; Alabia et al., 2020). This model explained well in the association between the PDO and the spatial and temporal abundance variabilities for O. bartramii and D. gigas and, again, demonstrated PDO-induced opposite habitat conditions between these two squids.

In our studies, by the cross-validation test of GAM, TPC₁ of O. bartramii was believed to be the primary TPC_i affecting the interannual abundance trend. Regarding D. gigas, TPC₁ was the second important factor (Table 2). This TPC_i was related to PDO according to the step regression analysis (Table 3). Moreover, the PDO affects the abundance of O. bartramii and D. gigas through its habitat environment (including temperature), in line with the statement of previous studies, and TPC₁ of O. bartramii and D. gigas accounts for the most proportion of variance of habitat temperature variabilities in PCA. It can be concluded that the PDO is an essential climatic factor affecting the interannual abundance fluctuations of squids through habitat temperature variabilities.

A limitation of our studies is that ONI and PDOI were used in the annual mean, only to conduct the preliminary analysis between habitat temperature variabilities and climate variabilities. It can be confirmed in the PDO situations since the warm or cold phases last for many years. Nevertheless, the monthly ONI data indicated that the events of El Niño and La Niña last for 5∼19 months in 2003∼2020. As a result, we might be unable to well examine the impacts of El Niño and La Niña in the relationship between TPC_i and the annual mean of ONI. From the step regressions results, TPC₄ of D. gigas was also a crucial TPC_i and related to ONI and ONI_–1 (Tables 2, 3) while not all sub-models in the cross-validation test identified TPC₄ as significant (85.62%, Table 2). None of the key TPC_i of O. bartramii could be explained by ONI and ONI_–1 (Table 3). The effects of El Niño and La Niña on these two squids’ abundance have also been widely documented (Ichii et al., 2002; Waluda et al., 2006; Chen et al., 2007; Yu et al., 2021a). Moreover, the long-term abundance fluctuation of D. gigas associated with the environmental changes by El Niño was reported at the Gulf of California (Robinson et al., 2016; Frawley et al., 2019). After the El Niño events in 2009, the landings of D. gigas declined due to the warmer waters, weak winter/spring winds, and extremely low Chl-a along the Gulf of California in the following years. From our analysis, a similar trend in decrease can be found in TA of D. gigas after the 2015∼2016 El Niño events (Figure 3B) and TPC₄ declined (Figure 6D). From the GAM results, TA increased with TPC₄ (Supplementary Figure 2B). But the TA of D. gigas kept steady during 2007∼2014 and TPC₄ increased after the El Niño events in 2009 (Figures 3B, 6D). In previous studies, only a 2-month lag effect between the SST anomaly at Niño 3.4 and Niño 1 + 2 regions was observed, as well as the SST at these two squids’ fishing (feeding) ground, which consequently affected their abundance (Yu et al., 2021a). This indicates the stock/group of O. bartramii and D. gigas we studied responses quickly for the El Niño and La Niña when habitat temperature conditions changed. Therefore, it is not clear whether the effects of El Niño and La Niña on the abundance fluctuations of D. gigas we studied in the waters off PEEZ are just below the annual scale or can be extended to interannual scale. It is suggested that further research needs to collect the data with a wider year range to evaluate it.

Additionally, TPC₅ of O. bartramii and D. gigas was also an important TPC_i affecting the interannual abundance trend based on the cross-validation test of GAM. However, it was not related well to PDO and PDO_–1 (Table 3, P = 0.07). Moreover, there is another key TPC_i (TPC₃ and TPC₄ for O. bartramii, Table 2) that could not be explained by the PDOI and ONI we examined. One possible reason is that other environmental and climate variabilities drive the variations of key TPC_i. These factors, which are independent of the El Niño and La Niña and PDO effects but related to habitat temperature variabilities, indirectly affect the interannual fluctuations of squid abundance. In other studies of cephalopod’s fishery, for example, common Chinese cuttlefish catches were demonstrated to be negatively affected by both Asian Monsoon Index (AMI) and winter SST in the East China Sea (Pang et al., 2018). Meanwhile, the SST had an association with Asian Monsoon (Li et al., 2008; Bao and Ren, 2014). November Arctic Oscillation had been reported to impact the subsequent SST at tropical Pacific Ocean (Chen et al., 2015, 2018). There is evidence that the Arctic Oscillation-related sea ice drives the variations of tropical Pacific SST (Rigor et al., 2002; Purich et al., 2016). These emphasized again the importance of combining other potential factors in future works.

Finally, we conclude that the PDO and El Niño and La Niña can impact the abundance fluctuations in the same way but play an opposite role for the suitable habitat between the O. bartramii and D. gigas from the results of Yu et al. (2021a,b) and our analysis. Nevertheless, it is not correct that their abundance had synchronous variations because of PDO or El Niño and La Niña. For example, no significant relationship could be found between these two squids if we viewed the abundance series from the same fishery data source on the annual and interannual scales (Figure 3). And there was no clear increase in the abundance fluctuations if we observe the abundance of D. gigas at the cold PDO phase within 2007∼2013 (Figure 3B and Supplementary Figure 1). From our studies, it can be confirmed that, at least, on the interannual scale, there is another potential key factor related to habitat temperature variabilities affecting the abundance fluctuations. The results from Yu et al. (2021a,b) were based on the seasonal fishery data (from September to November) on the fishing ground and resulted from similar life histories in the feeding season at the feeding ground. Even though the seasonal variations are crucial to the stock assessment and fishery management, studies of the overall, annual, and interannual abundance dynamics are essential for keeping a long-term sustainable fishing activity and supporting the studies of the ecosystem regime shift. Therefore, their complex life histories and related small- and large-scale environmental and climate variabilities should be comprehensively considered in future analysis to better understand the abundance dynamics of O. bartramii and D. gigas.