Impacts of climate change on the global potential distribution of two notorious invasive crayfishes

Zhixin Zhang, César Capinha, and Nisikawa Usio contributed equally to this work. 1Graduate School of Marine Science and Technology, Tokyo University of Marine Science and Technology, Tokyo, Japan 2Centro de Estudos Geográficos, Instituto de Geografia e Ordenamento do Território IGOT, Universidade de Lisboa, Lisboa, Portugal 3Institute of Nature and Environmental Technology, Kanazawa University, Kanazawa, Japan 4Cat Drop Foundation, Drachten, The Netherlands 5Key Laboratory of Animal Ecology and Conservation Biology, Institute of Zoology, Chinese Academy of Sciences, Beijing, China 6Instituto de Oceanografía y Cambio Global (IOCAG), Universidad de Las Palmas de Gran Canaria (ULPGC), Telde, Las Palmas, Spain 7Jiangsu Provincial Culture and Tourism Department, Nanjing, China

Climate change can greatly alter the distribution of suitable habitats to species including invasive species (Faleiro, Nemésio, & Loyola, 2018;Guisan, Thuiller, & Zimmermann, 2017), which may ultimately define the regions that will be susceptible to their negative impacts, including habitat competition, predation, hybridisation, and pathogen transmission (Huxel, 1999;Muhlfeld et al., 2014;Rosewarne et al., 2016;Schrimpf, Schmidt, & Schulz, 2014;Zhang, Yokota, & Strüssmann, 2019). Therefore, it is crucial to understand how future climates could alter the distribution of suitable habitats for invasive alien species, in order to effectively protect native biodiversity and ecosystem functioning. Species distribution models (SDMs), which aim at predicting habitat suitability by linking species distribution data with environmental data, have been extensively used to investigate the impacts of climate change on the potential distribution of a number of species including invasive alien species (Elith & Leathwick, 2009;Guisan et al., 2017;Guisan & Thuiller, 2005;Peterson, 2003;Thuiller, Guéguen, Renaud, Karger, & Zimmermann, 2019).
For the purpose of effectively preventing new invasions and mitigating adverse effects of the two invasive crayfish species on native species, several SDMs have been developed for predicting their current potential distribution Capinha, Leung, & Anastácio, 2011;Chucholl, 2016;Larson & Olden, 2012;Zeng, Low, & Yeo, 2016 these studies only focused on predicting habitat suitability under current climate conditions and did not consider climate change scenarios. Thus far, the impact of climate change on the potential distribution of the two targeted species was estimated by a single modelling algorithm (Capinha, Larson, Tricarico, Olden, & Gherardi, 2013;Liu, Guo, Ke, Wang, & Li, 2011) or considered limited geographical extents (Capinha et al., 2013;Gallardo & Aldridge, 2013).
For example, Capinha et al. (2013) (Elith & Graham, 2009;Pearson et al., 2006;Qiao, Soberón, & Peterson, 2015). Considering the highly invasive nature of the two crayfishes and their negative impacts on native ecosystems, it is of great importance to project their potential distribution under climate change scenarios at a global scale while accounting for variability in predictions from distinct modelling algorithms.
To reduce uncertainties in single SDM algorithms, the ensemble modelling approach has been proposed as a more reliable alternative given that it explicitly accounts for the variability of distinct single predictions (Araújo & New, 2007;Guisan et al., 2017;Thuiller, Georges, & Engler, 2014). The ensemble modelling approach can integrate predictions from different single SDMs thus can reduce prediction uncertainties and provide more accurate predictions than those from single SDMs (Araújo & New, 2007). Thuiller et al. (2019) have addressed the high levels of uncertainties in SDM studies associated with choices of SDM algorithms, general circulation models (GCMs), and representative concentration pathways (RCPs). Therefore, it is highly recommended to consider multiple SDM algorithms, GCMs, and RCPs in SDM studies (Thuiller et al., 2019). To examine the impacts of climate change on the global potential distribution of the two crayfishes, we (1) developed ensemble SDMs for the two invaders to (2) predict their current habitat suitability and (3) forecast their future potential distribution under future climate change scenarios derived from six GCMs under two RCPs. We hypothesise that the two crayfishes with different environmental tolerances may respond differently to future climate change.

