Habitat suitability modelling revealing environmental-driven abundance variability and geographical distribution shift of winter–spring cohort of neon flying squid Ommastrephes bartramii in the northwest Pacific Ocean

Abstract

To identify climate-related habitat changes and variations in abundance and distribution of Ommastrephes bartramii in the northwest Pacific Ocean, an habitat suitability index (HSI) model was developed in this study including sea surface temperature (SST), photosynthetically active radiation (PAR), and sea surface height anomaly (SSHA). The catch-per-unit-effort (CPUE) of O. bartramii gradually decreased between 2006 and 2015, and the latitudinal gravity centres (LATG) of the fishing effort shifted southward. Correlation analyses suggested that CPUE was positively related to SST, PAR, and the areas of suitable and optimal habitat, but negatively correlated with SSHA and the percentages of poor habitat. A significantly positive correlation was found between the LATG and the average latitude of the most preferred SST, PAR, and the average latitude of the area with an HSI between 0.9 and 1.0. From 2006 to 2015, the annual declined CPUE was highly consistent with the increase in areas of poor habitat and the reduction in areas of suitable and optimal habitats. The south-approaching LATG coincided with the southward migration pattern of the latitude of the HSI area in the range of 0.9–1.0. Moreover, compared to the El Niño events, the La Niña events, and normal climate condition yielded enlarged suitable habitat areas for O. bartramii, and the LATG moved further north.

Introduction

Large-scale climate variability, interacting with regional oceanic conditions, has produced profound impacts on the quality of habitats for fish, causing abundance variability, and extensive distribution shifts (Brander, 2010). Understanding how environmental conditions affect fish stocks is an important aspect of conserving and sustainably exploiting their resources (Brierley and Kingsford, 2009). For the short-lived ommastrephid squid, the interaction between environmental variability and squid species tends to be particularly prominent, due to the high sensitivity of squid to various environmental conditions (Anderson and Rodhouse, 2001). Numerous studies have reported that changes in the abundance and distribution of ommastrephid squid are strongly related to climate-driven environmental conditions in spawning and feeding grounds (Postuma and Gasalla, 2010; Robinson et al., 2013; Alabia, Dehara, et al., 2016; Kooij et al., 2016). Different ocean conditions may lead to various impacts on squid stocks (Sakurai et al., 2002). Therefore, assessing the impacts of varying ocean conditions on ecologically and economically important squid stocks is helpful to develop and implement appropriate management strategies.

The neon flying squid Ommastrephes bartramii is an oceanic squid extensively distributed in the subtropical and temperate waters of the world’s oceans (Yatsu et al., 1997). Ommastrephes bartramii plays a vital role in marine ecosystems, serving as a connecting link between top predators and low level prey (Watanabe et al., 2004). In the North Pacific, the O. bartramii population consists of two reproductive cohorts: an autumn cohort and a winter–spring cohort (Yatsu et al., 1998). Both cohorts have a large commercial value and attract international fishing vessels from China, Japan, Korea, and Chinese Taipei (Chen, 2010; Arkhipkin et al., 2015). The winter–spring cohort comprises two geographical stocks: the western stock and the central stock. The former is mainly targeted by Chinese squid-jigging fishing vessels, with annual catches fluctuating from 80 000 to 100 000 tons (Chen et al., 2008). Within its 1-year lifespan, the western winter–spring O. bartramii stock migrates from its spawning ground in subtropical waters to its feeding ground in subarctic waters (Bower and Ichii, 2005). The feeding ground of O. bartramii is located in the Kuroshio–Oyashio Current System (Chen et al., 2012), which is also characterized as important habitats for other epipelagic fish such as anchovy (Engraulis japonicus) (Komatsu et al., 2002) and Pacific saury Cololabis saira (Tian et al., 2003).

Ommastrephes bartramii is an environmental sensitive squid species showing interannual variability in abundance and distribution (Yu, Chen, Yi, et al., 2015; Ichii et al., 2017). Due to their different geographic distribution and life history, the autumn and winter–spring cohorts of O. bartramii may be influenced by different oceanographic processes but there are also some similarities between the two cohorts (Ichii et al., 2009; Nishikawa et al., 2014, 2015; Igarashi et al., 2017). For the autumn cohort of O. bartramii, Japanese researchers have conducted many studies to examine a possible link between annual variations in stock levels and oceanographic conditions (Ichii et al., 2011; Alabia et al., 2015a, b). They pointed out that changes in the abundance of the autumn cohort of O. bartramii can be largely attributed to the effects of the winter and summer positions of the Transition Zone Chlorophyll Front on the nursery and feeding grounds. Food availability is affected by climate variations and related to changes in the upper ocean structure, and is a critical factor in squid recruitment (Ichii et al., 2011). For the winter–spring cohort, most findings support the view that climate variability, such as El Niño and La Niña events, appears primarily to regulate the major oceanographic features across the feeding grounds of O. bartramii, and further influences their habitat distribution in ways that force abundance variations (Chen et al., 2007; Cao et al., 2009; Yu, Chen, Yi, Chen, et al., 2015). These studies show the importance of interannual climatic variability for this cohort relative to multi-decadal climatic events. The former event has yielded strong impacts on squid abundance and distribution and ultimately affected the annual fisheries catch. However, there has been little research evaluating habitat quality of O. bartramii in terms of biological response (e.g. abundance and distribution) to climatic and environmental conditions (e.g. ENSO, water temperature, and sea surface height).

