Integrating population genetics and species distribution modelling to guide conservation of the noble crayfish, Astacus astacus, in Croatia

The noble crayfish, Astacus astacus, is an indigenous European freshwater species. Its populations show significant declines caused by anthropogenic pressure on its habitats, climate change and the spread of invasive species. Diminishing populations’ trends and loss of genetic diversity highlight the need for effective conservation that will ensure their long-term survival. We combined population genetics and species distribution modelling (SDM) to reveal the impact of climate change and invasive species on the noble crayfish, and to guide future conservation programs of current populations. Our study showed that Croatian populations of A. astacus harbour an important part of species genetic diversity and represent significant genetic reservoir at the European level. The SDM results predicted substantial reductions of suitable habitats for A. astacus by the 2070; only 13% of its current potential distribution is projected to remain stable under pessimistic Representative Concentration Pathway (RCP 8.5) emission scenario. Moreover, most of the populations with high genetic diversity are located in the areas predicted to become unsuitable, and consequently have a high probability of being lost in the future. Further, SDM results also indicated considerable decrease of future habitat suitability for invasive crayfish species in Croatia, suggesting that climate change poses a major threat to already endangered A. astacus. The obtained results help in the identification of populations and areas with the highest conservation value which should be given the highest priority for protection. In order to preserve present diversity in areas that are predicted as suitable, we propose assisted migration and repopulation approaches, for enhancing populations’ size and saving maximum genetic variability. The result of our research emphasizes once again the benefits of multidisciplinary approach in the modern biodiversity conservation.

One of the greatest challenges faced by humanity is the mitigation of rapid biodiversity loss, associated with negative anthropogenic activities 1 . Indigenous crayfish species are among the most threatened animal taxa in European freshwaters where they are experiencing substantial population declines across their entire distribution ranges 2,3 . Thus, the requirement for appropriate conservation actions and policies are urgently needed 4 .
The noble crayfish, Astacus astacus, is a keystone species of high ecological, economic, and cultural importance in Europe 5 . It is an indigenous European freshwater species whose gene pool and wide current distribution have been shaped by geo-climatic events (i.e. Pleistocene glaciations) and anthropogenic impacts (i.e. translocations, pollution, habitat degradation). In Croatia, A. astacus is recorded in all three biogeographical regions: Continental, Alpine and Mediterranean. It is naturally distributed in the waterbodies of the Black Sea drainage, with a few recorded populations in the Adriatic Sea drainage that are of anthropogenic origin 6 . Largescale genetic analyses revealed that A. astacus encompasses several mitochondrial lineages that have separated and diversified during the Pleistocene glaciations in the western and southern Balkans 7-9 , as well as in the lower Danube basin 7 . Results of microsatellite analysis revealed a differentiation of northern European populations from central European populations, with the former exhibiting a lower genetic diversity 10 . Furthermore, Schrimpf et al. 7,11 and Laggis et al. 8 revealed that this species harbours the highest genetic diversity in south-eastern Europe, while, low genetic diversity was detected in central and northern Europe, resulting from founder effects due to postglacial re-colonization and frequent human translocations for economic reasons 7 . Climate change is also www.nature.com/scientificreports/ (1) To reveal genetic diversity and population structure of A. astacus from 17 localities in Croatia (Table 1), using mitochondrial DNA (mtDNA) and nuclear DNA (microsatellite) markers; (2) To assess potential suitable habitats for the current and future period under different climate change scenarios for endangered A. astacus as well as for two NICS (P. leniusculus and F. limosus), and to identify areas of their potential current and future distribution overlap in Croatia using SDM; (3) To combine genetic data from A. astacus with its potential future distribution areas, as well as with future potential distribution of both NICS in order to identify populations and areas of the highest conservation value and priority for protection.
We expect that a combination of SDM and genetic data will provide the required information needed to develop conservation programs for endangered A. astacus. Genetic characterisation will help identifying populations that should be given the highest priority in conservation, and which can also serve as suitable donor populations for possible repopulation and reintroduction programs not only in Croatia, but also in other European countries. Furthermore, we will be able to define areas and habitats that will be under the greatest pressure from NICS and climate change, as well as potential ark sites for this species long-time survival.

