Introduction

Humanity’s influence on the planet can lead to an extensive future scenarios of climate change (Moss et al. 2010). Depending on the projected future concentrations of greenhouse gases and consequent projected radiative forcing trajectories (in W/m2; Moss et al. 2008), these can be split between four representative concentration pathway (RCP) scenarios (Moss et al. 2010). One pathway predicts a rising of emissions until 2100 (RCP 8.5; Riahi et al. 2011), one which emissions peak before 2100 and then decline (RCP 2.6; Van Vuuren et al. 2011) and two stabilizing after 2100 without overshoot (of the 2 ℃ warming; RCP 4.5 and RCP 6.0; Masui et al. 2011; Thomson et al. 2011). Additionally, these pathways can be further split into different future timeframes: near- (2035), mid- (2050) and long term (2100 and beyond) (Lee et al. 2021). However as diverse these future scenarios are, they all predict that the physicochemical properties of the World’s Ocean will be (and are already) affected (Nazarenko et al. 2015; Fox-Kemper et al. 2021). These include changes in the surface water temperature, salinity, sea-ice extent, currents (Fox-Kemper et al. 2021), and pH (Lee et al. 2021).

Declining pH and oxygen, rising temperatures, and changes in current patterns are expected to produce major effects across marine life, both at physiological and ecological levels (Bijma et al. 2013; Sampaio et al. 2021). Rising temperatures affect both the physical conditions of habitats and physiology of organisms. Sea surface temperatures increase faster than other parts of the water column, leading to strong stratification of the upper layer of the ocean. This in turn will decrease the primary productivity of the ocean, as nutrients from the deeper parts of the water column will be trapped beneath the pycnocline (Bijma et al. 2013). In terms of physiology, organisms specialize on a range of temperatures, due to metabolic tradeoffs between the structure and function of their enzymes and cell membranes (Pörtner 2002; Bijma et al. 2013). Increase in temperatures beyond this range leads to deficiency of oxygen in blood (hypoxemia) and higher dependance on anaerobic metabolism (Bijma et al. 2013) which are exacerbated if low oxygen levels are present (Sampaio et al. 2021). However, cephalopods are a group of marine molluscs that seem to thrive in areas of the ocean where such metabolic restraints are already present (Rosa and Seibel 2010), and there is evidence that they may benefit from climate change (Coll et al. 2020).

Cephalopods fulfill an important role in numerous marine foodwebs and make up to a fourth of the global mollusk biomass (Rodhouse et al. 1996; Bar-On et al. 2018; de la Chesnais et al. 2019). Squids are particularly important for human consumption and they make up to 70% of the biomass harvested for cephalopods (Jereb and Roper 2010; Arkhipkin et al. 2015). Squid consist of the coastal and benthic spawning myopsids and oceanic oegopsids that spawn their eggs into the water column. The dominant squid families that are harvested belong to the Ommastrephidae and Loliginidae and, to a lesser extent, the Gonatidae and Thysanoteuthidae (Jereb and Roper 2010). These squid generally occur in highly productive areas (Jereb and Roper 2010) within a narrow range of temperatures and salinity. Changes in squid distribution have been observed, either sporadically due to changes in temperatures (due to el Niño/la Niña events; Field et al. 2007), seasonally due to ontogenic movements (Jereb and Roper 2010), or long-term movement of the species towards new habitats (Chen et al. 2006; Golikov et al. 2013; Arkhipkin et al. 2015).

