Species distribution models were created for North America, which we defined as Canada, the Caribbean, Central America, Mexico, and the United States. Due to limitations in available hydrology data (See Environmental Variables below), we excluded the U.S. state of Hawaii. The total surface area for the study region was 43,731,082 km². Faxonius rusticus and P. clarkii are both native to portions of North America but have expanded their ranges such that they are each invasive across large regions of the continent (USGS 2020).

North American distribution data for both species was collected from the U.S. Geological Survey Nonindigenous Aquatic Species Database (USGS 2020) and from the Global Biodiversity Information Facility (GBIF 2020). We collected and used data for both native and non-native ranges for two reasons. First, although habitat within species’ native ranges often contains high densities of occurrence records which can influence model predictions towards conditions in the native ranges, native range data has been shown to offer greater insights into which environmental variables contribute most strongly to the model’s predictions (Barbet-Massin et al. 2018). We addressed the issue of a high density of sampling points with rarefication (see below). Second, crayfish SDMs focusing on only non-native or native ranges have been found to under-predict occurrences in the other range type (Morehouse and Tobler 2013).

Duplicate records and data points that lacked latitude and longitude coordinates were removed. Data points with less than four decimal places for latitude and/or longitude were also removed because these records were not sufficiently accurate for an SDM resolution of 1 km². We imported these records into ArcGIS Pro™ version 2.5 (Esri Inc 2020) and spatially rarefied the data (Brown 2014) so that no records were closer than 10 km to any other record (Appendix 1: Fig. A1). Rarefication is commonly performed on presence-only data because it alleviates issues of overfitting and inflation of model performance variables – these issues result when spatially clustered records exist in areas possessing disproportionately large sampling efforts, which is a common issue for occurrence records retrieved from large databases (Veloz 2009; Boria et al. 2014; Hui 2023). The final data consisted of 904 F. rusticus and 953 P. clarkii records.

Yiwen et al. (2016) investigated which environmental variables are most highly correlated with distributions for a range of crayfish species. The study examined the distribution of eight crayfish species at the global extent and used step-wise regression to quantify the relationship between 34 environmental variables selected from climate, hydrological, and landcover datasets. The Yiwen et al. (2016) study similarly examined environmental variables using a priori selection but found more robust modeling occurred for our target species when using step-wise regression. A main aim of Yiwen et al. (2016) was to inform future efforts to model the location and spread of these species. Faxonius rusticus and P. clarkii were included in the Yiwen et al. (2016) study and we used their findings as the guide for variable selection in our models (Table 1). Variables determined most relevant by Yiwen et al. (2016) for F. rusticus and P. clarkii consisted of climatic data from both WorldClim (Fick and Hijmans 2017) and the Intergovernmental Panel on Climate Change (IPCC) (University of East Anglia Climatic Research Unit et al. 2020), topographic data derived from USGS Hydrosheds (Lehner et al. 2008), and vegetation data from USGS Landcover (Broxton et al. 2014). Per Yiwen et al. (2016), variables were not assessed or removed based on collinearity, to prevent potential violation of assumptions in their methodology (e.g., multivariate normality). All 29 variables were available in 30 arc-seconds resolution (1 km²) files, except for variables from IPCC, which we resampled from their dataset to 1 km² using the nearest-neighbor interpolation technique (Accadia et al. 2003) in ArcGIS Pro™. For the final models, we used 21 environmental variables for F. rusticus and 23 for P. clarkii (Table 1; additional variable descriptions available in Appendix 1: Table A1).

