Future distribution patterns of nine cuttlefish species under climate change

While most coastal communities are expected to, or have been, negatively impacted by climate change, cephalopods have generally thrived with shifting ocean conditions. However, whilst benefitting from the same physiological flexibility that characterizes cephalopods in general, cuttlefish have depth constraints imposed by the presence of a cuttlebone and are limited to specific locations by their particularly low vagility. To evaluate the potential effects of marine climate change on cuttlefish, Species Distribution Models (SDM) were applied to nine species of genus Sepiidae to assess potential changes to their future distribution (2050 and 2100), under four representative concentration pathway (RCP) scenarios (i.e., RCP 2.6, 4.5, 6.0, and 8.5; CMIP5). We show that future cuttlefish habitat suitability and distribution will potentially decrease. The species with the most extreme impacts, Doratosepion braggi (Verco, 1907), was observed to decline as much as 30.77% in average habitat suitability (from present 55.26% to 24.48% at RCP 8.5 in 2100), to Sepia officinalis Linnaeus, 1758 with a low maximum decrease of 1.64% in average habitat suitability (from present 59.62% to 57.98% at RCP 8.5 in 2100). Increases in habitat suitability were projected mostly at higher latitudes, while habitat decrease was predicted for the tropical regions and lower latitudinal limits of species’ distributions. As their habitats decrease in terms of habitat suitability, cuttlefish may not benefit from future changes in climate. Additionally, as potential “sea canaries” for coastal ecosystems, we may see many species and habitats from these systems affected by climate change, particularly in tropical regions.