Squid species have recently expanded their range poleward (e.g., D. gigas; Nigmatullin et al. 2001; Field et al. 2007; Zeidberg and Robison 2007; Arkhipkin et al. 2015). In the last sixty years, cephalopod species populations have also increased in abundance regionally (Doubleday et al. 2016) (e.g., L. forbesii; Chen et al. 2006). These population changes are thought to result from the adaptative capacity and phenotypic plasticity of cephalopods to cope with environmental change (Pecl and Jackson 2008; Hoving et al. 2013; Kooij et al. 2016; Liscovitch-Brauer et al. 2017; Jin et al. 2020). The overfishing of their competitors and predators has also been suggested as a cause of their population changes(Caddy and Rodhouse 1998; Field et al. 2007; Zeidberg and Robison 2007; Vecchione et al. 2009). Although the plasticity of some cephalopods allows them to cope and potentially benefit from climate change, responses are dependent on physiology, life cycle characteristics, and behavior, and hence are likely to be species specific (Rodhouse et al. 2014). Indeed, the controlled exposure to future climate conditions results in depressed metabolic rates and activity levels in squids (Rosa and Seibel 2008), along with malformations and stress in early life stages (Rosa et al. 2012) and during embryogenesis (Rosa et al. 2014). Additionally, recent studies applying species distribution models on coastal cephalopods (using Sea Surface Temperature [SST], Sea Surface Salinity, total Chlorophyll Mass Concentration at Surface, Dissolved Oxygen Concentration at Surface, and Ocean Surface pH; Boavida-Portugal et al. 2022) and the three main commercial European cephalopods (Octopus vulgaris, Sepia officinalis, and Loligo vulgaris; using Sea Bottom Temperature and SST; (Schickele et al. 2021)) suggest that they may face a decline under future climate change scenarios. Currently, we lack scientific information to properly predict the future responses of commercial squid.

To predict and project future distribution of squids, species distribution models (SDMs) have been used (Schickele et al. 2021; Boavida-Portugal et al. 2022; Borges et al. 2022). These models are based on the Hutchinsonian ecological niche definition (Hutchinson 1957), in which a niche is a “n-dimensional hypervolume” of the different abiotic and biotic variables a species lives in. In this context, SDMs employ a correlative approach to the characterization of the ecological niche, assuming that the current distribution of a species (or other target group of individuals) is a good indicator of the suitability of the underlying ecological conditions (Pearson 2007; Kearney and Porter 2009). In practice, it means that environmental data that cover the areas where the species has been observed, can be used to predict the species’ broader distribution and project the suitability of new habitats in space (e.g., project the invasive potential of species; Barbet-Massin et al. 2018; Sung et al. 2018; Borges et al. 2021) and/or time (e.g., project the potential impact of future climate change scenarios; Araújo and Rahbek 2006; Melo-Merino et al. 2020; Borges et al. 2022).

If the current increase in cephalopod populations is due to the increase of the effects of climate change, more drastic climate scenarios would mean an increase of cephalopod populations and an expansion of habitats range. To test this hypothesis, we used data collected from online repositories and a SDM framework to forecast the potential future distribution of the main 12 commercial cephalopod species (six oegopsids and six myopsids) under different climate change RCP scenarios (RCP 2.6, 4.5, 6.0 and 8.5) and timeframes (present, mid-century [year 2050] and long-term [year 2100]). With the overall habitat suitabilities from these forecasts, we can determine if cephalopods will benefit from worse climate change scenarios and expand their habitat range into a wider latitudinal amplitude.

Materials and methods

Data collection

We selected the largest decapodiform cephalopod fisheries species according to fishstatj for the year 2020 (FAO 2014) (Table 1). Species occurrence data were retrieved from Global Biodiversity Information Facility (GBIF.Org User 2021, 2022) and the Ocean Biodiversity Information System (OBIS; Grassle 2000). The search for occurrence points was limited to squid species by restricting the taxonomical search to “Decapodiformes”. Of these, Doryteuthis pealeii, Doryteuthis opalescens, Doryteuthis gahi, Loligo vulgaris, Loligo forbesii, Loligo reynaudii, Dosidicus gigas, Illex argentinus, Illex illecebrosus, Nototodarus sloanii, and Berryteuthis magister were selected for modeling (Todarodes pacificus data did not allow for proper model construction [49 occurrences that when modeled, did not cover its present distribution]). Moreover, given that only 8 occurrences were available for Doryteuthis gahi in these repositories, these were complemented with data from Xavier et al. (2016a). Each occurrence dataset was manually curated to remove any points in the wrong location (i.e., land, wrong ocean basin, outside of distribution range; Jereb and Roper 2010) (See Online Resource 1). To mitigate the effects of spatial sampling bias, a species-level thinning step was then implemented by randomly selecting points that were at least 20 km apart from each other—a distance that provided the best present forecasts when compared to the literature (Jereb and Roper 2010).

Table 1 Top 12 squid species by captured tonnage in 2020 (FAO 2014)
Full size table

