Impact of climate change on the distribution and habitat suitability of the world’s main commercial squids

Climate change is expected to have major negative effects on marine life across phylogenetic groups. Cephalopods, however, have life history characteristics that suggest they may benefit from certain climate change scenarios. Of all cephalopods, squids reach the greatest biomasses; as a result, they are of substantial importance for human and predator consumption. To test the hypothesis that the effects of climate change are beneficial for commercial squid, we used species distribution models on climate scenarios for the period between 2000 and 2014, as well as the years 2050 and 2100 (RCP [representative concentration pathway] 2.6, 4.5, 6.0, and 8.5; CMIP5). Our results suggest that consequences of climate change scenarios are species specific. In the North Pacific and Northwest Atlantic, squid’s habitat suitability may increase (from + 0.83% [Doryteuthis pealeii] to + 8.77% increase [Illex illecebrosus]), while it is predicted to decrease in other regions (from  − 1.03% [Doryteuthis opalescens] to − 15.04% decrease [Loligo reynaudii]). Increases in habitat suitability occurred mostly at higher latitudes (north of 50° N), while suitable habitat decrease was predicted for the tropical regions. These shifts in future habitat suitability were stronger under harsher emission scenarios. Starting in 2050 (with RCP scenarios 4.6, 6.0 and 8.5), as a result of warming of the Arctic, squid habitat may increase along both coasts of North America. In the Southern Hemisphere, squids may lose habitat with no poleward habitat alternatives to move into. Contrary to our hypothesis, these commercial squid do not stand to benefit from climate change. Since these squid are an important food source for marine megafauna and humans, it is imperative that climate change biogeographic impacts are considered for a sustainable management of this important group of molluscs.


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 Responsible Editor: R. Villanueva. greenhouse gases and consequent projected radiative forcing trajectories (in W/m 2 ; 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.

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).
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. 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)

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) 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 Fig. 3I), mostly affecting the northern, southern and western edge of the species distribution range (Fig. 2D).  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 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. 2012Rosa et al. , 2014Rodhouse 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 (

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 Page 11 of 16 129 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(Augustyn , 1991Augustyn 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 Fig. 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 (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(Rosa et al. , 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.