Introduction
Coastal systems are currently the marine ecosystems most affected by anthropogenic pressures (Halpern et al. 2008), due to their close proximity to urban centers (Hinrichsen 1996).These areas are valuable to humans, as they are highly productive (Nixon et al. 1986;Carr et al. 2006), maintaining high biodiversity levels (Kaiser 2020) and serving as nursery grounds for fisheries species, thus promoting the proliferation of heavy fishing efforts (Hinrichsen 1996;Chassot et al. 2010;Stock et al. 2017).Also, waste resulting from human activities (e.g., agricultural, industrial, and sewage waste) can lead to impacts such as eutrophication (Smith et al. 1999;Paerl and Otten 2013) and the associated high mortality rates in wildlife communities (Burkholder et al. 1992;Penela-Arenaz et al. 2009).Furthermore, anthropogenic climate change is poised to affect these areas through changes in temperature, local currents, chemistry, bathymetry, and geomorphology, which, in turn, are likely to impact biological communities (Harley et al. 2006;Cinner et al. 2012;Ranasinghe 2016).
Coastal benthic communities are under added high risk due to climate change (IPCC 2022).The shallow depths in which these communities live are under increased risk of stratification, as warming is expected to reduce vertical mixing (Hallett et al. 2018).These, in turn, translate into greater exposure to ocean warming (OW), oxygen loss, and ocean acidification over coastal regions (Sampaio et al. 2021;IPCC 159 Page 2 of 15 2022).Furthermore, shifts in precipitation regime, greater runoff and water input into these systems is expected in the future, with implications for the levels of salinity and primary productivity (Bindoff et al. 2019).Shell-bearing benthic and demersal coastal species, having limited dispersal capabilities, may be particularly affected by these combined phenomena.Currently, the distribution of coastal species has been either observed (Lu and Lee 2014;Pang et al. 2018) or projected (Boavida-Portugal et al. 2018, 2022;Borges et al. 2022) to shift poleward or into deeper waters (Perry et al. 2005) under climate change scenarios.Additionally, coastal communities are expected to face local extinctions (Albouy et al. 2013) along with structural and functional changes at the ecosystem scale (Osland et al. 2013).
Cephalopods (i.e., squid, octopuses, and cuttlefish) are considered climate change indicator species (Jackson and Domeier 2003).The coastal representatives of these molluscs are characterized by their life strategy of "live fast, die young" (Boyle and Rodhouse 2008), that allows them to adapt quickly to new thermal environments (Rodhouse 2013).They are also important sources of food and play an increasingly important role on both artisanal and commercial fishing (Caddy and Rodhouse 1998;Jereb and Roper 2005).Furthermore, cephalopods are important trophic mid-points on the food-web, between low trophic levels (Rodhouse and Nigmatullin 1996) and top predators (Clarke 1996;Croxall et al. 1996;Klages and Clarke 1996).Therefore, they meet most of the attributes for being good ecological indicators (Fulton et al. 2005): (i) short life cycles and fast growth rates; (ii) targeted by fisheries; and (iii) keystone species (Rosas-Luis et al. 2008); with populations reactive to environmental changes.Within the cephalopods, cuttlefish (family Sepiidae) are particularly dependent of the coastal systems.These demersal animals have developed an internal calcareous shell for buoyancy control.The ultrastructure of the cuttlebone imposes depth constraints that restrict particular species to particular depth ranges (between 300 and 600 m, devised from testing only a couple of species; Ward and Boletzky 1984;Sherrard 2000;Jereb and Roper 2005), limits the ability of this group to seek refuge across the vertical plane.Cuttlefish can venture into deeper habitat than what their cuttlebones allow (as specimens collected from depths in excess of 800 m show), be it by temporary flooding of the cuttlebone chambers, or partial collapse of these chambers (Birchall and Thomas 1983;Gower and Vincent 1996).However, these situations are likely to be extreme events in these animals' lives.Furthermore, their population dispersal is horizontally limited, as they have direct development with no planktonic phase (Jereb and Roper 2005;Xavier et al. 2016).Lastly, their eggs and early developmental stages are vulnerable to the effects of climate change, with cuttlebone development being affected by ocean acidification (OA; Dorey et al. 2013;Sigwart et al. 2016;Otjacques et al. 2020) and ocean warming (OW) leading to higher respiration rates, protein damage, membrane fluidity and organ malfunction, lower survival rates and to decreases in average body mass (Lacoue-Labarthe et al. 2009;Rosa et al. 2013;Sigwart et al. 2016).Hence, if cuttlefish closely follow the quality trend of the coastal systems with climate change, their populations can potentially decrease in the near-future (Harley et al. 2006).In fact, while overall coastal and oceanic cephalopod populations have increased over the last 60 years worldwide (Doubleday et al. 2016), both due to decrease in competition and predation from commercial fish species and to the ability of cephalopods to adapt to changing environments (André et al. 2010), populations of certain cuttlefish species have been in decline due to overfishing (Pang et al. 2018).
Species distribution models (SDMs), which allow for the estimation of species habitat suitability, have been used to predict and project cephalopod distributions (Schickele et al. 2021;Boavida-Portugal et al. 2022;Borges et al. 2022).If the presence of a species in a given location means that the area meets the environmental conditions for the species to survive (Pearson 2007;Kearney and Porter 2009), we can collect the environmental data from the areas where the species are present to predict the species' broader distribution and apply it in a broader region (e.g., project the invasive potential of species; Barbet-Massin et al. 2018;Sung et al. 2018;Borges et al. 2021) or timeline (e.g., project the potential impact of future climate change scenarios; Araújo and Rahbek 2006;Melo-Merino et al. 2020;Borges et al. 2021) and provide values of suitability of habitats across these new environments.Hence, if the ongoing cephalopod populational increase is due to climate change, future SDM projections with scenarios of worsening climate change should have higher habitat suitability for cephalopods.This should be especially evident in coastal ecosystems that are more exposed to these changes.In these areas, cuttlefish are the cephalopods more dependent on the environmental conditions, as they are limited in depth to the coastal area (due to their cuttlebone implosion depths) and have a benthopelagic lifestyle.
If cephalopods are benefiting from climate change, cuttlefish, as benthopelagic coastal representatives of this group should benefit as well, especially as the environment they occupy will be so severely and rapidly affected.To test this hypothesis, we used data collected from online repositories and SDM framework to forecast the potential future distribution of nine most conspicuous cuttlefishes under different climate change representative concentration pathway (RCP) scenarios (numbered according to the radiative forcing [W/ m 2 ; Moss et al. 2008] projections; RCP 2.6, 4.5, 6.0, and 8.5) and timeframes (present, year 2050, and 2100).With the overall habitat suitability from these forecasts, we can check if coastal cephalopods will benefit from worse climate Page 3 of 15 159 change scenarios and expand their habitat range into a wider latitudinal amplitude.

