Climate change is modifying oceanic biological response patterns through a variety of ways, such as changing habitats’ physical and chemical properties, and amplifying environmental fluctuations1,2,3,4. Thereby, it has been established as an important driver of changes to marine fisheries through altering primary production, food web interactions, life history, and distribution5. Myriad studies have examined the variety of species’ distributional responses to climate change and found that temperature tends to be among the largest limiting habitat factor for most poikilothermic organisms2,3. Thus, global warming is predicted on balance to increase fisheries catch potential in higher latitudes while decreasing in tropical regions due to the poleward shift in the geographical distribution of fish stocks1,6. Because of the unique biological characteristics of short-lived, semelparous, fast growing, highly fecund, gluttonous, and metabolically efficient Ommastrepes squids7,8, it is unclear whether or to what extent ocean warming will result in distributional changes.
Here we fill gaps in our knowledge regarding how squids respond to altered marine environmental variables driven by climate change in the Pacific Ocean. We used 1) Chinese squid-jigging fisheries data covering three stocks of two species (Ommastrephes bartramii in the Northwest Pacific Ocean, Dosidicus gigas in the offshore regions of Peru and Chile, and Dosidicus gigas in the equatorial waters of the Southeast Pacific Ocean, Fig. 1a, digitized by the National Data Center for Distant-water fisheries of China in Shanghai Ocean University; 2) marine environmental data, including sea surface temperature (SST), chlorophyll a (Chl), and sea surface salinity (SSS) downloaded from NOAA OceanWatch (https://oceanwatch.pifsc.noaa.gov/), and 3) Oceanic El Niño indices (ONI), which characterize ENSO events in the Pacific Ocean, downloaded from Climate Prediction Center (https://www.cpc.ncep.noaa.gov/data/indices/sstoi.indices) spanning from 2005 – 2018, to quantify the effects of climate change on the squids’ distribution. The Ocean Niño Index (ONI) values were divided into five levels (fONIs) including “Strong La Niña”, “Moderate La Niña”, “Normal”, “Moderate El Niño”, “Strong El Niño”20 During this period, strong La Niña and El Niño events both occurred twice (Fig. S1). A simulation study showed that the monthly aggregated catch in 0.5° longitude × 0.5° latitude cell can serve as a robust local abundance index to represent the dynamics of spatial distribution for squids9, thus, the original fisheries datasets were grouped into monthly 0.5° × 0.5° cells, and the environmental data were converted into 0.5° × 0.5° for each month to correspond to the spatial resolution of the fishery data. Finally, we modelled 6048 (Fig. 1B), 3334 (Fig. 1C), and 8878 (Fig. 1D) fished cells covering the three stocks across ten variables ‘Year’, ‘Month’, ‘Longitude’, ‘Latitude’, ‘SST’, ‘SSS’, ‘Chl’, ‘ONI’, ‘fONI’, ‘Effort’, and ‘Catch’ for each record.
Quantifying the preferred habitat distribution of squids
We developed a set of generalized additive mixed models (GAMMS) to quantify the preferred habitat distributions and to check whether poleward shifts in habitat distribution track climatic changes. Because squid is continuous (catch weight in metric tons), asymmetrical, with variances generally increasing with higher catches (Fig. S2), we used GAMMs with a gamma distribution and log-link. The models included Catch as response variable, Month, Longitude, Latitude, SST, SSS, Chl as fixed random effects, and fONI as random effect for each stock, respectively. Model structure was tested following a hierarchical approach with five types of smoother for each fixed variable21. The approach enabled us to examine factors that influenced the preferred habitat distribution of squids. We fitted the models using the ‘mgcv’ package version 1.810 in the R software environment. We determined the best-fitting model using AIC values (Table S1). The preferred habitat distribution during 2005 – 2018 were hindcasted by the final model for each squid (see Supplementary materials).
Generally, the most obvious pattern is that La Niña events (such as, 2010 and 2011) tend to increase suitable habitat compared to El Niño events (such as 2015) (Fig. S3 – S5) for the three stocks. However, annual variability and habitat preferences varied among stocks. For two stocks of D. gigas in the Southeast Pacific Ocean, we found sinusoidal annual distribution patterns: the suitable areas increase from April to November, and then decrease from December to the following April, whereas the unsuitable habitat areas have the opposite seasonal changes (Fig. S4, S5).
Environmental response curves of squids
Environmental response curves were created by plotting the scaled catch against the predictor of interest while holding all other predictors at their median values based on the distribution models for each ONI level and squid11. Most response curves were bell-shaped, or appeared as partially bell-shaped, limited by sampling data contrast. We found that the environmental response curves are both species and stock-dependent. Despite incomplete environmental sampling overlap, SST and Chl a response curve shapes were similar across D. gigas stocks (Fig. 2). Moreover, the estimated environmental response curves differed significantly across climate regimes. Scaled abundance indices for strong La Niña years were typically 2-3 fold higher than strong El Niño events across all stocks, and all environmental variables. Interestingly, O. bartramii response curves were most clustered (similar) across regimes, but increased dramatically in strong La Niña years.
Evaluating the dynamics of habitat for squids under climate change
The contour of annual suitable habitat area (SUA, the area with scaled catch > 0.6) was depicted to reflect the long-term dynamics in fishing grounds (Fig. 1E, 1F and 1G). The annual gravitational center positions (LonG – longitude of gravity, LatG – latitude of the gravity) were calculated12, then we fit linear models with LatG as the response variable and LonG as explanatory variable to observe whether a poleward shift occurred (Fig. 1H, 1I and 1J). The spatial response maps were made (similar to environmental response curves by changing longitude and latitude, simultaneously) to reflect basic distributional patterns under the five ONI levels (Fig. 1K, 1L and 1M).
The results show that O. bartramii annual SUAs are getting smaller and moving northwestward (Fig. 1E), the annual SUA are getting bigger with moving southeast for D. gigas in the equatorial waters (Fig. 1F), and the annual SUA are getting smaller with moving southwest for D. gigas in the offshore Peru and Chile (Fig. 1G). Similarly, the significant latitudinal moving occurred on three squid stocks in the Pacific Ocean (p < 0.05, Fig. 1H, 1I, and 1J). The spatial response maps demonstrate the preferred habitat distribution for each stock clearly, and La Niña event tend to produce more suitable habitat than El Niño event (Fig. 1K, 1L and 1M).