Results
Phylogenetic assignment of studied populations using mtDNA sequencing. Intraspecific phylogenetic relationships and haplotype relatedness within A. astacus were described by the Median-joining (MJ) networks (Fig. 1). Reconstruction based on concatenated mtDNA data indicated the existence of six previously reported genetic lineages undefined sensu Schrimpf et al. 7 and Laggis et al. 8 within A. astacus in Europe (Fig. 1). Both COI and 16S + COI MJ networks exhibited comparable results and based on the number of mutational steps could possibly indicate the presence of a new distinct lineage containing haplotypes from the two Croatian populations and several Slovenian populations (Lsh18/Hap51 and Lsh19/Hap61 in Fig. 1, Supplementary Table S1). Remaining novel concatenated haplotypes obtained from studied Croatian populations (Hap55-Hap60) were nested within formerly recognised mtDNA lineages. Precisely, haplotypes were recovered within two lineages, Lineages 2 and 4 sensu Schrimpf et al. 7 , with some populations harbouring crayfish with haplotypes from both lineages (populations JAN, MOT, OTU) ( Fig. 1, Supplementary Table S1). The most widespread were populations belonging to Lineage 4 sensu Schrimpf et al. 7 encompassing the whole A. astacus distribution range in Croatia, while Lineage 2 sensu Schrimpf et al. 7 was found only in a few populations (Supplementary Fig. S1 and Supplementary Table S1).
Population genetics. Genetic diversity. The final data set for the microsatellite analyses comprised 413 samples and 15 microsatellite loci; 269 successfully genotyped noble crayfish samples from 12 populations in this study, and reference data from five populations obtained in Gross et al. 9 . No evidence of linkage disequilibrium between pairs of loci tested over all populations was detected after Bonferroni correction (p = 0.0004), hence all   www.nature.com/scientificreports/ Population genetic structure was detected by the Bayesian clustering analysis implemented in the software STRU CTU RE. The Bayesian Assignment Test was applied in order to assign individuals into clusters. The Evanno method, as implemented in STRU CTU RE HARVESTER, revealed that the optimal number of clusters was two (ΔK = 2). Individuals were assigned to a certain cluster if their assignment probability was ≥ 0.8, where individuals with membership to a cluster below this threshold were considered to be admixed. Most individuals showed a high assignment to one genetic cluster. The cluster I included individuals from populations MAK, TOT, JAN and BED, whereas the cluster II comprised crayfish from populations MOT, BRE, BUR, ILO, BIJ, GLO, KIK, SLO, KUT, VEL, JAR, VUK and OTU (Fig. 2). In the populations OTU and BED evidence of admixture was observed in some individuals (Fig. 2). In addition, with the purpose of getting finer insight into genetic structure of A. astacus, we report the second most probable number of distinct genetic clusters, ΔK Structure in the distribution of genetic variation was also depicted by the principal coordinates analysis (PCoA) (Fig. 2), where the PCo1 axis accounted for 25.94%, while the PCo2 axis accounted for 16.84% of the variation in the data. The PCoA revealed the existence of two well separated distinct clusters, with indication of another cluster between them. These results were congruent with the results of STRU CTU RE (Fig. 2).
In order to reveal partitioning of genetic variance by AMOVA, populations were grouped according to their affiliation to the genetic clusters inferred by the Bayesian clustering analysis (Fig. 2). The results of the hierarchical genetic diversity analysis by AMOVA revealed that most of the genetic variation was represented among crayfish within populations (66.67% of variance) followed by variation among populations within clusters (27.01% of variance), while there was less variation between genetic clusters (6.31% of variance) (Supplementary Table S3).