To project suitable habitat into the future we used four climate change prediction scenarios from the IPCC’s 5^th Assessment Report (IPCC 2014) under the Hadley Global Environment Model 2 Earth System (HadGEM2-ES) circulation model (Collins et al. 2008). We chose the HadGEM2-ES because it proved reliable in the Coupled Model Intercomparison Project Phase 5 (CMIP5) centennial experiments, including multiple simulations with different greenhouse gas (GHG) emission predictions (Jones et al. 2011). HadGEM2-ES also possesses future climate predictions at 1 km² resolution which eliminated the need to downscale lower resolution files and provided us with more detailed future environmental variables. The IPCC’s 5^th Assessment Report contains a series of four representative concentration pathways (RCPs) of GHG emission scenarios (RCPs 2.6, 4.5, 6.0, 8.5). Scenarios range from a strong reversal of GHG emissions (RCP 2.6), to a constant increase of GHG emissions (RCP 8.5) between the years 2000 and 2100. We used RCP 4.5 and RCP 8.5 for future projections. The RCP 4.5 scenario is considered an intermediate outcome with a decline in GHG emissions beginning in 2040, while RCP 8.5 is representative of an unmitigated rise in emissions through the year 2100 (IPCC 2014). For both RCPs, we used the same occurrence data to generate SDMs for three time periods: 2020 (present-day), 2050 (average predicted climates for years 2041–2060), and 2070 (average predicted climates for years 2061–2080). In total, we produced five maps per species, representing three time periods (2020, 2050, and 2070) and two emission scenarios (RCP 4.5 and 8.5).

Current and future distributions were modeled and projected using MaxEnt version 3.4.3 (Phillips et al. 2020). MaxEnt uses presence-only species occurrence data and environmental variables to predict suitable habitat for the species (Phillips et al. 2006). We chose MaxEnt for generating SDMs because of its customizability and because it has been demonstrated to have high performance compared to other approaches when presence-only occurrence data is available (Elith et al. 2010; Merow et al. 2013). Models were created by randomly assigning 75% of the occurrence data to a training-set and the remaining 25% to the testing-set, which allows for validation and calibration of models (Elith et al. 2006). For each SDM, we used the 10-replicate approach to create 10 models and average them into a single model. This method is commonly employed to remove variation that may occur in individual models (Merow et al. 2013; West et al. 2016; Godefroid et al. 2020). We retained most default settings within MaxEnt but made two modifications to reduce potential for overfitting (Radosavljevic and Anderson 2013). First, we adjusted the regularization multiplier to a value of four, which is an increase from its default value of one. Regularization is a feature within MaxEnt that can be scaled to adjust how much predictive weight is given to occurrence points (i.e., a low regularization value produces tight and potentially overfitted models, while a high value produces a more generalized predictive model), and a Regularization value of four has been demonstrated to produce accurate models that are less prone to overfitting (Radosavljevic and Anderson 2013). Second, we implemented a 10^th percentile presence threshold, which omits regions with habitat suitability values less than the suitability values for the lowest 10% of species occurrence records. This threshold setting assumes the 10% of occurrences existing in least suitable habitat are not reflective of the species’ overall habitat range. This represents a more conservative (less permissive) modeling approach because it removes the influence of extreme values and outlier occurrences. We elected to use this method as the resulting models are less sensitive to extreme values than other threshold rules and it offered balance to the higher regularization value (Radosavljevic and Anderson 2013). For future projections, we matched present-day WorldClim variables with HadGEM2-ES future climate variables in MaxEnt, allowing the application to generate predictions based on changes in the climate variables. Other variables (e.g., land cover and hydrology data) were held constant between present and future projections, as future projection data does not exist for these variables. Future projections for P. clarkii were created without the frs (Frost days) variable, which had high model contribution (Table 2), but no future analogue (see Discussion). Pearson correlation showed frs to have a strong negative correlation (-0.91) with Bio6 – which did possess a future projection – so we used Bio6 to substitute frs in P. clarkii’s future projections.

