Global Patterns of Coastal Cephalopod Diversity Under Climate Change

Introduction

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.

Materials and Methods

Species and Climate Data

We obtained polygons of the extent of occurrence (range filling) for 161 coastal cephalopod species (79 cuttlefishes, 71 squids, and 11 octopus species; refer to the list of species in Supplementary Table 1) from the Food and Agriculture Organization of the United Nations (FAO) (Jereb and Roper, 2005, 2010; Jereb et al., 2016) and converted them to presence point data in a 1°×1° latitude/longitude grid using ArcGIS (ESRI, 2006). We found limitations in the number of octopus species available to use in this study, as there were several important species (especially in the Southern Ocean) that were not included, due to the lack of accurate distribution information. Also, from a methodological point of view, we decided to restrict our analyses to the continental shelf (200-m depth limit), as we assumed that surface conditions provided by the climatic models were not maintained at deeper depths. In addition, species associated with, but not restricted to, the continental shelf are more likely to be affected by climate change (Rosa et al., 2012b). To avoid statistical bias in ENM fitting, 5 species were excluded from the analyses, corresponding to those with fewer than 20 record points over the study area (Wisz et al., 2008).

The patterns of marine species distribution are strongly influenced by bathymetry (Dambach and Roedder, 2011), so to reduce false positives in the presence data, we refined the extent of occurrence maps by clipping off areas with depths falling outside the bathymetric range of the species (Jereb and Roper, 2005, 2010). The bathymetry of the ocean was obtained from ETOPO2 (2010) and resampled to a 1°×1° latitude/longitude grid.

We used 30-year averages of 5 climate variables (i.e., SST, sea surface salinity, total chlorophyll mass concentration at surface, dissolved oxygen concentration at surface, and ocean surface pH) from Earth System Models (ESM) developed for coupled model intercomparison project - phase 5 (CMIP5). There were 21 ESMs from 15 climate centers in 9 countries that modeled at least one of the variables analyzed (Supplementary Table 2). For each model and variable, we used the period 1976–2005 from the historical experiment, to establish the baseline period, and the period 2071–2100, to define our future scenarios. A stabilization scenario was used in this study (Representative Concentration Pathway, RCP4.5), with CO₂ concentrations projected to increase up to 650 ppm by 2100 (Vuuren et al., 2011). This scenario was chosen as it is the one that projected the rise in surface temperature by the end-century closer to the +1.5°C increase targeted by the Paris Agreement (United Nations, 2016), so we considered it the most realistic at present.

Climate data were publicly available from the World Climate Research Programme (WCRP).¹ SST and surface seawater salinity have monthly frequency while the other three variables have annual data. We only extracted the first layer (i.e., surface) for chlorophyll, dissolved oxygen, and pH. All variables were interpolated into a common 1°×1° grid prior to calculating multi-model means (for more on the methods and uncertainty associated with the ensemble, refer to Mora et al., 2013). We estimated the multi-model variability by calculating the SD of model means among ESM per variable and time period (Supplementary Figure 1). There are several methods to ensemble ESMs, but the average of several coupled climate models is usually found to agree better with observations than any single model (Tebaldi and Knutti, 2007).

To determine the extent of environmental differences between baseline and future climates, the Multivariate Environmental Similarity Surfaces (MESS) analysis was performed, as proposed by Elith et al. (2010). For each cell, the degree of similarity between the new environmental conditions and those in the baseline period was computed (negative values represent dissimilarity). As models are less reliable when predicting outside their domain (Barbosa et al., 2009), we have to carefully interpret the results for those areas. These calculations were performed using the modEvA R package (Supplementary Figure 2; Barbosa et al., 2014).

To match the resolutions of species and climate data, all datasets were resampled in ArcGIS (ESRI, 2006) to the 1° grid used for species. Data processing and statistical analyses were performed using R software (R Development Core Team, 2010).

Ecological Niche Models

