The effects of climate change on European distributions of four alien marine crab species

Widely recognized as a major threat to marine biodiversity, invasive species have become a fundamental global concern. With over 500 alien species identified in the Mediterranean alone, European seas are particularly susceptible to the potential ecological and economic threats of invasives. The rate of marine species introductions in the European Union (EU) continues to increase, with climate change facilitating their spread and impact. Crabs and other crustaceans are among the most successful groups of marine invasives and can have significant negative ecological and economic impacts where they become established. To assess the ecological and economic threats posed by these species and to develop monitoring, early response, and mitigation plans, it is important to be able to determine which areas are at highest risk of further range expansion, especially under expected climate scenarios. We studied the current and future distributions of four predatory brachyuran crabs that were previously identified as species of concern for European seas, namely Hemigrapsus sanguineus, Charybdis longicollis, Matuta victor, and Portunus segnis, under various climate change scenarios. Species distribution models were built using an ensemble modelling approach. The results show that the potential distributions for all species are much larger than the current known distributions. Under all predicted climate change scenarios, the climatic conditions for P. segnis, C. longicollis, and M. victor, in particular, are expected to improve in most of the Mediterranean Sea, resulting in an expansion of suitable habitat. The Adriatic and Aegean seas are of particular concern as results indicate that these seas are not only highly suitable under current climatic conditions but also will become more suitable under all climate scenarios. It is, therefore, important to further investigate potential impacts, to increase monitoring, and to explore possible management strategies for these seas in order to manage the invasion of these species and avoid future biodiversity and economic losses.

1 Introduction

Marine ecosystems are increasingly under threat from invasive species – those that have expanded into areas beyond their natural geographical ranges (Vermeij, 1996). Widely recognized as a key component of global environmental change and as one of the leading causes of worldwide biodiversity loss (Butchart et al., 2010; Clavero et al., 2009; Gentili et al., 2021; Jaureguiberry et al., 2022; Sala et al., 2000), invasive species are often transported across natural boundaries by anthropogenic activities, including shipping, aquaculture, ornamental trade, and tourism (Ojaveer et al., 2018), and can have significant detrimental ecological and economic impacts (Diagne et al., 2021; Grosholz, 2002; Jaureguiberry et al., 2022; Katsanevakis et al., 2016; Kouba et al., 2022; Otero et al., 2013; Tsiamis et al., 2020). Furthermore, the dynamics and distributions of marine invasive species can be significantly influenced by climate change, adding extra complexity to the challenges already posed by these organisms (Stachowicz et al., 2002). Notably, the alteration of temperature regimes in marine environments, specifically warmer waters, can enhance reproduction and growth rates of some species, providing them with a competitive advantage over native counterparts (Cockrell and Sorte, 2013; King et al., 2021; McKnight et al., 2021; Sorte et al., 2010; Witte et al., 2010). Moreover, due to changing temperature gradients, previously unsuitable habitats can become more favourable for these invasive organisms (Burrows et al., 2011; Witte et al., 2010). Additionally, increasing temperatures lead to ocean acidification, which can also affect invasion dynamics. In some cases, for example, invasive species have demonstrated a greater resilience to changing pH levels, further exacerbating their competitive edge (Parker et al., 2011).

Within the European Union (EU), marine invasive species have become a fundamental concern for policymakers, with the identification and mitigation of potentially harmful organisms (those with the potential to pose significant threat to local biodiversity, ecosystems, and economies) becoming a high priority (Roy et al., 2019; Tsiamis et al., 2020). Specifically, marine invasive species can outcompete natives, disrupt food chains, alter habitats, and transport pathogens and parasites (Innocenti et al., 2009; Schrimpf et al., 2014), ultimately leading to population declines, local extinctions, changes to ecosystem functioning (Blackburn et al., 2014; Gallardo et al., 2016; Guy-Haim et al., 2018), and economic losses for fisheries (Katsanevakis et al., 2018) and tourism (Van Beukering et al., 2014). Consequently, Tsiamis et al. (2020) identified a number of species that should be prioritized for risk assessment. This culminated in a list of 26 top-priority species, including mostly those that are still not fully established within the EU. Among these top-priority species, four crabs (Brachyura) were identified, namely Hemigrapsus sanguineus (De Haan, 1835), Charybdis longicollis (Leene, 1938), Matuta victor (Fabricius, 1781), and Portunus segnis (Forskål, 1775).