To validate models, we calculated the area under the curve (AUC) of receiver operating characteristic (ROC) plots. AUC output was generated within MaxEnt by evaluating the randomly selected 25% testing-set of occurrence data against the model created with the training data. AUC values and ROC plots are common tools for validating SDMs’ ability to accurately discriminate between localities where a species is present and localities where a species is absent (Elith et al. 2006), but additional validation methods are often advised, since MaxEnt, a “presence-only model”, relies on pseudo-absences (Lobo et al. 2008; Yackulic et al. 2012). We also calculated the continuous Boyce Index as a means of assessing the predictive abilities for our models, which can be used for absolute evaluations of presence-only SDMs (Hirzel et al. 2006). The continuous Boyce Index is a correlation value with a range between -1 and 1, where values below 0 indicate modeling errors, values near 0 indicate model prediction ability no better than random, and values near 1 indicate model predictions accurately capture and predict occurrences (Hirzel et al. 2006). We performed continuous Boyce Index calculations in R (version 4.1.2) using the “ecospat” package (Di Cola et al. 2017).

The final SDMs predicted habitat suitability on a scale from 0 and 1, with 0 being lowest probability and 1 being highest probability that a 1 km² cell contained suitable habitat for the species. Based on similar analyses performed in other studies, we organized the prediction outputs into four classes of suitability (Qin et al. 2017; Keliang et al. 2018): Least Potential (<0.2), Moderate Potential (0.2–0.4), Good Potential (0.4–0.6), and High Potential (0.6–1.0). We then used the reclassify tool in ArcGIS Pro™ to make presence/absence maps for habitat with Good or High Potential (>0.4) from each SDM, such that each pixel with a value of 1 represented 1 km² of habitat with a predicted Good or High Potential. Reclassification maps were compared between the 2020 (present day projections), 2050, and 2070 SDMs to evaluate how predicted ranges may change over time. Similar methodology for predicting habit suitability has been utilized to predict the presence of Dreissena in the Great Lakes (Zanatta et al. 2024). For both species and each climate change scenario, we subtracted the predicted Good and High Potential habitat in future projections from the Good and High Potential habitat in present day predictions to assess changes in predicted area (km²).

The MaxEnt variable analysis for model contribution, averaged from 10 replicates, showed the four highest contributing environmental variables for F. rusticus to be annual precipitation (Bio12, 36.9%), mean temperature of warmest quarter (Bio10, 22.9%), cloud cover (cld, 14.6%), and annual mean monthly max temperature (tmx, 10.9%), accounting for a collective 85.3% contribution to the SDM (Table 2). The four environmental variables offering the greatest contribution to P. clarkii’s SDM were frost days (frs, 71%), breakline emphasis (be_grd, 11%), minimum temperature of coldest month (Bio6, 4.3%), and annual mean monthly max temperature (tmx, 3.1%), accounting for a collective 89.4% contribution to the SDM (Table 2).

Response curves illustrate the effects of single variables on the predictive power of the models (Phillips et al. 2006). Using response curves (Fig. 1), we analyzed positive prediction thresholds (Probability of Presence >0.4) for each species’ top contributing environmental variables. For F. rusticus, annual precipitation (Bio12) ranged from 700 to 1,250 mm, mean temperature of the warmest quarter (Bio10) ranged from 20.5 to 28 °C, cloud cover (cld) ranged from 62 to 82%, and annual mean monthly maximum temperature (tmx) ranged from -8 to 7 °C. For P. clarkii, frost days (frs) ranged from 0 to 27 days, breakline emphasis (be_grd) ranged from -500 to 250 meters, minimum temperature of the coldest month (Bio6) ranged from -7 to 10 °C, and annual mean monthly maximum temperature (tmx) ranged from 2 to 21 °C.

Figure 1. 

Response curves for environmental variables that provided greatest contributions to species distribution models for F. rusticus (A–D) and P. clarkii (E–H). Response curves analyze the range of a single variable’s values on the predictive power of a model. Curves show the mean response of the 10 replicate MaxEnt runs with occurrences from predictive models represented as “X” at 0.0 on the Y-axis.

Mean (± standard deviation) AUC values from the 10 replicate present-day models were 0.945 ± 0.009 for F. rusticus and 0.935 ± 0.007) for P. clarkii. These values indicate both species’ SDMs have robust performance. The continuous Boyce Index correlation values supported this with values of 0.994 for F. rusticus and 0.973 for P. clarkii. We conclude that our present-day models accurately captured the current ranges of both species in North America (Figs 2, 3) (GBIF 2020; USGS 2020).