Data collection
We collected geolocated data on species occurrence from the Global Biodiversity Information Facility (GBIF.Org User 2022) and the Ocean Biodiversity Information System (Grassle 2000).On the online repositories, we used the search term 'Sepia' to limit our search to the genus (as formerly defined to contain most sepiid species prior to a recent revision by Lupše et al. (2023)), and with the downloaded occurrences (Grassle 2000; GBIF.Org User 2022) we limited to the species with more than 50 occurrences.Of the observations obtained, only nine species were selected for modeling.These nine species represent five out of the 12 sepiid genera recognized by Lupše et al. (2023).The generic names recognized by Lupše et al. (2023) have been used here (Abbreviations: Ac., Acanthosepion; As., Ascarosepion; De., Decorisepia; D., Doratosepion; R., Rhombosepion; S., Sepia).To curate our occurrence dataset for points in wrong locations (i.e.: located on land, wrong ocean basin, outside of distribution range; Jereb and Roper 2005;Pierce et al. 2010;Reid 2016;Hastie et al. 2016; for the locations of the point dataset, See Online Resource 1), by manually curating it in QGIS (QGIS.org2022).Lastly, a species-level thinning step (rarefaction) was implemented by randomly picking occurrences distanced 20 km or more from each other, as this reduces spatial sampling biases (Elith et al. 2011;Radosavljevic and Anderson 2014).The distance selected (20 km) for this last step was the one that made our models present-day predictions resemble best the literature (Jereb and Roper 2005;Pierce et al. 2010;Reid 2016;Hastie et al. 2016).In total, 15,695 observation points were retrieved from the repositories for the nine species (Table 1).Of these, 2463 were removed due to errors in location.From the 13,232 curated locations, only 3407 were used after rarefication (Table 1).
We selected temperature, salinity, current velocity, depth, and chlorophyll for the environmental variables for modeling.These were retrieved from the Bio-Oracle project (Tyberghein et al. 2012;Assis et al. 2018).The values for these layers correspond to the measurements performed at the average seabed depth (except for chlorophyll, where values were retrieved for the sea surface), with a resolution of 0.00833(3) degrees (latitude and longitude; projection = " + proj = longlat + datum = WGS84 + no_defs").For each environmental variable, 3 raster layers were collected, one for the longterm mean, another for maximum and the last for minimum values.These were retrieved for both the present-day, as well as all future climate scenarios (years 2050 and2100;RCP 2.6, 4.5, 6.0, and 8.5;CMIP5).To assess the collinearity between the environmental variables, we used Pearson's method (r).We removed variables with a collinearity > 0.7 (Naimi and Araújo 2016) to decrease potential biases in inference statistics of the models (Dormann et al. 2013).After this step, only the mean values were kept for most environmental samples, with minimum chlorophyll also maintained.

Species distribution models
To build our species distribution models, we made use of the MaxEnt algorithm (Phillips and Dudík 2008).A total of five model runs were done using all data (presence points plus 1 000 [randomly selected] background points) for each species, where the data was randomly partitioned and shuffled between calibration (75%) and validation (25%) datasets at each run.Subsequently, we used True Skill Statistic (TSS; Allouche et al. 2006) to assess the performance of each model.Models with a TSS lower than 70% were discarded.We calculated the weighted median of the remaining models to produce the ensemble (Araújo and New 2007) by weighting each of the model runs forecasts with the square difference of its TSS score with 50%.The ensemble model TSS as well as for environmental variable contribution to the model used the same computation.The previous step was performed for each of the nine climate scenarios and for each species, producing 81 Habitat suitability forecasts (9 species × 9 scenarios [present + 2100*four RCPs + 2050*four RCPs]).These forecasts were restricted to the distribution range of each species, both the latitudinal and longitudinal ranges, ocean basins, and depth range (Jereb and Roper 2005).The ensemble forecasts were outputted as raster layers.
The average habitat suitability of each species was calculated by weighting each grid point on the rasters with its corresponding area and dividing the sum of all values with the sum of the whole species range area.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.Furthermore, main contributing environmental variables were retrieved.We used R (R Core Team 2021) to run all analysis, with resource to the packages dismo (Hijmans et al. 2017) and sdm (Naimi and Araújo 2016).Maps were built using QGIS (QGIS.org2022).
159 Page 4 of 15 Overall, the sepiids examined showed decreases in future habitat suitability (Figs. 1, 2; all other scenarios can be seen on the Online Resource 2; distribution of habitat suitability for all scenarios and species on the Online Resource 3), with losses of habitat focused mainly on lower latitudes (Fig. 3).