Species distribution models (SDMs). Model performances.
We evaluated model performance using area under the receiver operating characteristic curve (AUC) 36 . All SDMs for all species had excellent performance following interpretations for AUC values given in the literature 37,38 , with AUC > 0.9, regardless of the method used (Supplementary Table S4). The current ensemble model for A. astacus had an AUC value of 0.998, while for the NICS (P. leniusculus and F. limosus) AUC values were 0.999 for both species.
Current and future Habitat suitability. Based on model projections under current environmental conditions, A. astacus habitat suitability values (ranging from 0, indicating areas of no or low suitability, to 1 indicating areas of the highest suitability) largely corresponded to current known distribution of this species in Croatia (Fig. 3a). Largest continuous suitable habitat for this species was projected into Continental Croatia, along and between the Drava and Sava Rivers, and along the Kupa River towards the south into Alpine Croatia, while smaller and more isolated areas of suitable habitat were predicted in the area of Mediterranean Croatia, where this species is not indigenous.
Current projections for NICS revealed highly suitable habitats for F. limosus in the easternmost part of Croatia corresponding to the regions along the Danube River and lower parts of the Sava River, and the small areas of suitable habitat were predicted along the middle part of the Sava River that could enable this species spreading towards the west of Croatia (Fig. 3b). For P. leniusculus, suitable habitats under current conditions were predicted www.nature.com/scientificreports/ in the Continental Croatia. The suitable habitats were anticipated along and between the Sava and Drava Rivers, as well as along the Kupa River, overlapping with habitats suitable for A. astacus (Fig. 3c). Main trends in projected future habitat suitability under two considered RCP scenarios were similar for all species; therefore, we only report and show results for the more extreme RCP 8.5 pessimistic scenario (Fig. 3), while results for mid-range RCP 4.5 scenario are in the Supplement ( Supplementary Fig. S3).
Future projections for A. astacus suggest considerable negative impact of climate change on habitat suitability of this endangered species in Croatia (Fig. 3d). In particular, future climate change projections forecasted severe reduction in suitable habitat by 2070 in the easternmost parts of the distribution in Croatia (along and between the Sava and Drava Rivers) and to some (lesser) extent in the western part along the Kupa River towards the Alpine Croatia. In addition, future maximum habitat suitability values did not exceed 0.62, compared to current maximum of 0.98. Overall, potential future distribution of A. astacus was predicted to shift towards north-west with some gain of suitable habitat predicted in the area of Slovenia (Figs. 3 and 4). Ensemble model projections suggested that 87% of the current suitable habitat will be lost by 2070 under pessimistic RCP 8.5 scenario and only 13% will remain suitable (Fig. 4). Under mid-range RCP 4.5 scenario 65% of the current suitable habitat is projected to be lost and 35% remains stable.
Although the projected future suitable areas for NICS were wider compared to current ones, we found a severe decrease in habitat suitability values for both NICS under future climate predictions which were the  Table 1 Fig. 3e,f). In most global circulation model (GCM) projections, maximum habitat suitability values were below the threshold maximizing the sum of sensitivity and specificity. Consequently, binary maps did not provide any suitable areas for NICS in the future, regardless of the RCP scenario. We therefore show and interpret only continuous future habitat suitability projections for NICS. Under future environmental conditions F. limosus is predicted to gain suitable habitats towards the west from its current distribution, along the Sava and Drava Rivers, although with very low probability, while suitable habitats for P. leniusculus are predicted to remain relatively similar to current ones, however with lower probability (Fig. 3e,f). Under current conditions we found an overlap between suitable habitats for A. astacus and P. leniusculus in the north-west Croatia along the Drava River and southern tributaries of the Sava River, which seem to be suitable for both species (Fig. 4a). Contrary, no overlap was detected between current suitable areas for A. astacus and F. limosus (Fig. 4a).
Overlapping genetic variation of A. astacus with projected changes between its current and future habitat suitability indicated that majority of the areas harbouring highly diverse A. astacus populations are expected to A. astacus    www.nature.com/scientificreports/ become unsuitable in the future (Figs. 3d and 4b). The only populations (sampled for microsatellites) remaining in areas predicted to remain suitable in 2070 under RCP 8.5 scenario are MOT, BUR, BED, JAR and TOT (Fig. 4b). Considering mtDNA, overlap of COI haplotypes and future habitat suitability indicated that 50% (8/16) of COI haplotypes recorded in Croatia may be lost (namely, Lsh1, Lsh4, Lsh12, Lsh15, Lsh18, Lsh22, Lsh23 and Lsh24) (Fig. 4c). Majority of those haplotypes are distributed only in the eastern part of the Continental Croatia (Fig. 4c).

