Climate change scenarios predict global sea surface temperature (SST) to rise throughout the 21st century (IPCC, 2013), and marine biota is expected to respond to this warming trend by shifting their geographical (Pinsky et al., 2013; Jones and Cheung, 2014) and bathymetric ranges (Dulvy et al., 2008). Coastal marine systems are currently exposed to warming at a much faster rate than many other ecosystems (Harley et al., 2006) and as several species already live close to their thermal tolerance limits (Rosa et al., 2014), an increase in extinction rates of many marine organisms is expected (Lasram et al., 2010).

Cephalopods are known to play an important ecological role in marine trophic webs (both as prey and as predators), and many species are commonly defined as “keystone” species due to their strong influence on ecosystem dynamics (Rosa et al., 2013a,b). Thus, changes in cephalopod abundance can have a mixed impact on marine communities and fisheries, contributing to changes in their predators and prey abundance (André et al., 2010). In addition, given their short life spans and rapid growth rates, cephalopods are expected to respond faster than other marine species to changes in environmental conditions, making them good indicators of environmental change (Pierce et al., 2010). It is also important to note that these mollusks are a significantly growing component of global fisheries, with landings increasing steadily from the 1950s to reach about four million tons annually over the last decade (Doubleday et al., 2016).

The productivity of marine fisheries is likely to be affected by the alteration of ocean conditions such as water temperature, ocean currents, and coastal upwelling, as a result of climate change (e.g., IPCC, 2014; Lam et al., 2016). Such changes in ocean conditions may affect primary productivity, species distribution, community, and food web structure, which have direct and indirect impacts on the goods and services provided by marine ecosystems, which will have direct implication for the welfare of human society (FAO, 2018).

A central question in studies exploring the effects of climate change on biodiversity is how will changes be measured and characterized. We assumed that species respond individually to environmental changes and model distributions of individual species one at a time (Guisan and Thuiller, 2005) or that the distribution of species can potentially be influenced by the distribution of other taxa, using community-level modeling strategies instead (Gotelli et al., 2010). Beta diversity describes the extent of compositional change in the community between sites and also helps to reveal the assembly mechanisms that drive these differences (Bishop et al., 2015). Most studies examine species temporal turnover (STT) (Hillebrand et al., 2010; Poloczanska et al., 2013; Cheung et al., 2015; Pecl et al., 2017; Lotze et al., 2019), but frameworks exist to assess changes in species turnover both in time and space (Almeida-Neto et al., 2011; Baselga, 2012). STT is a widely used metric to assess these changes in composition (e.g., Almeida-Neto et al., 2011; Baselga, 2012); however, as a measure of beta diversity equivalent to the Jaccard’s dissimilarity index (Anderson et al., 2011), it mixes two components in one metric: changes in assemblage composition caused by a process of species loss or gain (i.e., the nestedness component of beta diversity) and changes in assemblage composition caused by a process of species replacement (i.e., the pure turnover component of beta diversity). Following the study by Baselga (2010, 2012), Albouy et al. (2012a) proposed a strategy to fully comprehend the potential effects of climate change on species assemblages by analyzing changes in species richness and changes in species composition together and highlighted a bivariate mapping strategy to picture simultaneously the spatiotemporal trend of both processes.

Beyond species composition, another important issue is how to account for functional and phenotypic differences in multispecies assemblages. Quantifying the distribution of traits in a community or the relative magnitude of species similarities can give us a good measure of the assemblage functional diversity (Cadotte et al., 2011). Body size is considered a fundamental species trait and a good indicator of ecosystem functioning due to its relationship to several functional traits such as growth, reproduction, and mortality (Brown et al., 2004). Also, the body size is an easy and cheap way to translate several covarying traits into a single one (Woodward et al., 2005). Commercial fishing is known to constrain the body size distributions of marine populations, as most fishing gear is size-selective and preferentially targets large-bodied organisms (Myers and Worm, 2003). The increased temperatures associated with climate change are expected to disrupt large-scale patterns in body size distributions (Sheridan and Bickford, 2011; Cheung et al., 2012) and ecosystem functioning (Fisher et al., 2010). Still, the effect of climate change on marine populations is less studied than the effects of fishing, so the large-scale projections of climate-mediated changes in body size distribution are urgent.

