Environmental matching reveals non-uniform range-shift patterns in benthic marine Crustacea

Empirical and theoretical studies suggest that marine species respond to ocean warming by shifting ranges poleward and/or into deeper depths. However, future distributional patterns of deep-sea organisms, which comprise the largest ecosystem of Earth, remain poorly known. We explore potential horizontal range shifts of benthic shallow-water and deep-sea Crustacea due to climatic changes within the remainder of the century, and discuss the results in light of species-specific traits related to invasiveness. Using a maximum entropy approach, we estimated the direction and magnitude of distributional shifts for 94 species belonging to 12 orders of benthic marine crustaceans, projected to the years 2050 and 2100. Distance, direction, and species richness shifts between climate zones were estimated conservatively, by considering only areas suitable, non-extrapolative, and adjacent to the currently known distributions. Our hypothesis is that species will present poleward range-shifts, based on results of previous studies. Results reveal idiosyncratic and species-specific responses, with prevailing poleward shifts and a decline of species richness at mid-latitudes, while more frequent shifts between temperate to polar regions were recovered. Shallow-water species are expected to shift longer distances than deep-sea species. Net gain of suitability is slightly higher than the net loss for shallow-water species, while for deep-sea species, the net loss is higher than the gain in all scenarios. Our estimates can be viewed as a set of hypotheses for future analytical and empirical studies, and will be useful in planning and executing strategic interventions and developing conservation strategies.