Discussion
In our research, the fine scale phylogenetics, population genetics and species distribution modelling were used to explore genetic diversity and structure of A. astacus, as well as the impact of climate changes and invasive species on its populations. Analyses of genetic data coupled with species distribution models revealed the vulnerability of this keystone species to climate change.
The phylogenetic network based on mtDNA displayed intraspecific relationships within A. astacus consistent with the findings of previous studies [7][8][9] . Our results confirmed the existence of several genetic lineages, with the indication of novel divergent lineage containing haplotypes from Croatia and Slovenia (Hap51/Lsh18 and Hap61/ Lsh19). Phylogenetic analysis indicates that all Croatian haplotypes belong to two mtDNA lineages (Lineages 2 and 4 sensu Schrimpf et al. 7 ) that were also recorded in different countries across Europe. Astacus astacus exhibits lower mtDNA diversity and lower genetic structuring, without an obvious geographical pattern 7 , compared to other native European crayfish species 17,[39][40][41][42] . Precisely, the MJ network showed weak phylogeographic structure and high haplotype-sharing even between geographically distant populations. This finding is consistent with the results of previous studies 9,11 showing that the contemporary distribution and genetic structure of A. astacus were shaped through past geo-climatic events, strong anthropogenic influence on its habitat and frequent human mediated translocations that partly eroded their genetic structure. Similarities between distant A. astacus populations in several cases were explained by artificial stockings from different countries or populations 10,11 . Such a case was also observed in our study; crayfish from population JAR were used for aquaculture in the geographically distant hatchery Otočac, and consequently samples from both populations belonged to the same mtDNA lineage and shared the same haplotype (Lsh2). The genetic lineages of A. astacus diversified during the late Pliocene and throughout the Pleistocene, within the period between 1.7 and 0.5 mya 9 . Current A. astacus lineage distribution shows a divergence pattern congruent with the phenomena of insularity and isolation of multiple southern glacial refugia during repeated climatic pulses in the Pleistocene that produced a mosaic of lineages 43 .
Population genetic analyses on A. astacus across the sampled localities revealed high within-population genetic diversity and moderate differentiation among populations, that differed from the results of previous studies using the same 9 or different microsatellite loci 7,8,10,11,24 . Overall, we detected a high number of alleles, proportion of polymorphic loci (P), allelic richness (A R ) and observed heterozygosity (H O ) in the study area. Genetic diversity, expressed as the P, A R and H O was higher in populations ILO, BED, KUT, VEL, while the level of genetic diversity was lower in populations MAK, VUK, BUR. Reduced genetic diversity in the populations MAK and BUR could be explained by the fact that they represent introduced populations 44,45 . Overall genetic diversity across the sampled localities of A. astacus was high when compared to the results of Gross et al. 10 , Schrimpf et al. 7,11 , Laggis et al. 8 and Panicz et al. 24 that used different set of microsatellite loci. A considerable number of private alleles was found in the majority of populations suggesting the presence of the unique genetic variation. Besides, private alleles are considered important in the long-term response to selection and the survival of populations and species 46 . We found that two populations (SLO and VEL), with significant homozygote excess, are vulnerable to inbreeding which may reduce the populations' genetic diversity, and consequently lead to the loss of adaptive evolutionary potential of the species 47 . Furthermore, we analysed whether the recent bottleneck events influenced the observed genetic structure of the studied populations, and found that two populations did experience a recent bottleneck event (JAR and TOT). Bottlenecks in small remnant populations with limited gene flow could lead to low effective population sizes and cause fitness reductions across at least part of the species distribution. In Croatia, as elsewhere in Europe, A. astacus populations are mostly isolated by natural (i.e., watershed boundaries) or artificial (i.e., anthropologically influenced) barriers, and their distribution being frequently limited to small fragmented areas (geographical regions). Therefore, there is a reasonable concern that they may undergo significant declines in effective population size and that much of their genetic diversity might be lost.
The results of STRU CTU RE and PCoA indicated the presence of genetic structuring among A. astacus populations in Croatia by identifying two main genetic clusters. Moreover, they revealed the presence of admixed individuals/populations assigned to different genetic cluster reflecting contributions of different ancestral groups or artificial translocation. Furthermore, results indicated to some extant populations' structuring according to different river basins what is similar to that found by Gross et al. 9 . The pairwise F ST values and AMOVA indicated moderate to high levels of genetic differentiation among studied populations demonstrating isolated populations with limited gene flow. The most genetically differentiated populations were MAK, TOT, JAN and VUK when compared to other studied populations. Contrary to other populations, the Vuka River (population VUK) flows directly into the Danube River which may explain the high F ST . Lower F ST values obtained in this study reflected well geographical proximity, with the exception of populations OTU and KUT. The native range of A. astacus is restricted to the rivers of the Black Sea basin, whereas population OTU belongs to the Adriatic Sea basin. Thus, low F ST value obtained for this population pair could indicate anthropogenic translocation between those two populations. Likewise, high values of F ST for MAK populations could also be explained by artificial stockings from an unknown source. Moreover, it should be pointed out that MAK is recorded in an urban lake in Zagreb City, and that it may have been introduced from the Sava River where Karaman 48 recorded A. astacus. Therefore, it is possible that this population represents a remnant astacofauna formerly present in the Sava River, with unique genetics that no longer exists elsewhere. This study discovered a higher value of global F ST (= 0.319) compared

