Among demons and killers: current and future potential distribution of two hyper successful invasive gammarids

Biological invasions represent one of the main contemporary pressures facing freshwater ecosystems, and a better understanding of invasive species potential distributions is essential to prepare for future stressors. Crustacean invaders contribute significantly to global invasions with the Ponto-Caspian region being one of the primary donor areas for the Palearctic. The amphipods Dikerogammarus villosus and Dikerogammarus haemobaphes, popularly known as “killer” and “demon” shrimps, are emblematic of successful Ponto-Caspian invaders of European freshwaters. However, the geographical areas in which the abiotic environment is potentially suitable for them have not been investigated. To address this gap, current and future potential distributions were studied for the European Western Palearctic considering two scenarios and time periods (2050 and 2070) as well as the association between anthropogenic activities and individual species habitat suitability. Results show large areas of central-western Europe are currently suitable for both species and indicate some potential for range expansion within colder European areas. In particular, D. haemobaphes has the potential to expand its range further west and within southern parts of Europe. Scenarios of future climate change don’t provide evidence for further range expansion compared to the current conditions and suggest a reduction of range overlap within the most suitable areas. Results reveal lowland areas are at greatest risk of colonisation as well as a significant association with anthropogenic activities for both amphipods. The outcomes of the research could be used by resource managers for preparing and managing future changes of both species distributions and facilitate decision-making for monitoring and control.

emblematic of successful Ponto-Caspian invaders of European freshwaters. However, the geographical areas in which the abiotic environment is potentially suitable for them have not been investigated. To address this gap, current and future potential distributions were studied for the European Western Palearctic considering two scenarios and time periods (2050 and 2070) as well as the association between anthropogenic activities and individual species habitat suitability. Results show large areas of central-western Europe are currently suitable for both species and indicate some potential for range expansion within colder European areas. In particular, D. haemobaphes has the potential to expand its range further west and within southern parts of Europe. Scenarios of future climate change don't provide evidence for further range expansion compared to the current conditions and suggest a reduction of range overlap within the most suitable areas. Results reveal lowland areas are at greatest risk of colonisation as well as a significant association with anthropogenic activities for both amphipods. The outcomes of the research could be used by resource managers for preparing and managing future changes of both species distributions and facilitate decision-making for monitoring and control.

