Uncovering glacial footprints and identifying phylogeographic units in the freshwater crab Potamon elbursi Pretzmann, 1962 (Decapoda: Potamidae) based on mitochondrial data

) Uncovering glacial footprints and identifying phylogeographic units in the freshwater crab Potamon elbursi Pretzmann, 1962 (Decapoda


Introduction
Population structure is influenced by a combination of present environmental processes and past evolutionary history of species (Hewitt 1999).Studies on the genetic structure of natural populations will enhance the understanding of how microevolutionary forces have interacted throughout the history of species and can help to reconstruct important past evolutionary processes (Avise 2000).In the Quaternary period, the earth has been covered repeatedly by ice sheets, and climatic oscillations related to glacial cycles were followed by great changes in population structures and distributions of species (Petit et al. 2003;Hewitt 2004;Kohli et al. 2021).
As a result of the climatic oscillations of the Quaternary period, aquatic ecosystems have been dramatically influenced (Dudgeon 2010;Novico et al. 2022).In particular, obligate freshwater species and land-locked organisms including freshwater crabs have been hugely affected by the habitat changes driven by such climatic fluctuations (Cumberlidge et al. 2009;Monteiro et al. 2021;Biswas and Sarkar 2022;Reed et al. 2022).The fact that freshwater crabs have low dispersal abilities and limited capability to cross land bridges, present them philopatric (Ng 1988;Baird et al. 2021;Biswas and Praveen Karanth 2021;Vecchioni et al. 2021), so that they spend their entire life in a specific freshwater drainage system.These crabs do not have larval stages, and their fertilized eggs develop directly into juveniles, involving parental care (Ng 1988;Dai 1999).Therefore, freshwater crabs are appropriate models for the reconstruction of evolutionary processes in freshwater ecosystems, demography of populations, estimating the impact of human activities on freshwater habitats, and the effects of Quaternary ice ages on genetic diversity (Rodrıguez and Ng 1995;Shih et al. 2006;Daniels et al. 2006b;Keikhosravi and Schubart 2014a;Keikhosravi et al. 2015).Potamon elbursi Pretzmann (1976), is a freshwater crab distributed in northern Iran, specifically in two different drainages; the Caspian Sea and Namak Lake (excluding the Caspian Sea coastlines).This species was revalidated and re-described recently by Keikhosravi and Schubart (2014b).
The Caspian Sea coastal region is the mistiest area of Iran; the climate is warm temperate, with hot and humid summers and fairly warm and rainy winters.Otherwise, the south face of the Alborz Mountains (i.e., the Kavir Desert and Namak Lake) experience aridity, with quite cold winters and hot summers.Although the Kavir Desert is an endorheic basin, the Namak Lake is fed by rivers flowing down from the Alborz and Zagros mountains draining to the Caspian Sea (Farajzadeh and Matzarakis 2009).Recent studies showed that populations of different taxa in the northern region of the Alborz Mountains are highly influenced by climatic fluctuations associated with glaciations (Naderi et al. 2014;Parvizi et al. 2018).Recently, Keikhosravi et al. (2015) described the small-scale population genetic structure of P. elbursi (in downstream of Ghezelozan and upstream of Namak Lake drainages) and they found the lack of genetic partitioning between drainages and only limited restriction of gene flow among populations.Additionally, they suggested a population expansion model beginning about 1 Mya in a Pleistocene context.In contrast, a later geometric-morphometric study by Kalate et al. (2017) revealed a significant difference between the populations of the two drainage systems.
In light of these previous and diverging results, it appears necessary to conduct a more detailed phylogeographic study of P. elbursi, covering its entire distribution range and using a dataset composed of both newly generated and previously investigated sequences.The genetic diversity of P. elbursi populations should be determined for both drainages of its occurrence, the Caspian Sea and Namak Lake, for a better understanding of the degree of genetic differentiation among populations.It is also aimed to investigate the effects of the Pleistocene climatic fluctuations on P. elbursi in order to illustrate the demographic history of these populations along the southern face of the Alborz Mountains.

Sampling strategy
This study covers two drainage systems, the southern Caspian and Namak Lake basins.A total number of 70 specimens from seven stations were collected during fieldwork from 2011 to 2013 (aiming to gather 10 specimens from each station).One walking leg from each individual was removed before the crab was returned to the river.Each sample was preserved in absolute ethanol in a labeled tube and kept on ice during transport to the lab for later molecular study (detailed sampling sites and sample sizes are given in Table 1 and Fig. 1).

DNA extraction, amplification and sequencing
Total genomic DNA was isolated from 0.5 to 1 g of muscle tissue of walking leg by using the Puregene method (Gentra Systems).DNA of almost the entire mitochondrial gene Cox1 (cytochrome oxidase subunit I, ~1500 basepairs) was amplified with polymerase chain reaction (PCR) and the primer combination LCO1490 (5′GGTCA ACAAATCATAAAGATATTGG3′) ( Folmer et al. 1994) and COH16 (5′CATYWTTCTGCCATT TTAG A3′) (Schu bart 2009).The marker was chosen because in previous studies (i.e., Keikhosravi and Schubart 2014b;Keikhosravi et al. 2015) the same region was used, and the combination of datasets allows more robust conclusions.The PCR conditions and the thermocycler profile were as follows: initial denaturation step at 94°C for 4 min followed by 36 cycles including denaturation for 45 s at 94°C, annealing for 45 s at 50°C, extension for 90 s at 72°C; a final extension at 72°C for 5 min.Products were outsourced to Macrogen Asia (South Korea) using the primer COH16.
Sequences were manually edited with Chromas Lite 2.01

Statistical analysis
Sequences were obtained from collected specimens (70 sequences belonging to 7 populations).In addition, 61 sequences (GenBank accession number: LN833869-LN833879) from the previously published paper (Keikhosravi et al. 2015) were added to our dataset from six populations.Therefore, a total of 131 sequences with readable sequence length of at least 711 bp were used for the analyses.
A statistical parsimony network analysis was characterized by TCS ver 1.21 (Clement et al. 2000), to determine the extent of genetic divergence within populations and the relationship among the haplotypes.Distribution patterns of the divergent haplotypes across the studied area were also investigated.With this analysis, sequences are alienated into a network of closely related haplotype groups with joined branches with less than 95% probability for parsimonious connection.The main indices of population genetic diversity, including the number of haplotypes (h), haplotype diversity (Hd), nucleotide diversity (π), and number of segregating sites (s) were calculated using DNASP ver 5 (Librado and Rozas 2009) and Mean Pairwise Differences (MPDs) with the use of ARLEQUIN ver 3.5 (Excoffier and Lischer 2010) were obtained for comparison of the population genetic structure.Gene diversity that also known as haplotype diversity, illustrate the possibility that two arbitrarily sampled alleles are different, though nucleotide diversity is defined as the average number of nucleotide differences per site in pairwise comparisons among DNA sequences (Nei 1987).

Population genetic structure
Genetic differentiation among populations was quantified by performing a global test of differentiation among specimens and computing pairwise Φst (Raymond and Rousset 1995).These analyses take into account the presence of indels, which show that some haplotypes are differentiated.The analysis was performed in ARLEQUIN and significance values of estimated Φst were tested using 10000 permutations.Gene flow among populations was calculated with Nm using the formula: Nm = (1 / Φst1) / 4 (Slatkin 1993).
To investigate patterns of historical population structure and identify the neighboring demographic groups with the highest genetic differentiation, an analysis of molecular variance AMOVA (Excoffier et al. 1992) was performed with ARLEQUIN.The percentage of diversity in the subclasses of different hierarchies including among groups (inter-drainages), among populations within groups (within drainages) and within populations was estimated according to geographic proximity.
In order to estimate the occurrence of isolation by distance among populations, we calculated the correlation between genetic and geographic distances among the sampling sites and mean genetic distance value for each population, therefore Mantel test was performed using Alleles in Space ver1.0 (Miller 2005) with 100,000 permutations (kilometer as distance index).Mantel index varies between +1 and −1, where the more positive and significant values of the Mantel index show the more positive correlation (direct), and the more negative and significant values indicate the more negative (more inverse) correlation between the two genetic distance matrices and the geographical distance among the populations (Miller 2005).The genetic distance between all pairs of individuals based on Kimura two parameter and the Fst values between population was calculated and arranged in form of a genetic matrix.The connection between the two matrices was assessed using Gen Alex 6.5.

Demographic history
We used two statistical tests, TAJIMA'S D (Tajima 1989) and FU'S FS (Fu 1997), in order to test for the possibility of a past population expansion event with effects on the recent demographic structure of genetic diversity.The analyses were implemented in the program DNASP and pvalues of TAJIMA'S D and FU'S FS (were 0.1 and 0.05 respectively) were generated using 1,000 simulations under a model of selective neutrality.A significantly negative deviation from zero can be taken as a result of recent expansion, significant increase in population size, directed selection or purifying selection.While a significantly positive value can result from genetic drift, decrease in population size and balancing selection effect during the evolutionary history of the population (Ramos and Rozas 2002).Mismatch distribution analysis (MMD) with 10,000 bootstrap replicates between the haplotypic pairs (Sherry and Rogers 1994;Rogers 1995) were also conducted using ARLEQUIN.Parameters of the demographic model expansion such as tau (τ), theta 0 and theta 1 (1000 permutations and α = 0.05) were obtained from ARLEQUIN.Harpending's raggedness index (Harpending et al. 1993) and the sum of squared deviations (SSD) between the observed and expected mismatches for each population were calculated using ARLEQUIN to examine demographic changes.This measure quantifies the smoothness of the observed mismatch distribution, and a nonsignificant result indicates an expanding population (Harpending 1994).A significantly high value of Harpending's raggedness index (p < 0.05) or high estimates of SSD suggests a nonfitting deviation and thus a rejection of the null hypothesis of a population expansion model.The spatial expansion hypothesis (both raggedness index and SSD) was tested using a parametric bootstrap approach (500 replicates).
We investigated mtDNA effective population size change through time using coalescent-based Bayesian skyline plot (BSP) implemented in BEAST v.2.6.3 (Drummond et al. 2005;Bouckaert et al. 2019).We used JModeltest vr.2.1.7 (Darriba et al. 2012) to choose the best model of nucleotide substitution and considered a strict clock prior with a substitution rate of 2.33% per MY for COX1 as was previously used for other Potamid crabs (Jesse et al. 2011;Parvizi et al. 2018Parvizi et al. , 2019)).We ran 3 independent MCMC analysis and set the MCMC chain length to 40 million, sampling every 4000 steps.We then combined the results of log and tree files of independent runs using LogCombiner v.2.4.7 (Rambaut and Drummond 2017) and discarded the initial 25% of generations as burn-in.Finally, we checked the convergence of all parameters and produced a BSP in Tracer v.1.6(Rambaut et al. 2018).