Environmental layers were retrieved from the Bio-Oracle project (Tyberghein et al. 2012; Assis et al. 2018). The environmental marine variables selected pertained to temperature, salinity, current velocity and chlorophyll, retrieving environmental layers of their long-term mean, maximum and minimum values, at both the surface and the average seabed depth (except for chlorophyll, for which only surface layers were collected). Environmental layers were retrieved at a resolution of 0.008(3) degrees (latitude and longitude) for both present (long-term averages between 2000 and 2014) and all future climate scenarios (years 2050 and 2100; RCP 2.6, 4.5, 6.0 and 8.5; CMIP5). Collinearity between the environmental variables was assessed with Pearson’s method (r), removing all variables with a collinearity > 70% (Naimi and Araújo 2016; Schickele et al. 2021), as to decrease biases in inference statistics of the models (Dormann et al. 2013). With this step, the long-term mean values were kept for all environmental samples, with seabed salinity completely removed (highly correlated with its surface counterparts). Additionally, long-term minimum chlorophyll was kept.

Species distribution models

We used the MaxEnt algorithm (Phillips and Dudík 2008) for building the species distribution models. For each species, a total of 5 model runs were executed using all data (presence points plus 1 000 [randomly selected] background points; see total N of presence points in Table 2) randomly partitioned and shuffled between calibration (75%) and validation (25%) datasets for each run. Each model, in turn, was evaluated using the True Skill Statistic (TSS; Allouche et al. 2006), and models performing poorly (TSS < 70%) were discarded. To compute the ensemble of the remaining models for each species, scenario, and period, weighted mean consensus forecasts were used (Araújo and New 2007), by weighting each of the remaining model forecasts with the square difference between its TSS score and 50% (the same computation was used to calculate the ensemble model TSS; see final N of model runs in Table 2). The previous step was performed for each species across all the climate scenarios and timeframes, producing a total of 99 Habitat suitability forecasts (11 species × [4 scenarios × 2 future timeframes + present]). These forecasts were restricted to the distribution range of each species, both by latitudinal and longitudinal ranges, ocean basins, and continental shelves (if the species is coastal or associated with the continental shelf [defined as 1200 m bathymetry]; if the species is a loliginid [defined as 500 m bathymetry]; Depth retrieved from Bio-Oracle; Jereb and Roper 2010). The ensemble forecasts were outputted as raster layers, and the average habitat suitability for each species was calculated on these by weighting each grid point with its corresponding area and dividing the sum of all values with the area of the whole species range. To assess latitudinal change, we calculated the difference between the row mean habitat suitability of each future scenario and the mean of the “present-day” scenario at the same latitude.

Table 2 Species presence points retrieved for the analysis
Full size table

All analyses were performed using R (R Core Team 2021), with resource to the packages dismo (Hijmans et al. 2017), sdmpredictors (Bosch et al. 2017), and sdm (Naimi and Araújo 2016). Maps were built using QGIS (QGIS.org 2022).

Results

The TSS value for the ensemble models (Table 3) ranged from 85.89% (Dosidicus gigas) to 98.61% (Nototodarus sloanii). Future projections (2050 and 2100) under the RCP 8.5 climate scenario for oegopsids (Fig. 1) and myopsids (Fig. 2) show decreases in habitat suitability values for all species in the lower latitudes and increases at the higher latitudes relative to the present-day (Fig. 3; all other scenarios can be seen on the Online Resource 2). Furthermore, the global patterns of habitat suitability shifts differ between squid species (Fig. 4; Table 2), with three species showing an increase (Fig. 4C, E and H), six a decrease (Fig. 4A, D, F, I, J and K) in habitat suitability across all scenarios, while two others present mixed results (Fig. 4B and G), with the direction of the shift differing according to the RCP scenario (distribution of habitat suitability for all scenarios and species on the Online Resource 3).

Table 3 Squid species distribution models’ True Skill Statistic (TSS) and mean habitat suitability (%) under 2020 conditions and different future (years 2050 and 2100) representative concentration pathways (RCP)
Full size table
Fig. 1
figure 1

Map of the distribution of main commercial oegopsids at the present (year 2020), year 2050 representative concentration pathways (RCP) 8.5 and year 2100 RCP 8.5. Map source QGIS. A Dosidicus gigas; B Illex argentinus; C Berryteuthis magister; D Nototodarus sloanii; E Illex illecebrosus

Full size image
Fig. 2
figure 2

