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).

Figure 1. 

Sampling locations of Potamon elbursi.

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′-GGTCAACAAATCATAAAGATATTGG-3′) (Folmer et al. 1994) and COH16 (5′-CATYWTTCTGCCATTTTAGA-3′) (Schubart 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 (http://technelysium.com.au) and aligned using BIOEDIT (version 5.09; Hall 2001). Sequences were submitted to the European Molecular Biological Laboratories database (EMBL) under the accession numbers: OR708714–16, OR708722–23, OR709672, OR709686, OR709689, OR709737, OR709739, OR709740, OR710927–28, OR710947–48, OR710977, OR711012, OR711004.

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).

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 / Φst-1) / 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.

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 p-values 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₀ and theta₁ (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 non-significant 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 non-fitting 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. 2018, 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).

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 WorldClim 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.

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).

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 star-shaped 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. 2H4-D). 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.

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.

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).

We also found evidence of low population differentiation among some of the populations (Φst < 0.5). Specifically, Ghamchay and Arpachay showed low Φst values in multiple pairwise comparisons (Φst ranging from −0.007 to 0.46). Interestingly, Sepidrood (belonging to the Caspian Sea drainage system) was not differentiated from those populations belonging to the Namak Lake drainage, despite their considerable geographical distance (Φst ranging from −0.03 between Sepidrood and Jajrood to 0.28 between Sepidrood and Darake). Some populations within the same drainages could not be differentiate from each other, for example, two populations belonging to the Namak Lake, including Kinevars and Jajrood, were not significantly distinct even with a large geographical distance (Φst = 0.013, p < 0.05). Some populations belonging to the Caspian Sea basin included Shahrichay with Ghamchay (Φst = −0.007, p < 0.05) and Sariaghol (Φst = 0.01, p < 0.05), Ghamchay with Arpachay (Φst = 0.012, p < 0.05), and Arpachay with Sariaghol (Φst = −0.02, p < 0.5) and Shahrichay (Φst = −0.03, p < 0.5) were not differentiable (Table 3).

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.

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.

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.

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 (θ₀ = 2N₀μ and θ₁ = 2N₁μ), in which N₀ and N₁ indicate population size in the past and present respectively, and μ demonstrates the mutation rate in the mitochondrial sequence which is estimated to be 2.33% per million years (Schubart et al. 1998). Effective population size in the past was N₀ = 0.002 and in the present was N₁ = 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). 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).

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. 

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.

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 predictions 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).

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. 

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