To constrain algorithmic uncertainty associated with ENMs, we implemented an ensemble forecasting method (Araújo and New, 2007). Models were fitted using seven different modeling techniques implemented in BioEnsembles (Diniz-Filho et al., 2009): (1) BIOCLIM, (2) Genetic Algorithm for Rule-Set Prediction (GARP), (3) Generalized Linear Model (GLM), (4) Generalized Additive Model (GAM), (5) Multivariate Adaptive Regression Splines (MARS), (6) Maximum Entropy (MaxEnt), and (7) Neural Network (NNET).

For each species, data were randomly partitioned into calibration (75%) and validation (25%) dataset, the procedure was repeated five times, maintaining the observed prevalence of species in each partition, and models for each species were fit and evaluated using the True Skill Statistic (TSS) (Allouche et al., 2006). The TSS is defined as sensitivity + specificity −1 and ranges from −1 to +1, where +1 indicates perfect agreement, and values of zero or less indicate a performance no better than random. It is a simple and intuitive alternative to area under the curve (AUC) in measuring the accuracy of species distribution models (Shabani et al., 2018). Randomly selected pseudo-absence locations from the background data were used with equal weighting for presences to maintain prevalence (Barbet-Massin et al., 2012). Numerous techniques for pseudo-absence selection have been developed over the years, and the reports of their effectiveness vary. In this study, we used an ensemble of seven different algorithms, which have different theoretical approaches and data requirements, so random pseudo-absences are interpreted as real absences, pseudo-absences, or background depending on the algorithm. This is a simple and widely used approach in ENM studies, but methods that include environmental profiling in a previous step have been receiving positive feedbacks (e.g., Senay et al., 2013; Iturbide et al., 2015). Weighted median consensus forecasts were computed (Marmion et al., 2009), and models performing poorly (with TSS values ≤ 0.5) were excluded from the final ensemble (Supplementary Figure 3). The TSS scores were interpreted according to the accuracy classification scheme described by Landis and Koch (1977): TSS > 0.8 excellent; 0.6 <TSS < 0.8 good; 0.4 < TSS < 0.6 fair; 0.2 < TSS < 0.4 poor; and TSS < 0.2 no predictive ability. Consensus projections were built using 100% of the data, as data partitions have been shown to add significant uncertainty to forecasts (Araújo et al., 2009). We estimated the multi-model variability by calculating SD among the seven different modeling techniques included in the final ensemble. The SD was calculated over the presence-absence three-dimensional matrix used for the analyses in each time period (Supplementary Figure 4).

Using projected future climatic conditions, we estimated changes in the geographic location of environmental niches for each species. We imposed some limitations to dispersal, as we considered it unrealistic for a cephalopod species (or larvae) to move beyond important geographical barriers or disperse across major oceans to reach climatically suitable areas. Thus, we used an adapted version of the ecoregions suggested by Spalding et al. (2007) and only allowed species to move to an adjacent realm between time periods (Supplementary Figure 5 shows the realms used in this study). We also imposed bathymetric limitations so that species would not be allowed to colonize a cell if it fell outside the bathymetric range of species (Albouy et al., 2012b).

We calculated the potential distribution of each species for each time period and calculated the coastal cephalopod diversity by stacking individual distributional maps on top of each other and inferring species richness in each grid cell. We projected diversity for coastal cephalopods as a whole and individually for its three main groups, namely, cuttlefishes, squids, and octopuses. We then quantified the potential changes in cephalopod species richness as the difference between the future and the baseline period.

Then, we analyzed potential changes in cephalopod assemblage composition (species replacement vs. nestedness) between the two time periods. Using SST (as described in Albouy et al., 2012a) and its decomposition, we chose beta ratio (βratio) as a useful index to describe the relative contribution of each component (species replacement vs. nestedness) in the overall amount of STT. βratio values smaller than 0.5 indicate that species replacement is the main driver of STT, whereas values greater than 0.5 indicate that STT is mostly caused by nestedness (if βratio = 1, i.e., nestedness is the sole driver of STT; if βratio = 0, i.e., only replacement explains STT). In this study, we determined βratio as the ratio between the nestedness component of the Jaccard’s dissimilarity index (βjne) and Jaccard’s dissimilarity index (βjac) (Dobrovolski et al., 2011).