Map of the distribution of main commercial myopsids at the present (year 2020), year 2050 representative concentration pathways (RCP) 8.5 and year 2100 RCP 8.5. Map source QGIS. A Doryteuthis gahi; B Doryteuthis opalescens; C Doryteuthis pealeii; D Loligo reynaudii; E Loligo vulgaris; F Loligo forbesii

Full size image
Fig. 3
figure 3

Habitat suitability average shift at each latitude for main commercial squid. Lines for the different years are represented with continuous lines (2050) and dotted lines (2100). Representative concentration pathways (RCP) scenarios are represented by blue (RCP 2.6), green (RCP 4.5), orange (RCP 6.0) and red (RCP 8.5). A Dosidicus gigas; B Illex argentinus; C Berryteuthis magister; D Nototodarus sloanii; E Illex illecebrosus; F Doryteuthis gahi; G Doryteuthis opalescens; H Doryteuthis pealeii; I Loligo reynaudii; J Loligo vulgaris; K Loligo forbesii

Full size image
Fig. 4
figure 4

Mean habitat suitability of main commercial squid on their current distribution, in 2020, and future scenarios (years 2050 and 2100). Representative concentration pathways (RCP) scenarios are represented by blue (RCP 2.6), green (RCP 4.5), orange (RCP 6.0) and red (RCP 8.5). A Dosidicus gigas; B Illex argentinus; C Berryteuthis magister; D Nototodarus sloanii; E Illex illecebrosus; F Doryteuthis gahi; G Doryteuthis opalescens; H Doryteuthis pealeii; I Loligo reynaudii; J Loligo vulgaris; K Loligo forbesii

Full size image

Eastern Pacific

In the Pacific coast of North and South America, for Dosidicus gigas (mean habitat suitability [MHS]: 17.75%; TSS: 85.89%), our models predict an increase in habitat suitability for latitudes above 40° along with a decrease at lower latitudes, namely between 20°S and 20°N (Fig. 1A, 3A). Mean global net habitat suitability (MHS) decreases from the present to the forecasts of 2050 and 2100 (Fig. 4A; Table 2) for all RCP scenarios. By the middle of the century, the RCP scenarios 2.6 and 8.5 show a smaller decrease in MHS than RCPs 4.5 and 6.0. In 2100, MHS for all RCPs (except for RCP 8.5) increases relative to 2050 (Fig. 4A; Table 2), although with values still lower than the present. Latitudinal differences between RCPs should also be noted, with different patterns in net habitat suitability shift between 20°N and 40°N: RCP 8.5 shows an northward increase in suitability by 2050 and a decrease by 2100; RCP 2.6, in 2050, shows an increase up until 30°N and then a decrease, while in 2100 presents a peak increase around 20°N and then fluctuates around the present values; RCP 4.5 and 6.0 show mostly decreases around these latitudes (Fig. 3A).

North American coasts

On the North Pacific coast of North America, an increase in habitat suitability for Berryteuthis magister (Mean habitat suitability: 15.16%; TSS: 94.56%) in projected future values relative to the present. This increase is most evident by 2050, in the south Bering and Okhotsk Sea, while, by 2100, this increase is restricted to the waters close to the coast in Okhotsk and the whole Bering Sea (Fig. 1C). A similar pattern is observed for Doryteuthis opalescens (Mean habitat suitability: 44.55%; TSS: 97.60%; Figs. 2B, 4G), with both an increase in habitat suitability around the south Alaskan coast and a decrease on the California Peninsula coast being evident by 2050, with the southern decrease in habitat suitability reaching the south the Monterey Bay by 2100. Latitudinally, these habitat suitability shifts translate into a permanent decrease around 45°N and a permanent increase north of 47.5°N for B. magister (Fig. 3C). Meanwhile, habitat suitability for D. opalescens decreases mostly around the latitudes south of 40°N (with two peaks around 50°N; Fig. 3G) and increases north of 52.5°N. These habitat suitability shifts result in mean net habitat suitability gains for both species, in all future scenarios and timeframes for B. magister (jumping to 20.46% [RCP 6.0 in 2100; Figs. 1C and 4C; Table 2]) and for D. opalescens in 2050 (raising to 46.93% [RCP 8.5 in 2050; Figs. 2B and 4G; Table 2]), with trends for 2100 resulting in either a habitat suitability decline, with higher CO2 emissions (RCP scenarios 8.5 and 6.0) (drop to 42.75% [RCP 6.0 in 2100; Figs. 2B and 4G; Table 2]) or an increase, with lower emissions (46.18%[RCP 2.6 in 2100; Figs. 2B and 4G; Table 2]).