in
A. astacus populations in Greece and across the Balkan Peninsula, respectively. A pattern of isolated populations of freshwater species that contain high genetic diversity is characteristic for the Balkan Peninsula, that is recognised as one of the freshwater biodiversity hotspots 43,49,50 . Currently this region is characterised by fragmented and complex habitats with frequently no suitable surface water connections. Therefore, restricted dispersal and gene flow among populations probably led to genetic isolation of numerous freshwater species in this area, including crayfish. However, limited gene flow may lead to reduced effective population sizes, lower genetic diversity and increase the risk of local extinction, resulting in cascading effects through freshwater ecosystems 51 . Moreover, geographically isolated populations with low dispersal capabilities such as crayfish could experience problems in accommodating to ongoing climate changes due to limited possibilities for migration and a shift in their distribution toward more climate-suitable habitats. Sensitivity to climate change in freshwater taxa was proved to be higher than in terrestrial taxa 52 , and vulnerability of freshwater crayfish to climate change, as well as to NICS has been demonstrated in many studies 3,30,33 . To evaluate the impact of climate change and NICS on the endangered A. astacus, we performed SDM. The models were able to capture the known ranges of A. astacus and two NICS in Croatia. Our predictions are concordant with previous studies of A. astacus distribution in Croatia 6 ; majority of areas currently suitable for A. astacus are located in the area of Continental Croatia, including parts of Alpine Croatia and small isolated areas in Mediterranean Croatia where species was introduced 6 . Likewise, current projections for NICS, F. limosus and P. leniusculus, revealed highly suitable habitats corresponding to their present distribution in Croatia, but also encompassing areas for their potential spread. Our current projections suggested overlap between suitable habitats for A. astacus and P. leniusculus, a competitor which negatively affects A. astacus populations in the rivers of the continental part of Croatia through competitive exclusion and A. astaci transmission 5,15 . On the contrary, modelling the current potential distribution of F. limosus in Croatia did not detect any overlap between current suitable areas with A. astacus, as expected, since waterbodies of eastern Croatia are inhabited by P. leptodactylus 6 .
Overall, our future projections demonstrated that climate change may have major negative effects on the distribution of A. astacus by reducing the surface of climate-suitable areas available for this native European species. This result is in line with findings for other endangered aquatic species in Europe 54 . Change in thermal and precipitation regimes caused by global warming will probably lead to drastic range contractions of A. astacus. Consequently, this could drive population declines across the species distribution range in Croatia. This conclusion supports the alarming studies of Capinha et al. 34 and Hossain et al. 52 that predicted extreme loss of habitat suitability for freshwater crayfish due to climate change. Thus, our results indicate that climate change-driven habitat loss represents a greater threat to A. astacus than the potential future distribution of the two studied NICS. A similar scenario was found for Austropotamobius pallipes in relation to the invasive P. leniusculus ( 30,53,54 , see below).
Future SDM projections suggested that the suitable habitat for A. astacus will likely shift towards the northwest and practically disappear from the easternmost parts of Croatia due to the severe reduction (87% of currently suitable habitats) in habitat suitability by 2070. Furthermore, the most suitable areas for A. astacus in the future were forecasted to be in the western Croatian waterbodies, some of which are at high altitudes where Austropotamobius torrentium is currently recorded 6 . Even though the Alpine region and its freshwater ecosystems represent suitable habitats for most crayfish species 53 , these two indigenous species might compete for habitats and resources 29 . To overcome this, potential ark sites for A. astacus should be placed in the rivers and artificial lakes at lower altitudes in the Alpine region, as well as within gravel pits and oxbows alongside the Drava and Sava Rivers in the north-western part of Continental Croatia that were predicted as suitable in the future. These lower altitude waterbodies of the Alpine region would provide suitable ark sites for the crayfish from mtDNA Lineage 2 and/or Genetic cluster II, while crayfish from mtDNA Lineage 4 and/or Genetic cluster I could find refugia within gravel pits along the Drava and Sava Rivers in the north-western part of Continental Croatia. Keeping in mind that those suitable habitats are inaccessible to A. astacus due to natural dispersal barriers, human interventions would be needed. Assisted migration (AM) as an adaptation strategy for mitigating the projected effects of climate change on species is widely proposed 20,21 , especially for those with a life history features that prevents them from migrating to suitable habitats. However, it is a controversial topic among conservation biologists, with numerous identified risks. Arguments against AM include: risk of translocated species becoming invasive with associated negative biological, ecosystem and socioeconomic effects; spread of diseases and pathogens that can be transferred into new host species; removing individuals from existing populations increases the extinction risks facing those source populations 19,55 . In order to overcome those arguments, careful planning encompassing risk assessments, cost-benefit analyses, conducting AMs on a small scale, with robust monitoring that would enable prompt corrective actions to be taken if needed, along with political and public promotion, could insure successful AM implementations 21 .
Regardless of the RCP scenarios, our binary projections did not forecast any suitable areas for NICS in the future which should be interpreted with caution due to (a) the small number of available occurrences for NICS in the Croatian waterbodies used for SDMs; (b) models that do not account for human-mediated dispersal of NICS 56 ; (c) the naturalised climatic niches of NICS that can differ from their natives' climatic niches 28 ; and (d) underestimated potential range expansion in the future due to the known issue of non-equilibrium of NICS with the environment within the invaded range 57 . Nevertheless, our continuous future habitat suitability projections showed that, even though the projected future suitable areas for NICS were more extensive than the current ones, drastic decrease in habitat suitability values for both NICS were displayed under future climate predictions. Explicitly, potential areas where A. astacus would overlap and compete with NICS virtually disappeared by 2070 under both climate change scenarios of high-warming (RCP 8.5) and low-warming conditions (RCP 4.5). This result is consistent with the results of Préau et al. 30  www.nature.com/scientificreports/ and P. leniusculus in France based on SDMs, despite substantial ecological niche overlap between the two species. Likewise, Gallardo & Aldridge 54 found that both endangered A. pallipes and invasive P. leniusculus were predicted to be negatively affected by climate changes in Europe. However, the range contraction was predicted to be more dramatic for the invasive P. leniusculus, leading to decreased overlap and consequently competition between the two species in the future, particularly in our study area. A more recent study by Zhang et al. 58 confirmed that invasive P. leniusculus may lose a substantial portion of suitable habitat in Europe by 2070 in response to climate change. Moreover, Capinha et al. 34 studied the potential distribution of indigenous crayfish species and NICS in Europe and found that climate-suitable areas were predicted to decrease by nearly 70% for A. astacus, 42% for P. leniusculus, and 49% for F. limosus by 2080. However, their models predicted that overlap of suitable ranges for native European crayfishes and invasive crayfishes would increase in the future which is contrary to our results. This may be because south-eastern Europe seems to be less suitable for P. leniusculus under changing climatic conditions 54,58 . It is therefore crucial to continue the monitoring of NICS invasions in the future. Estimated reduction in habitat suitability by the end of this century indicates potential loss of a significant portion of the A. astacus genetic variability, especially in the eastern part of Continental Croatia that may potentially lose populations with high and unique genetic diversity. Minimising such possible losses in the future requires viable A. astacus populations to be established and maintained in ark sites/climate change refugia. Our results exposed an alarming need to prioritise conservation planning and management that will support existing populations and potentially establish new ones in the areas of stable habitat suitability that are expected to sustain A. astacus into the future. Species responses to climate change will depend on their distribution shifts to accommodate climate changes, and/or rely on the adaptation based on the standing genetic variation. Keeping in mind low dispersal abilities and isolated populations, we argue that assisted migration and population mixing approaches will be probably needed in the future to enhance the size and genetic diversity of remnant populations in order to maintain the long-term survival of the species 30,34,59 . Based on our results, we propose several donor populations for future restocking and reintroduction strategies. Namely, populations ILO, KUT, VEL, BAČ and BIJ, contain high and unique genetic diversity both at the mitochondrial and nuclear level, but they are predicted to be lost due to unsuitable habitats in the future. Dispersal as a fundamental behavioural mechanism is of great importance for adaptation and species' responses to rapidly changing climate 26 . Strong dispersal limitations, habitat discontinuities and limited gene flow have a major effect on the ability of crayfish populations to withstand climate changes. Thus, assisted migration in climate change refugia seems a logical solution for slowing down genetic diversity erosion, reducing genetic load and the detrimental consequences of inbreeding, but also allows variations in allele frequencies [60][61][62] .
The adoption of such approaches for conservation purposes has gained significant momentum over the last few decades; reintroduction of the A. astacus into restored waterbodies has become common practice, even though the genetic origin of stocking material has rarely been considered 10 . Thus, potential ark sites should represent areas that maintain the highest contemporary genetic diversity in the species and predicted climatesuitable habitats for the future. Astacus astacus relocation should be preceded with a careful assessment regarding potential negative consequences of assisted gene flow that can impact the success of relocated populations, particularly when populations exhibit local adaptation to factors other than climate 63,64 . Also, introgression between local and translocated populations could result in outbreeding depression 60,61 . Still, Bláha et al. 23 found no significant decline in genetic diversity between the source and translocated A. astacus populations after introduction. Furthermore, their study showed that even though the source populations did not possess high genetic diversity, their distinctiveness still made them suitable for conservation purposes. In addition, it is critical that climatically suitable sites outside A. astacus historical range for conservation purposes should be free from diseases, such as crayfish plague caused by oomycete Aphanomyces astaci.
In conclusion, our results suggest that securing the future of A. astacus will require significant interventions. This paper provides a baseline to guide these actions. Specifically, SDM combined with population genetics provided essential guidance for conservation actions aimed at safeguarding endangered A. astacus in Croatia by revealing genetic structure and identifying sites most suitable for protection and sites where climate change constitutes a threat. In addition, our study corroborates SDM as a valuable tool for conservation planning of threatened crayfish species by identifying areas within a species' distribution that may be vulnerable and suitable areas for assisted migration as shown in studies on European crayfish 34,52 , A. pallipes complex 3,29,30,53,65 , and A. torrentium 3 .