Figure 2. 

Present-day (2020) species distribution models for the non-native habitat of Faxonius rusticus.

Figure 3. 

Present-day (2020) species distribution models for the non-native habitat of Procambarus clarkii.

The present-day SDM for F. rusticus predicted areas of Good and High Potential habitat occurring predominantly in the U.S. Midwest and Northeast, with an area of particularly high habitat suitability surrounding the Great Lakes region. Smaller clusters of Good Potential habitat were identified along the Rocky Mountain range between the Western U.S. and Canada, the southern border of Ontario, and within Canada’s Gulf of St. Lawrence region (Fig. 2). The combined surface area of F. rusticus’ present-day Good and High Potential habitats is 2.57 × 10⁶ km².

Present-day areas of predicted Good and High Potential habitat for P. clarkii occur predominantly along U.S. coastal regions and in the Southeast. Notable regions with high levels of Good Potential habitat are the southern Great Lakes region, the U.S. Pacific Northwest, and north-central Mexico (Fig. 3). The combined surface area of P. clarkii’s present-day Good and High Potential habitats is 3.48 × 10⁶ km².

Reclassified cutoff maps (Figs 4, 5) derived from future projections for RCP 4.5 and RCP 8.5 (Appendix 1: Fig. A2) show changes in predicted range and areas of Good and High Potential habitat (Prediction >0.4) for F. rusticus and P. clarkii. For F. rusticus there is a contraction in areas of Good and High Potential habitat (Fig. 4) across all projections, with some northern range expansion into Canada’s west and Great Lakes region (most prominent in RCP 8.5 – 2070), and there is a predicted reduction of total Good and High Potential habitat in all projections compared to the present-day model (Table 3).

Figure 4. 

Reclassified cutoff maps for Faxonius rusticus. A non-white pixel indicates area of habitat with a predicted Good (gray) or High Potential (black) for suitability (prediction >0.4), and a value of 0 (white) indicates area classified below the 0.4 prediction threshold. Reclassified cutoff maps were made for the years 2050 and 2070 for both RCP 4.5 and RCP 8.5.

Figure 5. 

Reclassified cutoff maps for Procambarus clarkii. A non-white pixel indicates area of habitat with a predicted Good (gray) or High Potential (black) for suitability (prediction >0.4), and a value of 0 (white) indicates area classified below the 0.4 prediction threshold. Reclassified cutoff maps were made for the years 2050 and 2070 for both RCP 4.5 and RCP 8.5.

Reclassified cutoff maps for P. clarkii (Fig. 5) showed northern range expansion for Good and High Potential habitat into the Great Lakes region, the western coast of Canada, and the southwestern portion of Alaska. All projections show an increase in total Good and High potential habitat (Table 3).

For both crayfish species, our 2020 models identified numerous large regions that have the conditions of currently lacking species occurrence data but have Good or High Potential for habitat suitability (Figs 1, 2). Areas that met these two conditions for F. rusticus include segments of the central and western border between the U.S. and Canada, portions of the lower Great Lakes, the southern border of Ontario, Canada’s Gulf of St. Lawrence region, and the lower half of British Columbia. Areas meeting these conditions for P. clarkii include portions of the lower Great Lakes, major river systems and aquatic habitats throughout the U.S. Pacific Northwest including the Columbia River Gorge, western Mexico, and portions of British Columbia. These respective areas may represent locations where F. rusticus and P. clarkii are already established but have not yet been detected. Alternatively, they may represent suitable habitats where the two invasive crayfish have not yet become established.