Similar results are observed for the squid species present on the North Atlantic coast of North America: Illex illecebrosus (Mean habitat suitability: 41.35%; TSS: 93.88%) and Doryteuthis pealeii (Mean habitat suitability: 29.60%; TSS: 94.69%), with the former increasing to 50.12% (RCP 8.5 in 2100; Fig. 4E; Table 2) and the latter to 32.23% (RCP 6.0 in 2100; Fig. 4H; Table 2). For I. illecebrosus, decreases in suitability across all future scenarios are mostly observed south of 35°N, and between 40°N and 45°N (Fig. 3E; Table 2), corresponding to its southern range limit and the offshore waters of Nova Scotia (Fig. 1E). Meanwhile, for this species, increases are observed north of 50°N, especially for 2100 and RCPs 2.6 and 8.0 in 2050 (Fig. 3E), corresponding to the coasts of Greenland, Hudson Bay, and Baffin Island (Fig. 1E). For D. pealeii, latitudinal habitat suitability shifts mostly result in decreases south of 35°N and just south of 45°N (Fig. 3H) around the Gulf of Mexico and the offshore waters of Nova Scotia (Fig. 2C), while increases are observed north of 45°N (Fig. 4H; Table 2) on the Newfoundland and St. Lawrence waters (Fig. 2C).

Southern Hemisphere

All four squid species exclusively found on the southern hemisphere presented downward trends in their median net habitat suitability, with Illex argentinus (Mean habitat suitability: 47.81%; TSS: 94.62%) decreasing to 44.01% (RCP 8.5 in 2100; Fig. 4B; Table 2), Doryteuthis gahi (Mean habitat suitability: 35.76%; TSS: 94.20%) decreasing to 29.20% (RCP 8.5 in 2100; Fig. 4F; Table 2), Nototodarus sloanii (Mean habitat suitability: 9.76%; TSS: 98.61%) down to 7.98% (RCP 8.5 in 2100; Fig. 4D; Table 2), and Loligo reynaudii (Mean habitat suitability: 65.42%; TSS: 98.62%) decreasing to 41.80% (RCP 8.5 in 2100; Fig. 4I; Table 2). However, both I. argentinus, D.gahi, and L. reynaudii present similar (or higher) median net habitat suitability under the RCP 2.6 scenario, plus I. argentinus and L. reynaudii have the similar trends at RCPs 6.0 and 4.5. Decreases for the South American and New Zealand species center between 35°S and 45°S (Figs. 1B, D, 2A, and 3B, D, F); while, for L. reynaudii, decreases are observed north of 32°N under the RCP 8.5 by the year 2100 while, under the other scenarios, decreases appear at around 13°S, 30°S and south of 35°S (with a peak around 33°S; Fig. 3I), mostly affecting the northern, southern and western edge of the species distribution range (Fig. 2D).

European and Western Africa coast

Both species (Loligo vulgaris and Loligo forbesii; Median habitat suitability: 55.80% and 55.84%; TSS: 95.10% and 97.40%, respectively) showed lower values of global net habitat suitability under all future scenarios (Fig. 4J, K; Table 2). By the year 2050, RCP 6.0 resulted in the highest values and RCP 8.5 the lowest while, for 2100, RCP 2.6 has the highest value and lowest decrease, closely followed by RCP 4.5 and then RCP 6.0 and RCP 8.5. Decreases in habitat suitability for L. forbesii were observed on the whole Mediterranean Sea, Mauritanian coast and Atlantic islands, starting on the eastern Mediterranean and offshore Atlantic, moving towards the Tyrrhenian Sea and the northwestern African coast (Fig. 2F). L. vulgaris was similar, but in smaller scale, with the decrease in habitat suitability moving from the eastern Mediterranean, but only impacting the southern Adriatic and eastern half of the Aegean Sea, while the Atlantic coast lost a portion from Cape Verde to Ras Nouadhibou peninsula (Mauritania) (Fig. 2E). These changes were translated latitudinally in decreases in habitat suitability south of 45°N (Fig. 3J, K) and small increases around 58°N and north of 60°N.