Overall habitat suitability trends
For 2100, habitat suitability (HS) declines increase with RCP severity for all targeted species (Fig. 1A-I).For 2050, an overall decline in mean habitat suitability (mHS) over time was projected across all emission scenarios for most cuttlefish species (Fig. 2 and Table 2).The only exception was D. braggi, which represented the only species with a projected increase of mHS (1.65% increase; RCP 2.6 year 2050; Table 2).The species with the highest differences in mHS between present and future scenarios was D. braggi (i.e., biggest decrease: 30,77%; RCP 8.5, year 2100; Table 2), and S. officinalis Linnaeus, 1758 presented the lowest differences between present mHS and future scenarios (biggest decrease: 1.64%; RCP 8.5, year 2100; Table 2).
The overall trends in mHS varied between species and RCP scenarios under consideration.The mHS values from the projections for different RCP for the year 2100 deviated more from the predicted present values than the same values for 2050 (Fig. 2).The mHS for the 2050 RCPs' are very similar to each other, with the exception of Ascorosepion apama (Gray, 1849) (Fig. 2A), S. bertheloti (Fig. 2C), D. braggi (Fig. 2D) and Acanthosepion pharaonis (Ehrenberg, 1831) (Fig. 2I) and do not show a value order related to RCP.However, mHS for the 2100 RCPs are lower than the mHS for the predictions for the present, with lower values according to the severity of RCP (I.e.: 2.6 > 4.5 > 6.0 > 8.5) (Fig. 2).3I).For 2100, HS declines increase with RCP severity (Fig. 1I; Fig. 3I), peaking around 25ºS (up to ~ 50% HS decline), 10ºS (~ -40% HS), the equator (~ -30% HS) and 25ºN (~ -30% HS; Fig. 3I).On the other hand, for 2050, and specifically between the equator and 30ºN, RCP6.0 resulted in the lowest HS decline (~ -5% HS) while RCP2.6 resulted in the highest (~ -25% HS), being closely followed by RCP8.5 (~ − 25% HS; Fig. 3I), with HS declines still increasing with RCP severity for the remaining latitudinal range.This same pattern found between the equator and 30ºN for the year 2050 RCPs is also observed for the mHS (Fig. 2I).Increases in HS were observed north of 30ºN (Yellow Sea, north of the South China Sea and southern coast of Japan; Fig. 1, 3I), with this trend being exacerbated over time and with RCP severity (Fig. 3I).

African Coast
The HS for Decorisepia australis (Quoy and Gaimard, 1832) (mHS: 41.98%; TSS: 96.00%) and S. bertheloti (mHS: 54.86%; TSS: 94.75%) mostly decreased throughout their distribution around the African coast.Decorisepia australis experienced a decrease throughout its latitudinal range in HS, except for an increase between 5ºS and 10ºS in 2050 (mouth of the Congo River for RCPs 6.0 and 8.5 and at the south coast of the Horn of Africa for RCP 2.6; Figs. 1,  3B) where it experiences a 40% increase peak (RCP 6.0; Fig. 3B).Largest decreases were observed on all future  3B).Between 5ºN and 20ºN, strongest decreases in 2050 were observed on the RCP 2.6 projection (and RCP 6.0 showed the weakest decreases, to punctual increases) and in the year 2100, it changes to the RCPs 6.0 and 8.5.This same pattern is also observed for the mHS (Fig. 2B), where RCP 6.0 presented the highest (and RCP 2.6 the lowest) values for the year 2050.South of 5ºN, lower values are exacerbated over time and with RCP severity (Fig. 3B).Within the latitudes close to the equator (± 5º), low differences (< 5%) were observed between habitat suitability based on future projections.

Eastern Atlantic
On