The high AUC values and continuous Boyce Index calculated for both models strongly indicate that the models are robust and provide useful habitat predictions; however, presence-only SDMs are prone to placing higher prediction values near areas with a high density of occurrence records and lower prediction values in areas of few occurrence records (Veloz 2009; Radosavljevic and Anderson 2013). This is evident within the 2020 models, where both species’ SDMs contain areas of minimal presence prediction that overlap with a small number of isolated occurrence records. These spread-out and under-predicted occurrence records may be the result of humans unwittingly transporting crayfish from larger populations to new localities, as vectors like fishing, the pet trade, and the disposal of organisms from classroom aquariums have all been identified as drivers for invasive crayfish introductions (Keller and Lodge 2007; Larson and Olden 2008; DiStefano et al. 2016). Because presence-only models are inherently weighted towards spatially correlated data (Veloz 2009), these isolated records may indicate that there is a wider range of suitable habitat than accounted for within the models presented in the current study; however, there is a lack of sufficient occurrence data to verify this. These regions should be evaluated to see if sufficient sampling is currently being performed within them. This would offer greater context to the models, as well as determine if these areas possess larger populations of crayfish and lack sufficient sampling, or if these areas do not possess viable crayfish populations and represent habitat of low suitability. Despite these limitations, the 2020 models’ predictions for areas of higher suitability show general range agreement with previous studies that examined the distributions of these two crayfish species (Morehouse and Tobler 2013; Yiwen et al. 2016; Egly et al. 2018; Zhang et al. 2019).

While our models share general distribution trends with other SDMs developed for F. rusticus (Morehouse and Tobler 2013; Yiwen et al. 2016) and P. clarkii (Yiwen et al. 2016; Zhang et al. 2019), our models possess several distinctions from previous work. In Yiwen et al. (2016), present-day SDMs were generated for both species; however, these were generated at a coarser resolution of 2.5 arc-minutes (~4.5 km²) and models were not assessed for degrees of habitat suitability or future changes. Our models were generated at a finer resolution and utilized to assess general distribution trends over space and time, both for the present and under several climate change scenarios. For F. rusticus, Morehouse and Tobler (2013) assessed present and future distributions for the species but used variables solely from the World CLIM dataset and spatially examined only the United States portion of North America. In Morehouse and Tobler (2013), a single future projection was calculated by comparing present-day native and invasive distributions to infer where uninvaded, yet habitable, range may be, rather than informing models with future climate change variables. For P. clarkii, Zhang et al. (2019) utilized World CLIM environmental variables and a single measure of human influence in the modeling process, and maintained default MaxEnt settings, which resulted in a present-day model that shared general range agreement with ours but predicts a much broader area being highly suitable. This is likely a result of the coarser resolution (~10 km²) and lack of hydrological variables used in Zhang et al. (2019). The future projections in Zhang et al. (2019) examined a combination of future climatic models for RCP 2.6 and RCP 8.5. In their study, the HadGEM2-ES model for RCP 8.5 showed a decrease in suitable habitat for 2050 and 2080 projections. This is the opposite of our findings for P. clarkii’s suitable habitat projections. These differences could possibly occur because the present-day model in Zhang et al. (2019) used default settings and produced a broader distribution range, which may have experienced greater predicted habitat reductions when future environmental variables were applied. Alternatively, the addition of hydrology variables in our models may have driven our fundamentally different future projections.

For both species, temperature related variables made some of the largest contributions to the final SDMs (Table 2). MaxEnt’s response curves (Fig. 1) show that when these variables were analyzed in isolation, a rough bell-curve shape for increased presence probability forms across each variables’ range of values. For F. rusticus, the response curve for mean temperature of the warmest quarter shows an increase in F. rusticus presence probability from 14 °C, a peak between 22 and 27 °C, and a zero-presence probability at and above 35 °C. These ranges correspond with previous studies focused on relationships between F. rusticus and water temperatures (Layne et al. 1985; Mundahl and Benton 1990; Claussen et al. 2000). For example, in lab experiments looking at F. rusticus survival across a range of water temperatures Mundahl and Benton (1990) found a bell-shaped curve with a peak between 20 and 25 °C. Faxonius rusticus slow their metabolism in the winter and possess the ability to thermally acclimate to seasonal shifts in temperature, allowing them to regulate between water temperatures of 5 and 25 °C (Layne et al. 1985). Water temperatures outside of these ranges compromise F. rusticus locomotive abilities, hindering foraging and predator avoidance (Layne et al. 1985; Claussen et al. 2000). Although our metric (mean temperature of the warmest month) measures air temperature, water temperatures in streams have been shown to increase 0.6 to 0.8 °C for every 1 °C increase in air temperature (Morrill et al. 2005) and lakes have similarly displayed a 0.4 to 1 °C increase for every 1 °C increase in air temperature (Robertson and Ragotzkie 1990). These lab observations support the variable in our model and indicate mean temperature of the warmest quarter creates suitability thresholds for F. rusticus via changes in water temperature.