We also projected the changes in the distribution of mean body size of assemblages using the measurements of the maximum body size of species provided by FAO (Jereb and Roper, 2005, 2010; Jereb et al., 2016), by comparing current and projected future distributions of mean body size, under climate change. To limit the effect of extremely large body sizes and account for non-normal distributions, we used the logarithm of body size and applied the geometric rather than the arithmetic mean (Fisher et al., 2010). Body size analyses were averaged at the level of Exclusive Economic Zone (EEZ), so they can relate with the potential economic impacts of climate change on cephalopod fisheries [cephalopod global capture production per country (FAO, 2016) was used as a proxy for countries depend on the resource].

In addition, we quantified the potential effect of climate change on species range sizes by calculating the relative loss or gain (0.5 threshold in the probability of occurrence) of the sizes of the potential geographic ranges (measured as the number of cells occupied by a species) between the future and baseline period.

Results

Model Prediction Accuracy and Environmental Variable Importance

The predictive accuracy of the seven ENMs used in these analyses was classified from “fair” to “excellent” (according to the study by Landis and Koch, 1977 classification scheme), with a mean TSS criterion of 0.69 ± 0.08. The model with the lowest TSS was GARP (TSS = 0.45 ± 0.34), and the highest TSS was obtained with GAM, with TSS = 0.84 ± 0.09. The remaining ensembles performed at “excellent” levels with a TSS of 0.81 ± 0.11 (Supplementary Figure 3).

The variable responsible for the highest percentage (55%) of environmental space that was not observed under the baseline scenario was pH (i.e., ocean acidification). Chlorophyll was responsible for 21% of new environmental space, SST for 13%, oxygen for 8%, and sea surface salinity for 3% (Supplementary Figure 2).

Cephalopod Hot Spots and Projected Changes in Richness

From the latitudinal gradient analysis, we identified the zeniths of coastal cephalopod diversity for the baseline period at 8°N and 9.5°S latitude (in the Western/Central Indo-Pacific area), with a mean richness of 39 ± 15 species and 38 ± 13 species, respectively (Figure 1A). When looking at the major coastal cephalopod families, we identified a hot spot of cuttlefish richness around 12°N (19 ± 8 species), in the Bay of Bengal area (Figure 1C). Squid diversity peaks around 7.5°S (19 ± 4 species), in the central Indo-Pacific (Figure 1E). Finally, octopus diversity zeniths are identified around 40°N (3 ± 2 species) in the Mediterranean Sea and around 17°N (2 ± 1 species) on the wider Caribbean region (Figure 1G).

Figure 1. (A–H) Projected richness for the baseline and the end-of-century period, under IPCC AR5 RCP4.5 scenario, for coastal cephalopods and its main groups [cuttlefishes (n = 76), squids (n = 69), and octopuses (n = 11)].

Hot spots in the future scenario are projected to shift toward higher latitudes and present less diverse assemblages. Future coastal cephalopod richness zeniths are projected around 21°N latitude (Figure 1B; 13 ± 11 species) at the Bay of Bengal and the China Sea area and around 18°S (13 ± 6 species) in the Madagascar coast area. These latitudinal patterns are driven mostly by cuttlefish [Figure 1D; zeniths at 21°N (5 ± 6 species) and 19.5°S latitudes (4 ± 3 species)] and squid [zeniths at 20°N (8 ± 5 species) and 18°S latitudes (8 ± 3 species)] diversity, despite the high diversity of this last group all throughout the Indo-Pacific area [Figure 1F; zenith around 1°S (8 ± 5 species)]. Octopus hot spots are projected to slightly shift North but remaining in the same areas (the Mediterranean Sea and the Caribbean region), and still, they are projected to suffer a threefold decrease in richness (Figure 1H).