The habitat for pelagic fish species is closely related to environmental conditions in the oceans. Examining the relationships between fishery data and the environments around the fishing ground provides a good way to evaluate the habitat preference and detect the habitat hotspots for fish species. The habitat modelling approach is a way of modelling the distribution of a species and has provided useful tools for identifying the habitats of species and relating these to climate variability (Elith et al., 2010). A climate-driven habitat suitability index (HSI) model is quantitatively determined using the associations between a range of environmental parameters and a relative species abundance or occurrence, but not explicitly including ecological factors such as predator–prey interactions (Tanaka and Chen, 2016). Given that the HSI model can generate superior model performance and reliable predictions in ecological studies in comparison to other statistical models, the HSI model has been extensively applied to predict the spatial and temporal distribution of suitable habitats for many demersal and pelagic fish species (Chang et al., 2012; Chen, 2013; Li et al., 2017).

In this study, we developed an integrated HSI model to understand climate-related abundance variations and distribution shifts of the western winter–spring stock of O. bartramii in the northwest Pacific Ocean from 2006 to 2015. Based on the weighted arithmetic mean method, the HSI model was constructed using fishery data from 2006 to 2014 using three crucial environmental variables: sea surface temperature (SST), sea surface height anomaly (SSHA), and photosynthetically active radiation (PAR). The data from 2015 were applied to test and validate the model’s performance. The purposes of this study are (i) to characterize and identify the squid suitable habitat in relation to environmental conditions; (ii) to propose a possible mechanism causing the variability in squid abundance and distribution with respect to the physical and biological consequences of the El Niño and La Niña events.

Material and methods

Squid-jigging fishery data

This study used 10-year logbook data collected by the Chinese Squid-jigging Science and Technology Group of Shanghai Ocean University in the autumn (September to November) from 2006 to 2015. The total study area spanned between 36–48°N and 150–170°E in the northwest Pacific Ocean. Annual catches of O. bartramii exploited by Chinese squid-jigging fishing vessels accounted for more than 80% of the total squid catches in the study area. The data contained fishing time (year and month), fishing effort (fishing days), the location of fishing activity (latitude and longitude in degrees), and daily catch (tons). The fisheries dataset in this study were highly quality-controlled and included 90 438 original fishing samplings during 2006–2015. Previous studies have evaluated the effects of spatial and temporal scales on habitat modelling for O. bartramii in the Northwest Pacific Ocean, they suggested that the optimal temporal and spatial scales with the lowest coefficients of variation were month and 0.5° grid scale (Gong et al., 2014). Therefore, the fishery data in this study were aggregated into a spatial resolution grid of 0.5° latitude × 0.5° longitude and by month. The study area in Figure 1 showed the location of the Chinese squid-jigging fishing ground for the western winter–spring cohort of O. bartramii.

Figure 1.

Location of fishing grounds for western winter–spring cohort of neon flying squid O. bartramii in the Northwest Pacific Ocean.

Among the total catches of Chinese O. bartramii fisheries, the western stock of winter–spring cohort covered more than 95% and there was no bycatch species. Moreover, all the Chinese squid-jigging vessels were modified from the same type of inshore bottom trawler, and were equipped with almost identical engines (120 KW × 2), squid-attracting lamps (112 KW × 16) and fishing powers and all operated at night (Wang and Chen, 2005). Due to the homogeneous characteristics and operating behaviours of the Chinese squid-jigging vessels, catch-per-unit-effort (CPUE) could be considered as a potential index of squid abundance on the fishing ground (Chen et al., 2007; Cao et al., 2009). In this study, the monthly CPUE within a 0.5 × 0.5° fishing grid was calculated by the following equation (Cao et al., 2009):

$${\text{CPUE}}_{\mathit{ymij}}=\frac{\sum {\text{Catch}}_{\mathit{ymij}}}{\sum {\text{Effort}}_{\mathit{ymij}}},$$

(1)

where ${\text{CPUE}}_{\mathit{ymij}}$⁠, $\sum {\text{Catch}}_{\mathit{ymij}}$⁠, and $\sum {\text{Effort}}_{\mathit{ymij}}$ are the monthly nominal CPUE [tons (t)/day (d)], the total catch for all fishing vessels within a fishing grid, and the sum of fishing efforts of all fishing vessels within one fishing grid, respectively, at longitude i, latitude j in month m, year y.

Environmental data

The selection of biophysical environmental variables for habitat modelling was largely informed by factors that predominantly influence the abundance and distribution of O. bartramii stock. SST was regarded as a critical factor influencing squid physiological metabolism, and had a superior predictive power to detect the fishing ground of O. bartramii (Yu, Chen, Yi, et al., 2015). The SSHA field was coupled with stratified upper mixed layers. Previous studies suggested that changes in SSHA could affect O. bartramii distribution (Chen et al., 2010). In addition, PAR played an important role in controlling ocean productivity (Wu et al., 2011), which might have influences on food availability for O. bartramii. Therefore, SST, PAR, and SSHA were selected and used as the environmental predictors in the development of the HSI model in this study.

