Projecting future climate change impacts on the distribution of the ‘Octopus vulgaris species complex’

Introduction Historically considered to be a single cosmopolitan species, the so called Octopus vulgaris species complex (OVSC) is now recognized to be a group of (at least) six cryptic species: O. americanus (in the west Atlantic), O. vulgaris (in the northeast Atlantic and Mediterranean Sea), O. aff. vulgaris (in the region of South Africa), O. tetricus (southeastern Oceania), O. sinensis (northwestern Pacific), and O. djinda (western Australia). The potentially different environmental preferences of this highly cryptic species complex may result in distinct consequences under future environmental conditions. Methods The present study employed species distribution models (SDM) using MaxEnt to investigate potential changes in habitat suitability and geographical distribution of the OVSC in the future (i.e., 2050, and 2100), across four representative concentration pathway scenarios (RCP-2.6, 4.5, 6.0, and 8.5, CMIP5). Results Differential responses were observed in the OVSC species analyzed. Specifically, O. vulgaris and O. tetricus exhibited a severe loss in distribution across their predicted range; O. americanus exhibited projected extirpation close to the equator, with limited expansion towards the poles; O. aff. vulgaris was projected to lose half of its current distribution; O. sinensis exhibited moderate losses, with projected increases in northern areas; and finally, O. djinda exhibited limited losses to its distribution. Except for O. sinensis, increasing RCP severity exacerbated changes in mean habitat suitability and projected distribution gains and losses. Discussion Ultimately, this study provides information on the potential biogeographical effects of marine climate change on a key worldwide ecological and economic resource to further disentangle the effects over each OVSC species, with the goal of assisting toward the sustainable management of octopus species at the global scale.


