Spatial Habitat Shifts of Oceanic Cephalopod (Ommastrephes bartramii) in Oscillating Climate

Short- and long-term climate oscillations impact seascapes, and hence, marine ecosystem structure and dynamics. Here, we explored the spatio-temporal patterns of potential squid habitat in the western and central North Pacific across inter-decadal climate transitions, coincident with periods of persistent warming and cooling. Potential habitat distributions of Ommastrephes bartramii were derived from the outputs of multi-ensemble species distribution models, developed using the most influential environmental factors to squid distribution and occurrence data. Our analyses captured the underlying temporal trends in potential squid habitat in response to environmental changes transpiring at each climatic transition, regulated by phase shifts in Pacific decadal oscillation (PDO) from 1999–2013. The spatial differences in environmental conditions were apparent across transitions and presumably modulate the local changes in suitable squid habitat over time. Specifically, during a cold to warm PDO shift, decreases in the summer potential habitat (mean rate ± standard deviation: −0.04 ± 0.02 habitat suitability index (HSI)/yr) were observed along the southern edge of the subarctic frontal zone (162°E–172°W). Coincidentally, this area also exhibits a warming trend (mean temporal trend: 0.06 ± 0.21 °C/yr), accompanied with the prevalence of cold-core mesoscale eddies, west of the dateline (mean temporal trend in sea surface height: −0.19 ± 1.05 cm/yr). These conditions potentially generate less favorable foraging habitat for squid. However, a warm-to-cold PDO transition underpins a northward shift of suitable habitat and an eastward shift of regions exhibiting the highest rate of potential squid habitat loss (170–160°W; mean temporal trend: −0.05 ± 0.03 HSI/yr). Nonetheless, the emergence of the areas with increasingly suitable habitat regardless of climate transitions suggests the ecological importance of these regions as potential squid habitat hotspots and climatic refugia.