Discussion

According to our findings, climate change does not result in increased habitat suitability for commercially important squid species. Models for commercial squid around north America forecast that most species (except Doryteuthis opalescens in 2100 for RCP 6.0 and 8.5) may benefit from warming (borealization) of the Arctic (Polyakov et al. 2020). However, a decline in habitat suitability is expected for all other species, especially those on the Southern Hemisphere. Although the squid considered here fall into two taxonomic groups (Myopsida and Oegopsida) with different spawning strategies, these groups did not consistently differ in future habitat suitability changes.

Limitations

We found high resemblance between (ensemble) models’ predictions for the present-day distribution of the studied squid and their documented distributions (Jereb and Roper 2010). However, caution should be taken into the future projections here made, as there are limitations of the species distribution models (SDM), that may explain differences between (future) projections of habitat suitability and future observed species distribution. First, it should be noted that while we considered multiple variables in our analysis, we only had one biotic variable which was primary productivity (chlorophyll) with both present and future values. This productivity at the local habitat level, limits the resources available to sustain metazoan biomass, including squid. Second, prey biomass predictions and other species (predators and competitors) were not included in our models, potentially leading to an overprediction of habitat suitability. Squid prey, predators and competitors may limit squid distribution. Third, our models do not consider key abiotic factors. For example, dissolved oxygen is affected by climate change and may limit the distribution and habitat suitability of vertically migrating squid and their competitors (Rosa and Seibel 2008, 2010; Rodhouse et al. 2014). One of the species here studied, D. gigas, can cope with a wide range of dissolved oxygen and as a result, has been expanding along with the oxygen minimum zones (OMZ) in the Pacific (Zeidberg and Robison 2007; Rosa and Seibel 2010). Alongside, bottom waters oxygen levels limit the distribution of coastal squid where OMZs occur, as for example the L. reynaudii (Sauer et al. 2013; Van Der Vyver et al. 2016). Likewise, pH represents an abiotic factor worth considering, as previous climate change studies have highlighted the potentially negative effects of ocean acidification on squid physiology and behavior (Rosa and Seibel 2008; Rosa et al. 2012, 2014; Rodhouse et al. 2014) and hence, an important environmental determinant for their distribution. Fourth, occurrence data for these species do not consider the life stage of the animals observed. Squid perform ontogenic migration during their life cycle (Jereb and Roper 2010), which takes them to a wide range of environmental conditions. It is likely that our dataset is biased towards the adult phase, neither covering egg and paralarvae stages, the most sensitive stages to environmental changes (Oosthuizen et al. 2002; Rosa et al. 2012). Finally, many species selected in this study have either few occurrence points available online (below 100 points; T. pacificus, D. gahi and D. gigas) or have their sampling skewed to the waters of countries with extensive resources for surveys (particularly the U.S.A coastlines), with limited to no information in other areas of their range (D. pealeii, D. gahi, I. illecebrosus and D. gigas). Furthermore, due to the difficulty in distinguishing some of these squids from closely related species (e.g., D. pealeii vs. D. pleii, D. gahi vs. D. sanpaulensis, N. sloanii vs. N. gouldii, L. vulgaris vs. L. forbesii; Jereb and Roper 2010), many occurrence points may represent cases of misidentification, especially in areas where the distribution of two similar species overlaps. This issue may disproportionally affect D. pealeii, likely contributing to an underestimation of the habitat suitability for this species within the Caribbean basin.

Species-specific and regional changes in habitat suitability

The 11 squid species here studied did not change their distribution equally (as mentioned in the beginning of this discussion), and the consequences of climate change seemed species specific. However, within the same region, squid presented similar patterns of future habitat suitability change.

