Water Temperature at Different Depths Affects the Distribution of Neon Flying Squid (Ommastrephes bartramii) in the Northwest Pacific Ocean

Ommastrephes bartramii can vertically swim during its life, and previous studies have suggested the need to account for preferred habitat distribution influenced by water temperatures at different depths. To explore the impacts of deep-water temperature on O. bartramii spatial distribution, we constructed generalized additive models (GAMs) based on the Chinese squid-jigging fishery data and the Argo deep water temperature data during 2005–2018 in the Northwest Pacific Ocean to analyze the relationships between the local abundance of O. bartramii and deep-water temperatures. The results showed that the variables including surface water temperature (T0), water temperature at 30- and 100-m depths (T30 and T100), and water differences between T0 and T30 (D0–30) significantly affected the spatial distribution of O. bartramii. The suitable ranges of each variable are different, > 15.5°C for T0, 11–18°C for T30, < 6°C for T100, and 4–4.5°C for D0–30. The areas occupied by the suitable T30 seemed to reflect the outline of fishing ground, whereas the areas with suitable T100 were to indicate the high density of O. bartramii. The predicted suitable habitat area and high-density area for O. bartramii are also regulated by El Niño–Southern Oscillation (ENSO) events. We demonstrated how the estimates of O. bartramii spatial distribution would vary influenced by deep-water temperatures in the Northwest Pacific Ocean. This information may help develop an appropriate method for investigating the effects of deep-water temperature on species with vertical migration.


INTRODUCTION
Neon flying squid (Ommastrephes bartramii) is an oceanic cephalopod economic fish species widely distributed in the subtropical and temperate waters of the Northwest Pacific Ocean, with high economic value and special ecological status (Yu et al., 2013). O. bartramii has a 1-year life cycle and is divided into two main stocks, namely, the autumn cohort and the winter-spring cohort (Yatsu, 1992;Yatsu et al., 1997). These two stocks can be further subdivided into four stocks, namely, the central stock of the autumn cohort, the eastern stock of the autumn cohort, the western stock of the winter-spring cohort, and the central-eastern stock of the winter-spring cohort, based on their geographical range (Yatsu et al., 1998). Of these stocks, the western winterspring cohort has become an important target in the locations around 150-165 • E (Chen et al., 2008), and this cohort undertakes the migration from subtropical waters around southeast Japan to the subarctic boundary during the first half of summer and then northward toward subarctic domain in August to November, matures, and starts spawning migration in fall (Murata and Hayase, 1993).
The population structure (Yatsu et al., 1997), migration (Murata and Nakamura, 1998;Ichii et al., 2009), and distribution in relation to the marine environment of O. bartramii have been extensively studied (Murata et al., 1983;Yatsu et al., 2000;Ichii et al., 2004). Reports from these studies indicate that sea surface temperature (SST) is the dominant factor affecting the distribution of O. bartramii in the Northwest Pacific Ocean (Bower and Ichii, 2005). Previous studies also suggested that changes in the abundance of O. bartramii resources may be linked to mesoscale environmental changes (El Niño and La Niña events) . Furthermore, they found that La Niña events tend to yield high habitat suitability indices, while the El Niño years result in relatively lower habitat suitability indices . During the last two decades, La Niña and El Niño occurred several times (Figure 1). The monthly SST anomaly between 5 • N and 5 • S and 120 • E-170 • W indicated that strong La Niña event happened in 2007, while strong El Niño event happened in 2015 during the main fishing seasons for O. bartramii (Figure 1; Li et al., 2020). Saputra et al. (2018) insisted that La Niña and El Niño events can also affect deep-water temperature in the Pacific Ocean. Due to frequent diel vertical migration for O. bartramii (Murata and Nakamura, 1998), deep-water temperatures should play an important role in gathering O. bartramii to format fishing grounds (Yatsu and Watanaba, 1996). Deep-water temperature has widely been used to explore the distributions of tuna. For example, Song and Zhou (2010) used water temperature gradients to build a habitat model to analyze the distribution of bigeye tuna in the Indian Ocean. Deepwater temperature and dissolved oxygen in each water layer were utilized to build a habitat model for bigeye tuna in the Indian Ocean (Feng et al., 2007). Guo and Chen (2009) also used deep-water temperature to predict the distribution of skipjack tuna in the Western and Central Pacific Ocean. At present, the application of deep-water temperature to explore the distribution of O. bartramii has been given very less attention.
In this study, we developed generalized additive models (GAMs) to explore the relationships between water temperatures at different depths and hindcast the preferred habitat distributions of O. bartramii for revealing the mechanism of formatting fishing ground-driven by water temperatures in the Northwest Pacific Ocean.