Crustaceans are among the most commonly reported coastal marine invaders, with an estimated 73 brachyuran and crab-like anomuran decapods reported as being alien species (Brockerhoff and McLay, 2011), and more recently, 56 species of predatory brachyurans were identified as spreading outside of their natural ranges (Swart and Robinson, 2019). With high tolerances for abiotic stress, as well as genetic and phenotypic adaptations that increase invasion potential (Tepolt and Somero, 2014), predatory brachyurans have emerged as one of the most successful invasive taxonomic groups (Rato et al., 2021). These conspicuous invaders are considered particularly tolerant to global warming (Giomi and Pörtner, 2013) and ocean acidification (Wittmann and Pörtner, 2013). Furthermore, given reports of their ecological impacts, including predation on commercially important invertebrates (Brousseau et al., 2001; Tyrrell et al., 2006) and their potential to cause shifts in community structure, in other parts of the world (Le Roux et al., 1990), including North America (Ens et al., 2022; Kimbro et al., 2009) and southern Africa (Swart, 2017), this group should certainly be considered high priority when assessing their potential impacts in European waters.

One of the first Lessepsian immigrant crustaceans to be identified in the Mediterranean, P. segnis, underwent a long period of establishment, followed by rapid expansion, which is thought to have begun in 2001 (Castriota et al., 2022). High demand for this species led to the development of target fisheries, and they have high value in some markets (Mili et al., 2020); thus, the presence of P. segnis represents a valuable fishery resource within parts of its invasive range (Hamida et al., 2019; Liquete et al., 2013). However, as is also expected for the later-established C. longicollis (Innocenti et al., 2017; Stasolla et al., 2015), predation on fishing catches and damage to fishing nets have a negative impact on fishery economies (Özgül and Akyol, 2019; Tsirintanis et al., 2022). Furthermore, P. segnis is known to provide strong competition with at least two species of edible crabs, Eriphia verrucosa and Carcinus aestuarii (Ariani and Serra, 1969; Tsirintanis et al., 2022), and is a prolific bivalve predator (Giraldes et al., 2016), potentially offsetting their economic potential by reducing availability of native fishery species. Both P. segnis and C. longicollis are littoral and sub-littoral species are are expected to have major negative impacts on local biodiversity where, as generalist predators, they may impact mollusc, fish, and crustacean populations (Hamida et al., 2019; Stasolla et al., 2015).

