Lineage‐level distribution models lead to more realistic climate change predictions for a threatened crayfish

As climate change presents a major threat to biodiversity in the next decades, it is critical to assess its impact on species habitat suitability to inform biodiversity conservation. Species distribution models (SDMs) are a widely used tool to assess climate change impacts on species’ geographical distributions. As the name of these models suggests, the species level is the most commonly used taxonomic unit in SDMs. However, recently it has been demonstrated that SDMs considering taxonomic resolution below (or above) the species level can make more reliable predictions of biodiversity change when different populations exhibit local adaptation. Here, we tested this idea using the Japanese crayfish (Cambaroides japonicus), a threatened species encompassing two geographically structured and phylogenetically distinct genetic lineages.


| INTRODUC TI ON
Amidst the current global climate emergency (Ripple et al., 2019), our planet is changing rapidly in an unpredictable manner. According to the latest assessment by the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services (IPBES, 2019), climate change is one of the five most influential direct drivers of biodiversity decline, along with habitat loss, direct exploitation of organisms, pollution, and invasion of alien species. It has been demonstrated that climate change indiscriminately threatens species in virtually all ecosystems-terrestrial (Burrows et al., 2011;Pinsky et al., 2019), subterranean Mammola, Piano, et al., 2019), freshwater (Woodward et al., 2010) and marine (Burrows et al., 2011;Pinsky et al., 2019). Climate change can affect species by altering physiological performance, phenology, geographical distribution and species interactions, among others (review by Hughes, 2000;Walther et al., 2002). For instance, there is mounting evidence that species are rapidly rearranging their distributions along elevational and/or latitudinal gradients in response to changing climates (Burrows et al., 2011;Lenoir et al., 2020). To better protect and manage biological resources in an era of climate change, it is therefore fundamental to gain a better understanding of how climate change will redefine the geography of life (Román-Palacios & Wiens, 2020).
Species distribution models (SDMs) are powerful tools in this regard, as they can estimate species' habitat suitability by determining statistical relationships between species' range data and environmental predictors, and can also be used to forecast how suitable areas may vary under different climate change scenarios (Araújo et al., 2019;Booth et al., 2014;Elith & Leathwick, 2009;Guisan et al., 2017). To date, SDMs have often been used to investigate the potential impacts of climate change on target species of conservation concern (Araújo et al., 2019;Elith & Leathwick, 2009;Guisan & Thuiller, 2005;Zhang, Mammola, Liang, et al., 2020). For instance, Zhang, Mammola, Liang, et al. (2020) developed a SDM for the Critically Endangered Chinese giant salamander Andrias davidianus (Blanchard, 1871)-the world's largest extant amphibian-and found that this species might be extremely vulnerable to climate change and might lose more than two-thirds of its suitable habitat in the future.
As the term "species distribution model" suggests, the species level is the most commonly used taxonomic unit in SDMs (Qiao et al., 2017;Smith et al., 2019). One basic assumption underlying species-level SDMs is "niche conservatism" (Guisan et al., 2017;Peterson et al., 1999), positing that species should exhibit uniform responses to climate and retain similar niche characteristics over space and time. However, recent studies have emphasized that taxonomic resolution below (or above) the species level should be considered in SDMs for some systems (e.g. Collart, Hedenäs, Broennimann, Guisan, & Vanderpoorten, in press;Guisan et al., 2017;Pearman et al., 2010;Serra-Varela et al., 2015;Smith et al., 2019). Different populations of the same species inhabiting geographically and ecologically distinct environments may experience local adaptation and thus have niches that are divergent to some extent (e.g. . For such systems, modelling habitat suitability below the species level should lead to more accurate range estimates and climate change projections (Collart et al., in press;Pearman et al., 2010;Smith et al., 2019). In threatened species conservation, for example, climate change responses of phylogeographic lineages have recently been considered for some taxonomic groups from terrestrial and marine systems (e.g. Cacciapaglia & van Woesik, 2018;Collart et al., in press;D'Amen et al., 2013;Lecocq et al., 2016), but such studies have rarely been conducted in freshwater systems (but see Zhao et al., 2020).
Here, we examined how predictions of climate change responses can differ when taking into account intraspecific genetic heterogeneity for a threatened freshwater arthropod, the Japanese crayfish Cambaroides japonicus (De Haan 1841) (see Species description in Methods for details). Most freshwater crayfish have specialized habitats and limited mobility, which contribute to their vulnerability to climate change (Hossain et al., 2018). Cambaroides japonicus is largely sedentary and its intraspecific phylogeography has been resolved for its entire distribution. There are two distinct genetic lineages (eastern and western) in C. japonicus delineated by partial sequences of mitochondrial and nuclear DNA (Koizumi et al., 2012). On the basis of molecular phylogeographic analysis, Koizumi et al. (2012) inferred that the two genetic lineages have undergone distinct evolutionary histories over several million years in different geographical regions characterized by differential climatic factors and topography. Given these different histories over a substantial time-scale, it is likely that these lineages have experienced at least some local adaptation, which would lead to an expectation of differential responses to climate change. In this case, it is important to compare future range projections at the levels of species and lineage to assess the impacts  (Benito Garzón et al., 2019;Peterson et al., 2019;Smith et al., 2019).
In this study, we quantified realized niches (i.e. the portion of the fundamental niche currently used by the species; Guisan et al., 2017;Soberón & Nakamura, 2009), developed SDMs and made future predictions to examine how climate change might influence C. japonicus by constructing species-level versus genetic lineage-level models.
We sought to address the following hypotheses: (1) there will be divergence in realized niches between the eastern and western lineages of C. japonicus, and (2) there will be considerable differences between the future projections of the SDMs at the species and lineage levels. Our results address the uncertainty in climate change predictions for species with high intraspecific variation in genetic diversity, as well as highlight the importance of developing SDMs below the species level. Moreover, our results should provide a useful guide for designing conservation strategies for C. japonicus lineages.

| Species description
Cambaroides japonicus is designated a threatened species (Category II; Vulnerable) in the Ministry of the Environment, Japan's Red List (Ministry of the Environment, Japan, 2020) and data-deficient but with a decreasing population trend in the IUCN Red List (IUCN, 2020).
The species is the only crayfish species native to Japan, endemic to freshwater areas of Hokkaido and northern Honshu in northern Japan (Miyake, 1982). Natural populations of this species have declined substantially in recent decades. Multiple stressors including habitat loss, water quality degradation, and predation and competitive exclusion by invasive crayfish are thought to be major factors for the observed decline of C. japonicus (Kawai, 2003;Usio et al., 2001).
Cambaroides japonicus is also a popular pet for ornamental aquariums and is distributed through Internet auctions and aquarium shops (Kawai, 2003;N. Usio, personal observation). Nevertheless, there is no environmental law to protect this species except for the southernmost population in Akita Prefecture which was designated a national monument (Miyake, 1982). In addition to current pressures, another important yet largely neglected threat to C. japonicus might be climate change, given that the distributions of other crayfish species were predicted to be driven by temperature and other climatic conditions (e.g. Capinha et al., 2013;Gallardo & Aldridge, 2013;Zhang, Capinha, Usio, et al., 2020).

| Study area and species distribution data
Considering the known distribution of C. japonicus, we selected northern Japan (138°E-146°E, 40°N-46°N) as our study extent and restricted the analysis to freshwater areas ( Figure S1). We identified freshwater areas with MERIT Hydro (Yamazaki et al., 2019), a global hydrography dataset at 3 arc-second resolution (~90 m at the equator), and then resampled this layer to the spatial resolution of our environmental predictors (30 arc-seconds; ~1 km at the equator).
We obtained occurrence data for C. japonicus from published literature (Koizumi et al., 2012), data collections by the Bihoro Museum, Hokkaido, Japan (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019), and field surveys (N. Usio 2003, K. Tanaka 2007, S. Niwa 2000-2019, and N. Ichijo 2003-2012. From these sources, we collected a total of 497 presence records. Based on the phylogenetic results of Koizumi et al. (2012), which determined a clear spatial partition between the two lineages, we divided the study extent into two regions by drawing a boundary of straight lines to assign presence data to either the western or eastern lineage ( Figure S1). Note, however, that eight presence records at the boundary were unable to be assigned to any lineage due to an absence of phylogenetic information; these were excluded from lineage-level SDM analyses. Overall, we assigned 325 records to the western lineage and 164 records to the eastern lineage.
In accordance with standard practices for SDMs, we cleaned and spatially thinned the presence records using a 5 km thinning distance to avoid spatial sampling bias (Kramer-Schadt et al., 2013) (see Supporting Information for Presence data processing). After processing the presence data, we retained 113 records for the species-level SDM, 47 records for the eastern lineage SDM, and 60 records for the western lineage SDM ( Figure S1).

| Predictor variables
For modelling realized niches and potential distributions for C. japonicus, we initially considered 14 environmental predictor variables with a spatial resolution of 30 arc-seconds: eight bioclimatic predictors from CHELSA (Karger et al., 2017), three land use and land cover predictors (Li et al., 2016), and three topographic predictors (Amatulli et al., 2018;Hengl, 2018) (see Table 1 for details). We chose these variables based on our expert opinion that bioclimatic extremes, seasonality, and means likely influence the species' physiological performance; land use and land cover can have strong relationships with its habitat preferences; and water velocity and rain infiltration are closely related to topography and can affect habitat quality. The Japanese crayfish is typically found in small streams, ponds or mesic areas (where groundwater emerges) in broadleaf forests (Kawai, 2003;Usio, 2007). Therefore, we also generated a broadleaf forest variable using the EarthEnv dataset (Tuanmu & Jetz, 2014) by summing the layers for evergreen broadleaf trees, deciduous broadleaf trees and mixed/other trees (Table 1). We checked for collinearity among 15 predictors by calculating the pairwise Pearson's correlation coefficient (r) and retained one predictor when two or more were highly correlated (i.e. |r| > .70) (Zhang, Capinha, Karger, et al., 2020;. Given biological importance, collinearity analysis results, and data availability under present-day and future scenarios, we selected the following eight predictors for modelling: maximum temperature of warmest month, minimum temperature of coldest month, precipitation of wettest month, precipitation of driest month, slope, and percent cover of forest, water, and wetland ( Figure S2).
We addressed the inherent uncertainty in future climate projections by considering four different global circulation models (CESM1-BGC, CMCC-CM, MIROC5, MPI-ESM-MR) (Karger et al., 2017) and obtained future projections for land use and land cover from Li et al. (2016). We assumed that slope would remain unchanged in the future.

| Niche comparisons
We used n-dimensional hypervolumes to quantify the realized niche space for the eastern and western lineages. In brief, we performed a principal component analysis on the eight selected predictors and used the first four principal components to describe species niches as they accounted for more than 75% of the total variance ( Figure S3). We used the R package hypervolume (Blonder, 2019) to construct hypervolumes for the two lineages. We used the volume of each hypervolume as a measure of the size of the realized niche in multidimensional space (note that the volume of a hypervolume is a unitless measure). We then used the R package BAT  to evaluate niche overlap between the two lineages. In particular, we expressed differentiation between hypervolumes as the sum of two components: niche shifts (i.e. replacement of space between hypervolumes) and niche contraction/expansion (i.e. net difference between hypervolumes) (Carvalho & Cardoso, 2020; . Niche differentiation varies from 0 (complete overlap between hypervolumes) to 1 (complete separation between hypervolumes).

| Species distribution modelling
We performed an extensive SDM analysis incorporating explora- com/jamiemkass/ENMeval; Muscarella et al., 2014). We evaluated each model using spatial block cross-validation, which is recommended for assessing model transferability, or how well models extrapolate to new environments (Roberts et al., 2017). We used the "block" partition option to divide the study region into four regions containing an equal number of presences: three blocks were used for model training and the withheld block for model validation, and this process was repeated until all blocks were withheld. Using new partitioning options in ENMeval 1.9.0, we specified spatial block partitions with different orientations for each occurrence dataset to best ensure each block represented discrete and contiguous areas: longitudinally for the eastern lineage, latitudinally for the western lineage and a combination of both for the species level ( Figure S7).
We chose optimal settings from among the candidate models by applying sequential criteria on performance metrics (Kass, Anderson, et al., 2020;Radosavljevic & Anderson, 2014). We first selected models with the minimum average 10% omission rate, or the pro- characteristic curve for validation data. The AUC value ranges from 0 to 1, with higher values indicating better model discriminationthis metric is problematic for determining the absolute performance ability of presence-background SDMs, though it is acceptable to use for relative comparisons across models with the same data (Lobo et al., 2008).
We additionally performed post hoc procedures on the selected models via tuning to assess performance with different methods, determine variable importance and marginal response curves, and visualize the extent of extrapolation across spatial blocks. We first assessed the predictive abilities of the selected models using the continuous Boyce index calculated on the full dataset, which unlike AUC can be used for absolute evaluations of presence-background SDMs (Hirzel et al., 2006). The continuous Boyce index varies from -1 to 1, whereby values above 0 indicate model predictions consistent with distribution data, values of 0 indicate performance no better than random, and values below zero refer to incorrect model predictions (Hirzel et al., 2006). We then ran a series of null model simulations (n = 100) using the complexity settings of each selected model to determine if the empirical models predicted validation data better than models built with random occurrences (results for simulations with 1000 iterations, which do not change our conclusions, are available in Figure S9). To do so, we used a recently described null model approach for SDMs available in ENMeval 1.9.0 that eval-  (Elith et al., 2010) using the rmaxent package (Baumgartner & Wilson, 2020) to determine how environmentally similar training data was to validation data for each spatial block.
We plotted histograms to visualize the extent of extrapolation that occurred during cross-validation for each model.
We used the selected models and environmental variables to predict current and future habitat suitability. We made predictions using Maxent's cloglog transformation (bounded by 0 and 1) (Phillips et al., 2017) and additionally made binary range maps by thresholding the continuous predictions by the 10% omission values. All analyses were performed in R (R Core Team, 2020).

| Comparison of realized niches
The volume of the four-dimensional hypervolume for the western lineage (1,971.14) was larger than that of the eastern lineage (1,848.67), and niche overlap between the two lineages was intermediate (0.50) (Figure 1). Niche differentiation was mainly due to contraction/expansion (>90%), while niche shifts only contributed marginally (<10%). Niche differentiation between the two lineages was mostly observable along PC2 (mainly explained by precipitation of driest month, minimum temperature of coldest month, and maximum temperature of warmest month) (Figures 1 and S3).

F I G U R E 1
The four-dimensional hypervolumes for western and eastern lineages of Japanese crayfish. To visualize the shape and boundary of the hypervolumes in two dimensions, a random selection of 10,000 stochastic points for each hypervolume was used

| Best-performing SDMs and variable importance
The three Maxent models we selected had different settings that resulted in different levels of complexity. The western and eastern lineage models were relatively simple (L1.5 with six non-zero coefficients and L0.5 with seven non-zero coefficients, respectively), while the species-level model was more complex (LQH3.0 with 12 non-zero coefficients) (Table 2 and Figure S4). All three models had continuous Boyce index values well above 0 (western: 0.65, eastern: 0.66, species: 0.92), indicating that each had predictions consistent with the distribution of presence data (Hirzel et al., 2006) (Table 2).
For both 10% omission rate and validation AUC, the western lineage (both p < .05) and species-level (both p < .01) models performed significantly better than null models, and effect sizes for the specieslevel model were larger for both metrics (Table S1 and Figure S5).
However, although the selected eastern lineage model performed the best on withheld spatial blocks compared to other candidate models, its cross-validation performance was not significantly better than null models (p = .094 for AUC and p = .933 for omission rate; Table S1 and Figure S5). Thus, there is evidence that the transferability of the eastern lineage model may be poorer than that of the other models, resulting in higher uncertainty for its future projections. The environmental variables with the highest permutation importance differed among models: precipitation of wettest month was shared by all, while precipitation of driest month was shared by the lineage-level models, slope was shared by the western lineage and species models, and others were particular to each model (eastern: minimum temperature of coldest month, species: forest) (Table S2). For most environmental variables, the marginal response curves were either positive or relatively neutral, though precipitation of the wettest month was negative for all models ( Figure S6). Regarding the MESS analysis, the eastern lineage model showed considerable extrapolation across spatial blocks, whereas

| SDM projections
Lineage and species SDMs resulted in different suitability projec- The lineage-level and species-level SDMs resulted in largely different projections for future conditions, with reasonable agreement between GCMs (Table 3 and Figure 3). Taking the MPI-ESM-MR for RCP 8.5 in the 2050s as an example, according to the lineage-level model predictions, the western lineage is predicted to experience range restrictions (9.2% decrease overall) in north-central Hokkaido, southern Hokkaido and northeast Honshu, while the eastern lineage is predicted to slightly expand its range (8.9% increase overall) in two distant areas in east-central Hokkaido and eastern Hokkaido. According to the species-level model, large areas in Hokkaido and northern Honshu are predicted to become unsuitable (27.0% decrease overall, 16.0% decrease in east, 32.4% decrease in west) (Figure 3 and Table 3).

| D ISCUSS I ON
In this study, we quantified realized niches for two genetic line-

| Intraspecific variation in SDMs
Species distribution models have been frequently used in biodiversity assessments with the most common application being to estimate habitat suitability at the species level (e.g. Araújo et al., 2019;Collart et al., in press;Elith & Leathwick, 2009;Guisan & Thuiller, 2005;Guisan et al., 2017). The species-level SDM relies on the "niche conservatism" assumption and does not take into account intraspecific phyloge-

Incorporating local adaptation and intraspecific variation into
SDMs stems from the recognition that populations of species inhabiting largely different habitats over significant time-scales will often show adaptations to their respective local conditions, resulting in intraspecific niche variation (Collart et al., in press;Smith et al., 2019). While these local adaptations are known to result in clear morphological differences among populations in some organisms (e.g. DeWoody et al., 2015), crayfish are known to exhibit a wide range of body shape and colour differences among intraspecific populations, at least in part due to abiotic habitat characteristics (including light availability and water chemistry) and food resources therein (Holdich, 2012;Thacker et al., 1993). In crayfish, local adaptations to climate may be rather represented by intraspecific genetic differences that are only discernible via phylogenetic analysis, as in our study species (Koizumi et al., 2012).

| Conservation implications
Our findings have important implications for designing conservation strategies for C. japonicus. With the purpose of preserving genetic diversity for conservation efforts, the concept of the Evolutionarily Significant Unit, also termed "Unified Species Concept" (De Queiroz, 2007), was developed to define separately evolving lineages below the species level as the only reasonable conservation unit (Ryder, 1986 (Ikeda et al., 2016). Long-term population data can also help validate our SDM projections, which should be treated as hypotheses in need of testing.
In conclusion, this study represents a first attempt to explore climate change impacts on C. japonicus, which to date have remained largely unexplored for this species. As multiple nuances exist, further studies are required to test whether intraspecific variation in environmental requirements between the two lineages reflects different adaptation at the physiological levels. Future experimental studies and monitoring programmes would be needed to explore these possibilities.

ACK N OWLED G EM ENTS
We thank César Capinha (University of Lisbon, Portugal) and Geng Qin (

CO N FLI C T O F I NTE TR E S T
The authors declare there are no competing interests.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13225.

DATA AVA I L A B I L I T Y S TAT E M E N T
Presence records of the Japanese crayfish Cambaroides japonicus are presented in Figure S1. Environmental predictors can be retrieved from online databases (see details in Table 1