Introduction
Neon flying squid (Ommastrephes bartramii) is a large epipelagic cephalopod widely distributed across the subtropical and temperate oceans [1]. It is among the internationally (e.g. Japan, China, South Korea and Taiwan) exploited fishery resources of high commercial value in the North Pacific [2,3]. It has an averaged capture production amounting to 19,245 tonnes from 2001-2010 and is one of the major squid species comprising the world capture production of cephalopods for this decade [4]. Multinational industrial fishery started when the catches of the Japanese flying squid (Todarodes pacificus) experienced a dramatic downturn in the 1970s and commercial harvests of the neon flying squid in the North Pacific waters proceeded since then [3].
O. bartramii population in the North Pacific Ocean encompasses two seasonal cohorts: winterspring and autumn spawning groups [5]. Their spatial and temporal distributions are linked with the basin-wide oceanic circulation and their life history stages are highly responsive to the variability in ambient oceanographic regimes and concomitant changes in epipelagic environment [6][7][8]. Furthermore, both squid cohorts undergo extensive seasonal and meridional migrations between their spawning (subtropical region; 20-30°N) and foraging sites (Kuroshio-Oyashio transition and subarctic frontal zones; 30-50°N), within their brief 1-year life-span ( Figure 1A) [5,9,10]. Concurrent with its summer feeding migration, industrial Japanese fishing activities primarily occur in the offshore waters of the western and central North Pacific ( Figure 1B). In an ecological context, O. bartramii adopts significant roles in the epipelagic ecosystems owing to its mid-trophic position in oceanic food webs. Previous work on the feeding habits and diet of neon flying squid revealed that it is an opportunistic/generalist feeder that preys on crustaceans, small mesopelagic fishes and other cephalopods. It is also a major prey item for larger predators and is, hence, thought to create vital mid-level trophic links between animals of tertiary trophic level and apex predators [11][12][13]. Further highlighted are the major oceanographic (currents and frontal zones) and topographic features, along with the putative geographical extents of the annual migration between the spawning and foraging grounds for the autumn cohort of the neon flying squid [7]. (B) Spatial distributions of Japanese squid fishing effort (expressed as the total number of fishing vessel per 0.25° geographic grid) in summers between May and July 1999-2013. .
A growing number of research points to the significance of looking into regional and short-to mid-term patterns in species' habitat distributions, as interactions of environmental gradients vary across multiple spatio-temporal scales [14][15][16]. In particular, habitat distribution patterns in marine ecosystems are largely dictated by the species' sensitivity to diverse environmental and climatic cues.
As cephalopods are generally short-lived taxa, they are known to swiftly respond to drastic climate and environmental perturbations, thus, making them effective biological sentinel of climatic changes [17,18]. In light of the rapidly changing climate, sentinel species can provide valuable information on ecosystem function, distinguish ocean and human health-related risks and forecast future changes [19][20][21]. Taking these characteristics into account and building upon the results of our foregoing work [6,22,23], we hereby extended our present analyses of squid-environment interactions to implement a simple yet ecologically pertinent approach to quantify the changes in potential squid habitat in the North Pacific under the recent thermal regimes, marked by distinct phase transitions of the Pacific decadal oscillation (PDO) between 1999 and 2013. Here, we used the outputs from the multi-model habitat ensemble to infer the potential squid habitat over a 15-year period and quantify the rates of habitat change across different PDO-governed thermal regimes. In doing so, we were able to evaluate the spatial changes of potential squid habitat distributions in response to the decadal shifts in the dominant year-round pattern of sea surface temperature (SST) variability in the North Pacific Ocean.
Furthermore, our analyses, albeit on a shorter time-scale relative to earlier studies on climate refugia [24,25], provided some insights on potential sites that could act as open-ocean refugia under changing environmental conditions in the western and central North Pacific Ocean. In the literature, climate refugia are increasingly recognized for their importance in ecology and conservation studies, as they facilitate the persistence of biodiversity under changing climates [24,25]. While this concept has been broadly investigated in terrestrial ecosystems, identification of these sites in aquatic and marine systems could assist conservation efforts and spatial planning in a changing ocean climate. In the global ocean, climate refugia are generally characterized by areas where some or all variables influenced by climate variability exhibit the least changes and often overlap with oceanographic features that preclude climate change exposure [26]. Examples of these oceanographic features are seamounts and persistent frontal structures that create retention areas and habitat hotspots for many marine species [27][28][29]. In the present study, potential local refugia, referred to sites where the inferred squid habitat exhibited increasing trends regardless of the PDO phase transitions.
Overall, this study seeks to examine the spatial patterns of temporal trends in potential squid habitat in the North Pacific over the 15-year period, punctuated by persistent warm and cold phases of PDO-associated oceanographic regimes. Specifically, we aim to (1) estimate the annual changes in habitat suitability for the neon flying squid from late spring to summer (May-July) in the western and central North Pacific, across contrasting PDO phase transitions, and (2) identify climaticallysensitive zones to shifting contemporary patterns of oceanographic conditions and suitable climatic squid habitat over time.