Material and methods
Genetic diversity and population structure. Sampling and DNA extraction. We collected A. astacus samples across its entire distribution range in Croatia (Supplementary Table S1, Fig. 5). Specimens were collected by hand or baited traps in accordance with ethical standards and with permissions of local authorities. One pereiopod from each individual was sampled and stored in 96% ethanol at 4 °C. Genomic DNA was extracted from the pereiopod's muscle tissue using GenElute Mammalian Genomic DNA Miniprep kit (Sigma-Aldrich, St. Louis, MO) following the manufacturer's protocol and stored at −20 °C.
Phylogenetic assignment of studied populations using mtDNA. Samples used for phylogenetic reconstruction are reported in Supplementary Table S1. Mitochondrial 16S and COI genes were amplified and sequenced with universal primers 16Sar/16Sbr 66 and LCO-1490/HCO-2198 67 allowing comparison with previously published A. astacus sequences [7][8][9]68 . Polymerase chain reactions (PCR), purification and sequencing were performed according to Gross et al. 9 . Sequences were edited and aligned in Bioedit v. 7.2.5 69 . The final COI alignment did not contain any length variants or ambiguous sites and its final length was 623 bp, while 16S alignment contained 1 length variation and its final length was 475 bp. In order to perform phylogenetic analysis comparable with  Table S1). Median-joining (MJ) network approach 71 was used to visualise intraspecific relationships among haplotypes within A. astacus using PopArt 72 . Since A. astacus is characterised by low diversification in mtDNA 9 and within-species data sets have fewer characters for phylogenetic analysis which diminish the statistical power of traditional phylogenetic methods 73 , we used phylogenetic networks that are better suited for description of intraspecific evolutionary relationships. Two MJ networks were reconstructed in order to determine nonhierarchical phylogenetic relationships between A. astacus haplotypes. Median-joining network I comprised COI sequences (623 bp long) from Croatian populations obtained in this study and in the study by Gross et al. 9 . Median joining network II included concatenated 16S + COI sequences (825 bp long) obtained in this study and assembled with all available sequences at European level [7][8][9]11,68 . This approach enabled us to associate haplotypes obtained in the present study to the haplotypes obtained in previous research and indirectly to the lineages sensu Schrimpf et al. 7 and groups sensu Laggis et al. 8 .
Population genetics of studied populations using microsatellites. For microsatellite analyses we amplified 19 species-specific tetranucleotide repeat microsatellite loci developed by Gross et al. 74 , and following modified protocols and procedures as in Gross et al. 9 . Microsatellite loci were genotyped on Applied Biosystems 3500 XL Genetic Analyser (Life Technologies, USA) using internal GeneScan 600 LIZ Size Standard v2.0 (Life Technologies, USA). Genotypes were scored using GeneMapper v.5 software (Life Technologies, USA), and were double-checked manually by two experts. Since four loci had overlapping allele size ranges (Aast4_26, Aast4_47, Aast4_10 and Aast4_30) they were omitted from further data analyses which were performed using 15 microsatellite loci. Also, several samples from different populations with more than two non-amplified loci were omitted from further analysis. Microsatellite loci were tested for potential presence of genotyping errors due to null alleles, stutter peaks and large allele dropout using MICRO-CHECKER v.2.2.3 75 . Pairwise linkage disequilibrium between all pairs of loci was tested using Fisher's exact test in GENEPOP v. 4.7.2 76 . Null allele frequencies based on the expectation-maximization (EM) algorithm 77 and corrected F ST values using the ENA method were estimated using FreeNA 35 with a number of bootstrap replicates fixed to 10,000. The estimations of F ST , with and without null allele correction, were compared for each population using t-test in STATISTICA 13 (StatSoft. Inc).
Population genetics analyses were conducted with the microsatellite genotype data of 12 A. astacus populations obtained in this study that were supplemented with the microsatellite genotype data of five Croatian A.   76 . All probability tests were based on the Markov chain algorithm using 10,000 dememorization steps, 100 batches and 5000 iterations per batch. Significance levels were adjusted applying the Bonferroni correction to correct for the effect of multiple tests. Recent reductions in the effective population size using allele frequency data and potential signatures of recent bottlenecks were tested using the heterozygosity excess method implemented in BOTTLENECK v.1.2.02 81 under three different mutational models: infinite allele model (IAM), stepwise mutation model (SMM) and two-phase model (TPM). Significant deviations from mutational-drift equilibrium were tested using the Wilcoxon sign rank test with 10,000 simulations.
Population genetic differentiation and structure. Genetic differentiation between all population pairs was estimated through pairwise F ST values using FSTAT v.2.9.4 79 . Genetic structure among studied populations and assembling of individuals into groups (genetic clusters) was assessed using the Bayesian model-based clustering approach implemented in STRU CTU RE v.2.3.4. 82 . The conditions performed were 10 runs for each genetic cluster (K) between 1 and 17 using a 100,000 burn-in period followed by 100,000 MCMC iterations, under the admixture model, with correlated allelic frequencies. The number of optimal K was inferred using the protocol defined by Evanno et al. 83 as implemented in STRU CTU RE HARVESTER v. 0.6.93 84 . STRU CTU RE graphical results were plotted with CLUMPAK 85 . In addition, structure in the distribution of genetic variation was visualized by principal coordinates analysis (PCoA) using Nei's genetic distance in GenAlEx v. 6.51. Hierarchical analysis of molecular variance (AMOVA) was carried out using ARLEQUIN v. 3.5.1.2 86 in order to estimate partitioning of genetic variance among groups, among populations within groups and within population. Populations were grouped according to their affiliation to the genetic clusters inferred from STRU CTU RE; The cluster I included individuals from populations MAK, TOT, JAN and BED, and the cluster II comprised crayfish from populations MOT, BRE, BUR, ILO, BIJ, GLO, KIK, SLO, KUT, VEL, JAR, VUK and OTU (Fig. 2).