Introduction
Empirical and theoretical studies suggest that marine species respond to ocean warming by shifting ranges poleward and/or into deeper depths (Poloczanska et al. 2013;Saeedi et al. 2017;Pinsky et al. 2019;Lenoir et al. 2020;Chaudhary et al. 2021). Moreover, in association with human activities (e.g., aquaculture and maritime trade), changes in ocean temperatures have aided expansion and settlement of species into places that would otherwise be inaccessible. This leads to geographic shifts and unprecedented rates of non-native species settlement, representing a significant threat to ecosystem function (Molnar et al. 2008;Ojaveer et al. 2015;Chan et al. 2019). Deep-sea ecosystems, which have long been thought to be extremely stable in terms of physico-chemical conditions, have shown growing evidence that they also may experience abrupt changes and climate-driven temperature shifts (Danovaro et al. 2001(Danovaro et al. , 2004, putting the most extensive ecosystem on Earth and its largest reservoir of biomass in danger (Sweetman et al. 2017;Folkersen et al. 2018). However, studies on the potential effects of climatic changes and species range shifts on this marine ecosystem remain scarce, and the effects on ecosystem functioning and potential biodiversity gain or loss remain almost completely unknown (Danovaro et al. 2001;Hoegh-Guldberg and Poloczanska 2017).
Crustaceans are one of the most dominant taxa of both shallow-water and deep-sea communities (Saeedi et al. 2019a, b;Saeedi and Brandt 2020). They are also among the most successful groups of aquatic alien invaders in the ocean, with life history, behavioral and physiological traits that favor a high rate of establishment and success as invasive species (Karatayev et al. 2009;Hänfling et al. 2011). Their impacts are often substantial due to the complex trophic role of many of these species, leading to cascading effects throughout invaded ecosystems, resulting in fast spread and costly management and/or eradication measures (Thresher and Kuris 2004). Therefore, the early detection of climate-driven range shifts and prompt implementation of appropriate management strategies are crucial (Fogarty et al. 2017). The application of environmental matching approaches, of which a key concept is the linkage of species distributions and environmental conditions, can assist in the exploration and forecasting of species range shifts in response to changing environments Hänfling et al. 2011). Insights gained from such analyses can help elucidate current patterns and provide a first approximation of distributional patterns of marine biodiversity due to anthropogenic climate change (Cheung et al. 2016;Robinson et al. 2017).
Concerning marine crustacean fauna, so far, environmental matching has been seldom applied within the context of climate change and its potential impacts on future distribution patterns (Poloczanska et al. 2013;Melo-Merino et al. 2020;Lenoir et al. 2020). Limited by the high cost and logistical and technological limitations in sampling, the marine realm does not allow easy definition of the full spatial distribution of many marine species, and there are large taxonomic biases in the study of species/diversity spatial patterns (Reygondeau 2019). The few macroecological studies on marine crustaceans have focused on distributional patterns of economically relevant, well-studied species, and have been restricted to well-sampled geographical locations (i.e., North Sea and other coast of the Atlantic continental shelf; Neumann et al. 2013;Melo-Merino et al. 2020;Lenoir et al. 2020). Fortunately, in recent years, distributional, demographic, genetic and climatic data have become increasingly available digitally through global initiatives, such as the Ocean Biodiversity Information System (OBIS, www. obis. org). As a result, unprecedented possibilities to search for general patterns and mechanisms in the global distribution of marine biodiversity have been created (Wüest et al. 2019).
To fill the knowledge gap regarding future climatic change effects on the marine crustacean fauna, and to provide an initial assessment of climate-related horizontal distributional shifts for crustacean species, in this study we use a correlative maximum entropy approach. We estimate the direction and magnitude of potential climate-related range shifts of 94 shallow-water and deep-sea crustacean species representing 12 orders, to both 2050 and 2100, for two representative concentration pathway emission scenarios (RCP 2.6 and 8.5). In addition, we discuss our results in light of species-specific traits that could characterize invasiveness. We hypothesize that both shallow water and deep-sea crustacean species will follow trends seen in other marine organisms by shifting their distribution ranges poleward. Our estimations can be viewed as a set of hypotheses for future analytical and empirical studies, and instrumental in planning and executing strategic interventions, and in developing conservation strategies.

Species selection
Species were selected to represent the diversity within Crustacea, and based on availability of biological information and occurrence data representing distributions. In total, occurrence records for 94 globally distributed benthic crustacean species were assembled, 78 shallow-water and 16 deep-sea species. Benthic species were targeted in order to avoid problems related with the characterization of three-dimensional structure of marine space, as benthic species are less affected by the disparity between environmental predictor values and occurrence point recorded (Bentlage et al. 2013;Duffy and Chown 2017). To confirm their functional group as benthic species, we crosschecked species taxonomic classifications with information available at the World Register of Marine Species (WoRMS; www. marin espec ies. org). To classify the species according to depth, we used the maximum depth associated with records in the examined databases. Species were categorized into two groups: shallow-water (0-500 m) and deep-sea (> 500 m), as 500 m is the depth at which seasonal variation in physical parameters (e.g., temperature) as well as the influence of sunlight becomes minimal, following the classification by the World Register of Deep-Sea Species (WoRDSS; http:// www. marin espec ies. org/ deeps ea/; Glover 2019) (Fig. 1).

Distribution records
Occurrence records were obtained from the Ocean Biodiversity Information System (OBIS; www. obis. org) and Global Biodiversity Information Facility (GBIF; www. gbif. org), and complemented with occurrence data obtained from our previous deep-sea expeditions (Malyutina and Brandt 2013;Brandt and Malyutina 2015;Malyutina et al. 2018;Fig. 1). Records from before the year 2000 were excluded to avoid uncertainties related to geocoding errors. All species names were matched against WoRMS and synonyms reconciled. Distribution records were visually inspected for suitability, and dubious records were corrected (e.g., reversed latitude and longitude fields), or removed following the protocol in Cobos et al. (2018) (Fig. 1a).
In total, we obtained 51,852 occurrence records, which were rarefied through spatial thinning using a 10-km distance to avoid problems derived from spatial autocorrelation, resulting in a final count of 12,885 records used to calibrate and create final models (for number of occurrences used to calibrate each species model see Supporting Information Appendix S1, Table S1.1). To reduce sampling bias, thinning distance was chosen considering the spatial resolution of variables (~ 9.2 km at the equator), and the effect in geographic clustering and effective number of remaining points after exploring shorter and longer distances alternatives (i.e., 5 km, 20 km, 50 km and 100 km). After thinning, we selected species with unbiased distributions, i.e. with a good representation of their known distribution, and with a minimum of 30 records, which is above the minimum number necessary to obtain models with good predictive power in ecological niche modeling implementations (Papeş and Gaubert, 2007;Wisz et al. 2008;Shcheglovitova and Anderson 2013;Proosdij et al. 2016;Galante et al. 2018; Fig. 1).

Fig. 1
Schematic figure representing the steps to create environmental models for shallow-water and deepsea crustacean species. The main steps included: data capturing; ecological niche modeling; variability and uncertainty assessment. R script was used to perform all steps; it is publicly available at: https:// github. com/ msimo es123/ Crust aceaE NM Then, the complete set of occurrences of each species was randomly split in two subsets, one set for model training (75%), set aside for creating candidate models to be evaluated with testing data, to evaluate model accuracy (25%). We used a random partition of data as we did not have records that could be considered independent for testing. Partitioning data with spatial blocks has been suggested to potentially benefit model calibration when they need to be transferred (Muscarella et al. 2014). However, using spatial blocks could derive incomplete characterizations of species niches in our case, as records for many of the species in our study were distributed in regions with highly contrasting climatic conditions. After all filters were applied, 94 globally distributed crustacean species were assembled, 78 shallow-water and 16 deep-sea species, with counts of occurrences ranging from 30 to 800 records (see Supporting Information Appendix S1, Table S1.1). To clean and manipulate the data, all steps were performed using the statistical software R 3.5.1 (core Team, 2018) and packages 'spThin' (Aiello-Lammens et al. 2015), "raster" (Hijmans et al. 2015) and "rgdal" (Bivand et al. 2015).

Environmental data
Environmental data were obtained at 5 arcminutes (ca. 0.083°, ~ 9.2 km at the equator) spatial resolution extracted from Bio-ORACLE (http:// www. bio-oracle. org/; Assis et al. 2018), including maximum depth benthic layers, and bathymetry (i.e., depth) obtained from the General Bathymetric Chart of the Oceans (GEBCO; https:// www. gebco. net/). Limiting environmental factors for marine species include temperature, salinity, ice concentration and bathymetry (Belanger et al. 2012;Goldsmit et al. 2017). Thus, considering biological relevance and data availability, models included mean bottom temperature, mean salinity, and mean ice thickness, mean current velocity, and bathymetry. For deep-sea species, ice-thickness was not included due to lack of ice in the deep sea. Further, to avoid overfitting and inflation of model accuracy with overly dimensional environmental space and collinearity among variables, we performed Pearson's correlation coefficient (r) analysis to examine the cross-correlation of the variables, wherein correlation between variables greater than 0.6 were excluded. Finally, to test the statistical relevance of depth to the environmental models, we created two environmental sets, without depth, 'set 1', and including depth, 'set 2', to test the statistical relevance of depth to the environmental models (Fig. 1b).
Following Assis et al. (2018), environmental layers for present conditions were produced with climate data describing monthly averages for the period 2000-2014, and future layers were produced for 2040-2050 and 2090-2100 by averaging simulation results of three atmosphere-ocean general circulation models (AOGCMs: CCSM45, The Community Climate System Model 4; HadGEM2-ES5, Hadley Centre Global Environmental Model 2; MIROC55, Model for Interdisciplinary Research on Climate 5), in an attempt to reduce the uncertainties among diverse AOGCMs (for more information on climatic layers see Assis et al. 2018). The future layers used to project the environmental models included scenarios for both years 2050 and 2100, and two representative concentration pathway scenarios (RCP), RCP 2.6, presenting a peak-and-decline scenario ending on very low greenhouse gas concentration levels by the end of twenty-first century; and RCP 8.5, which shows high greenhouse gas concentration levels. For additional information on quality-control of climatic layers see Assis et al. (2018).

Ecological niche modeling
Ecological niche models were generated using the maximum entropy algorithm (Maxent; Phillips et al. 2006). Most crustacean species, due to their planktonic larvae, are able to disperse passively across long distances with ocean currents (Cabezas et al. 2010(Cabezas et al. , 2013Slusarczyk et al. 2019), hence, calibration areas included a buffer of two decimal degrees from the occurrences of each species (Barve et al. 2011), representing areas that populations of the species could have already had potential to reach due to their dispersal capabilities and, at the same time, restricted to areas adjacent to regions that had already reported records (Peterson et al. 2014). Seeking to reduce model overfitting, we first calibrated models by creating 126 candidate models (per each species), with parameterizations resulted from the combinations of nine regularization multipliers (β: 0.1-1.0 at intervals of 0.2, 2-5 at intervals of 1), seven feature classes representing combinations of linear, quadratic and product responses; and two distinct sets of variables, without depth and including depth; Fig. 1c). Selected model parameterizations were the ones that resulted in significant models (partial ROC with E = 5%, 500 iterations, and 50% of data for bootstrapping), omission rates lower than a previously defined error rate (E = 5%; Anderson et al. 2003), and the lowest AICc value per each species (Akaike information criteria; Warren and Seifert 2011), in that order. The mean AUC ratio of each final model is also estimated, representing the mean of all ratios between partial ROC, AUC values from each iteration and the AUC value representing a random model. The AUC ratio varies from 0 to 2: values greater than 1 indicate that model predictions are better than null expectations. These ratios are calculated using the partial ROC analysis ) and the p-values resulting from these analyses are calculated considering the number of times the iterated AUC ratios are below 1. Omission rates and significance were measured using the testing data to validate models created with 75% of the data (training set), AICc was measured on models created with all data. The chosen predictors and parameters settings were used to create final models with 10 replicates by bootstrap, logistic outputs, and transferred to globally, under current and future environmental scenarios (2050 and 2100). Final model projections were created allowing extrapolation and clamping in Maxent.
Finally, loss and gain of suitable habitats were calculated by comparing the geographic projections of niche models in current and future scenarios. The comparisons were categorized as follows: (1) if current and future areas were not suitable, these areas were defined as stable non-suitability; (2) when current and future areas were suitable, these were defined as stable suitable areas; (3) when current was suitable and future not suitable, loss of suitable areas was identified; and (4) when current was not suitable and future was suitable, gains of suitable areas were identified. After classifying areas of non-suitability, stability, gain and loss for all species, we summed the overlapping area of gains and of losses separately. To calculate the total area of gain and loss, binary maps were generated by considering only values above zero as areas of gain or loss. Then, the area was calculated in square kilometers, using ArcGIS version 10.6.

Model uncertainty
We measured model uncertainty based on the risk of strict extrapolation resulting from projections to non-analogous conditions, and the degree of variability resulted in final model projections. Model variation was assessed and represented geographically considering the amount of variance in model predictions that came from replicates and emission scenarios (RCPs), following (Cobos et al. 2019a). To assess the risks of strict extrapolation we used the mobility-oriented parity metric (MOP, Fig. 1e; Owens et al. 2013), which consists of measuring similarity between the closest 5% of the environmental conditions of the calibration area to each environmental condition in the area of transference (Supporting Information Appendix S2, Fig. S2.1). Areas of projection with values of similarity of zero compared to calibration areas indicate higher uncertainty as suitability in those regions, derives from model extrapolation only, and caution is required when interpreting likelihood of species presence in such areas (Alkishe et al. 2017). Binary maps of MOP results were generated considering only areas with zero similarity as strict extrapolation areas (Fig. 1e). Then, we excluded areas of extrapolation from binary final models, using binary MOPs, creating final binary models with no extrapolation ('post-MOP projections'; Fig. 1f). All modeling processes, uncertainty and variability analyses were performed using the 'kuenm' package (Cobos et al. 2019b) in R 3.5.1 (core Team, 2018). R script used to perform all steps is available at: https:// github. com/ msimo es123/ Crust aceaE NM.

Geographic centroids
To assess how suitable areas for the species shifted in latitude owing to climate change, first we calculated the geographic centroid of reduced suitable areas on post-MOP projections for current and future scenarios (Fig. 1f). For this, we created a concave hull based on species occurrences, of 200 km distance, and masked suitable areas with the resulted polygons to ensure that only relevant areas were considered in calculations (i.e., avoiding suitable areas in the opposite hemisphere when the species was only found in the northern hemisphere, or areas of strict extrapolation and outside the buffered concave hull; Fig. 1f).

Centroid shift
To calculate potential distances of centroid shift, we measured the distance in kilometers from current to future geographic centroid of suitable areas (henceforth, suitability centroids; Fig. 1f). To assess direction (i.e., northern hemisphere, increase of latitude; southern hemisphere, decrease of latitude) we calculated the differences in latitude of current centroids in comparison to future centroids, and to asses shifts between climate zones, we located the position of the centroid latitude of different time periods, within ranges of climate zone (i.e., tropical zone: 0°-23.5°, subtropical zone: 23.5°-40°, temperate zone: 40°-60°and polar region: 60°-90°; Fig. 1f). For constraint of suitable areas, estimation of suitability centroids and shifts, analyses were performed using "raster" (Hijmans et al. 2015) and "ellipsenm" (Cobos and Osorio-Olvera 2019) R packages. Schematic figure of the entire process of modeling and post-modeling can be found in Fig. 1.

Attributes characterizing invasiveness
When available, we assembled information on reproduction and feeding mechanisms, which are features that in combination with range shifts, can increase the potential of a crustacean species to become invasive (Hänfling et al. 2011). Here, we defined invasive species as a colonizer species that can establish populations outside their native distributional ranges and that have potential to spread and negatively affect native ecosystems or local human-mediated systems (Lockwood et al. 2007). These features are considered advantageous as a species maximizes the propagule pressure (i.e., 'introduction effort'), and range of resources (e.g., prey types) available to newly settled individuals (Hänfling et al. 2011). The information assembled was extracted from literature, the WoRMS database, and supplied by the experts listed in the acknowledgments. Within reproduction, brooder and freespawning modes (i.e., indirect reproduction) were considered, and within feeding mechanisms eight subcategories were considered: detritivores, herbivores, carnivores, parasites, scavengers, filter feeders, and omnivores.

Model statistics
In total, 11,844 candidate models were tested, 3723 were statistically significant models meeting omission rate criteria, 4194 were statistically significant models meeting AIC criteria, and 328 were statistically significant meeting 5% omission rate threshold and AIC criteria (mean AUC ratio of 1.29; Supporting Information Appendix S1, Table S1.2). Mean AUC ratios were higher than 1.24 for most taxa, indicating performance better than random (Supporting Information Appendix S1, Table S1.3). For each species, the selected parameters were the ones that created models that performed significantly better than random expectations (p-value < 0.01) and met the 5% omission rate threshold (i.e., false negatives, leaving out known distributional area). For parameters chosen for each species see Supporting Information Appendix S1, Table S1.3.

Model parameter settings
Between both environmental sets tested, the non-inclusion (set 1) of depth resulted in better models for ca. 29%, represented by 19 shallow-water and 4 deep-sea species; while the inclusion of depth (set 2) was selected to create the final models of ca. 71% species, representing 59 shallow-water and 12 deep-sea species. For models without depth (set 1), temperature and salinity were the most relevant variables for both shallow-water and deepsea groups (Supporting Information Appendix S1 Table S1.3; Appendix S2 Figure S2.3), with salinity only slightly more relevant than salinity for shallow-water species (Median contribution: shallow-water: 24%; deep-sea: 30.6%). In models constructed with the inclusion of depth (set 2), depth was recovered as the most relevant for both shallow-water and deep-sea groups, followed by temperature and salinity (Supporting Information Appendix S1 Table S1.3, Appendix S2, Figure S2.3). For 25% of the species, the statistically best final models were those with β = 0.2 (shallow-water: 20; deep-sea: 2), with β = 5 for only 5% of the species. The majority of best models contained linear-quadratic features (25%), then linear-quadratic-product (21%), closely followed by linear-product (17%), quadratic (11%), and linear (11%).

Patterns of climatic suitability
Under current scenarios, the sum of binary maps of post-MOP projections between all species of shallow-water and deep-sea species indicated the highest climatic suitability was concentrated between latitudes 70°N-45°N (Figs. 2a, c). General patterns for shallow-water species climatic suitability change for 2050 showed a consistent high suitability north of 45°N that dipped within the equatorial region. For deep-sea species, it was also shown that there was a high suitability north of 45°N, followed by a dip and an uniform increase until − 60°S, followed again by a dip (Fig. 2d; Supporting Information Appendix S1, Table S1.4, Appendix S2, Figure S2.4). Gains of suitable areas, resulting from differences between current scenarios and future projections, were mostly concentrated above 60°N (Supporting Information Appendix S1, Table S1.5, Appendix S2, and Figure S2.5). Losses of suitable areas for shallow-water species were higher above 60°N and steadily decreased until a minimum between − 20°S and − 30°S, while losses were mostly concentrated above 60°N for deep-sea species (Fig. 2; Supporting Information Appendix S1, Table S1.5, Appendix S2, Figure S2.6). Intensification of scenarios was observed for projections to 2100 for RCPs 2.6 and RCPs 8.5 (see vertical density graphs on Fig. 2; Supporting Information Appendix S2, Figure S2.4). Net gain of suitability was slightly higher than the net loss for shallow-water species, except for the projection 2100 RCP 8.5, while for Fig. 2 Environmental suitability under current climate scenario, and predicted to 2050 RCP 2.6 and RCP 8.5 emission scenarios for shallow-water and deep-sea Crustacea fauna. a-b Show trends for shallow-water species, c-d for deep-sea species. Vertical graphs on the right side summarize the density of suitability trends observed under current and 2050 RCP 2.6, RCP 8.5 and 2100 RCP 2.6 and RCP 8.5 emission scenarios deep-sea species, the net loss was higher than the gain in all scenarios (for projected gains and losses for each species see Supporting Information Appendix S1, Table S1.5).

Distance and direction of potential range-shifts
Shallow-water species were expected to shift an average of 431 km between suitability centroids by 2050 RCP 2.6 (RCP 8.5, ca. 620 km), and ca. 435 km by 2100 RCP 2.6 (RCP 8.5, ca. 1300 km) (Fig. 3a). Deep-sea species were expected to shift ca. 90 km between suitability centroids by 2050 RCP 2.6 (RCP 8.5, ca. 110 km) and ca. 130 km by 2100 RCP 2.6 (RCP 8.5, ca. 180 km) (Fig. 3a). Among the species studied, eight shallow-water and three deep-sea species were expected to have a shift between suitability centroids at least two times higher than the average by 2050 RCP 2.6 ( Fig. 3a), and six shallow-water and two deep-sea species were expected to shift between suitability centroids at least two times higher than the average by 2100 RCP 2.6 ( Fig. 3a; for distance of shifts projected for each species see Supporting Information Appendix S1, Table S1.6).
Concerning direction of potential distribution shifts, 82% of shallow-water species with current distributional centroid in the northern hemisphere, were projected to experience a shift in their environmental suitability centroids northwards (i.e., increase of centroid latitude) by 2050 RCP 2.6 (86% by RCP 8.5) and 2100 RCP 2.6 (90% by RCP 8.5). At the same time, all shallow-water species with distributional centroids in the southern hemisphere, were projected to experience a shift in their environmental suitability centroids southwards (i.e., decrease of latitudinal centroid), in all scenarios. On the other hand, for the deep-sea species, 80% of species with current distributional centroid in the northern hemisphere, and all species with current distributional centroid in the southern hemisphere, show potential distributional shift of environmental suitability centroids southwards, in all scenarios (Fig. 3b-e; for plots showing direction and magnitude of shifts projected for Fig. 3 a Circle bar graph indicating distance and latitude increase or decrease of distributional centroid positions of 78 shallow-water and 16 deep-sea crustacean species projected for the year 2050 and two representative concentration pathway emission scenarios (RCPs) 2.6 and 8.5. Blue and grey colors represent increase and decrease of latitudinal shift recovered for each species for 2050 RCP 2.6 and 8.5. Yellow stars indicate species (12 spp.) expected to present a distance shift above the average in 2050. Unit: kilometers. b-c, Maps show suitability centroid shifts between climatic zones, projected for shallow water species for b 2050 RCP 2.6 and RCP 8.5, c 2100 RCP 2.6 and RCP 8.5, and for deep-sea species d 2050 RCP 2.6 and RCP 8.5, (c) 2100 RCP 2.6 and RCP 8.5 1 3 2100 see Supporting Information Appendix S1, Table S1.6 and Appendix S2, Figure S2.7). Regarding shifts between climate zones (Fig. 3b-e), shallow-water species showed a projected net-increase of species richness in polar regions and decline of species in the temperate zones, with stability in remaining zones. Deep-sea species showed no changes of species richness in climate zones through time (see Supporting Information Appendix S2, Figure S2.7

Model projections and uncertainty
Differences recovered between RCPs presented higher variation in results of shallowwater than in deep-sea projections, with an average difference between RCPs of 2050 of ca. 0.7% (SD = ± 13%) and for 2100 ca. 20% (SD = ± 58%) for shallow water, and ca. 0.3% (SD = ± 1.9%) for 2050 and ca. 0% (SD = ± 4.4%) for 2100 for the deep sea. Mobility-oriented parity (MOP) results indicated that strict extrapolative areas for shallow-water models under current scenario were concentrated above 70°N. In 2050 projections, extrapolation areas abruptly reduced above 70°N and increased toward the mid-Atlantic Ocean, Indian Ocean, and Oceania, intensifying in these regions by 2100 (Supporting Information Appendix S2, Figures. S2.8a-e). For the deep-sea models, strict extrapolative areas were permanently recovered above 70°N (Supporting Information Appendix S2, Figure S2.8fj). For shallow-water and deep-sea species models, variability was higher when above ca. 60°N, and when derived from replicates, in comparison with variability derived from different RCP scenarios (Supporting Information Appendix S2, Figure S2.9).

Discussion
The use of environmental matching approaches, such as ecological niche modeling (ENM), has proven in recent years to be a powerful tool for conservation and has greatly increased our understanding of potential changes and impacts due to anthropogenic climate change (Evans 2012;Yates et al. 2018). However, the models used to estimate global change impacts on marine organisms have been considered limited, mostly focused on local studies, and varying substantially in complexity and underlying hypotheses (Planque et al. 2011). Additionally, ENMs of marine organisms have been constructed mostly for shallowwater species, while deep-sea species remain largely neglected, and scarcely represented in the literature (Vierod et al. 2014). As a consequence, based on patterns mostly observed on the study of shallow water species, major losses of diversity and biomass in the tropics accompanied by gains of suitability in poleward regions are expected, while future distributional patterns of deep-sea species remain poorly known (Oliveira and Scheffers 2018).
To our knowledge, despite representing only ca. 0.3% of the total diversity found within Crustacea (Giribet and Edgecombe 2012), our study is the largest effort to investigate the differential dynamics between shallow-water and deep-sea species on a global scale in the face of future climatic changes. We recovered idiosyncratic and species-specific responses, with the extent and direction of suitability shifts depending upon the niche characteristics of each species (VanDerWal et al. 2012), and areas lost and gained varying between taxa (Supporting Information Appendix S1, Tables S1.5). Nevertheless, although the predicted responses were species-specific, the dominant signal recovered showed poleward range shifts, with declines of species richness at midlatitudes, shifting from subtropics to temperate zone and temperate zone to polar regions. These results are strongly supported by trends reported in previous studies, indicating poleward shifts from warmer waters into cooler waters (e.g., Poloczanska et al. 2013;Neumann et al. 2013;Lenoir et al. 2020). Further, a large predominance of north-westerly range shifts was also observed, as has been similarly recovered from long-term observations of benthic invertebrates in the North Sea (Hiddink et al. 2015). Such patterns have been speculated to be related to the ability of species to expand their distributional leading edge (at the cold, polar boundary) and disperse and settle into areas that were previously too cold to inhabit (Sunday et al. 2012;Hiddink et al 2015). However, movements out of the tropics were not observed. For shallow-water species in particular, long distance range-shifts were estimated (e.g., T. paschalis), but no substantial shifts reaching cooler climate zones were predicted. This could be linked to the higher thermal limit and narrow thermal breath of tropical species, resulting in lower acclimation response to temperature changes, making them more vulnerable to warming than their temperate counterparts (Vinagre et al. 2019).
In our results, we also noticed a difference in magnitude of distance shifts between shallow-water and deep-sea benthic crustaceans. While relatively larger potential shift patterns were recovered for shallow-water species (up to 431 km by 2050 RCP 2.6, and ca. 435 km by 2100 RCP 2.6), shorter distances were forecast for deep-sea species (ca. 90 km by 2050 RCP 2.6 and 110 km by 2100 RCP 2.6). Values estimated for the shallow-water species were not so different in comparison with the average distances as described in previous studies for marine species (i.e. 72 km/decade; Poloczanska et al. 2013). However, for deep-sea benthic invertebrates, such estimations remain poorly explored. The contrast in magnitude between groups is linked to differences in thermal habitats, due to the strong links of marine species to temperature gradients (Sunday et al. 2012;Yasuhara and Danovaro 2016;Chaudhary et al. 2021). Temperatures in the deep sea are inherently more stable across the range of latitudes and through time and therefore less susceptible to superficial climatic changes; while in shallow waters, temperatures are warmer closer to the equator and decrease gradually toward the poles (Assis et al. 2018). Extrapolation analysis (MOP) -used to calculate environmental similarity from a given pixel in a transfer time/region to those within the calibration region (Supporting Information Appendix S2, Figure S2.3) -showed that extrapolation risks will decrease in northernmost areas for shallow-water species, while for the deepsea species these areas seemed to be permanently concentrated and increasing in the north, above 70°N (Supporting Information Appendix S2, Figure S2.6). This reinforces that shallow-waters are increasing in temperature in higher latitudes through time, while deep-sea shows higher stability and therefore, are less likely prone to movement of species toward polar regions.
The differences in magnitude of distance shifts might also be related to the unbalance between net gain and loss of climatic suitable areas observed in our results (i.e. net gain higher than loss for shallow-water species, and losses expected to be higher for deep-sea species; see Supporting Information Appendix S1, Table S1.5); while availability of suitable areas could allow for longer range shifts, a lack of suitable areas could lead to shorter shifts. Previous studies corroborate the responses recovered here, with longer shifts recovered for shallow-water species (e.g., Poloczanska et al. 2016;Sunday et al. 2012;Morley et al. 2018), but relatively small to non-significant contractions of future horizontal environmental suitability for deep-sea benthic species Basher and Costello 2016;Costello et al. 2017). However, sampling effort in the ocean is not evenly distributed (Webb et al. 2010;Boltovskoy and Correa 2017;Menegotto and Rangel 2018;Saeedi et al., 2019a), leading to a large bias in under-explored areas, undescribed species and species known only from single records that are not amenable to biogeographic analyses. Hence, such patterns might be influenced by the relative rarity of studies on marine species and/or their limited geographic scale, hindering the visualization of global patterns. For instance, our dataset presents a considerable latitudinal and bathymetric bias, as it includes mostly species with ranges in the Northern Hemisphere (see Supporting Information, Appendix S2, Figure S2.2), and much less than half of the taxa examined are deep-sea species (N = 16). Despite these biases, our results offer a robust hypothesis regarding future distribution patterns, and we speculate that an equal sampling in Southern Hemisphere species would result in similar pattern estimations, based on results from MOP analysis, and on the patterns recovered for other benthic invertebrates (e.g., Hiddink et al. 2015).
The selection of key environmental variables in which to explore the niche of species is a crucial step in model design (Melo-Merino et al. 2020). The inclusion of too many environmental dimensions can cause model overfitting (Anderson et al. 2003); thus, as good practice, various biological or statistical criteria are used to select the best set of environmental predictors. In models without depth (set 1), temperature and salinity were the most relevant variables, with salinity only slightly more relevant than salinity for shallowwater species. These results are not surprising given that temperature is a key environmental driver of performance (Cheung et al. 2010;Poloczanska et al., 2013;Saeedi et al., 2017), and that for animals in which osmoregulation is the main physiological mechanism maintaining hydromineral homeostasis, such as Crustaceans, salinity plays a fundamental limiting role (Thabet et al. 2017). Moreover, in models constructed with the inclusion of depth (63 shallow-water species; 12 deep-sea species; see Appendix S1 Table S1.3), depth was recovered as most relevant for both groups, followed by temperature and salinity. This aligns well with the findings of other studies (e.g., Pearman et al., 2020;Gonzalez-Mirelis et al. 2021), which have indicated that depth along with temperature and salinity are relevant for defining benthic species distributions. It is relevant to note that, previous studies have highlighted the relevance of oxygenation, pH and food supply (or POC flux) in the maintenance of the ecosystem functioning (Sweetman et al. 2017). However, as marine organisms are directly influenced by their environment through temperature, oxygen, and food availability, other factors, such as carbonate chemistry, are also important but have been less often implicated in shifting species distributions (Pinsky et al. 2019).
With the rearrangement of marine biota, invasion risks are also predicted to increase in the future (Cheung et al. 2016). Of all alien invasive species, crustaceans are often a significant group, as they are great colonizers, establishing outside their native distributional ranges with the potential to spread and affect both native ecosystems and local humanmediated systems (Engelkes and Mills 2011; Hänfling et al. 2011). Hence, uniting functional traits related to invasiveness to distribution model results could provide valuable insights and linkages between species distributions and their underlying biological mechanisms. Here, we focused on collecting features related to reproduction and feeding habits that could facilitate the process of colonization and settlement. Reproduction in particular is a critical step in the invasion process, while the successful establishment of invasive species may hinge on their ability to reproduce across wide-ranging environmental conditions. Among the species modeled in our study, amphipods and decapods -the most prominent examples of invasive crustacean orders (Hänfling et al. 2011) -were predominant and were among species with highest projected climate related range shifts. Within this order, free-spawning reproduction with large numbers of eggs within seasons and having planktonic larvae, aids propagule pressure, facilitating the rate of geographical spread of non-native populations (Lockwood et al. 2007). This could be one of the traits facilitating the wide range of decapod species such as Thor paschalis, distributed within the Indo-West Pacific Ocean, ranging from the Red Sea, Arabian Gulf and Madagascar eastwards to Australia (Anker & De Grave 2016). The species' exact distribution range remains unclear (Al-Kandari et al. 2020), and new records are often reported (e.g., Chace, 1997;Anker & De Grave 2016); and in our results, a trend of range expansion is also recovered, being one of the species with highest projected range shift (ca. 3715 km by 2050 RCP2.6; distance of shifts projected to each species see Supporting Information Appendix S1, Table S1.7).
The orders Cumacea (shrimp) and Podocopida (ostracods) follow in order of magnitude of projected range shifts in our study, and are also characterized by spawning, but eggs are brooded by females until they reach maturity (Watling 2005). This reproduction strategy might facilitate invasion success under non-favorable conditions, such as changes in seawater quality and environmental characteristics at spawning grounds, enabling offspring to survive under harsh environmental conditions (Hänfling et al. 2011). In particular, ostracods can possess different combinations of sexual reproduction, parthenogenesis, brood care, and egg-laying (Horne et al. 1998), favoring rapid spread into new environments. The order also shows flexibility in dispersal mechanisms, being passively dispersed by wind, birds, fish, water insects and drifting algae (Teeter 1973), with resistant double-walled eggs suited to withstand dissection during dispersal (Teeter 1973), and are usually small-sized (0.5-2 mm; Xing et al. 2018). The combination of these features is deemed advantageous for widespread distribution and aid of colonization of new environments, including in association with the gain of climatic suitability recovered in our results, and could confer high rates of migration and establishment in new environments. For example, Rabilimis septentrionalis (Brady, 1866) shows one of the largest range shifts projected in the order (ca. 590 km 2050 RCP 2.6; see Supporting Information Appendix S1, Table S1.7), and this species is known to have a large physiological breadth, i.e., high tolerances to seasonal salinity and temperature variations (Carter 1988), which could favor range expansion and settlement into new environments.
Regarding feeding habits, species in our study reflected well the diversity found within Crustacea, including omnivores, filter-feeders, scavengers, parasites and carnivores. Omnivore species have the highest threat to invade environments, by feeding on any trophic level in the invaded food web (van der Velde et al. 2009), and therefore adding pressure to the system due to predation of native species (e.g., Strecker et al. 2006). In fact, the large majority of Crustacea are in fact omnivorous, occasionally behaving as predators (Karatayev et al. 2009;Hänfling et al. 2011), and many benthic invertebrates classified as deposit-feeders are often also omnivores. In our results, species showing this particular feeding habit show great disparity in range shift estimates (i.e., 20 km to 3200 km; Supporting Information, Appendix S1, Table S1.7). Among them, Neomysis rayii had the largest projected shift, ca. 3200 km 2050 (RCP 2.6).
Besides reproduction and feeding habits, several additional traits should be considered for alien species to become invasive, e.g. thermoregulation, dispersal capability, and life history traits, which in concert with biotic interactions and evolutionary factors, can contribute to the success or failure of the establishment of a species in new environments (Hänfling et al. 2011). In particular, the complexity and broad-scale effects of climate change make it difficult to determine changes or distributional shifts a priori. Despite representing a small proportion of Crustacea diversity, the initial assessment presented here allowed the identification of general distributional patterns in the face of climatic changes and indicated features associated to invasiveness that could aid on the critical interpretation of ecological models. With the increase of anthropogenic disturbances such as fishing, mining, oil drilling, warming, and acidification in marine ecosystems, understanding the magnitude of changes in living marine ecosystems is essential for better management and conservation of their biological resources.

Conclusions
Our study provides a comprehensive assessment of patterns of climate change impacts on benthic marine crustaceans at the global scale, and results can be used as initial hypotheses in the development of future theoretical and empirical studies. Our results revealed idiosyncratic and species-specific responses, with prevailing poleward shifts and decline of species richness at midlatitudes, while more frequent shifts between temperate and polar regions were recovered. Shallow-water species are expected to shift longer distances than deep-sea species. The net gain of suitability was slightly higher than the net loss for shallow-water species, while for deep-sea species, the net loss was higher than the gain under all scenarios, which will potentially change the composition of assemblages, leading to reorganization and potentially also precipitate extinction of some species. Our estimations can be viewed as a set of hypotheses for future analytical and empirical studies, and should play an instrumental role in the planning and execution of strategic interventions and in the development of conservation strategies.