In this study, we provided the first forecast of global changes in coastal cephalopod species richness under projected climate change. Using an ensemble of ecological niche models (ENMs) (Diniz-Filho et al., 2009), we projected changes in habitat suitability for 161 coastal cephalopod species to an end-century stabilization scenario. We then examined how spatial and temporal components of coastal cephalopod assemblage diversity are projected to change toward the end of the century. Finally, we inferred the potential effects on body size distributions and their potential impacts in global cephalopod fisheries.

In this study, we provided a first attempt to understand the global patterns of cephalopod diversity within the neritic realm and explore changes projected to happen due to climate change. We have shown that, in the baseline scenario, the hot spot of coastal cephalopod richness is found in the Central Indo-Pacific region, particularly in the East China Sea and in the Eastern Philippines ecoregions (Figure 1, left panels). Given the limiting number of octopus occurrence data obtained for this study (10 species), it might seem that these coastal cephalopod hot spots are driven mainly by the high diversity of squids and cuttlefish; nevertheless, this may not be true as it is known that many octopus species are endemic in this region (Jereb et al., 2016). The Central Indo-Pacific region is described as a biodiversity hot spot for many marine taxa (Roberts et al., 2002; Tittensor et al., 2010), and several authors have suggested different hypotheses, based on particularly rich environmental conditions or historical geological events, that might have promoted speciation processes and/or refuge in this area (Renema et al., 2008; Cowman and Bellwood, 2013; Leprieur et al., 2016). Our projections are in line with the results of a recent study (Rosa et al., 2019) exploring for the first time the global patterns of species richness in coastal cephalopods, showing that despite the lack of data for some groups or geographic areas (e.g., octopuses), our models provide a consistent picture of the current global patterns.

Regarding changes in cephalopod diversity, we found that projected losses of habitat suitability for species were more important within the tropical areas whereas gains were greater toward the poles (Figure 2, right panels). These findings are consistent with studies that revealed poleward shifts in species distribution within the 21st century (Burrows et al., 2011; Poloczanska et al., 2013), as well as predictions of shifts in the 21st century (Pereira et al., 2010; Jones and Cheung, 2014). In the tropics, marine animals tend to have their critical thermal tolerances close to environmental temperature limits (Tewksbury et al., 2008), making them highly sensitive to warmer temperatures. In terrestrial organisms, the physiological adaptation to heat seems to be generally impaired (Araújo et al., 2013), although these patterns have not been fully explored in marine environments. If the pattern is true for marine organisms, then moving to cooler habitats at higher latitudes would constitute the more viable adaptive strategy. Another alternative is for animals to seek deeper colder water in response to environmental warming (Dulvy et al., 2008). However, moving toward the deep ocean might be unsuitable for coastal cephalopods, since most of them are highly dependent on the complexity and diversity of neritic habitats for reproduction (Boyle and Rodhouse, 2005).

The predicted range shifts in cephalopod distribution combined with the range contractions projected for most species will result in drastic changes in species composition. Our results point to nestedness being the main component of cephalopod temporal turnover, with species replacement only predicted to occur in a few areas (mainly in the higher latitudes; Figure 2, left panels). These changes in species composition will result in a less diverse assemblage, as they are projected to occur mainly due to the loss of species, but it is also vital to understand the degree of functional redundancy maintained and to better predict the consequences on ecosystem functioning and resilience (Albouy et al., 2012a). This is particularly important in coastal ecosystems since anthropogenic pressure may act in synergy with climate change intensifying local extinctions (Crain et al., 2009).

The distribution of maximum body size observed in the baseline scenario is consistent with the temperature-size rule (Atkinson, 1994), with larger animals occurring at higher latitudes (and lower temperatures). Rosa et al. (2012a) already reported this pattern for the same taxonomic group in the Atlantic Ocean, and in this study, we observed it across the globe. According to model projections, there is a tendency of a slight increase in mean body size toward areas with higher losses in projected richness (Figure 3), suggesting that these losses are affecting predominantly the smaller species (better represented in tropical latitudes). This trend can mean good news for fisheries in these areas, as this might potentially lead to a mitigation of the negative effects of climate change as measured by the total capture (ton). In contrast, for higher latitudes, the tendency seems to be the opposite. Projected changes are more favorable to countries at higher latitudes since their fisheries are expected to benefit from the predicted poleward shifts in species richness. Despite the projected gain in cephalopod diversity toward the poles, fisheries may need to adjust to the predicted reduction in mean body size.