By the end of the century, 69% of the continental shelf is predicted to experience some loss in adequate environmental niches, whereas only 12% is predicted to gain (Figure 2, right panels). Habitat loss for coastal cephalopods is predicted to occur mostly in the tropics, with peaks at 10.5°S and 8°N latitudes, with a mean loss of up to 39 ± 15 species. The gains in habitat are predicted only for the northern latitudes above 70°N but with much less intensity than losses (1 ± 1 species).

Figure 2. (A–H) Predicted changes in richness and composition between the baseline and the end-of-century period, under IPCC AR5 RCP4.5 scenario, for the coastal cephalopod assemblages. Changes in diversity are quantified using delta richness (DRS) and changes in composition using the βratio index.

Under the future climate change scenario, the potential geographic range sizes of coastal cephalopods are projected to decrease for 149 species (95%), of which 24 (15%) are projected to completely lose their suitable environmental space, and increase only for 7 species (5%). The same trend is true for the main groups within Cephalopoda, with cuttlefishes shrinking their potential range in 95% of the species (with 15% projected to completely lose their suitable habitat), squids in 96% of the species (of which 14% are projected to disappear), and octopuses in 100% of the species are projected to reduce geographic range (with 30% projected to completely lose adequate environment).

Projected Changes in Cephalopod Composition

Nestedness contributed more than replacement in explaining the temporal pattern of cephalopod turnover (Figure 2B, mean β ratio = 0.64 ± 0.41). It was also the key contributor to the total amount of STT in 48% of cells (mean β ratio = 0.95 ± 0.13). The replacement was the key contributor in only 32% of cells (mean β ratio = 0.13 ± 0.14). For the remaining cells in the continental shelf, the βratio was not calculated since there were no changes in predicted habitat suitability for all occurring species between periods. When looking at the latitudinal patterns of STT, the replacement is more intense in the northern latitudes, above 50°N, with squids being the main driver of this pattern. In all the remaining latitudes, nestedness is the main component of SST, with the highest values in the equatorial latitudes, driven both by the patterns of cuttlefishes and squids (Figure 2, left panels).

Current Patterns and Projected Changes in Body Size Distributions

Species body size patterns in the baseline scenario present a general trend of higher values toward higher latitudes, with a maximum mean body size of 5.91 ± 0.26 (log) cm around 53°S. Our models project decreases in maximum mean body size to occur in 37% of the globe, mainly in the higher latitudes above 50°N, and to increase in 44%, mostly in the intermediate latitudes (between 50°N and 30°N). Some of the areas projected to suffer a reduction in maximum mean body size are located near countries with higher dependency on the resource (e.g., Russia, USA, Chile, and Italy).

Discussion

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.

Figure 3. (A) Distribution of geometric mean body size (log) projected for the cephalopod group, in the baseline period. Lateral panel represents the mean body size predicted for the baseline (blue line) and for the future (red line) period. (B) Net differences in geometric mean body size predicted between the two periods. Inland shades of gray represent cephalopod global capture production per country (ton) [2016 data – FAO (2016)].

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.

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

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.

Funding

This study was supported by the Portuguese Foundation for Science and Technology (FCT) through a Ph.D. grant (SFRH/BD/51514/2011) to JB-P, project grants PTDC/BIA-BEC/103266/2008 and PTDC/MAR/0908066/2008 to RR, and PTDC/AAC-AMB/98163/2008 to MA and FG.

Conflict of Interest

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.

Publisher’s Note

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.

Acknowledgments

We would like to thank the Food and Agriculture Organization of the United Nations (FAO) for providing the spatial data of species and Dora Neto for the help with spatial data processing. We also thank Camilo Mora and Rollan Geronimo for providing the climatic data. MA also acknowledges the Integrated Program of IC&DT Call No 1/SAESCTN/ALENT-07-0224-FEDER-001755.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.740781/full#supplementary-material

Footnotes

1. ^ https://pcmdi.llnl.gov/mips/cmip5/

References