Introduction
Biological invasions are among the most important drivers of contemporary biodiversity decline (IPBES 2019) with no sign of reaching saturation in accumulation rates at a global scale (Seebens et al. 2017). The spread of invasive species can result in a wide range of consequences on biological communities with major effects on ecosystem structure and functioning, as well as the services they provide (e.g., Blackburn et al. 2014;Wainright et al. 2021). In freshwater ecosystems, biological invasions represent a multi-faceted challenge for numerous stakeholders including ecologists, natural resource managers and industry (Gallardo and Aldridge 2020;Guareschi and Wood 2019) which have been less widely considered compared to other ecosystems (Moorhouse and Macdonald 2015;Stephens et al. 2019). Since eradication and management of freshwater invaders remains challenging (Simberloff 2021), preventing introductions is generally considered the most cost-effective way to mitigate the spread of new invasive species and conserve biodiversity (Leung et al. 2002).
A common approach to support conservation decision making and increase the predictive power of biological invasion research relies on estimating areas which may be suitable for colonization through Species Distribution Models (SDMs also known as habitat suitability models or ecological niche models, see Jeschke and Strayer 2008;Jiménez-Valverde et al. 2011;Guisan et al. 2017). By combining records of species occurrences and environmental data, SDMs identify potentially suitable areas for invasive species range expansion and provide opportunities for preventing or mitigating future invasions of both terrestrial and aquatic taxa (e.g., Barbet-Massin et al. 2018;Polidori and Sánchez-Fernández 2020). Furthermore, this approach permits the exploration of potential distributions under different climatic scenarios and trajectories (depending on greenhouse gas emission trends) that could support biological invasion research to: (i) help inform early detection of new alien species (e.g., identifying areas with highly suitability that have currently not been invaded- Guareschi et al. 2013); (ii) identify potential refugia for threatened species (e.g., Gallardo and Aldridge 2013) and (iii) highlight areas prone of multiple invasions (e.g., accumulation of multiple alien species with shared bioclimatic suitability-Gallardo and Aldridge 2020).
The wide availability of global climatological datasets, models and future climate scenarios provides the opportunity for detailed geospatial analysis and has the potential to improve the performance of species distribution modelling (e.g., Title and Bemmels 2018;Mammola et al. 2021).
The occurrence and spread of numerous invasive species, within different ecosystems globally, is primarily limited by climate and modified by human activities (Gallardo et al. 2015). Aquatic ectothermic species, including many invertebrates, may be particularly sensitive to modified patterns of precipitation, temperature and anthropogenic land-use alteration that may influence habitat availability (e.g., climate-change losers or winners, see Domisch et al. 2013;Gallardo and Aldridge 2013).
Freshwaters are largely invaded by a non-random selection of taxa, with Crustacea providing numerous examples of successful invaders worldwide (Strayer 2010;Cuthbert et al. 2020;Mathers et al. 2020). In this research, we focussed on two emblematic invasive amphipods from the family Gammaridae: D. villosus (Sowinsky 1894) and D. haemobaphes (Eichwald 1841), commonly known as "killer shrimp" and "demon shrimp" respectively. Both species originate from the Ponto-Caspian region and are considered successful invasive species in numerous European countries (Bacela-Spychalska and Van Der Velde 2013; Rewicz et al. 2014). D. villosus (hereinafter Dv) has been recognized as one of the 100 worst alien species in Europe (Nentwig et al. 2018) with negative effects on macroinvertebrate assemblages ) and bio-behavioural traits (e.g., high fecundity, feeding habits, quick sexual maturation, capacity to drift and be easily transported) typical of an "ideal" invader (Rewicz et al. 2014;Taylor and Dunn 2017). D. haemobaphes (hereinafter Dh) has been less studied compared to Dv but has already been shown to have effects on widely used macroinvertebrate biomonitoring metrics (Guareschi et al. 2021a), leaf litter processing (Constable and Birkby 2016), and the introduction of pathogens into its invaded range (Bojko et al. 2018). Similar to Dv, it displays numerous features that promote its successful invasion including its flexible feeding habit (omnivorous and predacious) and high fecundity compared to native gammarids (Bacela-Spychalska and Van Der Velde 2013).
In this study we aim to: (i) identify and compare the environmental suitability of both species under current climatological conditions in Western Palearctic Europe; (ii) project the potential distribution of both species under future climate change scenarios for two time periods (2050 and 2070); (iii) explore the association between anthropogenic pressures (human footprint) and species suitability values obtained for current climate conditions for both species. Together the aims sought to help guide monitoring and the prediction of future distribution of both taxa, as well as their potential overlap.