The effects of ecological change on cephalopod populations driven by overexploitation of fishery resources are still to be fully understood. Yet, one might argue that under the combined effects of intense fishing pressure and climate change, fish are likely to be poor competitors in relation to cephalopods since the latter display faster growth, higher reproductive rates, short life cycles, and voracious opportunistic predatory habits (Rodhouse, 2008; Rosa et al., 2013a,b). Cephalopod biomass has not yet replaced fish biomass in the landings, but looking at the continuing growing trend (Doubleday et al., 2016), this hypothesis has to be considered. However, within the global climate change context, our results show that there might not be suitable environmental space for cephalopod species to do so. However, our models do not consider the potential for rapid acclimation and adaptation (Munday, 2014), which could give cephalopods time and evolutionary opportunities to adapt to future changes.

Given the increasing number of cross-factorial studies showing the deleterious interacting effects of ocean warming and acidification on the development and physiology of marine invertebrates (Portner, 2008; Findlay et al., 2010; Byrne and Przeslawski, 2013), such as cephalopods (Rosa et al., 2014), together with the availability of several new marine variables in the IPCC (2013), pH was included in our models. It proved to have a great impact on species distribution, as pH was the climatic factor responsible for unsuitable environmental space in 55% of the times, whereas SST was only in 13%. This led to more pessimistic results, in terms of range contraction due to loss of environmental space, when compared with other similar studies (e.g., Cheung et al., 2009; Poloczanska et al., 2013). Also, our assumption of limited dispersal and the fact that ENMs only take climatic variables into account, when characterizing the habitat suitability, should be taken into considerations when interpreting the results. Furthermore, ENMs predict potential niches, not the actual distributions (Peterson et al., 2011), so it is likely that many areas projected to be occupied in the present and in the future might actually not be. These assumptions and limitations call for careful interpretation of the projected changes in species richness and composition. Also, the underrepresentation of octopus species in this study, due to both bathymetric limitation imposed and lack of high-quality occurrence data for this group, makes the representativeness of the Southern Ocean incomplete.

Other important factor when interpreting the projections made in this study is taking into considerations that the projections made in the areas with new environmental space (Supplementary Figure 2) or high disagreement between models (higher SD in Supplementary Figure 4) call for ponderation. There are several sources of uncertainty when dealing with ENMs projections: the uncertainty introduced by the data, which in this case has to be taken into account as the representativeness of one of the groups is incomplete; the uncertainty linked with statistical algorithms, reduced in this study by the use of ensemble forecasting approach (Araújo and New, 2007) and exclusion of low-performance models, as evidenced in Supplementary Figures 3, 4 and finally the uncertainty related to climate projections, also addressed by the use of multi-model ensemble (Mora et al., 2013) and represented in Supplementary Figure 1. Despite these methodological limitations, the efforts to reduce and access uncertainty were made and can be observed that the disagreements between model projections are higher for baseline scenario than for future scenario, which derive from the consistency across projections made for the future, i.e., a loss of habitat suitability in the tropical area and a gain in the higher latitudes, pointing that the projected changes due to new climatic conditions are likely to occur. Nevertheless, the first-order geographical tendencies of these projections reflect major climatic drivers of change thus being likely to be ecologically meaningful (Garcia et al., 2015).

In this study, we used stacked multiple single-species ENMs to estimate community structure and species richness. Yet, it has been suggested that, because they do not account for species interactions, these models tend to overpredict species richness (e.g., Guisan and Rahbek, 2011; Calabrese et al., 2014). Joint species distribution models (JSDMs) are an extension of standard correlative ENMs that allow multiple species to be modeled simultaneously while accounting for species correlations not explained by available environmental predictors. By accounting for interactions between multiple species, the expectation is that JSDMs might allow for more accurate predictions and are being increasingly adopted in the literature (Wilkinson et al., 2021).

The global scale and complexity of climate change impacts, and the uncertainty in regional climate and earth system projections (Frölicher et al., 2016) calls for improved resolution of regional climate processes. Recent high-resolution climate models have allowed coastal processes to be better resolved and, in some cases, to reduce regional model biases (Saba et al., 2016), but more needs to be done in the terms of data collection and in the integration of biotic processes in more refined models. The use of these data in a “hybrid” mechanistic-empirical approach (as proposed by Robinson et al., 2011) could give us more accurate predictions of what will happen in the ocean of tomorrow.

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

JB-P and RR conceived the manuscript. JB-P and FG analyzed the data. JB-P wrote the first draft. All authors contributed to the final manuscript.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.