Discussion
Overall, the results obtained in the present study suggest a decrease in habitat suitability for the nine Sepiid species examined, with slightly different responses between the different species analyzed.However, the distribution of the loss in habitat suitability is not uniform, being most evident at lower latitudes, while some gains are projected on the higher latitudes (that do not compensate enough for the losses).These nine species represent five out of the 12 genera recognized in Lupše et al. (2023).
Despite the high resemblance of the ensemble models' projections for the present-day scenario, with the current descriptions of the distributions of the targeted cuttlefish species (Jereb and Roper 2005;Pierce et al. 2010;Reid 2016;Hastie et al. 2016), limitations inherent to species distribution models (SDM) and the data used may explain divergences between the habitat suitability values projected and to be observed over incoming decades.First, the absence of trophic interactions or prey biomass as variable predictors can lead to instances of overprediction in habitat suitability, as prey availability and the existence of predators and competitors can limit the distribution of cuttlefish.Second, given that dissolved oxygen and pH projections are not available at a global scale for the future scenarios analyzed, these variables were not integrated in the present analysis, limiting the predictive power of our analysis, as both these abiotic factors are important for cuttlefish.Specifically, dissolved oxygen acts as a limiting factor in the distribution of many cuttlefish species (Augustyn et al. 1995;Melzner et al. 2004) while lower pH conditions leads to hyper-calcification of the cuttlebone (Dorey et al. 2013;Sigwart et al. 2016;Otjacques et al. 2020) and thus, may impact the future vertical distribution of these animals.Third, many observations of these animals were made post-mortem, with identification based on their cuttlebones alone rather than whole animals.Floating to the surface after the death and decomposition of the organism, these structures can end up hundreds of kilometers from their original location (Reid 2016), distorting their real distribution.With a large portion of data points retrieved from GBIF resulting from citizen science (INaturalist Contributors 2023), where many of these observations are made along the seashore (Org User 2022), such method of detection is likely frequent.Fourth, for many species targeted in this study, occurrence data available online were either low in extent (less than 200 occurrence points; S. bertheloti, As. latimanus and D. braggi), or skewed to waters of countries with extensive resources for surveys (particularly European, South African and Australian coastlines), with other parts of the ranges of these species have either very few or no observations (De. australis, Ac. pharaonis, As. latimanus, S. bertheloti, R. orbignyanum, and R. elegans).Fifth, presence data collected for this study do not include the life stage of the individual observed.Cuttlefish may occupy different depth habitats at different stages of their life cycle.Some species appear to move from deeper to shallower waters to breed and lay eggs, so juveniles may have differing depth limitations and requirements than adults (Jereb and Roper 2005).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 (Dorey et al. 2013;Sigwart et al. 2016;Otjacques et al. 2020).Sixth, species distribution models do not take into account the life and biogeographical history of the species analyzed.Therefore, although we have limited our projections to within the current range of the species and contiguous areas where they can migrate, predictions and projections obtained by these models do not (and will not) fit perfectly with current and future distributions.Finally, due to the difficulty in distinguishing some of the more closely related cuttlefish species, many occurrence points may have been misattributed to a different species (the distributions of some cuttlefish species are not well known for this reason), especially in areas where two confounding species overlap their distributions (e.g.: Sepia officinalis vs. S. hierreda vs. S. bertheloti vs. R. elegans vs. De.australis vs. R. orbignyanum, Ac. pharaonis vs S. ramani; Jereb and Roper 2005).This is especially relevant for the eastern Atlantic cuttlefish, as this may have contributed to the underprediction of the models for these species on the western African coast.
Ascarosepion apama is the largest cuttlefish in the ocean, and occurs around southern Australia.Two lineages have been recognized to occur off south-eastern Australia, one eastern and one southern, separated by Bass strait, due to historical vicariance, events involving the opening and closing of the strait due to sea level changes (Kassahn et al. 2003).One population of the southern lineage forms a unique and only spawning ground known for cuttlefish, with sometimes thousands of individuals gathering at the north end of Spencer Gulf, South Australia (Hall and Hanlon 2002;Hall and Fowler 2003;Kassahn et al. 2003).According to our findings, this spawning ground is projected to exhibit a considerable decrease in habitat suitability starting in 2050 and getting worse in 2100 (for all RCPs except 2.6; worse for each higher RCP), which may lead to declines in population recruitment, at least for the Spencer Gulf population, which is in line with population trends previously observed (Prowse et al. 2015).Populations inhabiting Bass Strait and Spencer Gulf are projected to undergo through the greatest habitat suitability loss for the species.However, if the spawning ground at Spencer Gulf is unique to the animals living in the region, and not to the entire species, 159 Page 10 of 15 we may see a future shift southward towards the coastline of Tasmania, together with the isolation of the westernmost groups of animals inhabiting southwestern Australia (this last effect appearing in 2100 for scenarios higher than 2.6).Finally, our modeling for this species predicted what previous studies have suggested: that temperature represents a primary contributor for habitat suitability, as it affects egg and juvenile development (Hall and Fowler 2003).Around the same region, the habitat of Doratosepion braggi is projected to undergo a considerable decrease in suitability, relative to all the other species (except for RCP 2.6 in 2050; worse for each higher RCP in the year 2100).This species has a very limited depth range (30-86 m; Jereb and Roper 2005) off the south coast of Australia, limiting their ability to migrate to new habitat.In this context, this species is particularly sensitive to projected environmental change, especially temperature, as it explained more than 90% of the distribution in our models.The only two options for escaping from rising temperatures are moving down the narrow continental shelf and further south around the Tasmanian coast.This species will be facing the largest declines in the set of cuttlefish targeted by this study and be reduced to two separated regions of Australia as early as 2050 (southwestern tip and the Tasmanian coast) regardless of future climate change scenario (except for RCP 2.6 in 2050).However, the apparent limited depth range of D. braggi may reflect the limited benthic sampling effort carried out off the southern coast of Australia, as relatively few biological collections bellow 100 m have been undertaken in these waters.
Ascarosepion latimanus and, particularly, Acanthosepion pharaonis represent the main commercial cuttlefish species in the Indo-Pacific.Ascarosepion latimanus inhabits shallow tropical waters (no deeper than 30 m), relying on coral reefs and mangrove areas (Jereb and Roper 2005) for shelter and egg laying grounds.Both these ecosystems are projected to experience distribution shifts in the future climate change scenarios (Record et al. 2013;Kornder et al. 2018) which may affect recruitment success and survival of the early stages of this species.Furthermore, with temperature being the greatest contributing factor for the models of this species (~ 80%), followed by depth (~ 11%) and salinity (~ 6%), future increases in coastal water temperatures and precipitation (and terrestrial runoff) will impact the survivability of this species in enclosed and calm seas.More specifically, in the Indonesian archipelago, the Andaman Sea and the Gulf of Thailand, starting in 2050 (all RCPs) and becoming worse with RCP severity by 2100 (with signs of recovery with RCP 2.6).Acanthosepion pharaonis does not differ much from the previous species in terms of trends in our projections, as it is projected suitable habitat loss in the Indonesian archipelago, the Andaman Sea and the Gulf of Thailand as well, starting in 2050 (all RCPs) and becoming worse with RCP severity by 2100 (with signs of recovery with RCP 2.6).However, as it can inhabit deeper (up to 100 m) and colder (spawning takes place between 18ºC and 24ºC) waters (Jereb and Roper 2005), this may explain the difference between main contributing environmental variables to the distribution between the two co-occurring species (for Ac. pharaonis temperatures contributes more than 90%).Also, as the greatest contributing factor is the seabed temperature, greater losses in habitat suitability in 2050 for RCP 2.6 compared to RCP 6.0 can be explained by average higher (seabed) temperatures on RCP 2.6 and average lower (seabed) temperatures in RCP 6.0.This is especially relevant between the equator and 30ºN, where the species has a higher habitat suitability in 2020 and wider continental shelfs to inhabit.Additionally, we should consider that Ac. pharaonis is in fact a group of five cryptic species (Anderson et al. 2007(Anderson et al. , 2011)).If the projections that start in RCPs 2.6 and 8.5 at 2050 and worsen by 2100 (for RCPs greater than 2.6) are to materialize, the cryptic species from western Indian Ocean, northeastern Australia and Persian Gulf/Arabian Sea ('Iranian') could be endangered by the loss of habitat suitability.
Off the sub-Saharan coast of Africa, Decorisepia australis and Sepia bertheloti represent the main stocks of cuttlefish.According to our models for De.australis, this species seems to be unlikely to inhabit the coasts of eastern Africa north of the equator, as it has been described previously (Gulf of Aden and Red Sea; Jereb and Roper 2005).Also, this species is currently caught mostly in waters with temperatures between 9º (Roeleveld et al. 1993) and 11 ºC (Augustyn et al. 1995) and the main contributing environmental variable for the distribution was temperature (~ 70%; followed by minimum chlorophyl a at ~ 20%), which further supports the absence in the areas previously described.Even if it is present there, projections from our models suggest that the habitat suitability of this species will shift southward, with loss of suitability in the north and off-shore (Agulhas and Namibian; Lipinski et al. 1992;Roeleveld et al. 1993) banks of South Africa by the year 2100, habitats that are likely to represent spawning grounds for the species (Jereb and Roper 2005), and increase slightly in the south-eastern coast of Africa.A similar pattern for S. bertheloti is indicated by our models, where however, due to the narrowness of the continental shelf and hence, the bathymetry range of S. bertheloti is limited, which may explain the high variance of values observed in our observations.Even so, according to our projections, if RCP 8.5 is to be observed, the distribution of the later species could become restricted to the Mauritanian coast by the year 2100, with the Gulf of Guinea and southward populations declining or even eradicated.
The European common cuttlefish (S. officinalis) is by far the most well studied species of the genus (Jereb and Roper 2005;Guerra 2006).Its cuttlebone implodes at depths greater than 150 m (Ward and Boletzky 1984;Sherrard Page 11 of 15 159 2000) and temperature limits of the species are known to be between 10 ºC and 30 ºC, with animals entering a lethargic stage below 10 ºC (Richard 1971;Bettencourt 2000) and above 30 ºC, due to oxygen limitation (Melzner et al. 2006(Melzner et al. , 2007)).Fittingly, the main environmental contributions for this species distribution models were seabed temperature (> 98%) and depth (~ 2%).Additionally, salinity being a minor contributor (< 0.1%) fits with the fact that, in comparison to other cephalopods, this species is less sensitive to salinity, being able survive in coastal lagoons and estuaries with salinities as low as 18 ± 2 (Boletzky 1983;Guerra and Castro 1988).Projections for the next 300 years suggest that this species may move into the arctic and colonize the American continent (Xavier et al. 2016)-a first for the sepiids.Our results partially support the previous findings, as all our projections point towards an increase in habitat suitability in the northern regions, which can potentially be used by this species as stepping stones into the Arctic basin and American continent.Indeed, our models for this species entail the lowest habitat suitability decrease across the whole group of target species, seemingly benefitting from future ocean changes.This may be in part a consequence of the coverage of the distribution of this species by our model compared to the other cuttlefish studied here (see D. braggi above for comparison, a less extensively surveyed species).By 2100, the losses in the south and gains in the north increase with RCP severity.However, this same trend is not clear by 2050, which can be explained by the variance and error of the modeling being greater than the differences between RCPs environmental values up to this year.Also, we should consider that this species, along with R. elegans and R. orbignyanum, presented the lowest changes in habitat suitability, further contributing to the odd results projected for the middle of the century.
Rhombosepion elegans and Rhombosepion orbignyanum (Pérez-Losada et al. 1996;Sanjuan et al. 1996;Khromov et al. 1998) are distributed throughout the eastern Atlantic Ocean, from the North Sea to South Africa (Pierce et al. 2010;Hastie et al. 2016).They are small-sized cuttlefish with slender cuttlebones with closely packed septa and modified sutures (Ward 1991), that allow them to reach greater depths than most other cuttlefish (~ 450 m).Hence, along with S. officinalis, they represent one of the species with the lowest levels of habitat suitability decrease for futures scenarios, having greater environmental flexibility and being able to occupy many habitats throughout their latitudinal range and are not so restrained as other species by their cuttlebone to depth limitations.However, the fact that both R. elegans and R, orbignyanum are more stenohaline and stenothermal than S. officinalis may explain the difference observed between these species (Jereb et al. 2010).Despite this, our results showed a similar contribution from temperature for the three species (> 95%) with slightly different contributions from salinity (S. officinalis = 0.05% vs. R. elegans = 0.43% and R. orbignyanum = 0.68%).Both rhombosepids seem to be on a downward trend in the Mediterranean Sea, especially R. orbignyanum in 2100, where they represent the most abundant cephalopod species (Jereb et al. 2010).However, a poleward shift in habitat suitability is to be expected, with the colonization of the northern parts of the English and Irish islands and R. elegans even venturing to Iceland.As such, much like the case for S. officinalis, this pattern may lead in the future to the colonization of the Arctic and American continent (Xavier et al. 2016) at the expense of a loss of habitat at lower latitudes, namely in the Mediterranean and North Sea.Further, for the same reasons noted for S. officinalis, this northern shift in habitat suitability is exacerbated under higher RCPs by year 2100 but not by 2050.The high variance observed in the latitudinal trends of our models for both R. elegans and R. orbignyanum can be attributed, in part, to the narrowness of the continental shelf across their African range (as previously mentioned for S. bertheloti).As there are very few pixels for these areas (narrow shelf), changes in a few pixels can easily be translated into large variations in the latitudinal trends observed.
The trends exhibited in the present study suggest a potential loss of cephalopod diversity at lower latitudes, which are currently the areas with the highest diversity (Rosa et al. 2019).Additionally, if cuttlefish are neritic 'canaries in the coalmine', as are other cephalopods for the broader ocean (Jackson and Domeier 2003), loss of habitat suitability may also occur for other organisms in the near-tropical regions, thus suggesting a poleward shift of many coastal communities.Also, since cuttlefish are demersal predators, their presence in new areas will put them in contact with new prey and competitors.Due to the nature of buoyancy control afforded by the cuttlebone, they are less energetically expensive with respect to depth than are fish (Webber et al. 2000) and hence, they can easily become aggressive predators in new shallow water trophic webs, where predators undergo daily vertical migrations (Guerra 2006), competitively displacing other species from their trophic roles and/or suppressing prey populations.Additionally, they are prey for many marine top predators (Jereb and Roper 2005;Guerra 2006).Hatchlings and juveniles can be targeted by benthic and demersal fishes, while adults are targeted mostly by marine mammals (Jereb and Roper 2005;Guerra 2006).If habitat suitability is translated into abundance, predators that rely on these species may also be forced to migrate or face a populational decline (Arkhipkin et al. 2015).
In many of the areas where these species are located, they are reliable marine resources and important sources of income, with most of the areas affected by habitat suitability declines located around the tropics and on the southern side of the Mediterranean Sea.If habitat suitability is translated into abundance, with these species often sustaining artisanal 159 Page 12 of 15 and industrial fisheries throughout their distribution, a considerable source of income and food security may be under threat (Jereb and Roper 2005).In fact, many of the fishing grounds of these species will be affected by habitat suitability decreases: De. australis in the banks around south Africa; S. bertheloti throughout the African coast; R. elegans and R. orbignyanum in the Mediterranean, where it often represents the most abundant cephalopod fisheries species; Ac. pharaonis on the Gulf of Thailand, Persian Sea, and Indonesian archipelago, where it locally is the dominant cephalopod species fishery; As. latimanus in the Indonesian archipelago and South China sea (Jereb and Roper 2005).This may lead to hardship for local fishing communities, especially those in the global south.