Biological data
Information on the current global spatial distribution of both invasive gammarids were obtained from published studies (e.g., scientific articles, reports) and from the occurrences uploaded on GBIF (Global Biodiversity Information Facility, 04 October 2021, Dv data https:// doi. org/ 10. 15468/ dl. h5ee3p, Dh data https:// doi. org/ 10. 15468/ dl. wfqj28). A complete list of the references and sources are provided in SM1. We incorporated data from both the native range and invaded regions as recommended when modelling invasive species distributions (Jiménez-Valverde et al. 2011;Sánchez-Fernández et al. 2011). Records covered the periods 1928-2021 and 1956-2021 for Dv and Dh respectively. However, the majority of records were obtained between 1990-2021. After records with taxonomic uncertainties or doubtful localities were removed, the final datasets comprised 2,868 records for Dv and 3,166 for Dh and can be considered representative of the current known distribution of both species.
The native range was delineated based on recent literature for both species (Rewicz et al. 2015;Jażdżewska et al. 2020). Both species originate from the Ponto-Caspian biogeographical region. This area, comprising the Black, Caspian, Azov, and Aral seas, is recognized as a region of extensive radiation of crustacean species (Cristescu and Hebert 2005) and one of the primary sources of invasive gammarids to Eurasian waters (Cuthbert et al. 2020). The strong invasive characteristics of both species has resulted in research focussing on the species' invasive ranges (e.g., Guareschi et al. 2021a;Hellmann et al. 2017) with the unusual pattern of displaying much more records in their invaded regions compared to the presumed native range (e.g., Rewicz et al. 2015). To minimize this spatial bias during niche modelling, we used a spatial thinning approach to ensure all pointers were a minimum of 50 raster cells (~ 46.3 km) distance from one another. To do achieve this, we used the R function "reduceSpatialCorrelation" from the package "SDMworkshop" (https:// github. com/ BlasB enito/ SDMwo rkshop). We obtained 137 occurrences for Dv (16 from its native range) and 103 (7 from its native range) for Dh prior to starting the modelling process (points displayed in Fig. 1).

Bioclimatic and environmental variables
In the current study "environmental potential distribution" and "bioclimatic/environmental suitability" have been considered as synonymous and defined as the geographical representation of the abiotic environment (climate and altitude) which is suitable for the species (e.g., Peterson et al. 2011;Polidori and Sánchez-Fernández 2020). In the absence of specific knowledge about the bioclimatic factors most likely to determine the distributions of Dv and Dh we started the modelling processes with a large pool of biologically meaningful variables. Data on current climatic variables (maximum and minimum temperature °C; precipitation mm) and altitude (m. a.s.l.) were obtained from WorldClim, version 2.1 (http:// www. world clim. org, Fick and Hijmans 2017) with a 30 s resolution (approx. 1 km 2 ) for the period 1970-2000.
From these climatological variables, we created the 19 bioclimatic variables identified in the WorldClim dataset using the R functions "biovars" (R package "dismo", Hijmans et al. 2020). This pool of variables was complemented with an additional set of 18 biologically relevant climatic variables, considered to be important for ecological and physiological processes when determining species distributions using the function "layerCreation" (R package "envirem", Title and Bemmels 2018) at the same spatial resolution. To avoid collinearity among variables we performed a hierarchical cluster analysis, resulting in a dendrogram indicating the similarity among all variables (Dormann et al. 2013). The chosen distance-threshold to select variables in the cluster was set at 0.3, (i.e. less than 70% correlation) following Polidori et al. (2021). We then inspected the Variance Inflation Factor (VIF) and excluded variables with the highest multicollinearity (VIF > 3) on an individual basis using a stepwise procedure (Heiberger and Holland 2015). A complete list of variables, definitions and the selection process steps are outlined in SM2. Estimation of species potential distributions: modelling procedure To provide a complete overview of both species we estimated the potential distribution through two complementary algorithms commonly used in SDM and biological invasion studies: Maximum Entropy Algorithm (Maxent: Phillips et al. 2006) and Generalized Additive Models (GAM: Hastie and Tibshirani 1990).
To run the models, we used the R functions "maxent" from the package "dismo" (Hijmans et al. 2020) and "gam" from "mgcv" (Wood 2011). To reduce uncertainties associated with single SDM algorithms, an ensemble modelling approach (Araújo and New 2007) was used to obtain a final averaged spatial probability of habitat suitability. This is a widely used approach in contemporary species-distribution modelling (e.g., Bosso et al. 2022). We obtained the ensemble models by averaging values with the "calc" function (fun = mean) from the R package "raster" (Hijmans 2021). Suitability was defined as a measure of the match between the environmental conditions of locations with species records and those currently without, ranging from 0 (low probability) to 1 (high probability).
All models were calibrated with information based on the current presence of the species. Absence information is rarely available (especially for invertebrates), and unreliable in the case of invasive species that are probably still expanding their distributions. For this reason, we randomly selected 10,000 background points, a common practice in SDM, inside the area represented by the Minimum Convex Polygon which was calculated using species records plus a buffer zone of approx. 100 km surrounding the edges of the polygon (e.g., Velazco et al. 2020). We created a buffer zone and excluded sea surface, whenever it was included inside the fixed area, using the software QGIS (version 3.16.8). We used the K-fold cross validation method to evaluate the predictive performance of the models. This technique consists of randomly splitting the dataset into K groups (folds) of approximately equal size and uses a sequentially one-fold (testing dataset) to validate the remaining K − 1 groups. In our study we divided our occurrence dataset into 5 folds (K = 5). In all validation runs, model performance was assessed using the area under the curve (AUC) and Boyce index. The AUC method provides a single measure of model performance and ranges from 0 (random) to 1 (perfect discrimination), while the Boyce index varies from − 1 to 1, where positive values indicate consistent model prediction compared to the test dataset and negative values indicate a model of poor quality (Hirzel et al. 2006).
Variable importance for the final ensemble models was assessed with a random forest approach following a bootstrapping procedure with 1,000 resamples and 10,000 random points. We used the current prediction of the ensemble model as the response variable for each gammarid species, and the selected variables used to derive the SDMs as predictors. The Random Forest algorithm was performed using the function "randomForest" (R package "randomForest", Liaw and Wiener 2002) and the variable importance obtained using the function "importance" with parameter type = 1 (mean decrease in accuracy). Analyses were carried out using R v 4.1.1 software (R Core Team 2021).

Future climatic scenarios
Effects of future climatic conditions on the potential distributions of both species were derived using the IPCC5 climate projections (IPCC 2014) from the CCSM4 global climate model for 2050 and 2070. Two different Representative Concentration Pathways, which reflect two possible changes in future anthropogenic greenhouse gas emissions, were explored (intermediate-low RCP 4.5 and intermediate-high RCP 6.0 scenarios according to IPCC5). The first assumes that global annual greenhouse gas emissions (measured in CO 2 -equivalents) will peak before 2050, while the second considers a delayed peak towards the end of the century (IPCC 2014). We deliberately avoided the unlikely worst scenario (coded RCP 8.5) following the recommendations of Hausfather and Peters (2020).
We obtained predictor variables for future climatic scenarios from WorldClim v. 2.1 (Fick and Hijmans 2017) and we created the same type of climatic and environmental layers as described above for the current scenario. Future potential distributions were estimated for all SDMs methods and the final ensembled models plotted.
Finally, once we obtained the projection for the two models and the ensemble for all different scenarios (current and future), we converted the continuous suitability maps into presence/absence (binary) layers using the threshold of suitability ≥ 0.7 (as highly suitable see Fournier et al. 2017) in order to assess potential spatial overlap between the two species.

Anthropogenic pressures and species suitability values
The Global Human Footprint (WCS-CIESIN, 2005) was used as a proxy of anthropogenic and propagule pressures, and its association with the suitability values of both species examined via Spearman rank correlations considering the entire dataset (> 20 M of cells). The Global Human Footprint (i.e., human influence index normalized by biome and realm, hereinafter coded HF) represents a global dataset of 1-km grid cells, created from nine global data layers covering human population pressure, land use and infrastructure, and has recently been identified as being associated with invasive species distributions (Gallardo et al. 2015). For both species, we graphically inspected the trend of the relationship between HF and suitability values using a loess function on 1,000 random subsets (n = 10,000) of the entire dataset. For each random subset we fitted the loess function and then predicted the expected results on 500 points equally spaced within the HF range (0-100) by plotting the trend line. Both Pearson and Spearman correlations, between HF and suitability values, were used to check for consistency of the correlations derived for the entire dataset. We opted for this approach because loess computation becomes infeasible for large datasets. Different settings for the random subsampling procedure were tested and provided similar results.

Results
A final set of 7 variables (out of the initial 38, see SM2) were used for modelling: BIO02 (mean temperature diurnal range), BIO08 (mean temperature of wettest quarter), BIO15 (precipitation seasonality), Elevation (altitude above sea level), Continentality (average temperature of warmest month-average temperature of coldest month), Thornthwaite aridity Index (metric of the degree of water deficit below water need) and the embergerQ metric (Emberger's pluviothermic quotient).
SDMs performed well for both species, with ensemble models displaying average values of AUC and Boyce index of 0.81 and 0.75 for Dv and 0.84 and 0.80 for Dh. Almost 4% and 5% of the European Western Palearctic may currently provide highly suitable environmental conditions (values of suitability ≥ 0.7, red colour in the figures) for Dv and Dh, respectively. These values increased to 11.5% and 10% if a suitable value threshold of ≥ 0.5 was used.
Both species displayed preferences for lowland areas, with records of average altitude below 100 m a.s.l. (from − 7 to 629, mean altitude = 34 m a.s.l. for Dv with n = 2,868 occurrences; and from − 4 to 618, mean altitude = 63 m a.s.l. for Dh, n = 3,166 occurrences). Based on the results of the random forest analysis (details and plots in SM3) "elevation" was the most important variable for both species followed by BIO15 (precipitation seasonality) and BIO08 (mean temperature of wettest quarter) for Dv and BIO08 and Continentality for Dh.
Dv current and future potential distributions Ensemble model results indicated that central-western Europe (e.g., Germany, Belgium, The Netherlands) and the British Isles currently provide highly suitable environmental conditions for Dv (Fig. 1a, 2) as well as the native range areas around the Black Sea. Highly suitable areas were also predicted for France and North-East Italy (Fig. 1a, 2) where the species has only been recorded a few times thus far and for Sweden, Denmark and Iceland where the species has not yet been recorded. Mediterranean areas were not identified as being favourable for the species, except for a few restricted areas in northern Spain (e.g., the Catalonian region) and the Adriatic coast of Italy. Considering future climate scenarios, the amount of highly suitable areas (values ≥ 0.7) for Dv would decline (− 8.2% with RCP 4.5 and − 9.8% with RCP 6.0, at 2050) with some reduction in suitable areas especially in Italy, France and the Balkan region. In addition, there was no shift north in suitable areas (Fig. 3). Reductions were most marked for the 2070 projections (− 20.7% with RCP 4.5 and − 12.7% RCP 6.0, maps presented in SM4).

Dh current and future potential distributions
The current potential distribution of Dh includes highly suitable areas in central-western Europe and the British Isles. In addition, high levels of suitability were identified along the Atlantic coast of France, and Iceland (Fig. 1b, 2). In contrast to Dv, highly suitable areas in the western and southern parts of the European continent (e.g., coastal areas of Italy, Spain and Portugal) were identified, although similar to Dv high suitability was predicted to be maintained in NE Italy (Fig. 2).
The overall number of highly suitable cells for Dh were forecast to experience minor changes for the 2050 scenario (+ 1.3% with RCP 4.5 and − 1.8% with RCP 6.0) maintaining suitable habitat in Mediterranean areas and local contractions of suitable habitat primarily in Italy and France (a reduction in the most suitable habitat-red areas in Fig. 4). These patterns were most marked for the 2070 scenario projections (− 13, 2% with RCP 4.5 and − 3.2% RCP 6.0, maps available in SM4).

Species distributions overlaps and the role of anthropogenic pressures
The geographical overlap of highly suitability habitat for both species (i.e., values ≥ 0.7) under current climate conditions mainly occurs in the northern part of central-western Europe, British Isles and North-East Italy (45.8% of their cells, Fig. 2). The percentage of overlapping suitable areas decreased for both scenarios in 2050 (RCP 4.5: 40.2%-RCP 6.0: 40.9%) and 2070 (RCP 4.5: 33.2%-RCP 6.0: 36.5%, detailed maps presented in SM5). Significant correlations between the proxy of anthropogenic and propagule pressures (HF) and the current suitability values were obtained for Dv (Spearman rho = 0.33, p value ≤ 0.0001) and Dh (Spearman rho = 0.28, p value ≤ 0.0001) for the entire dataset. Further analysis using multiple randomized sub-sets of data confirmed these patterns using both Pearson and Spearman rank correlation approaches (results and plots presented in SM6).

Potential distributions in European Western Palearctic
The amphipods D. villosus and D. haemobaphes are emblematic examples of successful invasive species of European freshwaters, although little was known about their potential distributions at larger spatial scale. This research is the first attempt to project their distributions across the European Western Palearctic under different climate change scenarios and time-periods. Ensemble species distribution models confirmed that for both species, highly suitable environmental conditions exist in numerous geographic areas of Europe that have already been colonized. The results also indicated the potential to expand their ranges into new areas and identified some differences between the patterns for the two co-generic species. The species largely overlap in their current potential distributions in central-western Europe. However, there appear to be differences in some regions with the availability of suitable habitat for D. villosus being limited in Mediterranean areas while western and southern European regions still provide suitable habitat scores for Dh. This similarity may partially reflect the fact that "elevation" was by far the most important variable identified in both species' distribution models. Nevertheless, the two species differed in terms of the second most important variable, with "precipitation seasonality" being important for the distribution of Dv. Interestingly this variable has been listed among the ecologically meaningful variables in the study by Gallardo et al. (2012), while little is specifically known about Dh. Overall, rainfall patterns strongly influence the availability of aquatic habitats as well as the frequency and duration of droughts and floods (especially in rain-fed systems). As indicated by the predictions, both species seem to prefer areas with perennially flowing systems and avoid mountainous regions.
Examination of the current potential distributions of both Dv and Dh indicated that areas with limited records of invasion to date (e.g., France) and areas with no records thus far (e.g., Sweden, Denmark, Iceland, Ireland) support high environmental suitability. This highlights the potential for further geographical range expansion. The suitability of these areas, as well as some Mediterranean regions for Dh, indicates the need for local monitoring programs and taxonomic training for stakeholders to avoid the potential colonisation of these invasive species being overlooked.
It is important to recognise that outputs from SDMs are not without limitations and depend on data accuracy, the type of modelling technique, and niche stability of the species under investigation (e.g., Elith and Leathwick 2009). However, following the approach and recommendations of Sánchez-Fernández et al. (2011) andJiménez-Valverde et al. (2011), these constraints have been minimized. Modelling results would benefit from further data from the region where both species occur naturally, since they are currently poorly studied in their native range (Ponto-Caspian area, e.g., Crimea region) compared to their invaded ranges. Unfortunately, this regional information (e.g., based on field records) seems difficult to obtain in the short term due to local geo-political tensions (Sousa et al. 2022).
At the same time, the development of future SDMs would greatly benefit from further studies able to incorporate physiological information regarding thermal preferences and salinity tolerances, as well as drought resistance, of both amphipods based on data from their native and invaded ranges. After a series of short laboratory experiments (72 h) Dobrzycka-Krahel and Graca (2014) defined Dh as an eurytopic species, indicating that both high (e.g., 20 PSU) and low salinities (e.g., 0.1 PSU) reduced individual survival. Despite their euryhaline character, although Boets et al. (2010) suggested Dv displayed a preference for low conductivity waters, salinity may limit their colonization of mesohaline and polyhaline systems (e.g., saltwater intrusion into deltas of rivers) within bioclimatic suitable areas. Furthermore, the persistence and further spread of both species in semi-arid areas seems unlikely given their documented thermal tolerances (e.g., Dv critical water temperature varying from 26.2 ± 1.8 to 31 ± 1.4 °C- Wijnhoven et al. 2003).
Our results largely agreed with existing bioclimatic model projections for Dv in Great Britain (a minimum of 60% being bioclimatic suitability-Gallardo et al. (2012)) and recognises the potential for further range expansion of Dv into suitable areas in Ireland. Remarkably high levels of suitability were observed for specific areas of North-East Italy corresponding to the Po River delta, and the surroundings of the Venice lagoon for both Dv and Dh. Inland waters in this region are primarily located in a matrix of heavily exploited lowland landscapes ) and considering the recent trajectory of Dv range expansion in Italy (first recorded in 2003 with range expansion toward central Italy, Casellato et al. 2006;Catasti et al. 2017) specific monitoring programmes for early detection and range definition are recommended.
Future scenarios of climate change don't suggest an expansion of the environmental suitability for either species, indicating that they are not predicted to benefit from climate change. In fact, some potential contractions in range were predicted for both species for Italy, France and on the Balkan Peninsula. This was highlighted in both temporal projections (2050 and 2070) and emission scenarios (RCP 4.5, RCP 6.0), and could reflect the fact that "elevation" was held constant during the modelling. Different freshwater invaders may display different responses to climate change scenarios, for example Gallardo and Aldridge (2013) suggested a range contraction for the invasive crayfish Pacifastacus leniusculus (Dana, 1852) by 2050 in Europe, while predicted range expansion for the invasive bivalve Dreissena polymorpha (Pallas, 1771) within the same areas. Similarly, Zhang et al. (2020) predicted an expansion of climatically suitable areas in Europe for the red swamp crayfish, Procambarus clarkii (Girard, 1852), with concurrent contractions in other continents (e.g., North America and Asia). These results highlight the need to consider the context specific characteristics of each invasive species, and the difficulties of making wider generalizations regarding interactions between biological invasions and climate change, for many aquatic invertebrates.
To predict future biodiversity changes in the areas identified with high levels of suitability, for one or both invasive amphipods, is not straightforward. Moreover, the consequences may be difficult to forecast because of the plasticity of feeding behaviours that characterize the species (e.g., predator, omnivore, shredder, see Bovy et al. 2015;Worischka et al. 2018). Our findings indicate further geographical areas where both these invasive amphipods would encounter resident naïve taxa as well as co-occurrent non-native species. In this context, highly suitable regions for one or both species (e.g., NE of France and Italy, UK, Belgium) largely overlapped with those regions recently identified by Magliozzi et al. (2020) as at greatest risk of high cumulative impacts from freshwater invaders. The cumulative effects of multiple invaders may constrain aquatic biodiversity in rivers (Guareschi et al. 2021b) but also raise possibilities of new interactions among introduced species (e.g. competition) that require further investigation. Local scale features (e.g., habitat complexity and refugia potential) may also affect the establishment and successful range expansion of these species. For example, MacNeil (2019) found that the predatory activity of D. villosus on Crangonyx pseudogracilis Bousfield, 1958 and Asellus aquaticus (Linnaeus 1758) was influenced by water quality and habitat type. At the same time, in the geographical areas indicated as potentially highly suitable for both Dv and Dh (overlapping range), local factors and behaviour mechanisms may affect their co-existence but also the exclusion of one or the other, meaning that additive or synergic effects cannot be automatically assumed. The species share some environmental preferences: both have been described as predominately lithophilous (Borza et al. 2017a), with an affinity for hard substrates and cobbles (Clinton et al. 2018;Macneil and Platvoet 2013). However, Borza et al. (2017b), studying sites along the River Danube (East Europe), also reported differences in their current velocity preferences, with Dv being primarily associated with lentic/slow-flowing systems, but able to use refugia provided by obstacles in faster flowing systems. This difference may point to niche differentiation and may, at least partially, explain the co-existence patterns of both species in some waterbodies. Moreover, laboratory mesocosm experiments (with amphipods from Polish waters) highlighted the role of competitive interactions in influencing the dispersal of invasive gammarids (Kobak et al. 2016). D. villosus is considered a strong competitor within invertebrate communities (Jermacz et al. 2015) and its presence may unintentionally facilitate the dispersal of other invasive amphipods, like D. haemobaphes, which tend to avoid interactions with Dv and may therefore continue to migrate through hydrographic networks (Kobak et al. 2016).

Invasibility and anthropogenic pressures
Anthropogenic pressures on freshwaters, including dams, road infrastructure and urbanization gradients have been recognised as major factors associated with non-native fish richness in North American watersheds (Anas and Mandrak 2021). Similarly, field experiments in the Netherlands suggested that artificial structures may facilitate the spread of Dv (Macneil and Platvoet 2013). Our results also indicated a significant association between HF (a proxy of propagule pressures and human disturbance) and current environmental suitability values for both invasive amphipods, although they were only of moderate strength. Waterway utilisation (e.g., navigation) as well as local factors, such as local river properties and resident faunal community characteristics may also affect the successful spread and range expansion of both species. For instance, the removal of geographic barriers, (e.g., construction of canals which link previously unconnected waterbodies) and the increase in global container shipping trade, have potentially made it easier for invasive amphipods with flexible and opportunistic traits to colonise and invade new regions. At the continental scale, two inland migration corridors were identified as the main waterways utilised by Dh and Dv for range expansion in Western Europe: the southern (via the Danube River) and central corridor (via the Dnieper River), with further opportunities provided by the interconnection of river basins at the regional level via man-made canals (Bij de Vaate et al. 2002). This illustrates that the extensive network of artificial canals in the United Kingdom, coupled with inter-basin water transfers, may further promote the regional dispersal of D. haemobaphes and other invasive invertebrates (Gallardo and Aldridge 2018;Johns et al. 2018).
It has been suggested that low habitat specificity could be key in population establishment, and invasion success for some invertebrates (e.g., spiders: De Smedt & Van Keer 2022). Dv has been hypothesized to be highly adapted to anthropogenic modification of aquatic ecosystems (e.g., Danube River, Borza et al. 2017b), but with preferences for large rivers and canals with good chemical water quality (e.g., Belgium, Boets et al. 2014). Boets et al. (2014) highlighted the high dissolved oxygen preferences of D. villosus, and that improving chemical water quality may facilitate the colonization of additional watercourses. Likewise, Guareschi et al. (2022) recently illustrated that the expansion of D. haemobaphes in UK rivers was associated with the presence of other invaders in rivers with greater channel width and water depth compared to uninvaded control sites, but with no clear preferences for altered or simplified macroinvertebrate communities.

Final remarks and management considerations
Given the difficulty in controlling invasive amphipods (Wood et al. 2021), our maps can be best considered "Risk-Maps" for Dv and Dh in the European Western Palearctic region. These maps can be used to anticipate future trajectories of change for both species distributions and to provide resource managers with a powerful spatial and temporal basis to inform decision-making regarding the allocation of resources for monitoring and control of these amphipods.
The results of this research confirmed the suitability of large areas of Europe and indicated geographic range expansion into colder lowland European regions where the bioclimatic conditions appear highly suitable for both species (e.g., Ireland, Iceland, Denmark and Sweden). The results also suggest that the range of Dh may expand further, well beyond that it currently occupies, in parts of Italy, Spain and Portugal under current climate conditions, as well as expanding in areas already colonized (e.g., France). Three species from the genus Dikerogammarus, including Dv and Dh, have been already listed as potential alien fauna likely to expand towards Iberian freshwaters (Oliva-Paterna et al. 2020). Examination of environmental suitability conditions under future climate scenarios didn't indicate further expansions of their distribution compared to the current conditions. Overall, actions are recommended (e.g., survey campaigns) within aquatic ecosystems in areas with high bioclimatic suitability that have not been invaded thus far and taxonomic/biomonitoring training for different stakeholders (e.g., environment regulators and researchers). This is imperative given these species (as numerous "new" alien fauna) are either not recognised or barely considered in the most common macroinvertebrate taxonomic keys used in Europe. These results presented here should not be considered a single point in time: when new records are confirmed, they should be incorporated into the modelling procedure to update and improve estimates of Dv and Dh potential distributions. Periodically repeating the analysis using the best available information would help validate existing models and update the current outputs (e.g., Puchałka et al. 2021). Moreover, specific high-resolution risk maps and assessments (e.g., Bosso et al. 2017) are also recommended as components of regional environmental restoration projects, especially in areas with high bioclimatic suitability or in areas where the presence of these species is already suspected.
Data availability Data analysed during this study are included in this published article. Detailed information about coordinates and sources are available in SM1.

Conflict of interest
The authors declare that they have no potential conflict of interest in relation to the study in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.