Another temperature related variable, frost days, made the greatest contribution to the P. clarkii SDM (Table 2). The response curve for frost days sharply increased from a baseline at 0 days, peaks at 1 day, then decreases until the maximum of 31 days. This variable is related to temperature with an increase in frost days indicating more frequent cold temperatures (O’Donnell and Ignizio 2012). When comparing frost days to the third and fourth most contributing variables (min temperature of the coldest month and annual mean monthly max temperature), the three curves suggest that high and low temperature extremes decrease the presence probability of P. clarkii. This is consistent with literature identifying the optimal water temperature range for P. clarkii to be between 21 and 30 °C (Peruzza et al. 2015). Like F. rusticus, P. clarkii can thermally acclimate to winter temperatures; however, it is relatively more tolerant to heat and more limited by cold, experiencing mortality in water temperatures of 10 °C and lower but not incurring a slowed metabolism until 30 °C (Powell and Watts 2006).

All future projections displayed geographic shifts in habitat classified as either “Good Potential” or “High Potential” (Table 3), indicating a warming planet can influence the size and location of these species’ ranges. Faxonius rusticus’ future projections consistently showed a reduction of higher potential habitat and varied degrees of northern expansion (Fig. 4), implying the effects of climate change on North American habitat will strongly alter F. rusticus’ suitable range. Future projections for P. clarkii displayed major expansion of higher potential habitat (Table 3) with northward shifts in range also observed (Fig. 5). The increase of higher potential habitat for P. clarkii contradicts the findings in Zhang et al. (2019), which predicted future range reductions; however, these differences are likely due to different modeling techniques (see Models of current distribution). For both species, projections made for RCP 8.5 contained more northern expansion of higher potential habitat than projections made for RCP 4.5. RCP 8.5 represents an unmitigated increase in GHG emissions over the next 100 years, resulting in higher global temperatures than the RCP 4.5 scenario (IPCC 2014). Our results indicate that the higher temperatures and precipitation levels associated with RCP 8.5 will lead to larger geographic shifts of higher potential habitat for both species.

Native species of crayfish have been shown less able than non-native species to acclimate to higher temperatures. As climate change progresses in North America, native species will likely experience declines in suitable habitat (Capinha et al. 2013; Dyer et al. 2013; Hossain et al. 2018). While native habitat is predicted to decline, our results highlight an additional concern to native species by anticipating an increased overlap between invasive and native crayfish habitat as invasive ranges expand or shift, which has also been observed in projections of future European crayfish populations (Dyer et al. 2013). An increase of shared habitat will impose a suite of stressors on native species, as F. rusticus and P. clarkii are aggressive competitors, carriers of disease, and destructive to new habitats (Twardochleb et al. 2013; Larson et al. 2019; O’Shaughnessey and Keller 2019).