Fishery Data
Generally, each Chinese squid-jigging vessel on the main fishing ground (38 • -45 • N and 145 • -165 • E) was equipped with twin 120 kW engines, 112 kW squid attracting lights, and 16 squid-jigging machines and produced one record in a fixed position per day (Figure 2). The Chinese total annual catch accounting for 80% of global catches of O. bartramii in the Northwest Pacific Ocean (Wang et al., 2020). Thousands of fishing records in logbooks were collected and digited by the National Data Center for Distant-water Fisheries of China in the Shanghai Ocean University. The dataset consists of fishing dates (year and month), fishing position [latitude (Lat) and longitude (Lon)], catch (metric tons), and fishing effort (days fished per vessel). We split the main fishing ground (38 • -45 • N and 145 • -165 • E) for O. bartramii into 0.5 • × 0.5 • (latitude × longitude) cells. A total of 6,048 fishing records were observed by aggregating by date and location during 2005-2018.

Environmental Data
Water temperature at different depths on the fishing ground for O. bartramii was downloaded from China Argo Realtime Data Center. 1 For Argo temperature data, self-sustaining Lagrangian profiling buoys are deployed every 300 km in the global ocean, totaling around 3,000, to form a vast Argo global ocean observation network to obtain punctual, realtime, wide-area, high-resolution global ocean data. Its deepwater temperature data covers a depth range of 0-1,975 m, with 58 layers of water temperature data (Li et al., 2017). The maximum swimming depth of O. bartramii in the Northwest Pacific Ocean is less than 300 m (Tian et al., 2009a). Thus, water temperatures at different depths (T 0 , T 5 , T 10 , T 20 , T 30 , T 40 , T 50 , T 60 , T 70 , T 80 , T 90 , T 100 , T 110 , T 120 , T 130 , T 140 , T 150 , T 160 , T 170 , T 180 , T 190 , T 200 , T 220 , T 240 , T 260 , T 280 , and T 300 ) were converted into 0.5 • × 0.5 • for each month to correspond to the spatial resolution of fishery data using the "mean" method (Wang et al., 2015). For example, T-values in a layer for a 0.5 • × 0.5 • grid were averaged by 4 values of 0.25 • × 0.25 • grid. In addition, the temperature differences among different layers could generate physical barriers to confine the vertical swimming and then reshape the distributions for O. bartramii (Chen et al., 2012). We also calculated the temperature difference values between adjacent water layers (D 0

Selecting Environmental Variables and Developing Generalized Additive Models
The simulation study showed that the monthly aggregated catch in 0.5 • longitude × 0.5 • latitude cells can serve as a robust local abundance index to represent the dynamics of spatial distribution for O. bartramii in the Northwest Pacific Ocean (Wang et al., 2020). Because the aggregated Catch in metric tons is continuous and asymmetrical, with variances generally increasing with higher catches, the full GAM with a gamma distribution and log link included Catch as response variables and Month, Lon, Lat, T-values of 27 layers, and D-values of 6 intervals as explanatory variables were developed to select significant variables. The model was fitted by maximum likelihood with the "lme" function implemented in the R package "nlme" (Pinheiro et al., 2007). The collinearity between explanatory was inspected with the "check_collinearity" function in the R package "performance, " and then variables with variance inflation factor (VIF) > 2 were deleted from the full model (Lüdecke et al., 2021). Meanwhile, the significance of each variable was checked by the t-test.
where s is a spline smoother function; Tx and Ty are the temperatures at x(m) and y(m) depth; Dx-y is the gradient between x(m) and y(m) depths. The performance of the GAM with significant explanatory variables was evaluated via 100 times cross-validation, which was detailed by Wang et al. (2020). In each run, O. bartramii dataset was randomly divided into two subsets for use as training data (70%) and testing data (30%), respectively. The models fitted on the training dataset were used to predict abundance indices based on the testing data. Finally, we compared the predicted catch with the observed catch using the following simple linear regression model: Y=α+β×Y (2) where Y and Y are the predicted and observed abundance indexes of the testing data, respectively, α is the intercept, and β is the slope of the regression model. The median, α, β, R 2 , and Akaike Information Criterion (AIC) scores in cross-validation runs were selected as the model performance metrics.

Exploring the Relationships Between the Local Abundance and Water Temperatures
Sensitivity analysis can reveal the influence mechanisms of environmental variables according to plotting the predicted response variable against the predictor of interest by holding all other predictors at their mean or median values based on the established models (Elith et al., 2008). We quantified the relationships between the local abundance of O. bartramii and the water temperatures of selected depths. The synthetic datasets with interested variables were divided into 40 equal intervals between their minimum and maximum values, whereas other predictor factors set at their medium value were built to implement the analysis of water temperature sensitivity for O. bartramii. Finally, the plots of predicted local abundance against the interested variable demonstrated the changing relationships.

Integrated Environmental Preferred Habitat Distribution Maps
The predicted spatial catches using the optimal model distributed on the fishing ground for the months of July-October during 2005-2018 were depicted to reflect the dynamics of the preferred habitat. Moreover, we exhibited the habitat distribution maps overlapping with the water temperature at different layers for three distinctive scenarios (strong La Niña year in 2007, strong El Niño year in 2015, and relatively normal year in 2017) to explain how the local water temperature is driven by climate change to reshape the distribution of O. bartramii.

Model Performances
The selected significant variables checked by the VIF values and t-test were used to construct the final spatio-temporal distribution model for O. bartramii ( Table 1). The formula is as follows: The final model was validated by 100 cross-validations, with 1.47-1.62 for α and 0.015-0.017 for β (Figure 3). The average determination coefficient was more than 50% (R 2 = 0.53, AIC = 20,130, Figure 3). The performances of the model were enough for further analysis (Wang et al., 2010).

Sensitivity Analysis
According to the sensitivity analysis, the suitable ranges (scaled catch > 0.6) of each selected environmental variable were different. The suitable ranges were > 15.5 • C, 11-18 • C, and < 6 • C for T0, T30, and T100, respectively. Increasing temperatures were beneficial for O. bartramii at the surface, though unfavorable to habitat at 100-m depth ( Figure 4A). Additionally, this squid species preferred areas where the water temperature changed slightly (4-4.5 • C) between surface and 30-m depth ( Figure 4B).

Spatial Distribution of Suitable Habitat Area
Several contour lines were plotted to signify the areas of suitable or unsuitable habitat. The contour line of Scaled Catch = 0.8, Scaled Catch = 0.6, and Scaled Catch = 0.3 represents the "highdensity area" (HDA), the "suitable area" (SUA), and the "lowdensity area" (LDA), respectively. The SUA and LDA varied depending on the fishing season for O. bartramii. Generally, the SUA increased rapidly during July-August and then decreased gradually in September and October (Figure 5). In July, the LDA  was larger than SUA, and the SUA was taking shape whereas the LDA was narrowing down slowly ( Figure 5A). By August, the SUA reached its peak and occupied most of the fishing ground ( Figure 5B). In September, the SUA was moving northwest, and the LDA was expanding in the low latitude area of the fishing ground (38 • -42 • N, Figure 5C). In October, the SUA was getting smaller and fragmented, and the fishing ground was occupied by the LDA in September (Figure 5D). Moreover, the HDA is more

Suitable Habitat Overlapped With Water Temperatures
We also exhibited the preferred habitat areas overlapped with water temperatures at the surface, 30-m depth, and 100-m depth El Niño year in 2015 (Figure 6). Specifically, the suitable edge (the contour of 15.5 • C) of the surface layer could reflect the boundaries at the high latitude of the fishing ground, and the suitable ranges between 11 • C and 18 • C of the 30-m depth layer were close to the SUA of the fishing ground. Moreover, the areas with 0-6 • C of 100-m depth water layer essentially indicated the positions of the HDA of the fishing ground (Figure 6).

DISCUSSION
As a short-lived species, the distribution of O. bartramii is sensitive to the changes of ambient marine environmental variables, and water temperature is often considered for a FIGURE 3 | Boxplots of R 2 , α, β, and Akaike Information Criterion (AIC) over 100 cross-validations for the final model.  poikilothermal animal (Ichii et al., 2009;Wang et al., 2018). This study examined the potential variables of water temperatures at different depth layers and water differences between adjacent deep layers which could affect the preferences of selecting suitable habitat for O. bartramii in the Northwest Pacific Ocean. The relationships between four identified variables and the local abundance of O. bartramii are non-linear, and the complex combinations of these relationships determine the final spatial distributions of preferred habitat during fishing seasons. Previous researches have shown that the variable SST (T 0 in this study) plays a critical role in the variability of distribution and could be an important indicator of the HDAs for O. bartramii. Tian et al. (2009a,b) constructed an HSI model and GAM to estimate the optimal SST range with 14-21 • C. Alabia et al. (2015) suggested that the formation of suitable fishing grounds for the O. bartramii was mainly related to the SST in summer and winter based on the HSI model (11.6-18 • C in summer, 7-17 • C in winter). These findings are consistent with our results of suitable T 0 with > 15.5 • C.
Importantly, we had detected more critical water temperature indicators (the contour of 15.5 • C of the surface layer, 11-18 • C of 30-m depth layer, and 0-6 • C of 100-m depth layer) of suitable habitat area and HDA for O. bartramii. Especially, the combination of the steady horizontal water environment and the mild varied vertical water environment on the fishing ground could assemble squid easily (Figures 4B, 6). However, those cannot be validated directly from the perspective of O. bartramii biology in the Northwest Pacific Ocean due to very few studies existing at present.
To aid interpretability, we divided the O. bartramii fishing ground into 12 grids of 5 • longitude × 2 • latitude and then plotted the vertical profiles of mean water temperature for each grid in three distinctive years (2007, 2015, and 2017) using Argo temperature datasets (Figure 7). Intuitively, three rules of formatting abundant O. bartramii fishing ground in the Northwest Pacific Ocean could be concluded: (1) needing more lines of vertical temperature across horizontal suitable ranges; (2) needing more lines of vertical temperature mixed each other (steady horizontal water environment); (3) needing more lines of vertical temperature between surface and 30-m-depth layer parallel to oblique suitable water difference (slight varying vertical water environment). The first two are observed easily in La Niña year (2007) and normal year (2007), whereas the third only occurred in the HDAs (grids of F, G, and H in Figure 7) in August and September 2007. All three rarely happened together in El Niño year (2015) (Figure 7). Our interesting and primary explorations could be useful guidelines for future study.
In the Northwest Pacific Ocean, the known north-south seasonal migration during the life cycle may cause variations in local abundance on the fishing ground for the months of July-October for O. bartramii. This migration pattern could also be observed clearly according to the changing of SUA on the fishing ground. In July, O. bartramii starts entering the fishing ground for foraging until August with an increase in the abundance, and then, they leave from the fishing ground for spawning with the decrease in population during late October (Figure 5; Bower and Ichii, 2005).
The migration of the O. bartramii population is also related to changes in its availability of feeding. Hikaru et al. (2004) mentioned that O. bartramii was mainly distributed in the transition zone south of the subarctic boundary (TZ) during May-July and then migrated to the transition domain north of the subarctic boundary (TD) after July, which is corresponding to the change of SUA in our study. The vertical moving of preferred water layers for O. bartramii during May-October is directly decided by their diet changed from zooplankton to micronekton and fish during migration revealed by the stomach contents (Hikaru et al., 2004). El Niño and La Niña events have a significant impact on the primary productivity , which eventually altered the horizontal and vertical distribution of diet for O. bartramii in the Northwest Pacific ocean.
Empirically, the La Niña events tend to yield more SUAs, while the El Niño years resulted in relatively less SUAs (Tzeng et al., 2012;Syamsuddin et al., 2013;Yu et al., 2015). Moreover, we found that the La Niña events appeared to generate more HDAs for squid aggregation which could be determined by water temperatures at 100-m depth. The El Niño events fragmented the suitable habitat areas resulting in more LDAs and thus not conducive for O. bartramii growth.
In summary, the relationships between the local abundance of O. bartramii and water temperatures at different depths were estimated, and then, the preferred habitat distributions were investigated in the Northwest Pacific Ocean. The methods of the study may be useful for other marine species with the biological characteristics of vertical swimming. The results could provide basic information for sustainable management, such as the establishment of a marine protected area for O. bartramii in the Northwest Pacific Ocean.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
JW and YC carried out the experiment and wrote the manuscript. HL helped supervise the project and provided financial support. LL and JZ downloaded and processed the Argo data. XC and JW conceived of the presented idea. All authors contributed to the article and approved the submitted version.

FUNDING
This study was financially supported by the National Key R&D Program of China (2019YFD0901402 and 2019YFD0901401) and the National Natural Science Funding of China (NSFC41876141).