Species distribution modeling
For modeling analyses, records from 70 localities in Iran were analyzed.We used data from 19 bioclimatic variables (Bio01-19), which were obtained from the World-Clim database (https://www.worldclim.org).
The open modeller v. 1.0.7 was used to read the corresponding environmental values for each occurrence point from the input rasters (Munoz et al. 2009).To prevent repetition, bioclimatic variables with Pearson correlation coefficient < 0.75 were selected using SPSS v16.0.
The selected layers and distribution records were employed by MaxEnt 3.4.1 (Phillips et al. 2017) to construct the final model.Ninety percent of the data were used as training samples and 10% as test samples.The procedure was repeated fifteen times.The area under the curve (AUC) in receiver operating characteristic (ROC) was used for evaluating the model performance.The predicted models were run for the current scenarios and then projected for the past scenarios.The current and past models were evaluated using ENMTools to find the extent of niche overlapping between models.

Genetic diversity
A total of 131 sequences with readable sequence length of at least 711 bp of the Cox1 mtDNA belonging to 13 populations were aligned.The alignment contains 36 variable sites including 17 singletons and 19 parsimony informative sites.We found 25 different haplotypes with relatively high haplotype diversity (Hd = 0.868) and low nucleotide diversity (π = 0.0047), in contrast.Bijar population showed the highest (Hd = 0.76) and Taleghan population showed the lowest (Hd = 0.00) haplotype diversity.Accordingly, the least nucleotide diversity was related to Taleghan population with π = 0.00 (Table 2).

Haplotype network analysis
The statistical parsimony network revealed that the majority of haplotypes are "rare haplotypes" (20 out of 25 haplotypes).Four main haplogroups were estimated in the studied drainages.Haplotypes A, B, C and F contained samples from both drainages while haplotypes E and D included only samples belonging to the Caspian Sea drainage (Fig. 2).Sixty-one percent of haplotypes are found in the two largest haplogroups, group 1 (western populations) and group 3 (eastern populations), comprising 34 and 27 percent, respectively.Both haplogroups presented  a star shape topology, an indication of recent expansion.
Haplotype E, as the main haplotype in haplogroup 1, is shared by the populations of Sariaghol, Sharichay, Bijar, Ghamchy, Arpachay and Galerood, of which every population, excluding Bijar and Ghamchy, carries its private haplotypes giving the group a star-shaped topology.The main haplotype in haplogroup 3 is shared by all populations with rare haplotypes around it, thus presenting a starshaped topology as well.Haplogroup 1 includes all populations belonging to the Caspian drainage -except Molaali, which does not share any haplotypes with the group.All populations, except Sepidrood, that were associated with haplogroup 3 belong to the Namak Lake drainage.Sepidrood belongs to the Caspian Sea drainage, which is geographically distant from the haplogroup it partially contributed to (i.e., haplogroup 3; Namak Lake).
A similar pattern was conversely observed for Jajrood, but with a slight contribution with a geographically unrelated haplogroup.Although Taleghan and Sepidrood both belong to the Caspian basin and are geographically closely related, Speidrood slightly shared haplotypes with Taleghan.Sariaghol belongs to haplogroup 3, excluding the only individual which is shared with individuals included in haplogroup 4 (see Fig. 2H4D).Kinevars is partially shared with haplogroup 4, although it belongs to the Namak Lake drainage (see Fig. 2).Bijar also has partial contribution to haplogroup 1.Therefore, the populations of the two drainages partially share their haplotypes.Haplotype E (in haplogroup 1) takes almost a central position relative to the other haplotypes, which is an indication of it being the ancestral haplotype.

Genetic differentiation (Φst)
Mean pairwise genetic differentiation (Φst) among populations showed varying levels of genetic differentiation among populations, apart from their drainage origin (Φst ranging from 0.95 between Taleghan and Ghamchay to −0.007 between Shahrichay and Ghamchay; Table 3).
The Namak Lake drainage system populations were more genetically similar to each other compared to the populations within the Caspian Sea drainage system.The majority of the populations were characterized by Φst > 0.5, indicating substantial genetic differentiation, particularly among those populations that are distantly located from each other.Notably, Molaali and Taleghan (both belonging to the Caspian Sea drainage system) have the highest pairwise differentiation rates with other populations.Additionally, Jajrood, Kinevars and Darake, which belong to the Namak Lake drainage system, stood out as highly genetically differentiated from the rest of the populations (Φst > 0.5).

Genetic structure (AMOVA)
Three groups were defined based on geographical distance and drainage association as follows: (1) the Namak Lake Basin (Jajrood, Darake and Kinevars), (2) the Central Caspian Basin (Sepidrood, Molaali and Taleghan), and (3) the Western Caspian Basin (Arpachay, Ghamchay, Galerood, Shahrichay and Sariaghol) as well as the Southern Caspian Basin (Telvar and Bijar).In contrast to our expectations, the results showed that the percentage of variance among populations within these groups is greater than among the groups.Therefore, there is no distinct genetic structuring among the three geographic groups in P. elbursi (Table 4).However, distinct genetic structure among the three groups and between the two studied basins can be concluded, when we removed Sepidrood population from the AMOVA analysis.Since Sepidrood is not distinct from other Namak Lake populations, the AMOVA analysis was not able to show structuring of populations between the two basins and groups.

Mantel test
There was a positive and significant correlation between genetic and geographical distances (r = 0.468, p < 0.005; Fig. 3).Therefore, it can be argued that with increasing geographical distance among P. elbursi populations, the genetic distance is significantly increased.

Demographic analysis
Neutrality tests, including FU'S FS and TAJIMA'S D were applied to assess signatures of recent historical demographic events at level of mitochondrial lineages, identifying the effects of natural selection on the studied gene (Tajima 1989;Fu 1997).Although TAJIMA'S D test was not significant, FU'S FS test as a more powerful test (Ramos and Rozas 2002) was significantly negative indicating the effect of recent population expansion and a significant increase in the population size (Table 5).
Initially, the MMD analysis of total populations showed multimodal distribution, thus the populations were considered independently based on geographical distance and drainage association (see Figs 1 and 2).Multimodal distribution was also eliminated when the Caspian Sea basin was considered separately.The central and southern Caspian Sea populations showed a multimodal distribution, while the western Caspian Sea populations, along with the Namak Lake group, showed a unimodal distribution, indicating recent expansion in the effective population size (Table 6; Fig. 4).However, the hypothesis of recent population expansion cannot be ruled out for multimodal groups since the SSD and Harpending's raggedness value were not significant for all groups.In addition, these results are in line with the neutrality tests and network topology, showing star shapes with several rare haplotypes for the Western Caspian Sea and Namak Lake populations.In this analysis, the effective population size change in P. elbursi was obtained by calculating two consecutive values of the θ parameter (θ 0 = 2N 0 μ and θ 1 = 2N 1 μ), in which N 0 and N 1 indicate population size in the past and present respectively, and μ demonstrates the mutation rate in the mitochondrial sequence which

Galerood
Molaali Darake Jaj rood Taleghan is estimated to be 2.33% per million years (Schubart et al. 1998).Effective population size in the past was N 0 = 0.002 and in the present was N 1 = 2.304, which clearly indicates an increase in population size.In addition, the τ parameter (τ is the time from the beginning of the population expansion) was also calculated for the mitochondrial gene in this species by using the equation τ = 2μt.This calculation allowed us to estimate the starting point of population expansion in P. elbursi.The results showed that the beginning of population expansion in this species is estimated at 0.64 million years ago (τ = 3), which is in consistent with the MMD results indicating a recent demographic expansion of P. elbursi (Table 6; Fig. 4).

Sepidrood
The Bayesian skyline plot also showed an increase in the mtDNA effective population size of the studied samples which began at approximately 20,000 years before present (Fig. 5).

Species distribution modelling
According to the Pearson correlation coefficient and considering the habitat of the species, 9 environmental variables were chosen (Table 7).All models were run in 15 replicates for two periods of time (the current and last glacial) and the AUC values of all models were obtained with high accuracy values.Contribution performance of layers for each period and scenarios are presented in Table 7.Based on the AUC values, maps were predicted to be at a very good level (AUC > 0.950) (Fig. 6).According to these results (Table 7), Bio18 (mean diurnal range) and Bio10 (Mean temperature of warmest quarter) had the highest contribution in the current habitat suitability pre-  dictions for P. elbursi, respectively.However, the highly contributed variables for the last glacial period distribution were suggested to be Bio18 (mean diurnal range) and Bio19 (precipitation of coldest quarter).The suitability of habitats for this species changed from the last glacial period to the current time (Fig. 7).Based on the niche overlap analysis, habitat suitability among current and last glacial scenarios only overlapped about 39%.Therefore, the suitable areas for P. elbursi changed from the last glacial to the current period.In the last glacial period, this species distributed in north and northwest of Iran, but in present it is more restricted to the northwest.Comparing maps of the two periods shows that the northwestern populations are more dispersed at present than they were in the past, while the northern populations appear to have declined (Fig. 7).

Discussion
Population structure analyses of Potamon elbursi indicate a clear distinction among some local populations as well as between drainages when we excluded the Sepidrood population.The extent of separation is positively correlated to geographical distance, suggesting that regional adaptation, as a consequence of past glacial fluctuations or ancestral separation, may have an important role in this regard.In addition, the isolation of some populations, high haplotype diversity, low nucleotide diversity, and limited gene flow among some populations were observed.In line with the morphometric study by Kalate et al. (2017), high genetic differentiation can be expected among populations that can also be distinguished by geometric morphometric methods.Morphometric findings accurately showed morphological differentiation between the Caspian Sea and Namak Lake drainages, while a previous molecular study failed to show a distinction between the two basins (Keikhosravi et al. 2015).Explaining the causes of morphological differences among populations is not straightforward, but it can be argued that morphological variance is a result of interactions between environmental and genetic conditions (Daniels et al. 1998).Therefore, the apparent differences in crabs can be due to their evolutionary responses to various environmental conditions that populations are exposed to, such as temperature, available food and regional climate (Daniels et al. 1998;Torres et al. 2014;McLay 2015).Intrapopulation genetic diversity is an important parameter in determining the resilience of species against environmental stress.High levels of genetic diversity increase the ability of populations to respond to natural selection and thus the overall health of a population (Kalinowski 2005).Gene flow is one of the main drivers of genetic diversification, but higher gene flow leads to low genetic differentiation among populations.On the other hand, genetic drift can alternatively become an important cause of genetic differentiation, if Nm < 1 (Kalinowski 2005).In the present study, the amount of gene flow in  Precipitation of coldest quarter (Bio19) 5.9 5.2 P. elbursi was calculated as 0.1, indicating limited gene flow, which was also confirmed by Φst analysis among several populations.Direct development, the absence of larval stages, exclusive philopatric behavior (i.e., strategy to remain at birthplace) and parental care characterize freshwater crabs as potentially highly secluded organism with low gene flow (Cumberlidge et al. 2005;Daniels et al. 2006a, b;Yeo et al. 2007).In contrast, in marine crabs with extended larval planktonic periods, high levels of gene flow among populations are expected which can result in low genetic distinction among populations (Zhou et al. 2016).
It has been proven that climate change, particularly alternating episodes of glacial periods on earth, has been one of the major factors in the formation of the current phylogeographic patterns (Avise 2000).Recent studies (e.g., Ahmadzadeh et al. 2013;Ashrafzadeh et al. 2016;Parvizi et al. 2018) have shown significant possible impact of Pleistocene climate oscillations on some aquatic and terrestrial organisms in Iran, in areas sustained with mountain glaciers.The Alborz and Zagros mountains are among the main barriers for the geographical separation and climate variation in the Middle East (Kehl 2009), as well as the isolation and separation of the crab populations in these areas (Keikhosravi and Schubart 2014a;Parvizi et al. 2018).During glacial periods, populations colonized suitable refugia (tolerable conditions to survive) and species distribution ranges became condensed, therefore, the availability of such suitable habitats has played a key role in the survival of populations (Hartl and Clark 1997).After the last glacial period, the populations began to spread out from their refugia and expand their distribution range (Avise 2000;Hewitt 2000;Provan and Bennett 2008).Since the last glacial period had a slight impact on the southern face of the Alborz Mountains (Bobek 1963), several refugia emerged in these areas (Parvizi et al. 2018), providing suitable conditions for populations of freshwater crabs to survive during the last glacial period.While these populations contracted their range in refugia and became temporarily isolated, genetic drift and bottleneck events have affected their micro-evolutionary patterns.Along similar lines, glacial refugia were probably the reason for the primary isolation of the western and eastern clades of the Chinese endemic freshwater crab Sinopotamon acutum and local adaptations after dispersal (Fang et al. 2015).According to our results, some localities along the western Namak Lake had suitable conditions to play a role as microrefugia during the last glacial period.Since the temperature of lowlands decreases during glacial maxima, some populations of P. elbursi probably migrated to higher elevations in the western Alborz and the Namak Lake drainages.After glacial periods ended and warming conditions prevailed, the freshwater crabs increased activities and expanded their range.Davis et al. (2016) showed in a study on Anomura (hermit crabs, king crabs, porcelain crabs, mole crabs and squat lobsters) that the warming of the climate has a positive effect on the distribution of freshwater crustaceans.According to our species distribution modelling, daily alteration and temperature annual range are two factors which mostly influence the distribution of P. elbursi.These results suggest that populations of P. elbursi are able to disperse and expand their ranges during optimal conditions.
The result of niche modeling also showed that in the last glacial period P. elbursi was distributed in the north and northwest of Iran, whereas currently it is more restricted in the northwest and the northern distribution has declined.We assume that in the northern area the number of crab populations has reduced mainly due to extensive human activities (urban growth, agricultural expansion, deforestation, dam construction and water pollution) and local climatic fluctuations.In response to these habitat changes, populations of P. elbursi have possibly contracted their range to upstream areas of rivers in the mountains to find more favorable habitats.But in the northwest area, due to the existence of dispersal corridors and less human activity, these populations had the opportunity to increase their distribution range.Consequently, populations in the northwest are less isolated and share their haplotypes with nearby populations and genetic distance among them is minimal.Despite the philopatric nature of freshwater crabs, individuals are able to avoid droughts by migrating among rivers during the rainy seasons.If there is sufficient moisture, they can also enter the territory of neighboring populations by moving across the river (Ng et al. 2001;Cumberlidge et al. 2005;Daniels et al. 2006a, b;Cook et al. 2008;Poettinger and Schubart 2014;Keikhosravi and Schubart 2014a).Therefore, high gene flow, as shown in our results for some populations, can be excluded for the Namak Lake Basin more than for the Caspian Sea Basin.For example, the Taleghan population (belonging to the Caspian Sea drainage) is located in a humid area, and since the encountered conditions were optimal, their population have likely remained in their original habitat.While the results showed highly limited haplotype and nucleotide diversities in the Taleghan population (partially shared haplotypes with Jajrood and Sepidrood), it is also one of the populations with the highest genetic differentiation from other populations and is genetically isolated from all other populations by having private haplotypes.
Our results also show nonsignificant genetic differentiation among geographically distant populations.Sepidrood (belonging to the Caspian Basin), for example, is not distinct and shares haplotypes with Jajrood and Kinevars populations that belong to the Namak Lake Basin.However, Sepidrood is significantly distinct from geographically close populations like Taleghan, even though they belong to the same drainage.The reason for such indistinguishable population structure, as was suggested by Keikhosravi et al. (2015), may be due to the incomplete separation in lineages (Funk and Omland 2003) or human translocation between rivers and basins, which is common in crustaceans (Charmantier 1992;Noël and Guinot 2007;Jesse et al. 2009).Crabs may have been transferred coincidently with fish fingerlings between the two rivers over the past few decades as the most likely reason for such genetic homogeneity (Coad 1980;Keikhosravi et al. 2015).The population of Kinevars also shared some of its haplotypes with the Telvar and Bijar which belong to the southern Caspian Basin.Although these populations belong to different drainages, closeness of the headwaters can be another explanation for such homogeneity.Since Kinevars is located in an arid area and populations are faced with unfavorable conditions, they possibly cross the land bridges between headwaters during rainy seasons and thereby enter neighboring water system.Freshwater crabs have the potential to cross short terrestrial distances during periods of heavy rainfall, when headwaters of different populations originate from the same mountain regions (Daniels 1998(Daniels , 2006a, b;, b;Cook et al. 2008;Poettinger and Schubart 2014;Keikhosravi and Schubart 2014a;Schubart and Santl 2014).
The present results also show that P. elbursi populations have positive growth rate (effective population size in the past: 0.002, population size until now: 2.304).However, Keikhosravi and Schubart (2015) suggest that populations of this species are threatened by human activities, such as dam construction and habitant degradation.But due to life-history characteristics such as rapid growth, high fecundity, rapid sexual maturity, and comprehen-sive parental care (Cumberlidge et al. 2005) in a suitable condition, they can not only maintain their population stability, but also expand their ranges and occupy appropriate habitats.Although based on our SDM analysis P. elbursi must have partial sympatric distribution with its congener, P. ibericum (see Kalate et al. 2017), across the Caspian coast lines, it does not coexist with P. ibericum.This can be a result of interspecific competition between the two species.

Conclusion
Two main factors, including a lack of knowledge and destructive human development activities, have contributed to the rapid reduction in current biodiversity (Dobson 1996, Primack 2006;Thomas et al. 2022).Phylogeographic studies have become a fundamental procedure of any successful conservation attempt, as they allow to preserve the natural distribution of a species (Zaccara et al. 2004).Our phylogeographic study of P. elbursi clearly shows a high level of genetic differentiation in some populations and strong localization in others.Due to the increasing anthropogenic and climate change impacts on freshwater habitats, which have resulted in a reduction of 92-100% of suitable habitats for P. elbursi and five other potamids (Yousefi et al. 2022), conservation studies on P. elbursi and other species inhabiting freshwater ecosystems are crucial. 6.

Figure 2 .
Figure 2. The statistical parsimony network of Potamon elbursi was applied by TCS based on a 711 base pair alignment of the Cox1 gene for a total of 131 specimens.The size of the circles represents the frequency of the haplotype in the whole data set and each line illustrates one substitution and filled circle on lines representing missing intermediate haplotypes.A-F correspond to main haplotypes for the ease of citation.

Figure 3 .
Figure 3. Scatter plot of the Mantel test which was performed using Alleles In Space, assessing the relationship between genetic distance (Φst) and geographic distance among populations of Potamon elbursi.

Figure 4 .
Figure 4.The plot of mismatch distribution.Blue curved line represents the expected distribution under the sudden expansion model, and black peaked line represents the observed distribution.

Figure 5 .
Figure 5. Bayesian skyline plot showing mtDNA effective population size change in Potamon elbursi.The x axis shows time before present in million years (MY), going backward from left to right.The y axis shows the population size.Horizontal line represents the median parameter estimate with the 95% highest posterior density interval.

Figure 6 .
Figure 6.The receiver operating characteristic (ROC) curve and AUC of Potamon elbursi in Iran for (A) current, and (B) last glacial period (~ 22 ka BP) inferred from the MAXENT analysis.

Figure 7 .
Figure 7. Potential distribution of Potamon elbursi in Iran in the current and last glacial period (~ 22 ka BP) based on the MAXENT analysis.

Table 2 .
The main indices of genetic diversity of Potamon elbursi populations.Including the number of haplotypes (h), haplotype diversity (Hd), nucleotide diversity (π), and number of segregating sites (s).

Table 4 .
AMOVA result.Partitioning of molecular variance within and among sampling sites.Variance significance of components was tested by 1,000 permutations and p < 0.05.

Table 5 .
The neutrality test for two drainage systems of Potamon elbursi.The level of significance for D is p < 0.1 and for Fs is p < 0.05.

Table 3 .
Genetic distance Φst among Potamon elbursi population between two drainages.The level of significance is p < 0.01 and NS mean no significance.Boldface font indicates nonsignificant distance.

Table 7 .
Relative importance of bioclimatic variables included in the models.Boldface number indicates the highest contribution.