The Jumbo squid (D. gigas) in the Pacific East coast, is projected to expand poleward, which is in line with the documented patterns of migration towards higher latitudes during el Niño and recent slight expansion of their northern range (Field et al. 2007; Arkhipkin et al. 2015). This poleward shift is strongest in the northern hemisphere, where habitat suitability will potentially increase along the shores of Canada and USA. On the other hand, a decrease in suitability is projected for the Sea of Cortez and the areas from the Costa Rica Dome upwelling system towards the 140°W longitude, where this species is thought to follow the Pacific Equatorial current (Wormuth 1998; Jereb and Roper 2010). Indeed, across all future scenarios, this species may potentially be extirpated from tropical areas. Therefore, with these predicted extirpations and poleward expansion, future conditions may favor the split of this species into two different stocks without geographic continuity. Different size stock structures have been suggested in the past, either a 3 different sizes stocks (Nigmatullin et al. 2001) or a size continuum (Hoving et al. 2013). With a decrease in both area of habitat as well as its suitability, smaller animals could be benefitted (specially on the lower latitudes where habitat suitability loss concentrates), as these ecosystems carrying capacity will decrease, favoring smaller individuals. This has already been observed in this species (Hoving et al. 2013) and other squid in less favorable years (D. opalescens; Jackson and Domeier 2003; Zeidberg et al. 2006). The opposite (larger individuals) may also happen locally at higher latitudes, where the species is expanding into new habitat (Hoving et al. 2013). Finally, the gains in habitat suitability obtained at higher latitudes are less than the losses around the tropical regions. As a result, there is an overall loss of habitat for D. gigas under all future scenarios.

In the North Pacific, a northward increase in habitat suitability is expected for both B. magister and D. opalescens, with suitability decreasing in their southern limit. For B. magister, gains in the north will surpass the losses in the south, as there are large areas of low-depth high primary production continental shelves in the north, on the Pacific coast of Canada, the Bering Sea and Okhotsk Sea. Although the squid occurs between 0 and 1500 m, the highest densities of B. magister are generally found at depths greater than 200 m, near the bottom on the continental slope and in the mesopelagic zone (Jereb and Roper 2010) where they spawn. In this context, their depth use may limit their ability to take advantage of the habitat projected to become suitable (specifically on the northern parts of the Bering Sea). This may ultimately even lead to a population decrease, as spawning grounds on the southern limits of this species (mostly continental slopes) are negatively affected. D. opalescens is projected to experience a net loss of suitable habitat in 2100 for the two worst RCP scenarios, which is in line with the pattern observed during el Niño years, when the species’ abundance (and size) decreases throughout their home range (Jackson and Domeier 2003; Zeidberg et al. 2006).

Northward increase is experienced in squid species found on the northwest Atlantic Ocean, where projected squid habitat suitability increased poleward. Both I. illecebrosus and D. pealeii occur in this region, with I. illecebrosus preferentially using the colder northern waters and D. pealeii the warmer waters, in the South (Brodziak and Hendrickson 1999). Under future climate scenarios, both species are projected to shift their distribution poleward. Meanwhile, their northern limits of distribution are expected to increase in habitat suitability, where I. illecebrosus may gain new habitat on the continental platforms around Greenland, Iceland and Baffin Island. A habitat increase is also projected for D. pealeii around the great banks off Newfoundland, a pattern which has already been observed in warmer years (Dawe et al. 2007). These gains compensate for the losses projected in the southern distribution range of both species, with net gains increasing along with the severity of RCP scenarios. Should the borealization of the Arctic (i.e., increase in temperature, salinity and decrease in ice cover much like the surrounding temperate waters) continue for the next few hundred years, the models suggest that squids may be able to cross the Artic Ocean and settle in new ocean basins, as it has been hypothesized for the cuttlefish Sepia officinalis (Xavier et al. 2016b). In fact, this suggestion is supported by the recent evidence of squid migrating northwards around the North American continent (Arkhipkin et al. 2015; Burford et al. 2022).

In the southern hemisphere, climate change is likely to pose a considerable threat to squid (i.e., I. argentinus, D. gahi, N. sloanii, and L. reynaudii). With no nearby poleward continental platforms to colonize or facing other oceanographic features likely to prevent their expansion (e.g., Antarctic circumpolar current), these squid species stand at an oceanographic dead-end. Furthermore, the areas with the highest abundance of I. argentinus, D. gahi and N. sloanii (around 40°S; Haimovici et al. 1998; Jereb and Roper 2010) overlap with those projected to face the steepest decline in habitat suitability. The habitat of L. reynaudii is projected to decrease in size particularly offshore and in the southeast where currently two-thirds of the adult biomass is found (Augustyn 1989, 1991; Augustyn et al. 1993; Sauer et al. 2013). As a result, the population may be found closer to shore along the South African and Angolan coasts (Shaw et al. 2010). Additionally, if the OMZ associated with the Lüderitz upwelling cell increases in the future, this could further isolate the two different stocks present in South Africa and off the Cunene/Kunene river plum (Van Der Vyver et al. 2016). However, should milder RCP scenarios come to pass, L. reynaudii, D. gahi, and I. argentinus may rebound (or even prosper, in the case of the latter) by 2100, with models predicting the overall habitat suitability for that period to values like those currently observed.