The daily SST dataset was compiled from the National Oceanic and Atmospheric Administration (NOAA) Optimum Interpolation (OI) SST Version 2 with a 0.25 × 0.25° spatial resolution (https://www.ncdc.noaa.gov/oisst/data-access). The OI analysis was constructed by combining observations from different platforms (satellites, ships, buoys) on a regular global grid. To compensate for platform differences and sensor biases, this methodology included a bias adjustment for satellite and ship observations. Remotely sensed satellite data, the monthly Moderate Resolution Imaging Spectroradiometer PAR data were obtained from the NOAA OceanWatch dataset (http://oceanwatch.pifsc.noaa.gov/thredds/catalog.html). The spatial resolution of PAR was 0.05 × 0.05°. The daily gridded SSHA data with 0.25° grid resolution were obtained from the altimeters of the TOPEX, ERS, and JASON-1 satellites. The information was processed by the Archiving, Validation, and Interpretation of Satellite Oceanographic Data scientific team of the Collecte, Localisation, Satellite data centre (http://apdrc.soest.hawaii.edu/data/data.php). To examine the relationship between PAR and ocean productivity on the fishing ground of O. bartramii, the monthly ocean net primary productivity (NPP) data were also obtained from the high-accuracy standard products estimated from the Vertically Generalized Production Model developed by Behrenfeld and Falkowski (1997) (http://www.science.oregonstate.edu/ocean.productivity/index.php). All of the environmental data covered the study area and fishing years and were grouped on a monthly and 0.5 × 0.5° grid to match with the fishery data.

The El Niño and La Niña events were measured by the SST anomaly in the Niño 3.4 region above or below a threshold of ±0.5°C over at least five consecutive months. The definition was obtained from the NOAA Climate Prediction Center from website: http://www.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ensoyears.shtml. According to the definition of El Niño and La Niña events, there were three El Niño events (2006, 2009, and 2015) and three La Niña events (2007, 2010, and 2011) in the autumn period between 2006 and 2015, and the other years were ENSO-neutral climate years.

Developing an HSI model

HSI modelling, incorporated with one or more key environmental variables, can be used to create probability maps for identifying the availability of fish species (Vayghan et al., 2013). In this study, an integrated HSI model combined with three input variables (SST, PAR, and SSHA) was developed for O. bartramii stock. Before the habitat model was developed, the three parameters for collinearity were examined by SPSS software. To accurately predict the habitat suitability of O. bartramii, both CPUE and fishing effort were taken into account in relation to the biophysical environmental variables. CPUE was considered as an index of squid abundance and fishing effort was considered as an index of squid occurrence and/or availability, both have been successfully applied in the establishment of HSI models for fish species (Li et al., 2016).

The first step was to compute the probability of O. bartramii availability, also defined as a suitability index (SI) based on the associations between CPUE and environmental variables, and between the fishing effort and the variables by histogram graphs. Using the data from 2006 to 2014, the CPUE-based SI value was defined as the CPUE for each environmental class interval divided by the maximum CPUE for a given class interval in a specific month. The effort-based SI value was defined as the total fishing effort for each environmental class interval divided by the maximum fishing effort for a given class interval in a specific month. For each input variable (SST, PAR, and SSHA), the ultimate SI value was calculated by the average probability indices from CPUE and fishing effort. The equation was established as follows (Zainuddin et al., 2006):

$$\text{SI=}\frac{1}{2}({\text{SI}}_{i,\text{CPUE}}+{\text{SI}}_{i,\text{effort}}),$$

(1)

where SI_i,CPUE was the SI value calculated based on CPUE; SI_i,effort was the SI value calculated based on effort; i was month.

Then the calculated SI values for each environmental variable were used as the input data to fit the SI models with each class interval value. The SI model of each variable for each calendar month i (September, October, and November) was fitted by the following equations:

$${\text{SI}}_{i,\text{SST}}=\hspace{0.17em}\text{exp}[a{{(\text{SST}}_i-b)}^2],$$

(2)

$${\text{SI}}_{i,\text{PAR}}=\hspace{0.17em}\text{exp}[a{{(\text{PAR}}_i-b)}^2],$$

(3)

$${\text{SI}}_{i,\text{SSHA}}=\hspace{0.17em}\text{exp}[a{{(\text{SSHA}}_i-b)}^2],$$

(4)

where a and b are the estimated model parameters, which were solved with a least squares estimate to minimize the residuals between the observed and predicted SI values; SST, PAR, and SSHA are the class interval value in month i. Suitable environmental ranges for O. bartramii were identified as SI values above 0.6 based on the frequency distribution of fishing effort and CPUE in relation to the three environmental variables (Chang et al., 2013).

Finally, the SI models derived from each variable, for each month, were further combined into the integrated HSI model for that month. There are many commonly applied empirical models, such as the Arithmetic Mean Model (AMM), the Geometric Mean Model (GMM), the Continued Product Model (CPM), the Minimum Model (MINM), and the Maximum Model (MAXM) (Gong, et al., 2011). However, the AMM showed the best model performance and could accurately predict the suitable habitat for squid species (Chen et al., 2010; Yu et al., 2016). Moreover, different environmental variables might have different impacts with on squid habitat (Gong et al., 2012). Thus, the weighted AMM method was ultimately employed as the empirical HSI model to predict the habitat suitability of O. bartramii on the fishing ground in the northwest Pacific Ocean. The weighted AMM-based HSI model for each calendar month (September, October, and November) was developed using the following equation (Xue et al., 2017):

$$\text{HSI}=\frac{1}{{\sum }_{i=1}^n{w}_i}\sum _{i=1}^n(\mathrm{S}{\mathrm{I}}_i\times w_i),$$

(5)

where n is the number of environmental variables included in the HSI model; SI_i is SI values for each environmental variable i; w_i is the weight of variable i. In this study, ten different weighting scenarios were considered for w_SST, w_PAR, and w_SSHA (case1 [0 1 0]; case2 [0 0 1]; case3 [0.1 0.8 0.1]; case4 [0.1 0.1 0.8]; case5 [0.25 0.5 0.25]; case6 [0.25 0.25 0.5]; case7 [0.333 0.333 0.333]; case8 [0.5 0.25 0.25]; case9 [0.8 0.1 0.1]; case10 [1 0 0]) (Yi et al., 2017). The predicted HSI values ranged from 0 to 1. The areas with HSI ≥ 0.6, HSI ≥ 0.8, 0.2 < HSI < 0.6, and HSI ≤ 0.2 were regarded as suitable, optimal, normal, and poor habitat, respectively, for O. bartramii.

Model test and validation

The HSI values from September to November over 2006–2015 were predicted under different weighting scenarios. Under each HSI class interval ([0.0 0.1]; [0.1 0.2]; [0.2 0.3]; [0.3 0.4]; [0.4 0.5]; [0.5 0.6]; [0.6 0.7]; [0.7 0.8]; [0.8 0.9]; [0.9 1.0]), the parameters including the fishing grid squares with catch, CPUE, proportion of catch in each HSI class interval accounting for the total catches, and proportion of fishing effort in each HSI class interval accounting for the total fishing efforts were examined and compared for each scenario. Generally, a high HSI class interval corresponded to a high proportion of catch and fishing effort, as well as CPUE. Based on the HSI theory, the weighted AMM-based HSI model with the best prediction performance was chosen to predict the habitat suitability for O. bartramii. Moreover, the predictive ability of the HSI model with the best model performance was further evaluated by a cross-validation analysis. The fishery and environmental data during 2006–2014, representing 90% of all the data, were used for HSI development. The remaining 10% data (i.e. the data in 2015) were used to evaluate the HSI performance. The predicted maps of HSI in autumn 2015 were then overlaid with fishery data (fishing effort) to assess the model performance for predicting the potential suitable/optimal habitat of O. bartramii on the fishing ground. The goodness of fit of the linear relationship between HSI values and CPUE was also examined using R².

Linkage between climate-driven habitat variations and autumn abundance and distribution of O. bartramii

To examine the linkage between the climate-driven habitat variations and autumn abundance and distribution of O. bartramii, we first used the selected HSI model to predict the HSI values on the fishing ground from September to November during 2006–2015. The CPUE and latitudinal gravity centres of fishing effort (LATG) were used as proxies to indicate squid abundance and distribution, respectively, for O. bartramii stock. Statistical relationships between the monthly squid CPUE and six habitat variables (SST, PAR, SSHA, average HSI, suitable habitat areas, and optimal habitat areas) were examined. Similarly, the linear correlation coefficients were calculated between the autumn squid LATG and the following habitat variables: average latitude of the monthly most preferred SST, PAR, and SSHA, latitudinal gravity centres of HSI (LATG_HSI), and average latitude of HSI between 0.9 and 1.0, SI_SST between 0.9 and 1.0, SI_PAR between 0.9 and 1.0, and SI_SSHA between 0.9 and 1.0. The equations for calculating the LATG and LATG_HSI were in the form of Li et al. (2016):

$$\text{LATG}=\frac{\sum \left({\text{Latitude}}_{(i,m)}{\text{×Effort}}_{(i,m)}\right)}{\sum {\text{Effort}}_{(i,m)}},$$

(6)

$${\text{LATG}}_{\text{HSI}}=\frac{\sum \left({\text{Latitude}}_{(i,m)}{\text{×HSI}}_{(i,m)}\right)}{\sum {\text{HSI}}_{(i,m)}},$$

(7)

where Latitude_(i,m) is the latitude of the ith fishing unit in month m; Effort_(i,m) is the total fishing effort within the ith fishing unit in month m; HSI_(i,m) is the HSI value within the ith fishing unit in month m.

To better understand the process whereby ENSO influences squid abundance and distribution of O. bartramii in the northwest Pacific Ocean, the CPUE, LATG, and the area and distribution of suitable habitat were compared between composites of environmental and fishery data for the three El Niño (2006, 2009, and 2015) and three La Niña (2007, 2010, and 2011) conditions occurring in our study period. The oceanographic features of the fishing ground were also examined and compared among three different climate years (an El Niño year with extremely small suitable habitat, a normal year and a La Niña year with especially large suitable habitat) (Yu, Chen, Yi, Chen, et al., 2015).

Results

Variations in CPUE and LATG

Interannual and monthly variability were found in the September to November CPUE and LATG of O. bartramii during 2006–2015 (Figure 2). CPUE tended to be higher during 2008–2008 and relatively lower over 2009–2015. The LATG was likely to move southward from 2006 to 2015 with fluctuations. Both presented a negative trend from year to year, especially for LATG.

Figure 2.

Interannual variability of CPUE and LATG for O. bartramii from September to November during 2006–2015.

SI curves for each environmental variable

SI curves based on CPUE and fishing effort in relation to each environmental variable (SST, PAR, and SSHA) were established during September to November, as shown in Table 1 and Supplementary Data S1. Statistical fits for all these SI models were significant (p < 0.001) with high correlation coefficients and low Root Mean Square Errors (RMSEs). The environmental preferences of O. bartramii varied monthly, as shown from the SI curves (see Supplementary Data S1). In September, the inferred suitable ranges for SST, PAR, and SSHA, determined from the suitable SI values (SI ≥ 0.6), were 16–19°C, 25–29 E/m²/day, and −4 to 4 cm, respectively. In October, the suitable range of each variable was from 13 to 17°C, from 19 to 21 E/m²/day, and from −6 to 6 cm, for SST, PAR, and SSHA, respectively. However, the suitable SI values (SI ≥ 0.6) in November were associated with SST between 11 and 17°C, with PAR between 10 and 14 E/m²/day, and with SSHA between −6 and 6 cm.

Table 1.

Monthly fitted SI model of each environmental variable for O. bartramii in the Northwest Pacific Ocean.

Month	SI model	a	b	SSE	RMSE	R2	p
September	SISST	−0.148	17.282	0.050	0.004	0.970	<0.001
	SIPAR	−0.134	27.055	0.051	0.003	0.977	<0.001
	SISSHA	−0.030	−0.098	0.135	0.008	0.934	<0.001
October	SISST	−0.217	15.036	0.116	0.007	0.921	<0.001
	SIPAR	−0.293	19.722	0.075	0.007	0.939	<0.001
	SISSHA	−0.012	−0.600	0.059	0.003	0.979	<0.001
November	SISST	−0.155	12.894	0.039	0.003	0.977	<0.001
	SIPAR	−0.165	12.116	0.063	0.007	0.958	<0.001
	SISSHA	−0.020	−1.127	0.170	0.009	0.903	<0.001

SSE, sum of squares for error; RMSE, root mean square error.

HSI model validation

Comparing the number of fishing grids with catch, proportion of catch and fishing effort and CPUE in each HSI class interval under different weighting scenarios (see Supplementary Data S2), it was found that the case9 (weight of 0.8, 0.1, 0.1 corresponding to SST, PAR, and SSHA, respectively) accounted for the lowest catch and fishing effort in the poor habitat (HSI ≤ 0.2) and the largest catch and fishing effort in the suitable habitat (HSI ≥ 0.6). The above parameters under case9 were estimated and shown in Table 2. Results indicated that, with the HSI values increasing, all the above parameters generally showed a significantly increasing trend. The poor habitat, with HSI ≤ 0.2, attracted only 2.12% of catches and 2.41% of fishing efforts, indicating that these areas were not suitable for O. bartramii. The common habitat areas, with HSI values between 0.2 and 0.6, represented 20.73% of catches and 23.11% of fishing efforts. While in the suitable habitats with HSI ≥ 0.6, the highest percentages of catches and fishing efforts occurred and accounted for 77.15 and 74.48%, respectively (Table 2). These findings suggested that case9 was the best HSI model due to its good predictive performance. Comparing to PAR and SSHA, SST tended to be more important in the formation of suitable habitat for O. bartramii.

Table 2.

The parameters used for each HSI class interval correspond to fishing grids square with catch, proportion of catch, proportion of fishing effort, and CPUE.

HSI	Fishing grids square with catch	Proportion of catch (%)	Proportion of fishing effort (%)	CPUE (t/d)
0.0–0.1	94	0.48	0.65	1.45
0.1–0.2	105	1.64	1.76	1.51
0.2–0.3	122	2.16	3.31	1.61
0.3–0.4	123	6.17	5.34	1.87
0.4–0.5	132	4.01	5.75	1.81
0.5–0.6	148	8.39	8.71	2.10
0.6–0.7	192	11.05	9.74	2.10
0.7–0.8	271	19.78	17.96	2.18
0.8–0.9	285	22.27	21.19	2.09
0.9–1.0	301	24.04	25.58	2.22

The results from this table were determined from case9.

Furthermore, the integrated HSI models from case9 developed using data from 2006 to 2014 were cross-validated by the fishing efforts of the squid-jigging vessels from September to November in 2015, as shown in Supplementary Data S3. The predicted HSI values within the fishing ground were mapped and superimposed on the observed fishing efforts to examine the agreement of spatial distribution of the HSI values and actual fishery. It was observed that fishing locations in each month mostly distributed over suitable habitats at HSI ≥ 0.6. Higher fishing efforts generally occurred in the waters with higher HSI values. In addition, the goodness of fit of the relationship between the observed CPUE and the HSI (0–1 interval) was evaluated. The model was significant (p < 0.001) and could explain 89.5% of the variance (see Supplementary Data S4). These findings suggested that the weighted AMM-based HSI model from case9 in our study provided a reliable prediction of the suitable habitats for O. bartramii in the northwest Pacific Ocean.

Squid abundance and distribution in relation to habitat changes

Correlation analyses suggested that the autumn O. bartramii CPUEs during 2006–2015 were significantly positively related to SST, PAR, averaged HSI, percentage of suitable habitat (areas with HSI ≥ 0.6) and optimal habitat (areas with HSI ≥ 0.8) accounting for the whole fishing ground, but with a significant negative correlation with the SSHA and percentage of poor habitat (areas with HSI ≤ 0.2) (Table 3). These correlations implied that the O. bartramii CPUE was largely affected by the local environmental conditions on the fishing ground, especially the areas of the suitable and optimal habitats. Thus, the areas of poor, suitable, and optimal habitats from 2006 to 2015 were further examined as well as the averaged HSI (Figure 3). It was found that the monthly averaged HSI, and the areas of suitable and optimal habitat for O. bartramii, showed a dramatically decreasing trend, with a similar variability pattern for the CPUE. In contrast, the poor habitat areas for O. bartramii remarkably enlarged from 2006 to 2015, particularly in 2009 and 2015.

Figure 3.

(a) Monthly average HSI values on the fishing ground over 2006–2015; and the percentage of areas with (b) HSI ≤ 0.2, (c) HSI ≥ 0.6, (d) HSI ≥ 0.8 accounting for the fishing ground during autumn 2006–2015.

Table 3.

Correlations between different environmental variables and monthly squid abundance (indicated by CPUE) from September to November over 2006–2015.

Environmental variables	Correlation
	r	p
Average SST	0.361	0.025
Average PAR	0.367	0.023
Average SSHA	−0.446	0.007
Average HSI	0.433	0.008
Area with HSI ≥ 0.6	0.458	0.005
Area with HSI ≥ 0.8	0.476	0.004
Area with HSI ≤ 0.2	−0.476	0.004

Note that statistically significant to p < 0.05.

The O. bartramii LATG from September to November over 2006–2015 was significantly and positively correlated with the average latitude of the monthly most preferred SST and PAR, LATG_HSI (latitudinal gravity centres of HSI) and the average latitude of the area with HSI between 0.9 and 1.0 on the fishing ground. However, the LATG was not significantly correlated with the latitude of monthly preferred SSHA (Table 4). Furthermore, a high consistency was found between the south-approaching LATG from 2006 to 2015 and LATG_HSI as well as the latitude of the area with HSI in the range of 0.9–1.0 (Figure 4). This suggested that the distribution of O. bartramii was mainly dominated by the latitudinal location of suitable habitats, and the SST and PAR were more important in affecting the squid distribution than the SSHA.

Figure 4.

(a) Latitudinal gravity centres of HSI and (b) average latitude of HSI between 0.9 and 1.0 during autumn 2006–2015.

Table 4.

Correlations between different environmental variables and monthly squid distribution (indicated by LATG) from September to November over 2006–2015.

Environmental variables	Correlation
	r	p
Average latitude of the monthly most preferred SST	0.583	0.000
Average latitude of the monthly most preferred PAR	0.307	0.049
Average latitude of the monthly most preferred SSHA	0.205	0.139
LATGHSI	0.541	0.001
Average latitude of HSI between 0.9 and 1.0	0.583	0.000
Average latitude of SISST between 0.9 and 1.0	0.513	0.002
Average latitude of SIPAR between 0.9 and 1.0	0.578	0.000
Average latitude of SISSHA between 0.9 and 1.0	0.055	0.387

Note that statistically significant to p < 0.05.

Impacts of anomalous climatic events on squid habitats

To evaluate the impacts of ENSO variability on squid, the CPUE, LATG, and habitat conditions of O. bartramii were compared under the three different ENSO states (Figure 5). During the El Niño years, CPUEs were at a relatively low level, and the LATGs mainly distributed in the more southern waters. In contrast, the CPUEs were higher and the LATGs shifted to higher latitudes in the ENSO-neutral and La Niña years. Furthermore, a significant increase was found in the poor habitat areas, and the average latitude of HSI was in the range of 0.9–1.0 and located in the lower latitude under the El Niño state in comparison to an ENSO-neutral or La Niña state. In addition, it was found that the monthly mean SST and PAR anomalies during the El Niño years on the fishing ground of O. bartramii were lower than those in the ENSO-neutral years and La Niña years. However, the mean SSHA during the El Niño years was basically similar to the ENSO-neutral years, but was relatively higher than that in the La Niña years (Figure 6).

Figure 5.

The El Niño and La Niña signature in the (a) CPUE; (b) percentage of poor habitat on the fishing ground; (c) LATG; (d) average latitude of HSI between 0.9 and 1.0 for O. bartramii in the Northwest Pacific Ocean during autumn 2006–2015.

Figure 6.

Monthly mean SST, PAR, and SSH anomalies from September to November during El Niño years, ENSO-neutral years and La Niña years. The red dotted line indicated the average values for each variable under each climate condition.

To illustrate any differences in the habitat conditions of O. bartramii, as a response to the anomalous climatic events, we compared the habitat suitability and environmental conditions on the fishing ground in the northwest Pacific Ocean, focusing particularly on November 2010 (a strong La Niña year), 2008 (a normal climate year), and 2009 (a strong El Niño year). Figure 7 shows the spatial distribution of the habitat suitability. In 2010 and 2008, higher habitat suitability values were zonally distributed between 40 and 46°N. However in 2015, the suitable habitats became narrower, confined to the region between 40 and 43°N.

Figure 7.

Comparison of the predicted HSI values on the fishing ground of O. bartramii in November 2010 with a La Niña event, 2009 with an El Niño event and 2008 with normal climate condition.

The differences in November SST, PAR, and SSHA distribution among the three years were clearly presented in Figure 8. Positive SST and PAR difference values were widely distributed between 40 and 48°N, which coincides with the spatial location of the suitable habitat, indicating a higher SST and PAR in 2008 and 2010 compared to 2009. Moreover, the contour lines of the most preferred SST (13°C) and PAR (11 E/m²/day) for O. bartramii were located at a higher northern latitude in November 2008 and 2010 compared to 2009. Concerning the SSHA, however, the broad pattern of negative difference values shown in Figure 8 indicated that SSHA in 2008 and 2010 was lower than that in 2009 on the fishing ground of O. bartramii. In addition, we also examined the NPP on the fishing ground. Higher NPP values occurred on the fishing ground of O. bartramii in November 2008 and 2010, the NPP in 2009 was relatively low (Figure 9).

Figure 8.

Difference in SST, PAR, and SSHA between 2008 and 2009, and between 2010 and 2009 in November. The green contour lines indicated the preferred environmental variables for O. bartramii in 2010 (left panel) and 2008 (right panel); the yellow contour lines indicated the preferred environmental variables for O. bartramii in 2009.

Figure 9.

Comparison of the NPP on the fishing ground of O. bartramii in November 2010 with a La Niña event, 2009 with an El Niño event and 2008 with normal climate condition.

Discussion

This study reveals that the habitat suitability model developed for the western stock of the winter–spring cohort of O. bartramii in the northwest Pacific Ocean effectively captures this species’ variations under different climate conditions. The oceanographic-determined habitat changes estimated with our habitat suitability modelling are consistent with the observed variability and trends in effort and CPUE data for O. bartramii over the 2006–2015 periods. In this work, the AMM method was adopted to calculate the integrated HSI values. Evidence has shown that model results for habitat suitability from GMM, MINM, and CPM will be underestimated, and the HSI values may be very low. The MAXM is likely to overestimate the HSIs (Gong et al., 2011). The AMM refers to the combination of all the SIs into an HSI via the arithmetic mean. It assumes that each factor contributes equally to the model (Chen et al., 2010; Silva et al., 2016; Xue et al., 2017). However, fluctuations in any of the input factors in the model will produce large effects on the outputs. Therefore, variations in the environmental variables can be accurately reflected in the results from the AMM. Many studies support that the AMM can improve the understanding of the species–habitat relationships and have wide applications for fish such as swordfish Xiphias gladius (Chang et al., 2013), skipjack Katsuwonus pelamis (Yen et al., 2017), chub mackerel Scomber japonicus (Chen et al., 2009), as well as jumbo flying squid Dosidicus gigas (Yu et al., 2016). To better predict spatial–temporal distribution of O. bartramii habitat, we evaluated the impacts of differential weighting of environmental variables on HSI modelling. Our findings also suggest that the weighted AMM-based HSI model (weights of 0.8, 0.1, 0.1 corresponding to SST, PAR, and SSHA) in this study was the best model for evaluating the effects of climate-driven environmental changes on habitat suitability of O. bartramii, and for predicting the spatial dynamic of optimal habitats under different climate conditions.

In most fish habitat models, CPUE/effort is used as an abundance/occurrence index in developing the HSI model. However, for commercial fisheries, CPUE is not always a reliable abundance index (Bordalo-Machado, 2006). Fishermen tend to target fish based on their experience, and they know the areas where the fish are concentrated. Thus, over time, fishing vessels will generally operate in the most productive fishing grounds with locally high fish abundance, with a non-random distribution of fishing effort (Li et al., 2014). Consequently, CPUE tends to be biased to indicate an abundant fish population. Moreover, fishing effort could indicate fish occurrence, while it was also an indicator of fishing vessel distribution and could not necessarily represent real fish abundance. For squid-jigging fisheries, Tian et al. (2009) have compared the CPUE- and effort-based HSI models. They found that the CPUE-based HSI model tended to overestimate the ranges of optimal habitats and underestimate monthly variations in the spatial distribution of optimal habitats. Zainuddin et al. (2006) have suggested that considering both CPUE and fishing effort data into habitat hotspots modelling can accurately predict the fish habitat such as tuna. With those considerations, we chose a CPUE-effort coupled-based HSI model to define the optimal habitats for O. bartramii in the northwest Pacific Ocean.

It is clear that O. bartramii is not randomly distributed relative to environmental conditions. The suitable range of SST for O. bartramii varied monthly with much lower ranges in October and November. High fishing effort was found within specific ranges of each variable, which is generally consistent with results from previous studies (Fan et al., 2009; Tian et al., 2009). These favourable ranges yielded the highest probability of finding productive O. bartramii fishing grounds and represented suitable habitats with elevated squid abundance (i.e. high CPUE). SST, in this study, was regarded as a predictor of physiological condition of the O. bartramii habitat (Wang et al., 2010). In this study, a positive relationship was found between monthly mean CPUE and mean SST during 2006–2015. In general, O. bartramii tend to prefer warmer waters, and an increase in water surface temperature could favour squid physiology (Zheng and Chen, 2008). It has been suggested that, with their limited tolerance of an unfavourable SST range in the original habitat, O. bartramii will actively move to other regions in pursuit of an optimal habitat at comfortable temperatures (Chen et al., 2007). This will directly lead to the movement of the LATG to the most preferred SST in each month. Ommastrephes bartramii may also respond to oceanographic variables (e.g. SSHA) in an indirect way. The SSHA field is coupled with the dynamic of upper oceans and can be an effective indicator to predict squid abundance and distribution (Ichii et al., 2011; Gong et al., 2012). Our results showed a negative correlation between SSHA and CPUE, indicating that a high SSHA was unfavourable for squid aggregation. An abrupt increase in SSHA implied an increase in heat content and depth of the top of the thermocline or nutricline. This oceanographic process would reduce the enhancement of primary productivity, leading to the low availability of food for squid stocks (Ichii et al., 2011).

The expression PAR designates the spectral range of solar light that is useful to terrestrial plants in the process of photosynthesis (Meek et al., 1984). In the oceans, PAR is one of the most important factors driving the variations in ocean productivity (Wu et al., 2011), implying that PAR may have potential influence on the abundance and spatial distribution of squid species. Sanchez et al. (2008) related PAR to the potential habitat of the squid Loligo vulgaris in the northwestern Mediterranean using a geographical information system and generalized additive model. They found that PAR was a good predictor for detecting the habitat of L. vulgaris and the optimal PAR for L. vulgaris was high (43.9–52.9 E/m²/day). For O. bartramii in the northwest Pacific Ocean, a strong association was found between PAR and CPUE on the fishing grounds. High PAR could create high NPP in the habitat (Figures 8 and 9), implying enhanced food availability and good feeding conditions for O. bartramii. Therefore, high PAR is a benefit in generating a suitable habitat for O. bartramii and yields high squid abundance. It was further observed that the LATG migrated with the most preferred PAR because of the positive correlation between them, suggesting that PAR also played an important role in controlling the distribution of O. bartramii.

The far-reaching global impacts of El Niño and La Niña events have been recognized on terrestrial and marine resources (Holmgren et al., 2001). Fish responses to this mode of climate variability are likely to be species-specific and region-specific. For example, Zainuddin et al. (2004) examined the fishing ground distribution of albacore tuna (Thunnus alalunga) in the northwestern North Pacific. They clarified that La Niña events typically yielded an optimum combination of chlorophyll and temperature fronts, creating more productive fishing grounds. With regard to O. bartramii, Alabia, Saitoh, et al. (2016) revealed the potential impacts of ENSO-dominated environmental changes to the autumn cohort of O. bartramii in the central North Pacific. The Central Pacific (CP) El Niño/La Niña conditions in the succeeding summers would result in a dramatic reduction/enhancement of available habitats. However, the weaker and short-lived eastern Pacific (EP) El Niño in the autumn–winter periods would elevate potential habitats. It is concluded that large-scale changes occur in the oceanographic features of the spawning and feeding grounds of pelagic fish due to El Niño and La Niña conditions.

Outputs produced by the HSI modelling indicated substantial monthly and interannual variability in the proportion and distribution of suitable habitats for O. bartramii on the fishing ground between 2006 and 2015. The annually declining CPUE could be explained by the expansion of poor habitats and the contraction of suitable and optimal habitats from year to year. The yearly south-approaching LATG coincided with a southward migration pattern of the latitude covering the HSI area in the range of 0.9–1.0 spatially and temporally, suggesting that suitable habitats largely determine the O. bartramii distributions. Moreover, our models indicated that El Niño and La Niña events had strong influences on the winter–spring O. bartramii cohort in this study. The squid habitat changes were at least partially driven by the ENSO conditions. Compared to El Niño events, La Niña events tended to yield higher CPUE and enlarged suitable habitat areas for O. bartramii, and the LATG moved further north. The difference in abundance and distribution of O. bartramii under the two climate conditions were largely attributed to climate-induced massive environmental variations on the fishing ground. We combined all these changes together to infer a possible mechanism for the effects of El Niño and La Niña events on squid abundance and distribution: the El Niño events induce extensive areas of cool waters on the feeding ground. And worse, high SSHA and low PAR also cause more stratified upper mixed layer and inefficient photosynthesis. Thus, massive nutrient-depleted surface waters are advected into the feeding ground, which probably decreases the suitability of O. bartramii habitat and increases poor habitats. All of these effects would lead to declining squid abundance. However, the most preferred SST and PAR for O. bartramii would move southward, and squid would be likely to shift their distribution to southern regions due to their ability to track these favourable environmental conditions. Inversely, all of the changes are opposite under La Niña events.

In fact, the El Niño event is much more complex. There are two types of El Niño in the Pacific region. One is the CP El Niño and the other is the EP El Niño (Kug et al., 2009). These two different types of El Niño have various impacts on the same fish across the Pacific, such as the skipjack K. pelamis (Yen et al., 2017) and autumn O. bartramii (Alabia, Saitoh, et al., 2016). In this study, the impacts of different types of El Niño were not clear. Meanwhile, with only three warm and cold ENSO events in our study period, the ability of this study to evaluate ENSO impacts is highly limited by the small sample size and large degree of ENSO diversity. Other than the short-term ENSO event, more studies should be conducted to evaluate the effects of long-term climate change such as the Pacific Decadal Oscillation and state the differences between them.

In conclusion, combining the squid-jigging fisheries data with three crucial environmental variables, an integrated HSI model was developed to predict a suitable habitat for the western winter–spring O. bartramii cohort in the northwest Pacific Ocean. The relationship between habitat-driven trends and variations in abundance and distribution were revealed based on the modelling framework. To improve the understanding of the impacts of climate change on species abundance and distribution, and for better fisheries management, we advocate a wider application of this approach for other managed fish stocks.

Funding

This study was financially supported by the National Science Foundation of China (NSFC 41876141); the open fund of State Key Laboratory of Satellite Ocean Environment Dynamics, Second Institute of Oceanography (QNHX1818); and the Shanghai Universities First-Class Disciplines Project (Fisheries A).