First recorded along the French coast in the early 1990s (Dauvin and Dufossé, 2011) and later in the Mediterranean in 2001 (Schubart, 2003), H. sanguineus was in the earlier literature considered to be more herbivorous (McDermott, 1992), though more recently was characterized as a broad omnivore (Ledesma and O'Connor, 2001). Feeding on at least three commercially important bivalve species (Brousseau et al., 2001) within its invasive North American range, this species is also thought to be responsible for declining barnacle, mussel, polychaete, and macroalgae populations, indicating its potential impact on recruitment in these taxa (Tyrrell et al., 2006). It is also expected that H. sanguineus has the potential to cause substantive negative effects on sympatric populations of molluscs and crustaceans; however, its broad ecological and economic impacts are unclear (Epifanio, 2013), particularly in European waters. Lastly, Matuta victor, reportedly feeding predominantly on mysids, bivalves, fish, and crustaceans within its European invasive range (Innocenti et al., 2017), has recently been highlighted as a species of concern due to its expanding population in Turkish waters (Uysal et al., 2024). Despite few reports of its ecological and economic impacts, it is thought to be a threat to local ecosystems, particularly sandy beaches and coastal habitats (Uysal et al., 2024).

To better understand the potential impact of these four species, we aimed to predict the current and potential future distribution of each, under various climate scenarios in European waters, using ensemble species distribution models (SDMs). Species distribution modelling has become a cornerstone in predicting the suitable environmental niche of invasive species, which can be indicative of their potential spread. These models use climatic data and occurrence records to model the environmental niche in which species occur and project this to areas based on known or predicted environmental conditions (De Kort et al., 2020). These models provide insight into the potential range (expansion) of a species based on a region's environmental suitability (Zhang et al., 2019). The integration of ensemble modelling techniques can increase the reliability of these predictions (Ramirez-Reyes et al., 2021; Srivastava et al., 2019). Ensemble modelling involves combining outputs from diverse statistical and machine learning algorithms (such as artificial neural network, MaxEnt, and general linear models), thereby reducing the uncertainties inherent in single-model approaches (Araújo and New, 2007; Grenouillet et al., 2011). However, some studies show that ensemble modelling does not always improve predictions over single models but still provides a robust approach (Araújo and New, 2007; Parker, 2013). This synergistic approach leverages the strengths of various models while accounting for their weaknesses, thereby offering a more comprehensive depiction of a region's environmental suitability (Harris et al., 2024; Seni and Elder, 2010). This is particularly relevant when addressing invasive species' responses to changing climate patterns, making ensemble modelling an important tool for robust decision-making. These ensemble models have previously been applied to other invasive crab species in Europe, such as the Chinese mitten crab, Eriocheir sinensis (Zhang et al., 2019), but also in other regions such as the Bohai Sea in China, where the potential range expansion was modelled for the invasive two-spot swimming crab, Charybdis bimaculata (Zhang et al., 2024). By adopting this approach here, we aimed to provide robust predictions that facilitate scientists and policymakers in anticipating and responding to the potential impacts of each species, ultimately contributing to the preservation of native ecosystems and biodiversity.

2 Methods

2.1 Focal species

All four focal crab species (Brachyura) identified by Tsiamis et al. (2020) were used for our distribution modelling. The Asian shore crab, Hemigrapsus sanguineus, has already become established across the coasts of North America and Europe. It was first observed in Europe in 1999 in Le Havre, France, and has since spread along the Atlantic coast, all the way to the Baltic Sea (Karlsson et al., 2019). Charybdis longicollis is known for its agile swimming abilities, enabling it to traverse a variety of habitats. The species can reach very high densities within its non-native range, sometimes comprising up to 70 % of the benthic biomass in trawl catches, which is thanks to its adaptability and competitive nature (Deval, 2020; Galil, 1986). It was first observed in Europe near Türkiye in 1959 (Holthuis et al., 1961) and since then has expanded its range from Egypt to Cyprus (Galil, 2000). The moon crab, Matuta victor, originates from the Indo-Pacific region. Its preferred habitats include mangroves, intertidal zones, and sandy shores (Bom et al., 2020; Mohanty et al., 2019). Its adaptable nature and ability to tolerate a range of salinities have facilitated its spread beyond its original geographical boundaries. It was first observed in the Mediterranean Sea in 1902 in Port Said, Egypt (Global Biodiversity Information Facility (GBIF) Matuta victor occurrence). From there, it has spread to other areas in the Mediterranean Sea, including coastal areas of Israel, Lebanon, Türkiye, and Greece. The African blue swimming crab, Portunus segnis, predominantly inhabits shallow coastal habitats, including rocky shores, sandy bottoms, and seagrass beds (Spanier and Galil, 1991). This species was first recorded in the Mediterranean Sea in 1886, near Egypt, and plays a significant ecological role as both a predator and scavenger, influencing local food webs and ecosystem dynamics (Castriota et al., 2022). The three latter species are all thought to have arrived in European waters from the Red Sea through the development of the Suez Canal, often referred to as the Lessepsian migration.

2.2 Occurrence records

We developed distribution models for the four focal species using global occurrence records, including both the native and invasive ranges of each to ensure environmental requirements were fully characterized for all species (Broennimann and Guisan, 2008; Capinha et al., 2011; Zhang et al., 2020). However, our predictions of current and future environmental suitability focused on EU waters only, encompassing the region between 30 to 75° N and 30° W to 40° E. This area includes the Norwegian Sea, North Sea, Baltic Sea, Mediterranean Sea, Black Sea, and the eastern part of the Northern Atlantic Ocean. Species occurrence data were obtained from the online repository GBIF (GBIF, 2023a, b, c, d, GBIF Hemigrapsus sanguineus occurrence; GBIF Matuta victor occurrence; GBIF Portunus segnis occurrence; GBIF Charybdis longicollis occurrence) and with additional records from the published literature for Charybdis longicollis (Çiğşar et al., 2021; Doğdu et al., 2021; Firat et al., 2008; Galil and Kevrekidis, 2002; Innocenti and Galil, 2011; Lewinsohn and Holthuis, 1964; Özcan and Katağan, 2016; Yokes et al., 2007). Occurrence records that contained coordinates (in the coordinate system WGS 1984) that had three or fewer decimal places were excluded from the dataset to ensure that each record fell within the correct raster cell. Occurrence records were filtered using a 5×5 arcmin grid (corresponding to the grid size of the environmental data), keeping only a single occurrence record for each grid cell. Any points that were conspicuously distant from other points underwent further investigation and were removed from the dataset if there were indications of inaccuracies, in order to enhance data quality. Since true absence data were unavailable, we generated 10 000 randomly dispersed pseudo-absence points (Elith and Leathwick, 2009).

2.3 Environmental data

In addition to species occurrence data, species distribution models also require environmental data presented as spatial grids. Numerous environmental factors can impact species distribution, but using a limited number is often sufficient to predict the distribution of marine species (Belanger et al., 2012; Bosch et al., 2018; Goldsmit et al., 2017). Taking into account biological relevance and data availability for both the present and future scenarios, we selected eight variables. Specifically, we included six benthic variables at mean depth from Bio-ORACLE v2.2 which encompassed mean, minimum, and maximum temperature; the range and mean salinity; and mean current velocity (Assis et al., 2018). Mean, minimum, and maximum temperatures were selected because elevated temperatures are known to enhance growth and reproduction, but excessive temperature increases may result in reduced survival. By contrast, reduced temperatures may also result in decreased growth and reproduction (Thirukanthan et al., 2023). Mean salinity and its range were selected based on tolerance of many invasive crab species towards wide salinity ranges (Rato et al., 2021). Current velocity was included because this may relate to swimming behaviour for species, such as Portunus segnis (Luckenbach and Orth, 1992). Furthermore, we incorporated bathymetric data and distance-to-shore from the Global Marine Environment Datasets (Basher et al., 2018). Bathymetric data were included because they reflect the depth tolerance of benthic crab species (Bertini and Fransozo, 2004) and distance to shore as a proxy for land-based processes such as nutrient input and anthropogenic disturbances. All of these environmental data were at a spatial resolution of 5×5 arcmin. To reduce the effect of multicollinearity in predictors, we calculated a variance inflation factor for these variables but only for the locations of the presence and absence points. All variance inflation factors scored below 3, indicating no significant multicollinearity issues. Consequently, we retained all eight variables in our analysis.

Bio-ORACLE also provides spatial data for future projections based on four representative concentration pathway (RCP) emission scenarios. These RCPs represent distinct scenarios predicated on varying assumptions concerning societal, economic, and physical factors that could influence climate change. RCP2.6 is an overly optimistic scenario that postulates swift and sufficient societal action in response to climate change, with emissions peaking between 2010 and 2020. RCP4.5 anticipates emissions declining by 2045, RCP6.0 assumes that global emissions peak by 2080, and RCP8.5 projects a continuous rise in emissions throughout the 21st century. Presently, RCP4.5 stands as the most probable scenario. The Bio-ORACLE future projections draw from three atmosphere–ocean general circulation models (AOGCMs: CCSM4, HadGEM2-ES, MIROC5) from the Coupled Model Intercomparison Project (Assis et al., 2018). Employing various models mitigates model uncertainties and enhances the accuracy of these predictions. We incorporated all four RCP scenarios to forecast the future distribution of the four crab species of which predictions based on RCP4.5 and RCP8.5 will be presented in this paper and RCP2.6 and RCP6.0 in the Supplement (Sects. S2 and S5). For these future projections, we assumed that bathymetry and distance to the shore would remain unchanged throughout the projected period.

Figure 1 True skill statistic (TSS) and area under the curve (AUC) for the four modelled crab species. Each bar represents the mean of five repetitions for a specific modelling algorithm. Error bars represent the standard error of the mean. The dashed line displays the cutoff values. Models that scored below these cutoff values were not included in the ensemble model.

Figure 2 Averaged response curves based on all repetitions of all algorithms for Hemigrapsus sanguineus. On the x axis, the environmental variables' depth (a), distance to shore (b), current velocity (c), salinity range (d), mean salinity (e), minimum temperature (f), mean temperature (g), and maximum temperature (h) are given. Along the y axis, the modelled response (partial probability of occurrence) to these environmental variables is given.

2.4 Presence – background modelling

Species occurrence data, from both the native and non-native range of the four species, were used to fit the ensemble models. The data were divided into two groups: 80 % of the data points were allocated for model calibration and the remaining 20 % for model validation. To construct the ensemble models, we used the biomod2 package (Thuiller et al., 2024) for R, utilizing a range of algorithms including maximum entropy (MaxEnt), generalized linear model (GLM), general additive model (GAM), general boosting model (GBM), classification tree analysis (CTA), artificial neural network (ANN), surface range envelope (SRE), flexible discriminant analysis (FDA), and multiple additive regression splines (MARS). We ran five repetitions for each algorithm. Models were evaluated using the true skill statistic (TSS) and area under the curve (AUC), with relatively conservative cutoff values of 0.85 and 0.90, respectively (Liu et al., 2016). The final ensemble models comprised only those models that met both of these cutoff criteria. To calculate a weighted mean ensemble model, we considered the TSS values, assigning higher importance to models with higher scores compared to those with lower TSS values. We calculated averaged response curves for all modelled variables. In these response curves, a predicted response was plotted against a single environmental variable which was allowed to vary along its environmental gradient where all other modelled variables were kept constant. These plots are indicative of the effect of environmental variables in the ensemble model. The relative importance of environmental variables was determined using Pearson's correlation coefficient on the model's predicted values versus environmental variables. This metric helps assess the relation between a variable and the predicted habitat suitability. Values range from 0 to 1, where 0 signifies no correlation of the variable on a species distribution, and 1 indicates a very strong influence. Both response curves and variable importance were calculated using the biomod2 library.

3 Results

3.1 Hemigrapsus sanguineus

For H. sanguineus, all models except SRE performed above the thresholds for TSS and AUC (Fig. 1a and b). The response curves show an increase in environmental suitability for H. sanguineus after a maximum temperature of 12 °C (Fig. 2h), with values below that value strongly limiting H. sanguineus. The species also occurs close to shore with environmental suitability declining rapidly at greater distances from shore (Fig. 2b). Furthermore, environmental suitability increases with mean temperatures over 7 °C (Fig. 2g). All other environmental variables relate strongly to the known species distribution (Fig. 2a, c, d, e, f). Variable importance was highest for maximum temperature, distance to shore, and mean temperature, respectively (Table 1).

Table 1 The mean variable importance (± standard error) for the four modelled crab species based on all models that were included in the final ensemble model; n is the number of occurrence records that were used for each model.

Figure 3 Modelled distribution maps and European occurrence records for H. sanguineus. These maps include the current potential distribution (a), the occurrence records in European marine waters (b), the predicted distribution under RCP4.5 for the year 2050 (c), the change in probability of occurrence from the current situation to 2050 under RCP4.5 (d), the predicted distribution under RCP8.5 for the year 2050 (e), and the change in probability of occurrence from the current situation to 2050 under RCP8.5 (f).

Figure 4 Percent of the total area per region in which the habitat suitability ranged between 60 % and 100 %. Regions that had no habitat suitability within this range are not included in the figure. Blue bars represent the area under current climatic conditions, orange under RCP4.5, and red under RCP8.5. A table with more details on habitat suitability is provided in the Supplement.

According to occurrence records, the species is mostly found along the North Sea coast, the English Channel, and the Kattegat (Fig. 3d). Incidentally, the species was recorded in the western basin of the Mediterranean Sea, the northern part of the Adriatic Sea, and the Bay of Biscay. Under current climatic conditions, the species' potential suitable habitat in Europe is mostly concentrated along the Belgian, Danish, Dutch, and English coastal areas with some additional areas in the northern part of the Adriatic Sea, the northwestern coast of the Black Sea, and the Sea of Azov (Figs. 3a and 4). Under RCP4.5, the suitable habitat of the species could see a potential increase in the White Sea, Celtic Sea, and Irish Atlantic coast (Fig. 3b and c). The northern part of the Adriatic Sea will become less suitable for the species. Under RCP8.5, all areas that are currently suitable (except for the Adriatic Sea) will become more suitable (Fig. 3e and f). In addition, the Baltic Sea may become more suitable as well.

Figure 5 Averaged response curves based on all repetitions of all algorithms for Charybdis longicollis. On the x axis, the environmental variables' depth (a), distance to shore (b), current velocity (c), salinity range (d), mean salinity (e), minimum temperature (f), mean temperature (g), and maximum temperature (h) are given. Along the y axis, the modelled response (partial probability of occurrence) to these environmental variables is given.

3.2 Charybdis longicollis

For C. longicollis the MaxEnt and SRE models did perform below the AUC threshold, whereas the threshold for TSS was met by the ANN, CTA, FDA GLM, and some of the MARS models (Fig. 1c and d). From a minimum temperature of 8 °C the environmental suitability increases rapidly (Fig. 5f), while habitat is less suitable with minimum temperatures below 4 °C. The species also seems tolerant to fluctuations in salinity. They are often found in areas with salinity ranges of 5 PSS (Practical Salinity Scale) or more, being able to tolerate increased salinity (Fig. 5d and e). In addition, the species prefers high temperatures with suitability continuously increasing with maximum temperatures from 7 °C up to 32 °C (Fig. 5h). Other variables did not result in a strong response (Fig. 5a, b, c, g, h). Variable importance was highest for minimum temperature, salinity range, and maximum temperature, respectively (Table 1).

Figure 6 Modelled distribution maps and European occurrence records for C. longicollis. These maps include the current potential distribution (a), the occurrence records in European marine waters (b), the predicted distribution under RCP4.5 for the year 2050 (c), the change in probability of occurrence from the current situation to 2050 under RCP4.5 (d), the predicted distribution under RCP8.5 for the year 2050 (e), and the change in probability of occurrence from the current situation to 2050 under RCP8.5 (f).

The species has been recorded in the Persian Gulf, coastal Madagascar, the Red Sea, and in the eastern basin of the Mediterranean Sea along the coast of Türkiye, Israel, and Lebanon (Fig. 6d). Its current potential suitable habitat covers all coastal areas within the eastern basin and to a lesser degree the western basin and Atlantic coast of Portugal and Spain (Figs. 6a and 4). Under RCP4.5, this potential suitable habitat increases throughout the Mediterranean Sea but decreases slightly along the Atlantic coast of Portugal and Spain (Fig. 6b and c). Under RCP8.5, this increase is further exacerbated, making all coastal areas in the Mediterranean Sea highly suitable (Fig. 6e and f).

Figure 7 Averaged response curves based on all repetitions of all algorithms for Matuta victor. On the x axis, the environmental variables' depth (a), distance to shore (b), current velocity (c), salinity range (d), mean salinity (e), minimum temperature (f), mean temperature (g), and maximum temperature (h) are given. Along the y axis, the modelled response (partial probability of occurrence) to these environmental variables is given.

Figure 8 Modelled distribution maps and European occurrence records for M. victor. These maps include the current potential distribution (a), the occurrence records in European marine waters (b), the predicted distribution under RCP4.5 for the year 2050 (c), the change in probability of occurrence from the current situation to 2050 under RCP4.5 (d), the predicted distribution under RCP8.5 for the year 2050 (e), and the change in probability of occurrence from the current situation to 2050 under RCP8.5 (f).

3.3 Matuta victor

All M. victor models, except SRE, performed well above the thresholds for both the AUC and TSS (Fig. 1e and f). The response curves show that with mean temperatures below 5 °C the species is mostly absent and increasing strongly with higher temperatures (Fig. 7g). Distance to shore shows a rapid decline in probability of occurrence with increasing distances (Fig. 7b). The species also tolerates high maximum temperatures of up to 30 °C (Fig. 7h). Other variables did not result in a strong response (Fig. 7a, c, d, e, f). Variable importance was highest for mean temperature, distance to shore, and maximum temperature, respectively (Table 1).

The occurrence records showed that the species has so far been recorded in the eastern basin of the Mediterranean Sea along the coast of Cyprus, Türkiye, Greece, Lebanon, and Israel (Fig. 8d). The current potential suitable habitat encompasses most of the Mediterranean coastal areas where the eastern basin is more suitable than the western basin (Figs. 8a and 4). Under RCP4.5, most Mediterranean coastal areas become more suitable than they already are, especially the central parts surrounding Italy, Malta, and Corsica (Fig. 8b and c). Under RCP8.5, this increase in suitability is even further exacerbated (Fig. 8e and f).

Figure 9 Averaged response curves based on all repetitions of all algorithms for Portunus segnis. On the x axis, the environmental variables' depth (a), distance to shore (b), current velocity (c), salinity range (d), mean salinity (e), minimum temperature (f), mean temperature (g), and maximum temperature (h) are given. Along the y axis, the modelled response (partial probability of occurrence) to these environmental variables is given.

3.4 Portunus segnis

The MaxEnt and SRE models performed all below the AUC and TSS thresholds for P. segnis (Fig. 1g and h). All the GAM models also performed below the TSS threshold, as well as one MARS model. The response curves show that P. segnis tolerates increased salinity levels ranging from 35 to 40 PSS (Fig. 9e). Minimum and maximum temperature (Fig. 9f, h) also display a strong response in the species with minimum temperature above 15 °C and maximum temperatures above 20 °C, resulting in a strong response. In addition, the species is found at much closer distances to shore in comparison to the other species (Fig. 9b). Other variables did not result in a strong response (Fig. 9a, c, d, g). Mean salinity, distance to shore, and mean temperature had the highest variable importance.

Existing occurrence records from P. segnis show that the species is currently found in the eastern basin of the Mediterranean Sea along the coasts of Sicily, Malta, southern Greece (Crete and Rhodes), Türkiye, Cyprus, Israel, Libya, and Tunisia (Fig. 10d). The current potential suitable habitat encompasses all coastal areas in the Mediterranean Sea (Figs. 10a and 4). Under RCP4.5, the Adriatic Sea becomes more suitable as well as the English Channel (Fig. 10b and c). However, the probability of occurrence in the English Channel remains relatively low. Under RCP8.5, the suitability in the central Mediterranean Sea increases further as well as in the Adriatic Sea (Fig. 10e and f), while in the English Channel the suitability decreases.

Figure 10 Modelled distribution maps and European occurrence records for P. segnis. These maps include the current potential distribution (a), the occurrence records in European marine waters (b), the predicted distribution under RCP4.5 for the year 2050 (c), the change in probability of occurrence from the current situation to 2050 under RCP4.5 (d), the predicted distribution under RCP8.5 for the year 2050 (e), and the change in probability of occurrence from the current situation to 2050 under RCP8.5 (f).

4 Discussion

We aimed to estimate the current and future environmentally suitable regions, under present and expected climate conditions, of four priority brachyuran species identified as being of significant ecological and economic concern for European marine systems by Tsiamis et al. (2020). Our models showed that the potential suitable habitat of all four species studied is much larger than current known distributions. Under all predicted climate change scenarios, the climatic conditions for Charybdis longicollis, Matuta victor, and Portunus segnis will improve in most of the Mediterranean Sea, resulting in an expansion of suitable habitat. Conditions for Hemigrapsus sanguineus only improve under RCP8.5, a worst-case scenario in which no climate mitigation will take place.

The Adriatic and Aegean seas were predicted to experience the most significant changes of the focal species under all climate scenarios. This finding supports recent suggestions that these seas are at particularly high risk of biological invasions (Glamuzina et al., 2024; Katsanevakis et al., 2020). Climate models forecast larger increases in both salinity and temperature in these seas compared to other regions within the Mediterranean. The Adriatic–Ionian Bimodal Oscillating System (BiOS) plays a crucial role in altering salinity in the Adriatic Sea. Specifically, during the cyclonic phase, increased salinity occurs as a result of input from the Ionian Sea, while the anticyclonic phase brings reduced salinity from lower salinity input from the northern Adriatic (Civitarese et al., 2023). Further, BiOS is sensitive to changes in climatic variables, potentially affecting salinity in the Adriatic, Ionian, and Aegean seas (Gačić et al., 2014; Theocharis et al., 1999). In addition to the prediction that temperature and salinity will be the main drivers of changing species composition in the North Adriatic (Rizzi et al., 2016), nutrient input has also been highlighted as important, with flow from the Po River delta playing a significant role (Bongiorni et al., 2018; Nasi et al., 2020). Historical data indicate that salinity fluctuates seasonally in the Northern Adriatic, being higher in spring and summer and decreasing in autumn and winter due to river runoff (Appiotti et al., 2014). Furthermore, the Po River's discharge has been decreasing, leading to increased salinity in the Adriatic Sea – a trend expected to continue, based on climate change predictions (Appiotti et al., 2014). Thus, it is reasonable to expect that the predicted increase of suitable habitat for the modelled crab species within this region is accurate. In addition, the BiOS also affects the dispersal of biodiversity; in particular, it introduces Lessepsian invasive species into the Adriatic Sea (Duman et al., 2024; Mihanović et al., 2021), whereas the water circulation patterns within the Adriatic Sea facilitate the spread of invasive species internally (Kraus et al., 2019).

Overall, temperature and salinity resulted in the strongest response in the distribution of all four invasive crab species, as both factors significantly affect reproduction and larval development as well as dispersal in many crab species. For example, C. longicollis prefers slightly increased salinity levels, which affects the male–female ratio, resulting in a higher proportion of females (Yesılyurt et al., 2019), with similar observations reported for P. segnis (Tureli and Yesilyurt, 2018). Additionally, increased temperature and salinity are known to accelerate maturity in blue crabs, Callinectes sapidus, with smaller crabs maturing faster under these conditions (Fisher, 1999). This pattern has also been observed in P. segnis (Tureli and Yesilyurt, 2018). However, Giraldes et al. (2016) found that smaller P. segnis produce fewer and potentially less viable eggs compared to their larger conspecifics. Under lower salinity levels, osmosis can lead to swelling, which may explain the size differences but not the difference in egg numbers (Huang et al., 2022). P. segnis typically spawns during warmer months, but in tropical regions, reproduction occurs year-round (Safaie et al., 2013). Thus, temperature and salinity not only affect the occurrence of these crab species but may also enhance their reproductive success, thereby potentially adding to their invasive success and dispersal to nearby regions.

Distance to shore was an important predictor for M. victor. This is likely related to the species preference for shallower waters, occurring nearer to shore (Mohanty et al., 2019), but could also be related to anthropogenic disturbances that are more common in the nearshore. Disturbances such as eutrophication, physical disturbances, and pollution seem to be advantageous for M. victor as it is more tolerant than its competitors to such disturbances, although these findings are based on circumstantial evidence (Innocenti et al., 2017; Zviely et al., 2021). It is therefore possible that more enclosed areas that are exposed to high agricultural runoff, such as those in the Adriatic Sea, could be more susceptible to this species.

Our results accounted for various model and data issues that can significantly affect the accuracy and reliability of predictions made using species distribution modelling (SDM). Specifically, Hui (2023) describes several points that should be carefully considered when building and interpreting SDMs. SDMs should be fitted using only relevant variables that, from a biological/ecological perspective, can be expected to affect the distribution of the focal species. Furthermore, occurrence records should be filtered in order to remove clustering and inaccuracies. We accounted for these issues by selecting the most relevant variables and filtering the occurrence records. However, filtering resulted in low sample sizes for both C. longicollis and P. segnis, which resulted in poor performance for some of the applied algorithms. Our ensemble approach used relatively high threshold values, which resulted in the exclusion of these specific algorithms from the final models.

Although ensemble modelling is likely to improve the robustness of our predictions, it is not always an improvement over single-algorithm models (Araújo and New, 2007; Parker, 2013). Several studies have also shown that some single-algorithm models perform just as well, require less computation, and are easier to interpret (Kaky et al., 2020). MaxEnt, for example, is one of the most used single-algorithm models and often performs really well and sometimes better than ensemble models (Kaky et al., 2020). MaxEnt was used as part of our ensemble approach, although it was excluded for P. segnis and C. longicollis due to poor performance. MaxEnt performed poorly for these two species, possibly as a result of the lower sample size in combination with the global range that we modelled (van Proosdij et al., 2016).

Species distribution models model the environmental niche in which a species occurs (Thuiller et al., 2024). However, it should be noted that other factors may also impose limitations on a species' distribution. Such factors include predation pressure, food availability, and competition (Cannicci et al., 2018; Peer et al., 2018; Poggiato et al., 2025; Trainor et al., 2014; Yackulic et al., 2014). Invasive species are often opportunistic and can easily adapt to new circumstances (Simberloff, 2013). Thus, selection would likely favour invasive populations that are better adapted to the environmental conditions in their new range than the species would be in their native range (Vera-Escalona et al., 2023). Therefore, its environmental niche does not necessarily represent a species' tolerance to certain environmental factors. For example, Paz and Guarnizo (2020) compared the thermo-physiological tolerance limits of 50 reptile and amphibian species with niche estimates from SDMs. The differences that were observed between, for example, minimum temperatures were large for many species, indicating that inference about niche traits and tolerance limits from these models should be avoided. Predicting a species' distribution based on the niche where circumstances could be entirely different requires caution when interpreting such results, especially when sample numbers are low, as seen for P. segnis and C. longicollis in this study.

Adaptation to new niches either through ecological release, evolution, or favourable genotypes can explain niche shifts and incorporating this into SDMs can significantly improve models and provide more reliable prediction (Gallien et al., 2012; Rödder et al., 2009). In SDMs, this can be done by modelling a species' native and invasive range separately. Besides increasing the accuracy of the predictions, this will also provide insights into a species' adaptability to new environmental conditions, as demonstrated by Crickenberger (2016) when modelling the distribution of an introduced barnacle species in the USA. Here, the observed range retraction of Megabalanus coccopoma was modelled in its introduced range, using models trained on datasets containing only observation from the native, introduced, and global range. The models trained on the introduced range were only able to correctly predict the retraction, which can be indicative of a niche shift. The model fitted on global data was much more conservative and overestimated the range. Other SDM studies report similar discrepancies between the native and introduced niche, advocating for splitting data between native and introduced ranges when sample sizes are large enough and otherwise using a global dataset (Gallien et al., 2012; Rödder et al., 2009). Sample sizes for all our species were not sufficient to split between invasive and native range observations, especially in their invasive ranges, as this study focused on species that have not yet fully established in their invasive ranges (Tsiamis et al., 2020). Therefore, we used occurrence records from both the native and invasive ranges to fit the models, possibly resulting in a more conservative model (Crickenberger, 2016). As observations of these species will become more available in the future, it is advisable to compare these species' adaptations to these new environments.

For all species, the surface range envelope (SRE) algorithm performed well below the TSS and AUC thresholds. This algorithm defines the distribution based on the minimum and maximum values of the environmental variables. When a location falls within the ranges for all environmental variables, the area is considered suitable for that species. In the default settings of biomod2, this range is reduced to the 2.5th and 97.5th percentiles in order to reduce overpredictions caused by outliers. It is a relatively simple algorithm, but it is susceptible to overpredictions (Thuiller et al., 2024). The poor performance of SRE in biomod2 is also seen in other studies (Esmaeili and Eslami Barzoki, 2023; Liu et al., 2022; Tanaka et al., 2020). For example, Liu et al. (2022) developed species distribution models using biomod2 for another crab species (Portunus trituberculatus), also reporting SRE as the poorest-performing model. The dismo package uses a similar algorithm named BIOCLIM which calculates a range of percentiles while assuming that areas closer to the median values are more suitable (Hijmans et al., 2023). Although this is an improvement in comparison to the SRE algorithm in biomod2, it still underperforms in comparison to other algorithms (Elith et al., 2006). Hence, we suggest that the use of the SRE algorithm should be carefully considered when building future species distribution models, especially when using low TSS and AUC threshold values.

For predicting future environmental suitability, we used modelled climate data based on a combination of three global climatic models: the Community Climate System Model 4, the Hadley Centre Global Environmental Model 2, and the Model for Interdisciplinary Research on Climate 5 (Assis et al., 2018). These modelled data are based on four representative concentration pathways (RCPs). More recently, these data have been updated based on the Shared Socioeconomic Pathways (SSPs), including less earth system models and more variables that may be relevant for species distribution modelling (Assis et al., 2024). These new data include possibly relevant variables (e.g. pH and dissolved oxygen) that we were unable to use for our study and could have resulted in more accurate models. In addition, substrate type is an important variable that could help model species at a higher accuracy and scale. De La Barra et al. (2020) showed that this variable is very important in predicting the distribution of other portunid crabs. Nevertheless, substrate data were not available and difficult to attain on a large spatial scale (De La Barra et al., 2020). This does not invalidate our models, but it only allows for interpreting the predicted environmental suitability on a larger regional scale.

In conclusion, this study highlights the growing concern of marine invasive species in the EU, particularly under changing climatic conditions and demonstrates that a species' suitable habitat can expand significantly. Our research focused on four introduced crab species in the European marine environment, showing that suitable habitats for most of these species will increase. Habitat suitability for M. victor and C. longicollis should be expected to further increase more than others, especially in the Adriatic and Aegean seas due to predicted increases in temperatures and salinity. The effects of rising temperatures and salinity may also positively affect the reproductive success of these crab species, potentially exacerbating their invasive potential. Thus, we also highlight the potential importance of further study and the exploration of effective monitoring and risk screening strategies (Stasolla et al., 2021), such as early detection in high-risk areas (e.g. the use of crab condos; Hewitt and McDonald, 2013) for these two species to better understand the potential ecological implications of predicted range expansions.