Table 1
Species presence points retrieved for the analysis Raw' refers to presence points collected from the repositories; 'Curated' is the number of points that were attributed to correct locations; 'Rarefied' is the number of points used for the analysis after picking only the points 20 kms far apart a Generic placement uncertain; b Taxon recognized as invalid in Lupše et al. (2023) a Sepia brevimana Steenstrup, 1875 Sepia hierredda Rang, 1835 Sepia officinalis Linnaeus, 1758 5551 1202 a Sepia prabahari Neethiselvan and Venkataramani, 2002 b Sepia ramani (Neethiselvan, 2001) a Sepia vietnamica (Khromov, 1987) a Sepia vossi (Khromov, 1996) a Sepia zanzibarica (Pfeffer, 1884) Total 15,923 13,232 Page 5 of 15 159

Indo-Pacific, including northern Australia
North of 30ºN and south of 25ºS, HS values for the future projections increased with latitude, time, and RCP severity (Fig.3F).For Ac. pharaonis, projected decreases of HS were observed between latitudes 30ºN and 20ºN (southern coast of Oman and the Persian Gulf; Figs.1, 3I), between the latitudes 20ºN and the equator (Andaman Sea, Gulf of Thailand, Ganges delta and the west coast of India; Figs.1, 3I) and south of the equator/10ºS (western coast of Australia, and the Timor and Arafura Seas; Figs.1, Australian coastSimilar projections were observed for As.apama (mHS: 46.31%; TSS: 98.14%) and D. braggi (mHS: 55.26%; TSS: 99.60%), namely showcasing a marked decline in habitat suitability (HS) near the northern limits of their distribution.Moreover, a decline in HS was also observed for Bass Strait and the coast of southern Australia (Fig.1A, D), corresponding to negative peaks in habitat suitability at 32.5ºS (all RCPs), 36ºS (all RCPs except 2.6) and 40ºS (all RCPs) (Fig.3A, D), as well as north of 30ºS for As.apama.Increases of HS (of half the amplitude of the negative peaks) appear on the southern coast of Tasmania, south of 42ºS (all RCPs) (Figs.1, 3A, D).Furthermore, in 2050 RCPs' mHS for these species followed a different pattern than in 2100 (which had lower values with severity of RCP), where RCP 4.5 had lower values than 6.0 (or even lower than 8.5 for D. braggi).

Table 2
Species mean habitat suitability (%) under different climate change scenarios

Table 3
Species SDMs' TSS and environmental variable contribution (%) to the model