Lastly, the habitat suitability for northeastern Atlantic loliginid squids (L. vulgaris and L. forbesii; Hastie et al. 2016) is projected to decrease substantially, primarily around the Mediterranean basin and off the coast of northeast Africa, aligning with previous projections for L. vulgaris under climate change (Schickele et al. 2021). L. forbesii may be extirpated from the Mediterranean Sea, the northeast coast of Africa and the Atlantic islands, following the temperature-driven northerly abundance shift registered for this species since the 1990’s (Chen et al. 2006). This trend is also observed, albeit to a lesser degree, for L. vulgaris in the south Adriatic, east Mediterranean, and Mauritanian upwelling—the areas where this species is most abundant (Cunha et al. 1995).

Potential foodweb and exploitation effects

Assuming that habitat suitability is likely to translate into squid abundance, many squid predators (like marine mammals and seabirds; Clarke 1996; Croxall et al. 1996; Klages and Clarke 1996) dependent on the studied squid species may face in the future difficulties to find enough prey biomass on their current feeding grounds, having to either migrate to new areas (poleward), change target prey, decrease in number/biomass, or a combination of these options. In regards of fisheries, we may witness a change in the location of economically relevant fishing grounds for many of the commercial squid studied here. Additionally, a decrease in abundance (and therefore available stock) for several of these species is likely (assuming that habitat suitability is translate into squid abundance). The two most important squid fisheries worldwide (D. gigas and I. argentinus) can be hard hit, as they could experience both effects (range and abundance change) and their fisheries are currently unregulated on the high seas (Arkhipkin et al. 2022). Other locally important squid (N. sloanii, genus Loligo and Doryteuthis opalescens) will experience the same effect. If local extirpations are to happen, many fishing communities may face economic hardship (Downey et al. 2010; Arkhipkin et al. 2015), having to look for other marine resources to compensate for the loss of a profitable fishery.

Future studies

In the future, studies of SDMs in the pelagic environment should strive to produce three-dimensional projections of habitat suitability (Duffy and Chown 2017; Aspillaga et al. 2019), namely using present-day environmental data from repositories such as Copernicus marine services (Le Traon et al. 2019), which already provides marine information with depth resolution, and in terms of presence occurrence points, further effort should be put on providing depth information. Also, as different life stages of squid (i.e., eggs, paralarvae, juveniles and adults) require different environmental conditions (and therefore migrate during their lifetimes; (Jereb and Roper 2010)), different models should be made for each stage (which also implies further data collection on their life stage on the occurrence points). This is especially relevant within squid, as there are two different egg laying strategies for myopsids and oegopsids. Myopsids fix their eggs into the substrate, and therefore their eggs are more sessile and benthic (and dependent on seabed environmental conditions) while oegopsids release their eggs into the water column or even transport and care for them (Jereb and Roper 2010). The environmental stability provided by the water masses where oegopsid eggs are released and parental care that some squid provide can save many of these species from the nefarious effects that warming and acidification have on early life stages in contrast with myopsid’s eggs that can be exposed to different water masses that pass the seabed where they are settled. For myopsid embryos, lowered survivability (from normal 92–96% to 47% with + 2 ℃; Rosa et al. 2014) and increased occurrence of developmental abnormalities (Rosa et al. 2012, 2014) are expected with future warming, and hatchlings, having to deal with an increased metabolic rate and higher energy demands (however, this last point could be also applied to the oegopsid paralarvae; Rosa et al. 2012). Additionally, higher temporal resolution (i.e., Seasons or months) should be obtained for the environmental rasters. The decade average rasters used not only mask all the variation happening during the year cycle (during which squid migrate and change locations to follow the most suitable habitat; Jereb and Roper 2010), but also provide wrong information when collecting environmental data from the occurrence points during the model construction. Finally, this kind of studies in the future should cover a wider variety of animals, not only because there are others that have key roles in the ecosystem, but also to use as layers for other species distribution models.