Species distribution models (SDMs). Species occurrence data.
We compiled all known presence-only occurrences of A. astacus and the two NICS (P. leniusculus and F. limosus) from across Croatia from our own published and unpublished field sampling 6 . This resulted in a total of 174 occurrences for A. astacus, 22 for F. limosus and 17 for P. leniusculus (Fig. 3).
Environmental data. We initially considered 22 environmental variables from various sources and databases describing climate, topography and forest cover of the study area ( Table 3). The 19 bioclimatic variables were obtained from the WorldClim 1.4 database 87 , altitude and slope were derived from a digital elevation model from the NASA Shuttle Radar Topography Mission (SRTM) elevation data (https:// www2. jpl. nasa. gov/ srtm), while the variable percentage of forest cover in 1 km 2 was calculated from the Corine Land Cover 2018 dataset (https:// land. coper nicus. eu/ pan-europ ean/ corine-land-cover). All environmental variables were used at a spatial resolution of ~ 1 km 2 . Predictor variables for SDMs of A. astacus and the two NICS were then selected based www.nature.com/scientificreports/ on our expert knowledge about their ecological relevance for the target species (potentially influencing species' physiology and life history), excluding highly correlated ones based on variance inflation factor, VIF < 10 (usdm R package; 88 ). Thus, the final predictor set for A. astacus included ten, and for NICS nine environmental variables (see Table 3).
Modelling procedure. To assess the potential current and future habitat suitability of A. astacus and two NICS (P. leniusculus and F. limosus), we developed SDMs using an ensemble approach implemented in R package BIOMOD2 ver. 3.3-7 89,90 . For each species we applied three different modelling methods (Random Forest-RF, Generalized Boosted Model-GBM and Maximum Entropy-Maxent) with ten replicates for each method (a total of 30 models for each species). Occurrences were combined with 10,000 random pseudo-absences drawn across the study area for methods that require absences 91 . In each run, 70% of the occurrences were used for model calibration, and the remaining 30% were used for model evaluation using AUC 36 .
To build the current ensemble model we used only highly reliable models with AUC > 0.9 37 and obtained this ensemble as an AUC weighted average. The obtained current ensemble model was then projected under both current and future environmental conditions to obtain potential habitat suitability maps for each species. For future projections we used two RCP scenarios (mid-range emission scenario RCP 4.5 and pessimistic scenario RCP 8.5) and four global circulation models (GCMs) suitable for Europe 92 (CCSM4, MIROC5, MPI-ESM-LR and HadGEM2-CC) for the 2070-time period (average for 2061-2080). Since future projections for variables slope, altitude and forest cover were not available, we kept them as constant in our future projections, assuming that they will not change for our study area during the considered time period. The available data on forest cover change in Croatia during the last decades and current forest management structure and practices provide confidence that at least for our study area, forest cover may remain stable in the future 93 [https:// forest. eea. europa. eu]. An ensemble future projection for each RCP scenario was obtained by taking an average of the four GCM projections. Multiple RCPs and GCMs were used to address the associated uncertainties arising from different climate change predictions 94 .
To obtain binary presence/absence maps helpful in model interpretation and for calculating changes in habitat suitability for A. astacus, we applied a threshold maximizing the sum of sensitivity and specificity 95 to ensemble current and future continuous habitat suitability maps.
Finally, to estimate the effects of climate change on genetic diversity and structure of the focal species, we overlapped genetic data of A. astacus with its potential current and future suitable areas, as well as with future potential distribution of both NICS.