| Predictor variables
The average values of 19 environmental predictor variables from 1979 to 2013 were downloaded from CHELSA version 1.2 (http:// chelsa-clima te.org) at a resolution of 30 arc-seconds (Karger et al., 2017). Gallardo et al. (2015) addressed the importance of human footprint on species distribution. Therefore, the global human influence index (HII, from Socioeconomic Data and Applications Center, http://sedac.ciesin.colum bia.edu) at a resolution of 30 arc-seconds was also included into the models. The HII and 19 environmental predictor variables were resampled to 5 arc-minutes by bilinear interpolation. A multicollinearity analysis was performed for the 20 predictor variables, and only the ecologically most meaningful and non-correlated variables (Pearson's correlation coefficient |r| < 0.7) were used for further analysis (Dormann et al., 2013).

| Modelling procedure
We used 10 modelling algorithms (artificial neural network, classification tree analysis, flexible discriminant analysis, generalised additive model, generalised boosting model, generalised linear model, multiple adaptive regression splines, maximum entropy, random forest, and surface range envelop) with their default settings in the biomod2 package (version 3.3.7) (Thuiller et al., 2014) in R (R Development Core Team, 2014). A total of 10,000 global pseudoabsence records were randomly generated, and a five-fold crossvalidation approach with 10 repetitions was adopted to evaluate predictive performance of 10 algorithms. Model performance was assessed by two metrics: the true skill statistics (TSS) and the area under the receiver operating characteristic curve (AUC). Algorithms with TSS over 0.85 and AUC over 0.90 for both species were selected to develop committee averaging ensemble models (Costa, Muelbert, Vieira, & Castello, 2015;Thuiller et al., 2014). Committee averaging is simply the average of binary predictions from the distinct algorithms; this procedure has the advantage of indicating both the suitability of conditions and the level of agreement among predictions of distinct algorithms. Relative contribution of each predictor variable and response curves of important variables were determined (Guisan et al., 2017;Thuiller et al., 2014). All data were used to predict habitat suitability for the two species under current climate conditions. Potential distributions of the two crayfish under different GCMs climate change scenarios were predicted and future range size changes under different GCMs were calculated.

| Model accuracy and contributions of predictor variables
Based on multicollinearity analysis results, 10 predictor variables were selected: annual mean temperature, mean diurnal range, temperature seasonality, maximum temperature of warmest month, annual precipitation, precipitation of driest month, precipitation seasonality, precipitation of warmest quarter, precipitation of coldest quarter, and HII (Table 1, Figure S1). AUC and TSS results suggested that in both species, the random forest algorithm performed best and surface range envelop performed worst (Table 2). For P. leniusculus, all algorithms except the surface range envelop showed high predictive performance; while for P. clarkii, only four algorithms, namely generalised additive model, generalised boosting model, generalised linear model, and random forest had high predictive performance (Table 2). Therefore, these four models were used to develop ensemble SDMs for both crayfishes. The ensemble SDMs for both species showed high predictive accuracy ( Table 2).
Results of relative contribution of each predictor variable suggest that the two species have different environmental requirements.
Annual mean temperature, temperature seasonality, and maximum temperature of warmest month are the three most important predictor variables determining the distribution of P. leniusculus, while annual mean temperature, HII, and temperature seasonality contribute most to the potential distribution of P. clarkii (Table 1). Although different predictor variables contributed differently to the distribution of the two crayfish species, annual mean temperature is consistently most important. The response curves of the occurrence probability of the two species against annual mean temperature were presented in Figure 1, which clearly showed that the two crayfishes have different thermal preferences: the optimal annual mean temperature is about 5-12°C for P. leniusculus and 12-23°C for P. clarkii.

| Potential distribution under current and future climate conditions
The predicted potential global distribution of P. leniusculus under current climate conditions is presented in Figure 2. Occurrence records of P. leniusculus are mainly distributed in the west coast of North America, Europe, and the north part of Japan. The ensemble model also predicted these areas to be suitable for this species (Figure 2).
In addition, the ensemble model predictions showed that the potential suitable habitat for P. leniusculus in Europe is much larger than its current known range ( Figure 2) (Kouba et al., 2014). In addition to these known distribution areas, other areas including the east coast of North America, South America, Australia, and New Zealand were also predicted to be suitable for P. leniusculus (Figure 2 Note: For the two species, each algorithm was run 10 times and data are expressed as mean ± SE. The asterisks (*) represent algorithms selected to develop ensemble species distribution models for the two species. P. clarkii has neither been reported in Australia nor in New Zealand, our ensemble model suggests that these areas are also suitable for this species (Figure 3).
Results of the range size change of both P. leniusculus and P. clarkii varied among different GCMs, and habitat of P. leniusculus is consistently predicted to contract in future (Table 3). For instance, under RCP2.6 scenario in 2050s, the range size change of P. leniusculus varied from −45.3% (HadGEM2-ES) to −11.5% (BCC-CSM1-1), whereas in P. clarkii the range size change was from -17.5% (MRI-CGCM3) to 8.0% (IPSL-CM5A-LR) ( Table 3). Predictions based on mean values of the six GCMs suggested that P. leniusculus will lose suitable habitat in future while P. clarkii will expand its ranges (Table 3); besides, the two species respond differently to future climate change depending on the region (Figures 4 and 5,  (Table 4). In the future, Africa will still be unsuitable for P. leniusculus and the potential distribution of P. clarkii is predicted to increase with the exception of RCP8.5 scenario during the 2070s (  (Figure 4). Habitat of P. clarkii in Europe is predicted to expand further north and eastward under future climates ( Figure 5).

| D ISCUSS I ON
In this study, ensemble species distribution models were developed, respectively, for P. leniusculus and P. clarkii to predict their global habitat suitability under current and future climatic conditions.
Both ensemble SDMs exhibited excellent predictive capacities and models showed that the two invasive crayfish species have different environmental requirements and will potentially display different responses to climate change.

| Important predictor variables
Among the selected 10 predictor variables, annual mean temperature is expected to play an important role in determining the distributions of the two crayfish species. Response curves of the relationships between occurrence probabilities of the two species and annual mean temperature indicate that P. clarkii has a preference for higher temperatures than P. leniusculus. This result is supported by the fact that P. clarkii is considered a warm water species whereas P. leniusculus a cool water species (Paglianti & Gherardi, 2004;Usio et al., 2006). Besides, previous studies have reported that the preferred temperature for P. clarkii is 21-30 °C (Espina & Herrera, 1993;Suko, 1958), while the optimal temperature for the growth of P. leniusculus is about 15-20°C (Ackefors, 1999). Other studies have found that temperature can affect a variety of physiological processes including growth, which ultimately influences survival of the two species (Harlıoğlu, 2009;Kozák et al., 2009;Paglianti & Gherardi, 2004 River Basin in the south-eastern USA using maximum entropy modelling approach. Contrary to our results, they reported that geology, slope, and streamflow contributed more to crayfish distributions than water temperature. In our study, we initially also considered two topographic variables (elevation and slope), but the two variables contributed little to crayfish distribution (Table   S1). Given this result, together with potential altitudinal shifts of species in response to climate change, we excluded elevation and slope from our analyses. Previous studies have demonstrated that

SDM projections can be influenced by a variety of factors includ-
ing the scale of the study area (Barve et al., 2011;VanDerWal, Shoo, Graham, & Williams, 2009). Given that Krause et al. (2019) focused on a single river basin in the south-eastern U.S.A., we suspect that climatic variation within this region might be more or less uniform, which may lead to the low contribution of water temperature to crayfish distributions.

| Potential distributions under current and future climate scenarios
The two ensemble SDMs suggest that the present global suitable habitat for P. leniusculus and P. clarkii is much larger than their current known distribution ranges. Our results highlight extents of agreements with predictions from previous SDM studies on the two species Capinha et al., 2013;Gallardo & Aldridge, 2013;Larson & Olden, 2012;Liu et al., 2011), but there are also noteworthy differences in several regions. with suitable habitat shifting mainly north-eastward. Based on the study by Gallardo & Aldridge (2013), range contraction is expected to be up to 32% by 2050. Habitat contraction and the north-eastward shift were also found for P. leniusculus in Europe in our study.
Our predictions suggest that future climate change will lead to about 23-30% loss of habitat in Europe. These differences might be due to differences in modelling techniques and selection of current and future climate conditions. For P. clarkii, Capinha et al., (2013) show patterns of range expansion in Europe that are very similar to the ones identified here, and which corresponds mainly to an expansion of suitable climates in north-eastern regions, while most of currently suitable areas will remain stable.
Our study forecasted the distributions of P. leniusculus and P. clarkii under future changing climate conditions at a global scale using an ensemble model approach, which has three major advantages over previous projections made for these species. First, we obtained current climate data from CHELSA, which represents average climate data for 1979-2013 by a quasi-mechanistic statistical downscaling method (Karger et al., 2017). al., 2005), it has been frequently used in SDM studies including recent SDMs for crayfishes (see Gallardo & Aldridge, 2013;Liu et al., 2011). However, the current climate conditions from WorldClim represent the average climate data for 1970-2000(Fick & Hijmans, 2017, which failed to properly match with temporal range of occurrence records (since 2000) of the two crayfishes. In addition, Bobrowski and Schickhoff (2017) suggest that species distribution models based on predictor variables from CHELSA showed higher predictive performance than models using WorldClim climate data.
Second, we adopted an ensemble forecasting approach to construct the SDMs for the two species. There are a number of modelling approaches available for SDM studies including classification, regression, and machine learning methods. Previous SDMs on P. leniusculus and P. clarkii were mainly developed by single modelling algorithms, such as the artificial neural network approach  and the maximum entropy modelling method (Liu et al., 2011). Previous studies have revealed that different single modelling algorithms have different predictive ability and can produce largely variable results (Elith & Graham, 2009;Pearson et al. 2006;Qiao et al., 2015). This inter-model variability in predictive accuracy was also found in our study. We should notice that the maximum entropy modelling approach might be the most widely used SDM algorithm, plausibly because of its easy-to-use features. In our study, this algorithm was not considered for P. clarkii due to its lower predictive performance. To reduce the model-based uncertainty and produce reliable predictions, the ensemble modelling technique has been proposed and frequently used, which combines prediction results of multiple modelling algorithms thus can reduce the uncertainties in model predictions (Araújo & New, 2007;Guisan et al., 2017;Thuiller et al., 2014Thuiller et al., , 2019.
Third, a range of future climate projections from six GCMs were utilised to forecast future distributions of the two crayfish species.
It has been revealed that high uncertainties in future projections of species distributions would arise from differences among GCMs (Goberville, Beaugrand, Hautekèete, Piquot, & Luczak, 2015;Tang et al., 2018;Thuiller et al., 2019). Therefore, it is recommended that used as future climates which is believed to reduce uncertainties associated with different GCMs (Yan et al., 2017).
Previous studies have demonstrated that future climate change could result in extensive loss of coastal wetlands due to sea-level rise (Hinkel et al., 2018;Schuerch et al., 2018;Spencer et al., 2016).
Therefore, we should notice that our SDM projections may overestimate future habitat availability of the two crayfish species in the coastal areas. In addition, niche modelling cannot consider thermal flexibility of some species. Thermal flexibility has also been reported in crayfishes (Haubrock et al., 2019;Veselý, Buřič, & Kouba, 2015); for instance, Veselý et al. (2015) reported from laboratory experi-   (Larson et al., 2017) and P. clarkii (Tréguier et al., 2014). In areas where the dispersal of invasive crayfish species seems unavoidable, mitigation management strategies to conserve vulnerable wildlife may also be planned prior to the arrival of the invaders. For example, both P leniusculus and P. clarkii have been reported to displace native crayfish through predation, competitive exclusion or transmission of diseases in Japan and Europe. In such areas, captive populations of native crayfish may be established before extirpation by crayfish invaders.

| Management implications
Regarding Australia and New Zealand, these areas are considered suitable for both invasive crayfish species, but they have not been reported yet. These two countries are taking biosecurity extremely seriously and have strict regulations regarding importation of species (see the Biosecurity Act 2015 and the Biosecurity 2025 Direction Statement for New Zealand's biosecurity system for details). Existing biosecurity measures should be kept in place, and the ban on aquarium and ornamental trades of these invasive species in Australia and New Zealand should be continued.
However, there are also concerns that harvested crayfish may transmit crayfish plague to native crayfish species (Souty-Grosset et al., 2016); in addition, similar to what has been verified for the invasive Chinese mitten crab in the River Thames, economic exploitation might result in intentional releases of individuals into new watersheds for the sake of financial benefits (Clark, 2011). As a result, the derived economic activity and societal pressure might run counter population control or the eradication of the species in southern Europe.
Considering the potentially high costs and risks associated to crayfish eradication programmes, more efforts may be directed towards the prevention of species introductions into new regions. Our results highlight areas that are highly susceptible to the establishment of two highly problematic crayfishes under current and future conditions. These results could help in the development of strategies aiming to prevent further invasions.

ACK N OWLED G M ENTS
We would like to thank Dr Xiong Zhang (The University of

CO N FLI C T S O F I NTE R E S T
Authors have no conflict of interest to declare.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.