Pacific Decadal Oscillation (PDO)-Based Climatic Phase Shifts
During the duration of study, climate transitions were identified based on the patterns of Pacific decadal oscillation (PDO) from January to December, 1999-2013 ( Figure 2A). The PDO was initially introduced by [30] as the leading empirical orthogonal function (EOF) of the North Pacific (20°-70°N) SST monthly averaged anomalies (SSTAs). The anomalies were computed as the departures from the climatological annual cycle after the removal of global mean SSTs. Since its introduction, PDO has been linked to the other components of the climate system [31] and to its profound impacts on marine ecosystems and fisheries [30,32,33]. Here, monthly PDO data were downloaded from the Joint Institute for the Study of Atmosphere and Ocean (JISAO) (http://research.jisao.washington.edu/pdo/PDO.latest.txt, date accessed: 20 June 2017). The monthly time-series data were annually-averaged and plotted to visually identify PDO patterns. The annual averages of PDO index were used for the analyses as the neon flying squid has a one-year life-span [5,11], and thus, this temporal resolution better captures the climate conditions experienced by the species throughout its lifetime. Similarly, monthly sea surface temperature (SST) anomaly data computed using a 30-year base period  were from the International Research Institute for Climate and Society (IRI) (iridl.ldeo.columbia.edu/SOURCES/.NOAA/.NCEP/.EMC/.CMB/.GLOBAL/.Reyn_SmithOIv2/.mont hly/.ssta, date accessed: 20 June 2017). The data were used to examine the corresponding spatial SSTA patterns at each respective PDO phase. Due to the long-lived and year-long SST-based climate variability pattern associated with PDO signals [31,34], we were able to capture two PDO phase transitions within our study period (Figure 2A). Each of the periods within a given transition was represented by a basin-wide warming or cooling ( Figure 2B

Analytical Framework
Our analyses were implemented following the framework proposed in Figure 3. Here, multimodel habitat forecasts were generated using the monthly-pooled (May-July) information of squid occurrences and set of most influential environmental factors (i.e. sea surface temperature, SST and sea surface height, SSH) to squid habitat distribution. These variables are related to the thermal and mesoscale environments, respectively, that significantly impact the summer feeding areas for squid and other pelagic species. Annual seasonal averages of multi-model habitat predictions (May-July, 1999-2013) were computed and temporal trends for each climatic transition between two consecutive periods of PDO-specific phase were then calculated. These calculations were also carried out for SST and SSH to examine the rates at which the thermal and mesoscale environments were changing across PDO phase transitions. Finally, comparison of spatial patterns between temporal trends for potential squid habitat and oceanographic conditions was undertaken to identify areas that are particularly sensitive to climatic changes.

Squid Occurrences and Environmental Data
The monthly-compiled squid jigging vessel positions (Longitude, Latitude) from daily fishing logs (point data) from May to July 1999-2013 (corresponding to Japanese commercial jigging months in summer, Figure 1B) were provided by the Aomori Prefectural Industrial Technology Research Center (APITRC). The environmental variables used for the habitat model construction comprised the monthly-averaged sea surface temperature (SST) computed from the daily optimally-interpolated Advanced Very High Resolution Radiometer-Optimally interpolated (AVHRR-OI) SST (ftp.podaacftp.jpl.nasa.gov/allData/ghrsst/data/L4/GLOB/NCDC/AVHRR_OI/; date accessed: 04 February 2016) and monthly-averaged sea surface height (SSH) calculated from daily-resolved altimeter measurements from the Archiving, Validation and Interpretation of Satellite Oceanographic data (AVISO) website (ftp.aviso.oceanobs.com; date accessed: 04 February 2016) between May and July 1999-2013. The selection of environmental variables used in the analyses were supported by the results of earlier studies, suggesting that these factors were most influential to squid habitat distributions [22,23,35]. All the environmental data were resampled to a common spatial resolution (0.25° × 0.25°) prior to habitat model construction and subsequent analyses. Finally, using a bicubic interpolation [36], the grid at each fishing point location was also interpolated (0.25° × 0.25°) to extract the corresponding environmental values at the given fishing position.

Multi-Model Ensemble of Potential Squid Habitat
Using a set of environmental factors and squid occurrences, a total of 10 single-algorithm models available within the BIOMOD2 r package [37] were trained and fitted to predict the species habitat distributions in the study area. BIOMOD2 invokes regression and machine-learning based presenceonly models to compute for the habitat suitability index (HSI), where values close to or equal to 1 suggest areas of favorable species habitat. The parametrizations and settings used for the development of the single-algorithm models are summarized in Table 1.
Since the pixels with no squid jigging activity cannot be classified as true absences, a total of 500 pseudo-absences for each squid jigging observation compiled daily were randomly selected across all environmental layers. These data were then pooled together across each period and final pseudoabsences used for subsequent single-algorithm model runs were randomly subsampled from the pooled data using a 1:10 presence to pseudo-absence ratio [38]. The random sub-sampling of pseudoabsences was carried out to make model simulations computationally feasible yet within a moderate sampling prevalence [39]. The subsequent selection of the single-algorithm models used to generate a weighted mean multi-model ensemble was based on an arbitrary threshold for the true skill statistic (TSS). As there is no known consensus yet as to the most appropriate evaluation metric for ensemble forecasts [40], we used the TSS for model selection (TSS ≥0.80; Table 2). Unlike the other performance metrics (area under the curve (AUC) and Cohen's Kappa, hereafter Kappa), TSS is independent of the model prevalence and thus, provides a more robust measure [41]. Here, the weights for ensemble model were proportional to the computed TSS of best-performing models [37]. Finally, the compiled occurrences at each given period were randomly divided into training (70%) and testing (30%) sets, for fitting and evaluating the models, respectively.

Computing for Temporal Trends at each Climatic Transition
The temporal trend for potential squid habitat at each PDO phase transition between successive periods Figure 2), correspond to the slope of the cell-wise linear regressions between HSI and year, as response and independent variables, respectively. This calculation was also performed for environmental parameters (SST and SSH) against year at each given phase transition. Doing so enables us to quantify the magnitude and direction of habitat and environmental changes across space under different PDO phase transitions. All the analyses were performed in R version 3.6.0 [42] and all the mapping routines were carried out using general mapping tools (GMT) version 4.5.15 [36].
Moreover, to identify potential local refugia, we classified the temporal trends for HSI, SST and SSH in areas with the highest predicted squid habitat (165°E-160°W and 39-46°N) into equal quintiles, resulting into five categories of habitat and oceanographic changes at each transition (Table  3): large decrease, small decrease, largely unchanged, small increase, and large increase [26]. Subsequently, we mapped and examined the spatial overlap of small and large increases in HSI with, classified maps of SST and SSH at each PDO phase transition.

Squid Environmental Preferences
Preferred environmental ranges defining potential squid habitat were obtained from the multimodel simulations (from the best performing single-algorithm models; Table 2) and response curves for each of the environmental factors are shown in Figure 4A-B. The thermal conditions and mesoscale environment of squid habitat (HSI ≥ 0.5) were situated within the SST ( Figure 4A) and SSH ( Figure 4B) ranges from 12-17 °C and 40-120 cm, respectively. Period-specific spatial extents of optimal isotherms ( Figure 4C-E) showed northward advance under the cold PDO phases (P1 and P3; Figure 4C,E) relative to its alternate phase (P2; Figure 4D). In contrast, the optimal isolines for SSH showed minimal differences in their geographic extent across PDO periods ( Figure 4F-H). Nonetheless, the total area of pixels within the optimal SST and SSH limits was higher in the cold than in warm PDO phases (Table 4).

Spatial Patterns of Period-Specific Squid Habitat Forecasts
Potential squid habitat predictions obtained from the multi-model habitat ensemble highlighted differences in spatial distribution across PDO periods ( Figure 5A-C, left panels). In particular, periods of cold PDO phases covered slightly broader, albeit patchy zones ( Figure 5A,C left panels) of potential squid habitat relative to those predicted during the warm PDO phase ( Figure 5B, left panel). Moreover, the zonal and meridional HSI frequencies ( Figure 5A-C, right panels) further showed a northward advance of habitat suitability frequency peaks by as much as 1°N (~110 km) between the cold (P1 and P3) and warm (P2) PDO phases. Zonal pattern-wise, the frequency peak of habitat suitability also revealed a westward displacement in the latter period, concomitant with a spatial reduction in the total potential squid habitat area (Table 5). Overlain on the spatial maps are contours delineating the regions with moderate to high squid habitat suitability (HSI ≥ 0.5). Latitudinal bin corresponding to the HSI frequency peaks for each period were also labeled (values in brackets). Table 5. Summary statistics of the spatial extent of inferred potential squid habitat across the different PDO periods. The total area with moderate to high potential squid habitat (HSI ≥ 0.5), overall habitat suitability (HSI) average and standard deviation (SD) were computed for each PDO-specific periods from 1999-2013.

Temporal Trends in Squid Habitat and Environmental Conditions
Spatial distributions of temporal trends computed across contrasting PDO phase transitions for potential squid habitat and oceanographic conditions are shown in Figure 6A-C. Salient changes in trend patterns of potential squid habitat suitability emerged between the two phase transitions ( Figure 6A). Highest rates of habitat loss (negative trend) and gain (positive trend) were located at the southern edge and northern part of the subarctic frontal zone (SAFZ, 40-45°N), respectively. In particular, during the first transition (T1; 1999-2006), positive trends in squid habitat suitability (0.03 ± 0.02 HSI/yr; mean annual rate ± 1 standard deviation) were observed along the northern part between 172°E-170°W and southeastern edge between 170-160°W of the SAFZ ( Figure 6A, left panel). However, squid habitat losses were particularly pronounced at the southern border of the SAFZ between 162°E and 172°W (−0.04 ± 0.02 HSI/yr), further coinciding with the highest spatial density of the squid-fishing effort. During the second transition (T2; 2002-2013), the largest signal of habitat decline (−0.05 ± 0.03 HSI/yr) shifted east between 170° and 160°W ( Figure 6A, right panel), with an average trend significantly higher than T1 (p < 0.001; unpaired t-test). Areas with positive HSI trend expanded east of the SAFZ, generating a continuous zone of squid habitat accretion. While the zonal extent of high squid fishing effort density remained considerably constant between the two transitions, squid fishing activity during T2 exhibited meridional expansion towards the northern fringe of the SAFZ (43-45°N). This concomitant spatial change in squid-fishing pattern suggested that squid jigging fleets were tracking the captured north-directed shifts in potential squid habitat quite reasonably during this climate phase transition. Figure 6. Spatial maps of the temporal trends computed for (A) potential squid habitat, (B) thermal environment, and (C) mesoscale eddy field, during the PDO phase transitions from a cold to warm (T1, left panels) and warm to cold (T2, right panels) phase. The spatial signals of the decrease and increase in the climatic habitat suitability, temperature and mesoscale activity for a given phase transition are represented by blue and red shades, respectively. The solid contours correspond to the highest spatial density of squid fishing effort for each climatic transition, gray polygons represent the topographic features and the region within the broken lines correspond to the spatial domain of the subarctic frontal zone (SAFZ).
Temporal trends in thermal and mesoscale environment similarly showed pronounced meridional and zonal patterns ( Figure 6B-C). For instance, the latitudinal patterns in temporal trend for SST showed an extensive cooling signal (−0.07 ± 0.17 °C/yr) during the first transition (T1), further flanking a warming water parcel situated between 170°W and 170°E ( Figure 6B, left panel). During the second transition (T2), however, thermal trend reversed toward a general warming pattern (0.09 ± 0.09 °C/yr), with highest warming rate east of the basin (175°-160°W; Figure 6B, right panel). Moreover, temporal SSH trends also revealed characteristic spatial patterns, capturing the ambient mesoscale environment at each climatic transition. For instance, meridional patterns showed consistent elevated temporal SSH trends between 175°E-160°W and within the latitudinal boundaries of SAFZ (40-45°N) across all transitions ( Figure 6C, left and right panels). Overall, the mesoscale environment in the west and central North Pacific (140°E-160°W and 30-50°N) was characterized by the distinct overturn from cold-core (−0.39 ± 2.00 cm/yr) to warm-core (0.22 ± 1.08 cm/yr) dominated eddy field during the first (T1) and second (T2) transition, respectively.
Furthermore, the examination of temporal trends in potential squid habitat and oceanographic environment across PDO phase transitions showed salient spatial patterns, possibly regulating the transition-specific habitat responses. In general, the zonal habitat shifts showed differences in relative importance of the thermal and mesoscale forcing under contrasting PDO transitions. In particular, changes in habitat trends in the central (162°E-172°W; 40-44°N) North Pacific were likely modulated by transition-specific changes in both mesoscale and thermal conditions. Contrariwise, in the eastern North Pacific (170-160°W; 40-43°N), habitat trend changes appears to be responding more to thermal cues than mesoscale forcing. However, areas of increasing habitat irrespective of climatic transition have also been identified close to topographic features (i.e. emperor seamount chain) and northern part of the quasi-permanent temperate frontal zone (170°E-160°W; 43-45°N; SAFZ).
Further classification of the spatial overlap of regions classified with small and large increases in potential squid habitat with categorized changes in temporal trends of SST ( Figure 7A) and SSH ( Figure 7B) across phase transitions, showed that, indeed, areas sitting atop or close to the Emperor seamounts exhibited persistent habitat increases regardless of the PDO transition. Similarly, the patches of squid habitat increases proximal to the northern part of SAFZ (175°E-170°W; 43-45°N) were characterized by small increases to largely unchanged trends in thermal environment ( Figure  7A) and mesoscale eddy field ( Figure 7B) across the PDO phase transitions (left and right panels).

Discussion
While there are quite a vast number of studies looking at the patterns of suitable habitats for neon flying squid using the species distribution modeling approach [22,[43][44][45], the framework of analyses implemented in the present study provided both qualitative and quantitative means of elucidating the effects of the dominant modes of climatic variability in the North Pacific Ocean at a decadal time-scale. Our results underpinned the differences in the temporal trends for suitable squid habitat across climatic transitions and, more importantly, afforded a quantitative estimate of the spatial changes in the squid habitat suitability over time. This allowed the characterization of sites that were particularly exposed and vulnerable to climatic oscillations while highlighting climateresilient zones, exhibiting minimal changes in physical oceanographic properties over time.
Previous studies have described the spatial and temporal patterns of potential habitat of neon flying squid in response to dominant modes of high-frequency climate variability such as the El Niño Southern Oscillation (ENSO) [6,46,47]. Moreover, the impacts of the basin-wide regime shifts based on PDO patterns on relative abundances of autumn [7,8] and winter [35] cohorts of neon flying squid in the northwestern Pacific have been previously explored. Our analyses, however, focused on temporal habitat and environmental gradients in summer foraging zones of the autumn cohort (central Pacific; 160°W-160°E) under the contemporary PDO phase transitions, hence, affording insights into its distributional response to shifting climate and setting this work apart from previous studies. While the cold PDO phases (P1 and P3) revealed almost similar spatial extents of potential feeding habitats relative to the warm phase (P2), the averaged magnitude of squid habitat suitability was significantly higher in the latter (Table 3; p < 0.001; unpaired t-test). Such spatial and temporal patterns in potential squid habitat were likely modulated by the contrast in temporal trends in thermal and mesoscale features across a given PDO transition. Our previous analyses have shown a strong association of squid to warm-core eddies as important foraging features [22,23]. Henceforth, the prevalence of cold-core mesoscale circulation during T1 ( Figure 6C, left panel) potentially explains the decline in the squid habitat around the south margin of the SAFZ. Similarly, our results reinforced earlier findings on the important roles of bathymetric and non-transient frontal structures in creating local, yet persistent, zones of climate-resilient foraging grounds for migratory taxa in the western and central North Pacific Ocean [23,48]. For instance, in oceanic environments where climatic parameters often exhibit minimal variability over wide spatial extents, the presence of topographic features can decouple climatic processes at local scales from basin-wide conditions [26,29]. In fact, this appears to be the case in our study area, as shown by an increasing squid habitat suitability trends (habitat gain) in waters proximal to seafloor features and temperate frontal zone, regardless of the climate transition ( Figure 6A, left and right panels; Figure 7A-B). This suggests the potential of these sites as local refugia for squid and other pelagic resources [49,50] under changing environmental and climatic conditions in the open waters of the North Pacific Ocean.
In contrast, the climate-driven spatial potential habitat shifts across PDO phase transitions were particularly high between 39° and 47°N, with notable temporal trend reversal along the east-west direction. These results highlight and reassert the spatially-complex squid habitat responses not only to high-frequency modes of climate oscillation, but to low-frequency climate variability in the central and western North Pacific as well. In general, PDO impacts the basin-wide circulation patterns in the North Pacific basin, and ensuing changes in ocean currents have consequences on the heat, material and nutrient transport [51] which, in turn, could modulate the rates of biological production [52]. Further, zonal habitat shifts in the area within 162°E-172°W and 40-44°N suggest that these regions are most exposed to climate-driven changes in both thermal and mesoscale environments; thus, are potentially vulnerable to these impacts. The circulation dynamics in the transition zone (140°E-180; 35-40°N) between major frontal zones (i.e. sub-tropical and subarctic fronts) are highly variable due to the convergence of cold (Oyashio) and warm (Kuroshio) boundary currents [53,54]. The associated mesoscale eddies and flow instabilities that propagate northwestward of 40°N generate productive, albeit ephemeral feeding zones for numerous pelagic species [55][56][57]. In this region, the increasing trends in squid habitat during a mesoscale environment dominated by warm-core eddies (T2) suggest the higher affinity of squid to these features relative to the otherwise more productive cold-core rings. A similar response was recently observed from the tagging experiment of pelagic shark species in the Gulf Stream region [58]. The results from this work suggest that sharks use the core of anticyclonic ocean eddies as channels to forage at depths otherwise, unreachable due to the thermal physiological constraints. This could explain the reported affinity of neon flying squid (this study) and tuna species [55] to warm-core eddies in the pelagic waters of the Kuroshio region.
Moreover, in the temperate frontal zone (SAFZ; 170°E-160°W and 43-45°N), a region of increasing habitat suitability concurrent with warming waters were observed across all PDO phase transitions. In this area, the transition zone chlorophyll-a (TZCF) front were shown to assume its northernmost position during summer within our study region [6,23]. The TCZF is a permanent gradient in the sea surface color, persisting zonally across the North Pacific basin, with high and low chlorophyll to the north and south of the front, respectively [59,60]. Studies have suggested that this frontal feature could be important as a forage and migration habitat for pelagic species [61]. In the present analyses, periods of the extensive and sustained warming in these areas across the PDO phase transitions were likely to augment the potential squid habitat. This is particularly so, as the thermal conditions could become physiologically favorable (i.e. within the species' thermal tolerance), allowing the squid to move further north and closer to the boundary of the TCZF. This similar distributional response of neon flying squid to elevated temperature in the central North Pacific was also noted earlier [6,47].
Finally, the analyses implemented in this work offered, to some extent, insights that substantiate the growing consensus on the significant effects of environmental and climatic variability on species habitat distributions across the marine and aquatic systems [62][63][64][65]. These impacts have important consequences, in terms of ecosystem functioning and provision of its services [66]. Furthermore, the rapidly accelerating pace of climate change warrants the exploration of means to monitor its ecological effects (e.g. using biological sentinels) [67] and consequently to create conservation and adaptation strategies (e.g. establishment of climate-smart protected areas and identification of potential climate refugia) [26,68]. Here, we have explored these aspects briefly, nonetheless, it is worth noting that our present analyses are subject to potential caveats and highlight areas of future research that can be done to reinforce our current findings. These caveats stem largely from the limited temporal span of available data (i.e. species occurrence records) that allow the robust detection of long-term trends in climate oscillations (basin-wide climate teleconnections, e.g. PDO) and resulting climate-driven ecological impacts. Thus, future work should aim at addressing these key issues through the incorporation of more recent data when made available.

Conclusions
Climate is a major driver in the spatial redistribution of marine species through its effects on the environment and on species' optimal habitat. In the North Pacific, apart from high-frequency modes of climate oscillations, the low-frequency decadal-scale climate variability such as PDO significantly impacts the summer pelagic habitat of large oceanic squid. The effects ensued from the basin-wide changes in the ocean temperature and circulation, which in turn alters the foraging conditions for squid and translates to spatial shifts in the distribution of its potential habitat. Furthermore, the squid's distributional responses exhibited zonal and meridional patterns, highlighting the relative sensitivity and exposure of certain areas to climate-modulated changes. The emergence of squid habitat patches and persistent zones of increasing habitat suitability across the PDO phase transitions ( Figure 7A-B) underpins the importance of these sites as potential local climate refugia under anticipated threats of abruptly transpiring climate-driven environmental changes.