The maps showing future projections of Good and High Potential habitat (Figs 4, 5) predict northern expansion of both invasive crayfish species into areas that are currently predicted to be uninhabitable. For F. rusticus, our results predict an expansion of suitable habitat into Canada and small portions of Alaska. Climate change is expected to have an augmented increase in the temperature of this region due to a positive feedback loop from reduced snow cover (Christensen et al. 2007). Canada contains many uninvaded freshwater bodies and concern has already been raised about their degradation from the future spread of F. rusticus (Jansen et al. 2009; Phillips et al. 2009) – our results serve to support this concern. For P. clarkii, our results predict northern range shifts to be most prominent in the U.S. Midwest and the western coast of Canada. An important region of expansion is predicted to be the Great Lakes, which comprise 84% of North America’s fresh surface freshwater supply (EPA 2019). Over the last decade, this region has experienced an increase of P. clarkii establishments and the potential for further spread has been cited as a threat to the biodiversity, habitat quality, and economic services provided by the Great Lakes (Peters et al. 2014; Bunk and Egeren 2016; Egly et al. 2018; Smith et al. 2018; O’Shaughnessey and Keller 2019). The future projections for intermediate reductions (RCP 4.5) and no reductions (RCP 8.5) to GHG emissions both suggest spread and establishment of P. clarkii will continue in the Great Lakes region, resulting in overall decreases to native biodiversity and habitat quality (Twardochleb et al. 2013; DiStefano et al. 2016).

SDMs created with environmental variables possess some limitations because many types of abiotic and biotic relationships that may influence a species’ distribution are not accounted for (e.g., predation, disease, pollution, etc.), especially human vectors, which play a major role in shaping the spread of modeled species (Hui 2023). The complex relationships among these variables, especially on a scale as large as North America, make them difficult to include in models (Kissling et al. 2012). Also, the continuation of human-caused introductions, along with climatic shifts in habitat, lend to a non-equilibrium in both species’ predicted ranges where suitable habitat is likely to change (Hui 2023). As a result, we believe our present-day models best serve as only a snapshot of predictability in habitat suitability based on the occurrence records available at the time of modeling.

Another limitation of SDMs is that presence-only models are prone to being heavily weighted towards areas with high densities of occurrence records and large databases often contain records from areas of high sampling (Veloz 2009; Hui 2023). We took precautions to correct for and minimize any biases by spatially rarefying occurrence data and increasing MaxEnt’s regularization multiplier (Brown 2014; Radosavljevic and Anderson 2013); however, it is not possible to be certain that these adjustments completely addressed the issues. Despite these limitations, the high AUC and continuous Boyce Index values for our models support the accuracy and utility of our habitat predictions.

Projecting species distribution data with future climatic variables is a popular analysis method, but it is susceptible to accuracy issues from potentially exaggerating underlying model biases or not considering potential adaptations made by species (Santini et al. 2020). For example, a European population of P. clarkii displayed adaptations in the timing of its life stages, allowing it to acclimate to water temperatures below typical levels of suitability (Chucholl 2011). In addition, future projections may be hampered by the variables that did not possess future analogs (e.g., IPCC and HydroSHEDS), as they may not capture the full influence of changes those variables may make to predicted range. Variables without future analogs likely increase predictive power of present-day models but may cause a reduction in predicted suitable habitat amongst future projections. In one instance, we corrected for this potential reduction by replacing a highly influential variable without a future analog (frs – frost days) with a highly correlated future variable analogue (Bio6 – minimum temperature of coldest month). Even with this correction, we believe our future models may under-predict total suitable habitat. For these reasons, we acknowledge limitations in how we can and should interpret our models and future projections.

The SDMs made in this study should not be considered high-precision maps detailing habitat suitability down to a specific 1 km² resolution but rather viewed as tools to assess general spatial trends in habitat between snapshots of the present-day and general predictions for the future. For example, the future projection models predict a growth of habitat with higher suitability for P. clarkii throughout most of the Great Lakes region (Fig. 5). While P. clarkii may not become fully established in every km² of this predicted area, the general trend of expansion into new habitats is meaningful for future research and conservation efforts. Understanding where and when habitat shifts may occur is important for monitoring and preventing negative impacts from aquatic invasive species in North America. For both species of crayfish, the models provide valuable insight into present-day regions with high suitability that requires more stringent sampling efforts and regions that will likely become prone to new invasions in the future.