Introduction
Cryptic species complexes represent a particularly challenging topic in marine ecology. The term 'species complex' is employed when referring to a group of closely related species, historically classified as a single taxon, given the difficulty in their distinction through traditional morphology-based taxonomic methods (Knowlton, 1993;Bickford et al., 2007). Despite their typically high phenotypical proximity, which ultimately prevents their accurate discrimination, considerable genetic divergence can be observed (Espıńdola et al., 2016). Cryptic diversity is relatively common in the Tree of Life, with cryptic species occurring almost homogenously among metazoans (Pfenninger and Schwenk, 2007) and protozoa (Bass et al., 2007). Moreover, cryptic species can be found in all major terrestrial and aquatic taxa, across the entire latitudinal profile (Beheregaray and Caccone, 2007). Over the past decades, techniques such as molecular sequencing, DNA barcoding, and mitochondrial DNA (mtDNA) analysis have enabled the disentangling of such instances, revealing a vast array of cryptic species complexes (Hebert et al., 2003;Bickford et al., 2007). The cephalopod class is no exception, with well-known species complexes, including those of the big-fin reef squid, Sepiotheuthis cf. lessoniana (Lesson 1930) (Cheng et al., 2014), the pharaoh cuttlefish, Sepia pharaonis (Ehrenberg, 1831) (Farhadi and Anderson, 2021), the squid Lolliguncula sp. (Steenstrup, 1881) (Sales et al., 2014;Costa et al., 2021), Ommastrephes squids (Fernańdez-Álvarez et al., 2020); and the common octopus, Octopus vulgaris (Cuvier 1797).
Historically, O. vulgaris was considered a cosmopolitan species, inhabiting shallow benthonic waters in disjunct populations spread around tropical waters worldwide (Robson, 1932;Roper et al., 1984). Initially described as a single species with a rather broad range (occurring along the eastern North Atlantic and the Mediterranean Sea, in the Americas and the eastern South Atlantic, Oceania, and the Northwest Pacific Ocean), recent genetic and molecular analyzes have enabled researchers to discriminate genetically distinct species within the O. vulgaris complex (Söller et al., 2000;Leite et al., 2008;Amor et al., 2015;Gleadall, 2016). At the same time, the existence of geographically and genetically isolated populations also contributed to the inclusion of 'Type' species in the complex . To this day, the OVSC is known to be comprised of (at least) six Octopus species, including: i) O. vulgaris -previously named O. vulgaris sensu strictu and which occurs in the North Atlantic Ocean, from North Africa to the English Channel, and in the Mediterranean Sea (Mangold and Hochberg, 1991;Mangold, 1998;Rosa et al., in press); ii) Octopus americanus -occurring along the Atlantic continental shelves of the Americas (Avendaño et al., 2020); iii) Octopus aff. Vulgarisoccurring in the coastlines of South Africa and Madagascar (Oosthuizen and Smale, 2003;Oosthuizen et al., 2004); iv) Octopus sinensis, occurring in the north-western Pacific Ocean, in the region of Japan and the Eastern China Sea (Reid and Wilson, 2015;Gleadall, 2016); v) Octopus tetricus, occurring in the south-eastern Australian coastline and northern New Zealand (Guzik et al., 2005;Amor et al., 2017); and vi) Octopus djinda Amor, -formerly Octopus cf. tetricus, occurring in the south-western coastline of Australia (Amor and Hart, 2021). This species complex is highly valued worldwide, representing one of the world's most important cephalopod fisheries (Balguerıás et al., 2000;Sauer et al., 2021). Along with its high value, extensive exploitation is poised to negatively impact this group, with some populations exhibiting signs of depletion . Indeed, octopus catches have increased in the past decades (FAO, 2021), despite global cephalopod fisheries declining since 2014 (FAO, 2020). This increase in fishery pressures, together with poor taxonomic resolution and catch under-reporting at regional and even species level could potentially lead to overexploitation of octopus populations Sauer et al., 2021), as has been the case for certain species of the complex [e.g., O. sinensis (Gleadall, 2016)]. Likewise, small scale coastal fisheries are also becoming unsustainable in the long term, for example in Southern Europe, due to the growing interest in Octopus fisheries allied to depleting finfish stocks and rising prices (Pita et al., 2021). The difficulty implicit to species discrimination within the complex leads to increased difficulty when compiling data on each of these species. At the same time, population identification and stock discrimination preclude an accurate assessment of the status of the OVSC around the world (Sauer et al., 2021), raising concerns over the real impact of anthropogenic pressures on each of the complex's species. Additionally, within such pressures and given the potentially distinct habitat requirements between species, there is a particularly pressing need to disentangle the potential effects of anthropogenic climate change over these species.
Rising carbon dioxide (CO 2 ) concentrations in the atmosphere over the past centuries (Shukla et al., 2019;Masson-Delmotte et al., 2021) are responsible for significant changes to the earth's climate and the global ocean system (Pecl et al., 2017). Indeed, oceanic waters have been warming over the past decades due to the ocean's ability to absorb 90% of excess atmospheric heat trapped by greenhouse gases (Zanna et al., 2019). At the same time, increased absorption of approximately 30% of the excess atmospheric CO 2 has also been changing ocean chemistry, leading to decreasing pH levels (Gobler and Baumann, 2016). These changes in ocean chemistry and temperature have further contributed to increased ocean stratification and slowed and disrupted current patterns, leading to the expansion of areas with low oxygen (Gobler and Baumann, 2016). These phenomena are both consequences and drivers of ocean climate change and are set to yield a vast array of negative impacts on marine ecosystems worldwide . Among cephalopods, some lines of evidence indicate that this group may benefit from the impending shift in ocean conditions (Doubleday et al., 2016), although their responses are likely to be complex (Pecl and Jackson, 2008;Rodhouse et al., 2014). Their high phenotypic flexibility and consequent environmental plasticity (Liscovitch-Brauer et al., 2017), together with their 'live fast and die young' lifestyle (O'Dor and Webber, 1991) likely represent an advantage over other marine taxa when confronted with the various challenges of a changing climate. Indeed, the synergistic effects of warming and the overfishing of their predators and competitors, have been linked to the increased fitness of cephalopod species (Rodhouse et al., 2014;Doubleday et al., 2016). On the other hand, several experimental studies have shown potential deleterious effects of climate change for a wide array of cephalopod species, including octopuses (Repolho et al., 2014;Rosa et al., 2019), squid (Rosa and Seibel, 2008;Rosa et al., 2012;, and cuttlefishes (Rosa et al., 2013;Moura et al., 2019;Otjacques et al., 2020). At the same time, distribution shifts associated with climate forcing have already been observed for several octopuses (Ramos et al., 2014;Arreguıń-Sańchez, 2019;Ponce-Maŕquez et al., 2020), and other cephalopod groups (Golikov et al., 2013;Alabia et al., 2020;Oesterwind et al., 2022), suggesting a broad redistribution of cephalopod species over the coming decades as they seek suitable environmental conditions (Xavier et al., 2016;Rosa et al., 2019).
Given the relatively high potential for future climate change to induce differential changes to marine species, there is a need to accurately describe a species' ecological and geographical distribution. With the combination of biogeographic knowledge and modelling in climate change impact studies, it is possible to attempt the prediction of the future impacts on biodiversity and ecosystem health (Hannah et al., 2002). Species distribution models (SDMs) present a very useful tool in this regard. Indeed, its development and use have increased steeply over the past decades (Zimmermann et al., 2010), since these models offer a good framework for predicting changes in species distributions, across vast geographical spaces and regarding large species assemblages (Elith et al., 2006). These models take georeferenced occurrence data and environmental predictors for a defined geographical extent and establish the relationship between a species' occurrence and environmental conditions, defining a species' ecological niche and allowing researchers to perform projections on the potential changes in a species distribution ranges across time and space (Miller, 2010). While these models are bounded by a set of assumptions and limitations that must be taken into consideration during interpretation (Araujo et al., 2005;Araujo and Guisan, 2006;Heikkinen et al., 2006;Fitzpatrick and Hargrove, 2009), SDMs incorporating future climate change predictions are considered a very effective way to address some of the questions regarding climate change effects on biodiversity (Sinclair et al., 2010).
To this day, a growing body of literature has employed SDM frameworks with cephalopod species (Puerta et al., 2015;Alabia et al., 2016;Xavier et al., 2016;Boavida-Portugal et al., 2022). However, research employing this modelling approach to project potential biogeographical impacts of marine climate change on octopuses remains very scarce (Hermosilla et al., 2011;Lima et al., 2020;Ángeles-Gonzaĺez et al., 2021;Schickele et al., 2021;Boavida-Portugal et al., 2022), with few studies focusing on species of the OVSC and at this scale. In this context, the present study aims to evaluate the potential biogeographical impacts of future oceanic climate change on the distribution of the six OVSC species, by implementing an SDM workflow using MaxEnt modelling to predict present-time habitat suitability and species occurrence distribution and project these into two future periods (i.e., 2040/2050 and 2090/2100), across four Representative Concentration Pathway scenarios (RCP; 4.5,6.0,and 8.5,CMIP5). With this approach, the present study also aims to aid the correct management of octopus fisheries worldwide, particularly with regards to the six species studied. By identifying the most threatened areas of each species' distribution, it is possible to better inform policy makers and economic agencies to prevent the collapse of local or regional populations.

Collection and curation of occurrence data
To obtain the necessary geo-referenced occurrence data on the species of interest, the Global Biodiversity Information Facility (GBIF) [GBIF.org (24 May 2022)] (GBIF, 2022) and the Ocean Biodiversity Information Facility (OBIS) [OBIS.org (24 May 2022)] databases were surveyed. As a way of limiting the data retrieved in each database to only valid occurrences, a specific set of filters were used, positively selecting for those occurrences obtained through 'Human Observation' and 'Preserved Specimen', and filtering out duplicate observations, improper datum conversion points, missing (NA) values in either Latitude or Longitude, as well as rounded latitude/ longitude coordinates. When dealing with species complexes, there is the risk of misidentification of specimens with the original morphospecies. Online databases such as GBIF and OBIS, contain significant contributions from citizen-science databases (e.g., iNaturalist) and non-directed scientific surveys, which undergo different degrees of scientific validation (iNaturalist, 2022), but are still susceptible to frequent cases of species misidentification, mainly when dealing with highly cryptic species complexes. In this sense, since most retrieved occurrences were labeled as Octopus vulgaris, one of the most recent illustrations of the potential distribution of each species of the complex (Avendaño et al., 2020) was used to re-label each occurrence point. In the case of Octopus djinda, occurrence data from Amor and Hart (2021) was also included in the species' dataset, due to the relative lack of geo-referenced occurrences in either online database. For the present analysis, occurrences of O. sinensis from the Kermadec islands (Reid and Wilson, 2015) were excluded, and the extent of O. sinensis was limited to the Northwestern Pacific Ocean.
The compiled dataset was then curated, restricting all occurrences to the continental shelf area [i.e., from the surface to a depth of 200 m; (Laruelle et al., 2018)] and removing potential occurrences referenced on land. Since GBIF and OBIS do not include comprehensive records for the depth of each occurrence in most species, the occurrence data was converted into a spatial polygon object which was then used to extract a vector of the depth values at each occurrence's coordinates (using the function 'extract' from the package 'raster'), from a bathymetry raster layer obtained from Ocean Climate Layers for Marine Spatial Ecology (MARSPEC) (Sbrocco and Barber, 2013). The spatial polygon object was then converted back into a data frame and merged with the depth vector, and each species data frame was subset to exclude depths greater than 200 m. To remove data occurring on land, a second clipping of the occurrence data was performed by erasing all points outside a shapefile of the world's ocean bodies (downloaded from Natural Earth Data, https://www.naturalearthdata.com/). These restriction procedures were performed because SDMs must ideally restrict model calibration to accessible areas (Peterson and Soberoń, 2012). The curated occurrence dataset, together with the plotted occurrences for each species and the script used for data curation are supplied in the Supplementary Material (see Curated_species_datasets.zip). Dataset curation and the following analyses were all performed in R studio software (v.4.1.2) (R Studio Team, 2022). All scripts are available in the Supplementary Material.

Predictor variables
The predictor variables used in this study included one topographic variablei.e., bathymetryand three oceanographic variables: temperature, salinity, and current velocity. Concerning the oceanographic variables, both surface and maximum depth (benthic) layers were chosen, since octopuses are primarily benthonic, but can still be influenced by surface conditions. The choice of the oceanographic variables was primarily based on the availability of environmental predictors projected for the future periods (i.e., 2040-2050 and 2090-2100), and Representative Concentration Pathway scenarios (i.e., RCP2.6, 4.5, 6.0, and 8.5). The mean, max, min, and range for each predictor were obtained from Bio-ORACLE, which offers global geophysical, biotic, and climate layers at a common spatial resolution (5 arcmins) and a uniform landmask (Tyberghein et al., 2012;Assis et al., 2018). Bathymetry, in turn, was retrieved from MARSPEC at a resolution of 5 arcmins (Sbrocco and Barber, 2013) and used as a predictor variable, to incorporate the water column height into the spatial analysis.

Modelling
The species distribution models were built for MaxEnt modelling, using the 'megaSDM' package (Shipley et al., 2022). For this purpose, each predictor variable layer had to be reprojected to an equal-area projection (i.e., specifically the cylindrical equalarea projection -"+proj=cea +lat_ts=0 +lon_0 = 0 +x_0 = 0 +y_0 = 0 +datum=WGS84 +no_defs"), because conventional non-equal area projections have grids which vary in their area the further away from the equator, and MaxEnt randomly samples cells from the available geographic space, assuming cells of equal area in the entire extent of each predictor (Elith et al., 2011). To detect and prevent collinearity in the predictor variable list (i.e., a strong correlation between two or more variables), the function 'vifcor' from the package 'usdm' was used (Naimi, 2015). Specifically, the function tests the variables against each other until it finds a pair of variables with a maximum linear correlation greater than a previously specified threshold [i.e., in this case, 0.7 (Cohen et al., 2003)], excluding the variable with a greater variance inflation factor value, and repeating this process for all variable combinations. This was performed in two separate analyses, i.e., first to the surface predictors, and then to those from the benthos, to prevent instability in parameter estimation. The results from the 'vifcor' function with the remaining variable list for the surface and benthonic predictor stacks are present in the Supplementary material (see collinearity.rar), together with the correlation plots obtained using the 'ENMTools' package (Warren et al., 2021).
The SDM analysis followed the workflow described in the 'megaSDM' package documentation (Shipley et al., 2022). In short, the training area, where the occurrence and background points are located, and the study area, where the model will be projected, were defined using the functions 'TrainStudyEnv' and 'PredictEnv', which take specific raster stacks and manipulates them to standardize the present and future period's environmental data input, re-projecting, clipping, and resampling raster predictors when necessary (Shipley et al., 2022). The package then takes and manipulates the occurrence data and employs a series of measures to mitigate the inherent bias, typical of collected or downloaded occurrence data (Phillips et al., 2009;Boakes et al., 2010) and which decreases the overall accuracy of SDMs (Phillips et al., 2009;Beck et al., 2013;Varela et al., 2014). To mitigate environmental and spatial biases within the occurrence data, the package performs environmental filtering of the occurrence data by dividing the environmental values at each point into a pre-determined number of bins (n = 25 in the present study), and then selecting one point from each unique combination of bins, obtaining a subset of occurrence points that is filtered by the environment (Varela et al., 2014;Castellanos et al., 2019;Shipley et al., 2022). This method allows the removal of clustered or oversampled records, while still maintaining the range of environments in which a species was found (Varela et al., 2014). The number of occurrence records remaining after the environmental filtering process is present in Table 1.
Since the species occurrence data frames only feature presence data, there was the need to generate background points for each species. Background points are artificially created species occurrence points, which describe the environmental conditions of the training area, and allow the incorporation of true data or pseudo-absence data into the models (Peterson and Soberoń, 2012). The 'megaSDM' package uses a 'combined' method to generate the background points, which employs both random and spatially constrained sampling, each weighted by a user-defined threshold (i.e., 50% each for the present study) (Shipley et al., 2022). As such, half of the user-set number of background points (i.e., n = 1000 per species) was randomly sampled from the entire study area (Barbet-Massin et al., 2012), while the remaining half were sampled from within shapefiles buffered around each true occurrence point (with the radius of each buffer being proportional to the 95% quantile of the distance to the nearest neighbor). Using a combined method to generate background points reduces potential model overestimation of environmental suitability that typically arises in regions of greater occurrence point densities, usually in more easily sampled areas (i.e., spatial bias), and which random background point generation by itself does not consider (Lobo et al., 2010;Kramer-Schadt et al., 2013). At the same time, this method also reduces the susceptibility of extreme extrapolation errors and overfitting induced by spatially constrained methods, by assigning a percentage of the background points to be sampled randomly from the entire environmental extent (Radosavljevic and Anderson, 2014). During generation, the background points were also environmentally filtered, similar to the occurrence data, creating an even spread across available environmental space, while retaining its spatial weighting (Shipley et al., 2022).
The present habitat suitability and distribution of each species were estimated using the MaxEnt modelling technique, which employs maximum entropy methods and is particularly resilient when dealing with presence-only species records (Elith et al., 2010). This method also has a predictive performance that is consistently competitive with other high-performing methods (Elith et al., 2006;Feng et al., 2019). The 'megaSDM' package employs MaxEnt modelling using replication and the subsequent ensemble of the replicate model predictionsin the present case, a replicate number of 5 was used per species, with each model run using a different subset of occurrence points (Shipley et al., 2022). The area under the curve (AUC) value was employed to evaluate each model replicate and the final ensemble model. This was performed by taking the AUC value of each replicate and comparing it to the AUC of a null model, where multiple replicates of occurrence points were placed at random in the training area (Raes and ter Steege, 2007;Shipley et al., 2022). The 'MaxentProj' function removed all models with validation AUC values lower than a specified threshold (i.e., in this case, 0.70). The same function then projected all models onto the specified future environments, across all RCP scenarios, and created the ensemble of all replicate maps using the median value of each pixel, thus reducing the effect of outliers (Araujo and New, 2007;Shipley et al., 2022). The evaluation plots and tables for each species' models are supplied in the supplementary material (see Evaluation.zip). From the ensembles, the habitat suitability maps were obtained. To produce binary maps of probability of occurrence (0 or 1), a threshold value was applied to the continuous habitat suitability maps. This threshold consisted of the mean model TSS criteria of model evaluation, also named 'maximum test sensitivity and specificity' logistic threshold, which maximizes the specificity and sensitivity of the receiver operating curve (i.e., the ROC curve), and is particularly effective in presence-only data (Liu et al., 2005). All the model projections, including the habitat suitability ensembles and the binary maps, are supplied in the supplementary material (see Projections.zip).

Post-analysis
The post-analysis followed a similar framework to the one previously described by Borges et al. (2022). Also, since the OVSC species do not exhibit a global distribution, instead having their areas of accepted occurrence in specific continents or regions of the globe, the post-analysis was performed in subsets of the total extent of the model predictions, encompassing only the currently accepted areas of distribution. First, changes to each species latitudinal distribution (within their currently accepted range) were assessed by plotting each species latitudinal centroid (i.e., the arithmetic mean latitude for the species-occupied pixels) for each time and RCP scenario. Second, the latitudinal trends in mean habitat suitability were obtained, by converting each ensemble into matrix form and calculating each row's mean value. To obtain the projected changes in the latitudinal mean habitat suitability, the resulting present-day vector was subtracted to the 2050 and 2100 vectors, and the respective outputs were plotted along the latitudinal (y) axis. Afterwards, the binary maps of probability of occurrence were used to create a visualization of unidirectional range shifts (i.e., expansion and/ or contraction) and transitory fluctuations (i.e., range contraction followed by expansion, and vice-versa), per scenario, per species (Early and Sax, 2011) using the function 'createTimeMaps' (Shipley et al., 2022). The Time Maps for each species and RCP scenario are supplied in the Supplementary Material (see Projections.zip). Data regarding this analysis is present in the supplementary material (see Post_analysis.zip).

Variable contribution
The variable contribution for each species' ensemble model is presented in Table 2 (see Supplementary material for further details). Overall, Temperature (either benthic or surface levels) contributed the most for the determination of habitat suitability and species occurrence. Except for O. vulgariswhere benthic salinity maximum provided the highest contribution -all OVSC species exhibited a temperature-related environmental layer as their most contributing variable. In the top four most contributing variables, temperature layers dominated, with bathymetry and salinity occurred isolated (salinity mainly in O. americanus). Variables associated with Current Velocity, in general, were the least contributing variables (apart from O.

General patterns of habitat suitability
All six species of the OVSC exhibited a shift in their centroid of latitudinal distribution over time (Figure 1), despite relatively different responses between species and scenarios. Indeed, while O. vulgaris did not exhibit a clear trend over time and RCP severitywith the centroid barely moving until 2050 and with either a shift north (RCP-4.5 and RCP-6.0) or south (RCP-8.5) depending on the scenario ( Figure 1A) -the other species exhibited fairly clear trends. Specifically, the centroid of distribution for O. americanus exhibited a considerable northward shift for all RCP scenarios, with the magnitude of this shift increasing along RCP severity, particularly between 2050 and 2100 ( Figure 1B). For O. aff. vulgaris ( Figure 1C), there is a clear southward shift that is exacerbated over time in RCP-6.0 and RCP-8.5, while for the two less severe scenarios the centroid expansion contracts until the end of the century, relative to 2050. Akin to O. americanus, O. sinensis also exhibited a clear northward shift of its centroid until the end of the century ( Figure 1D), in a trend that is in general exacerbated with RCP scenario severityexcept for RCP-6.0, where the northward shift is smaller when compared to RCP-4.5. Similar to O. aff. vulgaris, O. tetricus exhibited the same southward shift in its centroid, with a considerable southward shift until 2050 which then contracts northward in RCP-2.6 and RCP-4.5 -although not returning to the present-day latitudeand the successive shift southward for the more severe RCP scenarios ( Figure 1E). The O. djinda ( Figure 1F) also exhibited a similar pattern with a shift southward.

Latitudinal trends in habitat suitability and occurrence distribution
In this section, the results of mean habitat suitability change across latitude and the plotted difference in projected occurrence distribution regarding present-day are presented in a species-byspecies fashion. In terms of global trends, O. vulgaris was projected to undergo considerable reductions in habitat suitability and consequently local extirpations in its projected future distribution, which were exacerbated with the increasing severity of the RCP scenario. Specifically, mean latitudinal habitat suitability is projected to decrease over most of its latitudinal distribution (in the northern hemisphere) until the end of the century (Figure 2), mainly along the North African and Mediterranean range ( Figure 3).
Changes in habitat suitability for this species are relatively similar across RCP scenarios until 2050 -the largest decreases in suitability occur generally over the same areas but are larger for the Mediterranean and North Sea latitudes (Figure 2, left). For the year 2100, there is a considerable decrease in latitudinal suitability for most of the species' range ( Figure 2, right), with the most extreme scenarios (i.e., RCP-6.0 and RCP-8.5) featuring the largest declines. Gains in mean habitat suitability were projected for a small latitudinal band in the northernmost latitudes, which increased between 2050 and 2100 for RCPs 4.5 to 8.5, moving northward in the two most severe scenarios ( Figure 2). The changes in mean habitat suitability induced the Projected changes to the centroid of latitudinal distribution (i.e., the mean latitude of the occupied pixels) calculated for 2050 (red circles) and 2100 (blue circles) and Representative Concentration Pathway scenario (RCP-2.6, 4.5, 6.0, and 8.  Changes in mean latitudinal habitat suitability for Octopus vulgaris between the present-day and 2050 (left side) and 2100 (right side), for each Representative Concentration Pathway scenario (RCP-2.6, 4.5, 6.0, 8.5; CMIP5), relative to the present-day baseline (i.e., the black dashed vertical line). Data to the right of the dashed line represents a relative increase in mean habitat suitability (red) and data to the left represents a relative decrease (blue). decrease in the specie's areas of occurrence ( Figure 3A). In RCP-4.5, mid-range latitudes exhibit greater losses between 2050 and 2100, compared to the previous scenario. In terms of projected distribution changes ( Figure 3B), there is a considerable increase in the areas of species extirpationagain, mainly along the coastline of Mauritania, in the Mediterranean Sea, and in the North Sea. Projections in RCP-6.0 emphasize the areas of distribution loss, mainly in the Mediterranean and the North seas. In the former, extirpation of the species until the end of the century expands in the Adriatic Sea and around Italy, while also appearing in the Balearic Sea ( Figure 3C). For the most extreme scenario, RCP-8.5, there was a severe increase in the areas of species distribution loss over most of Macaronesia and in the Mediterranean Sea ( Figure 3D). Northward, the species is projected to lose considerable areas of its distribution in the Sea of Ireland, while some potential exists for expansion in the northernmost limits of its distribution.
In the case of O. americanus, the ensemble models predict an overall decrease in mean habitat suitability in tropical and subtemperate areas near the equator until 2050, with slight increases in the higher latitudes of both the Northern and Southern hemispheres (Figure 4). This pattern is exacerbated by the increasing severity of the RCP scenario. Specifically, for RCP-2.6, the ensemble model predicted decreasing mean habitat suitability until 2050 in lower latitudes ( Figure 4A), while in higher latitudes there are predicted increases until the mid of the century as well, although with localized decreases between 2050 and 2100 in some latitudes ( Figure 4B). In terms of distribution, this led to projected localized extirpation by 2050 in the tropical and subtropical areas between Cuba and the Yucatań Peninsula Time maps of the projected changes in species occurrence distribution for Octopus vulgaris, showing unidirectional range shifts in occurrence distribution [i.e., range expansion (in red and orange) or contraction (in dark and light blue)] as well as transitory fluctuations (i.e., range contraction followed by expansion (pink), and vice-versa (purple)] across time (i.e., 2050 and 2100) and RCP scenarios (RCP-2.6, 4.5, 6.0, 8.5; CMIP5).
in the Gulf of Mexico, southward over the coastline of Brazil ( Figure 5A). North of the Yucatań Peninsula, there is some projected localized distribution loss by 2050, together with areas of transitory fluctuation of contraction by 2100 following expansion until 2050. In the northernmost regions of its distribution, O. americanus is also projected to undergo a considerable gain in distribution.
In RCP-4.5, the tropical and subtemperate latitudes exhibited an increased loss in mean habitat suitability by mid-century ( Figure 4C), with some areas also showing further losses between the present-day and 2100 (Figure 4D). At the same time, mean habitat suitability until 2100 in the northernmost latitudes increases further than in the previous scenario. This trend translated into an expansion of the areas of extirpation by 2050 from the region Brazil to the seas between Cuba and the Bahamas ( Figure 5B), and in the Yucatań Peninsula. In this area of the Gulf of Mexico, there is also projected distribution loss by the end of the century in the seas to the west of the Yucatań Changes in mean latitudinal habitat suitability for Octopus americanus between the present-day and 2050 (left side) and 2100 (right side), for each Representative Concentration Pathway scenario (RCP-2.6, 4.5, 6.0, 8.5; CMIP5), relative to the present-day baseline (i.e., the black dashed vertical line). Data to the right of the dashed line represents a relative increase in mean habitat suitability (red) and data to the left represents a relative decrease (blue).
Peninsula, and in the region of the Bahamas and the continental shelf near the United States. In RCP-6.0, there was a considerable relative increase in mean habitat suitability in the northernmost latitudes of this species' distribution ( Figures 4E,  F), together with some localized increases in higher latitudes in the southern hemisphere. However, most of the species distribution extent in the tropical and subtemperate bands is again considerably reduced. This led to an increased decline in O. americanus distribution in these areas, particularly in the Greater Antilles and in the Gulf of Mexico, where the models projected considerable losses between 2050 and 2100 ( Figure 4C). In higher latitudes, the increase in mean habitat suitability leads to increased areas of potential expansion in the northern and southern limits of the species distribution. Finally, RCP-8.5, the changes in mean habitat suitability were particularly exacerbated regarding the previous scenarios ( Figures 4G, H) leading to a potential expansion of the areas of distribution loss from the eastern Brazilian coastline, along the entire Central American shores by 2050, and northward along the coastline of most of the US. At the same time, expansion in the offshore areas of Canada, in the Gulf of Saint Lawrence and to the east of Newfoundland, and southward in Uruguay and Argentina, suggests an exacerbation of this shift towards the poles ( Figure 5D).
Regarding O. aff. vulgaris, there was a considerable decrease in mean latitudinal habitat suitability relative to the present-day, projected for most scenarios ( Figure 6). Indeed, until 2050, and save for localized suitability gains in the southernmost latitudes in RCP-2.6 and RCP-8.5 (both of a very small scale compared to the decreases), most of the latitudinal extent for this species is projected to undergo a considerable loss in suitability (Figure 6, left). Until 2100, however, RCP-2.6 and RCP-4.5 exhibit decreases in habitat suitability across most of the species' distribution which are relatively smaller compared to 2050. This relative increase in mean habitat suitability between 2050 and 2100 is also followed by a localized gain in suitability relative to the present-day in the southernmost latitudes, for both scenarios (Figure 6, right). For RCP-6.0 and RCP-8.5, however, the overall trend of decreasing habitat suitability regarding the present-day continued (Figure 6, right).
In terms of distribution projections (Figure 7), an oscillating decrease in mean habitat suitability in RCP-2.6 is observed, with a likely extirpation of this species from its northernmost areas of distribution -from the Quirimbas islands to the east of the Maputo Bay on the east-African shores, and sparsely along the western coastline ( Figure 7A). Considerable loss by 2050 is also observed along northern Madagascan waters. In these regions, there are also large extensions of areas projected to undergo transitory fluctuations, disappearing in 2050 but re-emerging by 2100, due to the oscillating values of mean habitat suitability. For RCP-4.5, the ensemble model projects a very different situation between the southern Indian Ocean and southern Atlantic Ocean African shores ( Figure 7B). Specifically, the models project some areas of fluctuation, with the species potentially disappearing in 2050 and remerging by 2100 in the northern region of Namibia. In the western areas of its distribution, offshore Namibia and South Africa, this kind of transitory fluctuations also occurs, alongside projected expansions into greater depths by the end of the century. In the southern Indian Ocean, the model projects O. aff. vulgaris to potentially Time maps of the projected changes in species occurrence distribution for Octopus americanus, showing unidirectional range shifts in occurrence distribution [i.e., range expansion (in red and orange) or contraction (in dark and light blue)] as well as transitory fluctuations (i.e., range contraction followed by expansion (pink), and vice-versa (purple)] across time (i.e., 2050 and 2100) and RCP scenarios (RCP-2.6, 4.5, 6.0, 8.5; CMIP5). disappear by 2050 from most of the Mozambican shores and from Madagascar's northern and southwestern coastlines. In RCP-6.0, this pattern of loss becomes more evident. The progressive decrease in mean habitat suitability until the end of the century further exacerbated the projected loss in Madagascar and Mozambican waters ( Figure 7C), with new areas of extirpation between 2050 and 2100 appearing in the southern parts of these regions. Despite some localized distribution restrictions, no considerable changes occur over the southern South African continental shelf. On the south Atlantic shores, there are sparse projected losses in the distribution in the Namibian coastline, and the area of projected fluctuation also seen in RCP-4.5 is projected to occur southward, in the western South African shelf. RCP-8.5 exhibited a more severe scenario, with the species being completely extirpated from the coastal areas of Madagascar ( Figure 7D) -despite a small region in the southeast resisting until the mid of the century but eventually disappearing by 2100. Changes in mean latitudinal habitat suitability for Octopus aff. vulgaris between the present-day and 2050 (left side) and 2100 (right side), for each Representative Concentration Pathway scenario (RCP-2.6, 4.5, 6.0, 8.5; CMIP5), relative to the present-day baseline (i.e., the black dashed vertical line). Data to the right of the dashed line represents a relative increase in mean habitat suitability (red) and data to the left represents a relative decrease (blue).
The situation is the same for the waters of Mozambique, but the species remains relatively undisturbed in South African and Namibian waters, despite localized losses in the latter region.
For O. sinensis, the models predict a different response of this species since no progressive exacerbation of the observed trends is seen with increasing RCP severity (Figure 8). In RCP-2.6, mean habitat suitability increases until the mid-and end of the century for most of its latitudinal distribution in the northwestern Pacific Ocean, with relatively small decreases being observed only at the lowest latitudes of its extent ( Figures 8A, B). In this scenario, there are no considerable changes to the species distribution, with losses being virtually non-existent, and projected expansions localized to areas near the western coastlines of Korea in the Yellow Sea, and around Hokkaido and Sacalina islands in the Sea of Japan ( Figure 9A). In contrast, RCP-4.5 exhibits a different trend, with mean habitat suitability decreasing considerably in the southern latitudes, where present-day suitability is predicted to be highest ( Figures 8C, D). In the areas of higher suitability, a temporary increase in habitat suitability until 2050 is observed ( Figure 8C) but quickly followed by a decrease below present-day values until the end of the century ( Figure 8D). In terms of distribution, the projected map is dotted with areas of projected loss in distribution between 2050 and 2100, mainly in the outer areas of the predicted distribution in the East China Sea, in the Yellow Sea, and across the coastlines of Japan ( Figure 9B).
Due to the oscillating habitat suitability in higher latitudes, areas of temporary emergence followed by extirpation are observed in the eastern regions of the Yellow Sea and dotted alongside the northern areas of the Sea of Japan. Also in these regions, there are instances of projected expansions, although restricted, localized in the Yellow Sea and mainly on Hokkaido Island. RCP-6.0 presents a very similar situation to the previous scenario in terms of habitat suitability and projected distribution Time maps of the projected changes in species occurrence distribution for Octopus aff. vulgaris, showing unidirectional range shifts in occurrence distribution [i.e., range expansion (in red and orange) or contraction (in dark and light blue)] as well as transitory fluctuations (i.e., range contraction followed by expansion (pink), and vice-versa (purple)] across time (i.e., 2050 and 2100) and RCP scenarios (RCP-2.6, 4.5, 6.0, 8.5; CMIP5). changes ( Figures 8E, F, 9C). However, for RCP-8.5, the ensemble model predicted somewhat similar response to RCP-2.6. In this scenario, there is an exacerbation of the increase in mean habitat suitability until 2100 in the northern regions ( Figures 8G, H). To the south, however, mean suitability decreases progressively until the end of the century. This leads to a relatively small but present restriction in the species' distribution at greater depths in the East China Sea ( Figure 9D). However, considerable expansion is observed in the north and eastern region of the Yellow Sea, on the coastlines of Korea by 2050, and on the shores of China until the end of the century, as well as in the Sea of Japan.
Regarding O. tetricus, this species exhibited an overall trend of decreasing habitat suitability over its northern extent of distribution in Australia, with gains in its southernmost limits ( Figure 10). This pattern was consistent across RCP scenarios, and the projected changes stayed relatively stable between 2050 and 2100, despite a slight overall increase between 2050 and 2100 for RCP-2.6, and an exacerbation of mean habitat suitability decrease for RCP-8.5 for the same period. Changes in mean latitudinal habitat suitability for Octopus sinensis between the present-day and 2050 (left side) and 2100 (right side), for each Representative Concentration Pathway scenario (RCP-2.6, 4.5, 6.0, 8.5; CMIP5), relative to the present-day baseline (i.e., the black dashed vertical line). Data to the right of the dashed line represents a relative increase in mean habitat suitability (red) and data to the left represents a relative decrease (blue).
In terms of projected distribution, the oscillating pattern of mean habitat suitability in the north under RCP-2.6 ( Figures 10A, B) led to extirpation until 2050 in the areas near Shark Bay in the west, and a temporary potential extirpation in eastern Australia, near the Coral Sea ( Figure 11A), and permanent extirpation in the areas north of this region. In the South of Australia, slight increases in mean habitat suitability led to expanding distribution towards greater depths in the Great Australian Bay and Tasmania. In the westernmost area of this species occurrence, specifically in the Bass Strait, the models also predicted a fluctuating pattern of expansion in 2050, followed by extirpation in 2100, due to oscillating mean habitat suitability in this region. Under RCP-4.5, the deeper western and southern coastlines of Australia are projected to undergo localized extirpation ( Figure 11B), as well as the higher latitudes of the species' distribution in both western and eastern Australian coastlines. Also, on the southern Australian coastline, there is an expansion of the transitory fluctuation patterns compared to RCP-2.6, evidencing an expansion followed by extirpation until the end of the century in this region. Lastly, in New Zealand, there is an overall expansion of the distribution towards greater depths all around its distribution. Projections for RCP-6.0 show a similar pattern as RCP-4.5 for both eastern and western Australian coastlines and for the Great Australian Bay, projecting losses in distribution until 2050 and 2100 for these regions ( Figure 11C). However, increasing mean habitat suitability in the southern latitudes until 2100 ( Figure 10F) led to a potential expansion of the species in the southwestern region of Australia (i.e., near the Bass Strait), while transitory fluctuations expand in New Zealandin the Northwestern Time maps of the projected changes in species occurrence distribution for Octopus sinensis, showing unidirectional range shifts in occurrence distribution [i.e., range expansion (in red and orange) or contraction (in dark and light blue)] as well as transitory fluctuations (i.e., range contraction followed by expansion (pink), and vice-versa (purple)] across time (i.e., 2050 and 2100) and RCP scenarios 4.5,6.0,8.5;CMIP5). coastlines of New Zealand, fluctuations between contraction in 2050 and expansion until 2100 are projected from North Taranaki Bight until Auckland, while expansion followed by contraction is projected in greater depths in the Southeastern coastlines. RCP-8.5 was the scenario with the most exacerbated responses, following the considerable decrease in mean habitat suitability between the present-day and the year 2100 ( Figures 10G, H). Specifically, the models projected a considerable loss in distribution along most of the eastern, western, and southern Australian coastlines ( Figure 11D), apart from the region of Tasmania. The Great Australian Bay area also exhibited a considerable area of transitory fluctuation, with a temporary expansion in the area in 2050, followed by its extirpation until 2100.
Regarding O. djinda, the ensemble models projected decreasing mean habitat suitability across most of its Changes in mean latitudinal habitat suitability for Octopus tetricus between the present-day and 2050 (left side) and 2100 (right side), for each Representative Concentration Pathway scenario (RCP-2.6, 4.5, 6.0, 8.5; CMIP5), relative to the present-day baseline (i.e., the black dashed vertical line). Data to the right of the dashed line represents a relative increase in mean habitat suitability (red) and data to the left represents a relative decrease (blue). distribution (Figure 12), in a pattern that is exacerbated between 2050 and 2100 for most RCP scenarios. RCP-2.6 exhibited the smallest magnitude of habitat suitability loss, relative to the other RCP scenarios, and the smallest temporal change (Figures 12A,  B). This scenario also exhibited a relative increase until 2050, followed by a decrease until the end of the century (albeit still above the present-day baseline). These patterns led to relatively small changes in the species distribution, with some projected losses until 2100 occurring in the northernmost range of its distributionin the western Australian coastline ( Figure 13A), and sparse gains together with a large span of momentary expansion followed by contraction along the southern coastline.
In RCP-4.5, the loss in mean habitat suitability increases until the end of the century (Figures 12C, D), while the southernmost latitudes exhibit a considerable increase in suitability until the end of the century. This led to a larger area of projected losses in distribution, relative to RCP-2.6, in the western Australian coastline, together with increased areas of expansion in the South ( Figure 13B). For RCP-6.0, the projections are very similar to the previous scenario, although the projected extirpation in the western Australian coastline is greater between 2050 and 2100, with further expansion on the southern shores during the same period ( Figure 13C), due to a relative exacerbation of the changes in mean habitat suitability relative to RCP-4.5 (Figures 12E, F). Finally, in RCP-8.5, despite projected habitat suitability loss being considerably exacerbated ( Figures 12G, H) leading to extirpation across the entire western Australian coastline ( Figure 13D), there is considerably lesser expansion projected for the southern areas of its distribution (e.g., mainly localized around the Great Australian Bay).

FIGURE 11
Time maps of the projected changes in species occurrence distribution for Octopus tetricus, showing unidirectional range shifts in occurrence distribution [i.e., range expansion (in red and orange) or contraction (in dark and light blue)] as well as transitory fluctuations (i.e., range contraction followed by expansion (pink), and vice-versa (purple)] across time (i.e., 2050 and 2100) and RCP scenarios (RCP-2.6, 4.5, 6.0, 8.5; CMIP5).

Discussion
In terms of present-day distribution of each OVSC species, the ensemble models were able to make relatively accurate predictions, with small exceptions that are discussed below. Regarding future projections, the biogeographical response to climate change-associated changes in environmental variables was different across OVSC species, hinting at potentially different long-term consequences and management requirements for each species within the so-called OVSC and historically referred to as Octopus vulgaris. However, a general trend of decreasing mean habitat suitability with increasing RCP severity, with associated losses in the current distribution range, was observed across species.
For O. vulgaris, the ensemble models predicted a distribution very similar to that presented in Rosa et al. (in press), with the species occurring along west Africa and in the Macaronesia region (Adam, 1962;Goncalves, 1991;Sańchez et al., 2015); in Changes in mean latitudinal habitat suitability for Octopus djinda between the present-day and 2050 (left side) and 2100 (right side), for each Representative Concentration Pathway scenario (RCP-2.6, 4.5, 6.0, 8.5; CMIP5), relative to the present-day baseline (i.e., the black dashed vertical line). Data to the right of the dashed line represents a relative increase in mean habitat suitability (red) and data to the left represents a relative decrease (blue).
the Mediterranean Sea and along the Iberian Peninsula (Rees, 1956;Sousa Reis, 1985) and northward into the Northeast Atlantic Ocean and the English Channel (Massy, 1928;Rees, 1950;De Luca et al., 2014). However, its occurrence was also predicted in Iceland and in the North Sea, where it currently does not occur (Goud et al., 2019;Oesterwind et al., 2022). Notwithstanding, previous SDM analyses performed on O. vulgaris in the same region, by Schickele et al. (2021), have led to similar predictions over the eastern North Sea, despite slightly different results. In terms of future changes, the present study projected a decrease in habitat suitability for the Mediterranean and the eastern Atlantic, exacerbated with scenario severity, but did not predict increasing habitat suitability and consequent northward distribution gains in the Baltic Sea (Schickele et al., 2021). These differences likely stem from the environmental predictors and geographical restrictions included in the analysis. Indeed, to avoid over-parametrization, Schickele et al. (2021) only employed sea surface salinity and mean annual sea bottom temperature (including mean, range, and variance), and in the case of Octopus vulgaris, the analysis excluded salinity altogether. The present study, however, employed bathymetry and current velocity as environmental predictors, but also used both max depth (i.e., sea bottom) and surface layers for current velocity, temperature, and salinity. This approach aimed to accurately describe the octopuses' habitats, while still avoiding overparametrization. Notwithstanding, both studies projected increasing environmental pressure in the southern areas of this species' distribution range, which could lead to severe impacts on populations at this trailing edge and potential socioeconomic repercussions, given the relative importance of this resource for regional fisheries (Sańchez et al., 2015;Sauer et al., 2021). Indeed, the present study projected severe contractions in the Mediterranean Sea and along the Iberian and French coastlines, with minimal expansion. Together, these results point to a Time maps of the projected changes in species occurrence distribution for Octopus djinda, showing unidirectional range shifts in occurrence distribution [i.e., range expansion (in red and orange) or contraction (in dark and light blue)] as well as transitory fluctuations (i.e., range contraction followed by expansion (pink), and vice-versa (purple)] across time (i.e., 2050 and 2100) and RCP scenarios (RCP-2.6, 4.5, 6.0, 8.5; CMIP5). potential severe decrease in habitat suitability towards the end of the century for this species, with severity increasing alongside the emissions implicit to each RCP scenario.
Regarding O. americanus, its predicted distribution was in line with its current accepted distribution (Avendaño et al., 2020). Indeed, the ensemble model predicted its occurrence from the north of Argentina, along the continental shelf of Yucatan, and into the northwest coast of the United States and Canada. In terms of projections, this species exhibited a clear pattern of decreasing suitability in tropical and subtropical waters (mainly in the regions of the Yucatan Peninsula, northern Brazil, and towards Florida) and increases in higher latitudes (offshore Argentina and alongside the north-western US and Canada), both exacerbated with RCP scenario severity. This is a slightly different response to that projected for O. vulgaris, which saw its available environmental niche decrease considerably with time and RCP scenario. It is worth noting that other American octopuses have been analysed with an SDM framework, where similar trends were observed. For instance, O. insularis (Leite et al., 2008), which occurs in tropical and subtropical waters of the Central and South American continents, sharing most of its distribution with part of O. americanus', is projected to undergo an increase in suitable niche space in the tropical Atlantic, potentially expanding into the temperate northern Atlantic, temperate South America, and temperate South Africa (Lima et al., 2020). In the present-study, however, O. americanus is projected to undergo a severe range contraction in lower latitudes, in the tropical and subtropical Atlantic, suggesting that suitable niche space decreases in these regions. Alongside this severe range contraction at lower latitudes, the increase in habitat suitability at higher latitudes, particularly in the northern hemisphere, presents a potential poleward shift, suggesting the possibility of a split in the distribution of O. americanus between North and South America which could isolate both populations and lead to genetic divergence in the long term. In fact, O. americanus is mainly a subtropical and temperate species (O'Brien et al., 2021), linked to waters between the 18 and 25°C depending on the life stage (Bastos, 2018;Ángeles-Gonzaĺez et al., 2020). Juveniles and adults of this species typically occur in cooler waters up to 200 m in depth (Avendaño et al., 2020), and shallow waters associated with upwelling systems (O'Brien et al., 2021), which limits their occurrence in tropical areas. Therefore, the prediction of O. americanus in warmer shallow waters in the tropics (e.g., in the Caribbean region) is likely an overprediction introduced by the occurrence points utilized. As will be discussed below, one limitation of the use of online databases is the potential for misidentification in highly morphologically similar and cooccurring species [e.g., as is the case between O. americanus and O. insularis (O'Brien et al., 2021)]. Nevertheless, the observed projected trend of reducing mean habitat suitability in equatorial and low-latitude tropical areas is likely to be accurate, with increasing oceanic temperatures leading further constraining the species in these areas. Yet, it is also important to note that one may also argue that adult individuals may also move into deeper, cooler waters, together with the projected regional displacement; alongside, the weakening of regional upwelling systems, such as in Brazil, associated with climate change could also lead to further contraction in coastal areas (de Souza et al., 2020).
For O. aff. vulgaris, the predicted present-day distribution was (again) quite similar to the accepted distribution range, with the species occurring in the coastlines of Namibia in the Southeastern Atlantic Ocean, and in the Indian Ocean shores of Mozambique and Madagascar (Oosthuizen and Smale, 2003;Oosthuizen et al., 2004). Future projections predict a severe pressure on the eastern range of its distribution, in the Indian Ocean, with decreasing habitat suitability leading to losses in northern latitudes. Indeed, over time and for all scenarios except RCP-2.6, this species is projected to potentially be extirpated from most of Madagascar and Mozambican shores. This could lead to a reshuffling in local food webs and a significant impact on subsistence and recreational fisheries in the eastern regions of South Africa, where most of this species' exploitation occurs (Robertson et al., 1997;Oosthuizen and Smale, 2003).
The O. djinda has only recently been shown to indeed be a distinct species than O. tetricus, occurring along the southwest Australian coast from Shark Bay to near Cape Le Grand (Amor and Hart, 2021; in press a). The ensemble model was able to accurately predict its present-day distribution but encountered issues of overpredictionmainly in the southern Australian shores, in disjunct areas along the Great Australian Bight. This species is considered not to occur in this region, which has been associated with sharp drops in sea surface temperature that prevent paralarvae dispersal and settlement, maintaining the allopatric distributions between the now O. djinda and O. tetricus to the east (Amor et al., 2014). Regarding the areas of present-day distribution, the future projections exhibited the same trend of decreasing mean habitat suitability in lower latitudes, with a (limited) increase in southern areas. In this context, the species is projected to lose some of its northern distribution, potentially expanding into the Great Australian Bight if conditions become more suitable over time. This pattern was consistent and exacerbated with increasing RCP severity. However, from RCP6.0 to RCP8.5, the areas of projected loss outpace those of potential expansion. Overall, the present results suggest a potential for southward expansion of O. djinda towards the shallow areas of the Great Australian Bay, since this species typically occurs in temperatures between 17-25°C and depths of up to 80 m (Amor and Hart, 2021). This could potentially lead to increased chances of overlap with other octopus species (e.g., Octopus kaurna, Stranks 1990) and top-down pressure in these habitats. At the same time, since O. djnda supports a highly productive fishery (Moltschaniwskyj and Hall, in press a), local fishermen communities may also be affected by these shifts.
O. sinensis was accurately predicted to occur over most of its present-day distribution. Indeed, this species is known to occur in the South China Sea (Gleadall, 2016), from the region of Taiwan (Lü et al., 2013;Reid and Wilson, 2015), northward into Japan in southern Hokkaido and the East China Sea (Sauer et al., 2021). Our projections translated into two positive (RCP-2.6 and 8.5) and two negative scenarios (RCP-4.5 and 6.0), with the least and most extreme emission scenarios both resulting in a considerable expansion of O. sinensis towards the most northern regions of its distribution (on the Sea of Japan and the Yellow Sea). This poleward expansion to newly suitable habitats would be beneficial to the species. On the other hand, the two middle-emission scenarios exhibited a considerable loss of habitat suitability and distribution, over most of its present latitudinal extent. Given that this species is of major economic importance in the southern waters of China and Japan (Sauer et al., 2021), this may lead to a loss of fishery grounds of this species. Indeed, in the 20 th century, O. sinensis was already overexploited in certain regions of Japan (Hamabe et al., 1976) and, despite current regulations, these synergistic pressures may undermine the sustainable use of this species in the region.
Lastly, O. tetricus also exhibited considerable overprediction of its area of present-day distribution. Indeed, this species was predicted to occur along most of the southeastern, southern, and southwestern shores of Australia, as well as most of eastern and northwestern New Zealand. However, its present distribution only encompasses the east coast of the Australian mainland, in shallow waters across the Tasmanian Sea, and in northern New Zealand (Amor et al., 2014;Amor and Hart, 2021;Moltschaniwskyj et al., in press b). Overall, until the end of the century, the ensemble models projected localized decreases in mean habitat suitability in northern latitudes, much as O. djinda, with consequent losses in distribution in these regions. For RCP-2.6, the species is inclusively projected to expand into the Great Australian Bight, together with increases into further depths in New Zealand and Tasmania. If this emissions scenario occurs, it is possible that this species could clash with other competitor species and further increase top-down pressure in these novel habitats. However, for the most extreme scenarios, together with shifts towards greater depths, the species is projected to be particularly pressured in its northern distribution latitudes in mainland Australia, with the Australian Bight becoming unsuitable for colonization (with the models projecting severe loss by 2100 for this area), and small expansions in the area west of the Bass Strait. Warming is known to significantly affect O. tetricus egg development speed at 25°C -the upper limit of their temperature range (Spreitzenbarth and Andrew, 2021). As such, with increasing temperatures over time, this species is particularly susceptible to negative effects imposed by warming, which would lead to a contraction of the species range in mainland Australia, pushing the species into being further pressured by other anthropogenic stressors, despite its recent poleward expansion into Tasmanian waters (Ramos et al., 2014). Indeed, in New Zealand, egg development for this species has been shown to cease above 21°C (Anderson, 1994), which further emphasizes the threat of warming in terms of reproductive fitness.
This study's results highlight the looming threat of marine climate change to ectotherm species, and namely in octopods. Indeed, temperature was one of the main contributing variables (either surface, benthic, or both) to determining habitat suitability in the ensemble models. The four RCP scenarios employed in this analysis encompass the wide gradient of possible climatic futures, mainly in terms of global temperature rise. The first scenario, RCP-2.6, projects the global temperature rise to stay below 2°C by the end of the century, requiring a global CO 2 emission decline until zero by the year 2100. The intermediate scenario RCP-4.5 projects peaking CO 2 emissions by the middle of the century, resulting in a 2-3°C increase in global temperatures. RCP-6.0 is a high emissions greenhouse gas scenario, resulting in a potential temperature increase of between 3 to 4°C. Lastly, the "worst case scenario" RCP-8.5 projects temperatures to increase by over 4°C due to continued increasing emissions during the 21 st century (IPCC, 2014;Schwalm et al., 2020). In the present analysis, temperature was one of the most contributing variables for the majority of the OVSC species, emphasizing its potential future impacts to these species' worldwide distribution. Indeed, temperature (and thermal stress) has been shown to promote a decrease in egg incubation time, hatchling size, weight at first sexual maturity, and average generation time (Andre et al., 2010). In this sense, it is likely that habitats featuring the greatest thermal pressure over time will significantly decrease in terms of habitat suitability, leading to potentially deleterious effects on octopus populations and promoting poleward shifts whenever possible. In this context, decreasing habitat suitability and potential distribution contraction in tropical and subtropical areas were to be expected, since the thermal tolerance of marine animals in these regions is closer to the environmental temperature limits in their habitats (Tewksbury et al., 2008;Nguyen et al., 2011;. Salinity was also particularly relevant, being the most contributing variable for O. vulgaris and the second most to O. americanus. Indeed, salinity gradients are of a high importance for the O. vulgaris complex (Hermosilla et al., 2011;Moreno et al., 2014;Iglesias et al., 2016), for instance in determining the preferred areas of recruitment. Changes to the salinity profile could, then, further condition the habitats that are suitable for these species. Likewise, with temperature increase, temperate and subtemperate areas become more suitable for these species, opening more habitats at increasing latitudes for these species to move into [i.e., poleward shifts (Burrows et al., 2011;Poloczanska et al., 2013). The present study revealed a potential poleward shift in O. americanus and poleward expansion in O. sinensis, while also projecting severe range contractions in lower latitudes for O. tetricus and O. aff. vulgaris, and to a lesser extent for O. djinda. This is in line with previous research (Xavier et al., 2016;Schickele et al., 2021;Boavida-Portugal et al., 2022), and contributes to the notion of cephalopod borealization (Xavier et al., 2016). This type of change in the distribution of species which are both important predators and preys can lead to significant ecosystem-level impacts which must be addressed. First, areas which are projected to undergo octopus range contraction could very well suffer a decrease in top-down control of benthic communities, while projected poleward shifts could increase this top-down control on lower trophic levels, in higher latitude areas (Allen, 1971;Colleter et al., 2012;Schickele et al., 2021). However, the present models do not account for adaptation, and recent abundance increases in commercially important cephalopods at a global scale suggest the existence of potential benefits under an already changing environment (Doubleday et al., 2016). It is still uncertain if this pattern of abundance increase will be maintained under further, more prolonged, and severe environmental changes, with distinct species bound to have different responses. Notwithstanding, suitable habitats are projected to be shifting for most of the OVSC species, as have been for coastal cephalopods (Boavida-Portugal et al., 2022), which emphasizes the need for further research, namely featuring field validation of said projections to clarify the potential road ahead for this marine group.
As with all modelling projections, there is the need to address potential caveats of the SDM analysis. First, all models were subject to some degree of over or under prediction in all species. The first issue is common, since SDMs assume that a species will completely occupy areas predicted to be climatically suitable, ignoring potential real-world limitations such as unaccounted environmental predictors, the existence of geographical barriers preventing dispersion, etc. (Araujo et al., 2005). Underprediction, however, is linked to the nature of the occurrence point data, in terms of sample size and spatial biases. Occurrence data retrieved from online databases, such as OBIS and GBIF, are unlikely to represent a complete sample of a species known distribution, since data in these databases is highly dependent on sampling effort, which is frequently geographically skewed (Beck et al., 2013). Also, regarding o c c u r r e n c e d a t a f r o m o n l i n e d a t a b a s e s , s p e c i e s misidentification may introduce some bias in the data retrieved. Notwithstanding, this issue is not likely to be prevalent in most species, as there is no co-occurrence with other highly similar species. This issue may be more pressing in the case of O. americanus however, the overall trends projected for the future are not invalidated. Finally, the relatively low sample size for O. djinda and O. sinensis (post-environmental filtering sample sizes of 19 and 28 respectively; Table 1), could also be introducing uncertainty in the predictions, since sample sizes under 50 to 100 occurrence points suffer from relative decreases in predictive accuracy, depending on the algorithm and modelling conditions used (Stockwell and Peterson, 2002). However, all model runs for each species were maintained by the 'megaSDM' package during evaluation, suggesting reliable performance. Another issue refers to the assumption of climatic niche conservation between native and non-native ranges, which could potentially introduce biases and lead to inaccurate predictions (Pearman et al., 2008). As was already mentioned, the high phenotypic flexibility and strong adaptation potential of cephalopods (Boyle et al., 1996;Liscovitch-Brauer et al., 2017) suggest that Octopus species could adapt to future climate change-related conditions. Notwithstanding, this issue is minimized by the reduced geographical extent of the present study's projections and the fact that predictions were primarily made within each species present-day distribution. In this context, future research should aim to increase sample size whenever possible, as well as increasing the ecological relevance of environmental predictors used. In this sense, the inclusion of dissolved oxygen, pH, and the distribution and ecological relationships between the OVSC species and potential predators and competitors would contribute towards more ecologically relevant models (Warton et al., 2015;Abdulwahab et al., 2022). Despite these limitations, the overall tendencies of the present study's projections reflect major climatic drivers of change, and thus are bound to be ecologically meaningful (Garcia et al., 2016).

Final remarks
The present study contributes to the body of research on the biogeographical impacts of climate change on cephalopod species, particularly regarding the impacts on species complexes which are often treated as single species. In this regard, the results of this study highlight the potentially differential responses exhibited by each of the OVSC species. While results suggest increased pressure and potential extirpation at lower latitudes for some species (e.g., O. americanus and O. aff vulgaris), widespread distribution contraction were also identified (e.g., O. vulgaris and O. tetricus), and the potential for poleward distribution shifts as well (e.g., O. americanus and O. sinensis). These different responses require that research conducted on the O. vulgaris complex be adapted to include the different environmental requirements of each species, as well as the different abiotic pressures they will be subjected to in the future. This is of paramount importance for the correct evaluation and management of the fisheries stocks of the OVSC species, and the sustainability of regional octopus' fisheries. The present study aimed to contributed towards fisheries management and conservation, by providing a set of future projections where the identification of potential areas of increase/decrease in habitat suitability and the consequent potential changes in each species distribution towards the end of the century are possible. By identifying such areas, it is possible to complement potential stock evaluation and stock management to predict increased environmental stress which could further increase anthropogenic pressure already present through fisheries exploitation.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

Author contributions
RR and FB conceptualized the study; FB was responsible for data gathering and curation; FB conducted the SDM analysis and post-analysis processing; FB wrote the original draft of the manuscript; FB, MG, CS, JP, and RR reviewed and edited the manuscript. All authors contributed to the article and approved the submitted version.

Funding
This study had the support of national funds through Fundacão para a Ciencia e Tecnologia (FCT), under the project LA/P/0069/2020 granted to the Associate Laboratory ARNET; the MARE grant MARE UIDB/04292/2020; and the