Skip to main content

Advertisement

Microevolution of the noble crayfish (Astacus astacus) in the Southern Balkan Peninsula

Abstract

Background

The noble crayfish (Astacus astacus) displays a complex historical and contemporary genetic status in Europe. The species divergence has been shaped by geological events (i.e. Pleistocene glaciations) and humanly induced impacts (i.e. translocations, pollution, etc.) on its populations due to species commercial value and its niche degradation. Until now, limited genetic information has been procured for the Balkan area and especially for the southernmost distribution of this species (i.e. Greece). It is well known that the rich habitat diversity of the Balkan Peninsula offers suitable conditions for genetically diversified populations. Thus, the present manuscript revisits the phylogenetic relationships of the noble crayfish in Europe and identifies the genetic make-up and the biogeographical patterns of the species in its southern range limit.

Results

Mitochondrial markers (i.e. COI and 16S) were used in order to elucidate the genetic structure and diversity of the noble crayfish in Europe. Two of the six European haplotypic lineages, were found exclusively in Greece. These two lineages exhibited greater haplotypic richness when compared with the rest four (of “Central European” origin) while they showed high genetic diversity. Divergence time analysis identified that the majority of this divergence was captured through Pleistocene, suggesting a southern glacial refugium (Greece, southern Balkans). Furthermore, six microsatellite markers were used in order to define the factors affecting the genetic structure and demographic history of the species in Greece. The population structure analysis revealed six to nine genetic clusters and eight putative genetic barriers. Evidence of bottleneck effects in the last ~5000 years (due to climatic and geological events and human activities) is also afforded. Findings from several other research fields (e.g. life sciences, geology or even archaeology) have been utilized to perceive the genetic make-up of the noble crayfish.

Conclusions

The southernmost part of Balkans has played a major role as a glacial refugium for A. astacus. Such refugia have served as centres of expansion to northern regions. Recent history of the noble crayfish in southern Balkans reveals the influence of environmental (climate, geology and/or topology) and anthropogenic factors.

Background

Quaternary climatic oscillations have influenced the flora (e.g. [1,2,3]) and fauna (e.g. [4,5,6]) of Europe. Pleistocene glaciations had substantial impact on the flora and fauna directly via major biogeographic events (e.g. displacement) and habitat alterations, while indirectly, via fluctuations in environmental conditions [7,8,9]. The distribution and structure of obligate freshwater organisms frequently reflects historic, geological processes (such as tectonic activity, sea level change and glaciation) due to their dependency to the aquatic environment [10]. The long-term survival of many species depended on refugia and their current distribution, genetic structure and diversity reflects such historical processes [7,8,9, 11]. In Europe, the Balkan Peninsula is considered one of the major glacial refugia for many species [8]. It, also, played an important role in the colonization of eastern and western parts of Europe [7]. The majority of the European temperate fishes seem to have colonized the continent from the Black Sea through the Danube and Dnieper rivers [8].

More specifically, Greece with its complex geographic landscape [12], coupled with its geographic location (southernmost part of the Balkan Peninsula), offered suitable conditions for many species during glaciations (e.g. [13, 14]). In this area, the geological processes (tectonic and seismic activity) as well as climatic conditions forged a complex geographic landscape with rich habitat diversity [12]. Consequently, several restricted alluvial reaches [12] and small climatically stable areas [15] were formed, influencing the flora [1, 2] and fauna [13, 14]. During the mid- and late Holocene, the anthropogenic influence was intensified, generating the modification - to some extent - of the natural environment [16, 17].

An example of obligate freshwater species that could serve as a “model organism” for tackling historical events and biogeographical processes is the noble crayfish (Astacus astacus). It is a well-established European freshwater crustacean and its distribution expands from Norway (North end) to Greece (South end) and from France (West end) to Russia (East end) [18]. It plays an important role on the freshwater ecosystems, with its wide habitat usage and biological features [19]. The species is harboured in a plethora of inland freshwater habitats (streams, rivers, lakes, ponds and reservoirs [19]) and it is an opportunistic, omnivorous feeder [2, 3]. Sexual maturity is reached between 2 and 5 years of age [19, 20]. Females may be reproductively inactive for long periods (from one to several years), depending on the ambient water temperature [21]. The estimated life expectancy varies between 13 and 20 years [20], with the upper limit considered as uncommon [22], anecdotal [23] or even doubted [20]. Cukerzis (1988, cited in Keller [24]) estimated the probable longevity of the noble crayfish in Central Europe to be 7-8 years. The species has high oxygen demands [20], is sensitive to organic pollution [25] and is considered as a water quality indicator [26].

The economic and cultural value of the noble crayfish is well known in the Central and Northern Europe [19]. At the same time, several humanly induced interventions (e.g. habitat degradation, introduction of foreign species and crayfish plague) affect negatively the noble crayfish ([23] and references therein). Despite their wide distribution, the IUCN estimated a decline rate of the species between 50% and 70% [23]. Thus, the “IUCN Red List of Threatened Species” classifies the noble crayfish as a vulnerable species [23], while the “Red book of threatened fauna of Greece” classifies it at the unknown status [27]. It is worth noting that international treaties and conventions have attributed special protection to the noble crayfish (Bern Convention Appendix III, EU Habitats Directive Appendix V and directive 92/43/EEC). In order to compensate the decline or loss of populations, reintroduction of the species has been widely applied in Europe as a management tool [28]. However, variation of regulations between European countries and regional authorities based mainly on cultural traditions [19] complicates the preservation of the species.

Translocations may affect noble crayfish by influencing its range expansion, gene flow and gene pool, as identified by several studies using various molecular markers [16, 17]. It is important to notice that the majority of those translocations have disregarded the genetic structure of noble crayfish populations. The necessity of recording the genetic pool prior translocations has been recently pointed out [16, 17, 29]. To our knowledge, only one study has used microsatellite markers as monitoring tool of already translocated noble crayfish populations in the Czech Republic; the latter research was focused on the genetic variation of two translocated populations to evaluate the success of the managerial protocol applied in the last decade [30].

In the past, increasing attention was paid upon the genetic structure of the noble crayfish populations using allozymes ([31] and references therein), ISSR-PCR [32], RAPD-PCR [33] and microsatellites within rDNA-ITS1 ([34] and references therein). During the last years, three large-scale genetic analyses have been published, diagnosing the genetic diversity of the species between large geographical areas [16, 17, 29]. Microsatellite analysis differentiated Northern European (Estonia, Finland and Sweden) populations from Central European (Germany and Czech Republic) populations, with the former exhibiting lower genetic variation [17]. Furthermore, Schrimpf et al. [16, 29] revealed higher genetic diversity in the Black Sea basin populations compared to those of Southern Baltic, North and Adriatic Seas.

To our knowledge, there is currently no molecular study that encompassed samples from the southernmost limit of the species distribution (i.e. Greece). It should be noted, there are a few studies incorporating Greek crayfish to their genetic analysis, but they all refer to the stone crayfish, Austropotamobius torrentium [5, 35]. The few, sporadically published, studies about the noble crayfish in Greece are reporting on its geographical distribution and/or general information ([36] and references therein).

The lack of genetic information on the noble crayfish of Greece, coupled with the particular geographic features of the region, initiated the present study. This study largely focuses on Astacus astacus populations in the southernmost Balkan Peninsula, which has not been included in previous assessments [16, 29] although it is considered as a potentially important glacial refugium. Mitochondrial (COI and 16S) haplotypes were used to further investigate the phylogenetic relationships of the noble crayfish in Europe. Knowing that mitochondrial DNA may have little practical value in population/conservation genetic studies of widespread organisms [37], microsatellite markers were used in order to detect potential genetic differentiation and recent demographic events that may have shaped the population structure of the noble crayfish in its southernmost distribution range (i.e. Greece). The genetic make-up and the biogeographical patterns of the Greek populations are analyzed to identify factors affecting its current structure and distribution. To accomplish the above goals (i.e. phylogeography of A. astacus in Europe and the genetic structure of the species in its southernmost distribution range) extensive sampling in the continental Greece was performed in an exhaustive manner. The genetic structure of the species is meant to be utilized as a baseline for management and conservation activities.

Methods

Sample collection and DNA extraction

Two hundred eighty four potential noble crayfish sites have been investigated in continental Greece. The selection of the potential sites was inferred from the bibliography ([36] and references therein), information from local residents and habitat requirements of the species (high oxygen demand and availability of shelters [20]). In order to maximise our chances to capture crayfish in areas where sporadic observations were reported from local residents, a wide search of the river catchment was performed. For instance, in Pinios river catchment (Thessaly district) a total of 53 potential sampling sites were visited and explored. Similarly, in Kalamas river catchment (Epirus district), 7 sampling sites were searched for crayfish. The thorough examination of the continental Greece resulted in the sampling of noble crayfish in 21 sites (from a total of 284 explored sites).

Most of the noble crayfish samples were collected by hand during the day. When the conditions were not favorable, LiNi traps [38] baited with meat were used during the night. The only exception was one individual in KEF site (see Table 1), which was captured while moving <2 m from the riverbank fleeing from it. In all potential sites sampled by hand, an area between 50 and 250 m was thoroughly examined (following the upward direction of the river flow) and lasted between 30 min and 1 h 30 min (sampling was performed by at least two researchers). Sampling was repeated in sites where the number of sampled individuals was less than 20 or if the presence of noble crayfish had been reported previously. Individuals were transferred to the laboratory in an isotherm with ice packs (individuals from different sites were stored in separate plastic bags). Fishing license for research purposes was retrieved from the Greek authorities. In total, 284 noble crayfish were collected from 21 sites (Table 1) and stored at −20 °C. Genomic DNA was extracted from pereiopods muscle tissue (≈ 10 mg) following the protocol of Estoup et al. [39].

Table 1 Information on sampling collection, where sites, abbreviation, habitat type, geographical coordinate, collection method, year of sampling (Year), sample size, number of mitochondrial sequences (mtDNA), number of microsatellite genotyped data (nDNA) and status are given

Mitochondrial amplification and sequencing

Phylogenetic analysis of the noble crayfish was carried out using two mitochondrial markers: 1) a partial 16S sequence using 1471 and 1472 primers [40], and 2) a partial PCR-amplified sequence of the Cytochrome Oxidase subunit I (COI) using the universal LCO1490 and HCO2198 primers [41]. PCR reactions were adapted from Maniatsi et al. for COI [42] and Baxevanis et al. for 16S [43], with the following modifications: 0.25 μl KAPA Taq DNA Polymerase (KAPA Biosystems, South Africa), 2.5 μl 10× KAPA Taq Buffer A (KAPA Biosystems, South Africa), and 0.75 (for 16S) and 1.5 (for COI) mM MgCl2. The PCR thermal profile for COI followed Trontelj et al. [5]. For 16S the annealing temperature was slightly modified (45 °C) from Crandall and Fitzpatrick [40]. Sequencing reactions (both directions) were electrophoresed on a PRISM 3730×l DNA analyzer (Applied Biosystems, Foster City, USA) and prepared by Macrogen Inc. (Seoul, South Korea).

Mitochondrial analysis

The dataset used for the analysis comprised 37 randomly chosen noble crayfish (from a total of 284 individuals) collected from 21 sites of Greece (one or two individuals per site; Table 1) and sequences retrieved from the GenBank (19 for 16S and 32 for COI [6, 29, 35]; accession numbers DQ320033, KX370092, KF888279 to KF888295 for 16S and AY667146, KX369672, KF888296 to KF888325 for COI; Additional file 1). All produced sequences during this study have been deposited in GenBank (accession numbers KY048193 to KY048202 for 16S and KY067207 to KY067228 for COI; Additional file 1). Identity of amplified regions was confirmed using the BLAST searches. Sequences (16S and COI) were viewed in Bioedit v. 7.2.5 [44] and aligned in ClustalX v. 2.1 [45]. COI sequences were screened for pseudogenes following the procedure described in Buhay [46]. For each mitochondrial marker (COI and 16S) the heterogeneity of nucleotide frequencies, substitution and saturation and presence of phylogenetic signal were checked in DAMBE v. 5.5.29 [47]. The substitution model was determined in PartitionFinder v. 2.1.1 [48] using the Bayesian Information Criterion (BIC). Sequences (16S and COI) derived from the same individual were combined for further analysis (Additional file 1). “Concatenation” approach was followed since: 1) the two loci are linked on the mitochondria and maternally inherited [49], 2) the two genes have different resolving power [50], 3) the concatenation method can perform equally or better than methods that attempt to account for sources of error introduced by incomplete lineage sorting [51], and 4) the selection of one mitochondrial region can influence the results obtained [52].

Phylogenetic inference was performed using a strict clock model and a coalescent tree prior in Beast v. 1.8 [53]. The Markov Chain Monte Carlo (MCMC) analysis run comprised 108 generations, sampled every 104 generations. In order to determine convergence and appropriate burn-in, the Effective Sample Size (ESS) values, densities plots and trace logs for each parameter were visualized in Tracer v. 1.5 [54]. The first 106 generations (10%) were discarded as burn-in. ESS values for all parameters were >486, larger than the threshold value of 200 identified by Tracer v. 1.5. The best fit tree was found using the Maximum clade credibility tree option in the Target tree type implemented in TreeAnnotator v. 1.8. Graphical representation of the relation of the number of haplotypes and samples with the phylogroups, were implemented in R v. 3.1.2 [55]. Pairwise within and between genetic distances of the haplotype groups for 16S and COI (Kimura 2-parameter and p-distance model), were computed in MEGA v. 6.06 [56].

Molecular dating was performed in Beast v. 1.8, using the same parameters as in the phylogenetic inference analysis (see also section “Results” for the best fit models generated for COI and 16S). To our knowledge, there are no fossil records and geophysical events in order to date the phylogroups of Astacus astacus. In order to estimate divergence time, different mutation rates were taken into account: 1) mitochondrial Arthropod substitution rate (2.3% pairwise sequence divergence per million years [57]), and 2) Decapod substitution rates for 16S (0.65-0.88% pairwise sequence divergence per million years) and COI (1.66-2.33% pairwise sequence divergence per million years) genes [58]. Both approaches have been previously used in crayfish [6, 35].

Microsatellite genotyping

The analysis comprised 278 samples from 16 sites [i.e. from the 21 sampled sites, 5 sites were excluded from the subsequent analysis due to low sampling size (n < 6); Table 1]. Six microsatellite loci (Aas8, Aas766, Aas1198, Aas2489, Aas3040 and Aas3950) were used [59, 60], with the forward primer IRD800-labelled. The PCR reaction and amplification were modified from those of Kõiv et al. [59, 60] (Additional file 2). Reagents were the same as described in mitochondrial section. In all PCR reactions 0.5 μl 10× Bovine Serum Albumin (BSA) was added.

Microsatellite analysis was performed in a semi-automated Li-COR® 4200 DNA Analyzer (Li-COR, Nebraska, USA), using 6.5% acrylamide sequencing gels (Li-COR®KB Plus™). Data was genotyped in SagaGT software (Li-COR, Nebraska, USA). Allele scoring was facilitated by using the same individuals for each locus as a reference sample between gels, in combination with the molecular genetic marker (50-350 bp, Size Standard IRDye™ 800 Li-COR®). Scoring accuracy was increased by independent genotyping of the samples by two researchers.

Microsatellite analysis

Genotyping errors were assessed using Micro-Checker v. 2.2.3 [61] (95% CI and 104 repetitions). Mean number of alleles per locus, allele frequencies and observed, expected and unbiased expected Nei’s heterozygosity [62] were calculated using Genetix v. 4.05 [63]. Genepop v. 4.2.2 [64] was used to test linkage disequilibrium, conformity to Hardy-Weinberg Equilibrium for each locus and to infer genetic differentiation for all pairs of sites. For all probability tests, Markov chain method was applied using 104 generations, 20 batches and 5000 iterations per batch. Significance was assessed following the sequential Bonferroni procedure [65]. Allelic richness, number of alleles and Weir and Cockerham’s estimators (FIT, FIS and FST) were computed in FSTAT v. 2.9.3.2 [66]. RSTCalc v. 2.2 [67] was used to calculate RST (104 permutations).

Genetic structure was inferred by the Bayesian clustering method implemented in Structure v. 2.3 [68]. The conditions performed were 20 runs for each genetic cluster (K) between 1 and 16 using a 105 burn-in period followed by 2 x 106 MCMC iterations, under an admixture model with independent allele frequencies. The number of K was determined via Structure Harvester [69]. Label Clump v. 1.1.2 [70] and Distruct v. 1.1 [71] were also employed to merge the results and generate the graphical representation of clusters, respectively. Discriminant analysis of principal components (DAPC) [72] in R package Adegenet v. 1.4-2 [73] was also used to assess the number of clusters and the relationships between populations. The number of K and the optimal number of components to be retained was assessed using the find.clusters (BIC) and optim.a.score functions, respectively.

Monmonier’s maximum difference algorithm implemented in Barrier v. 2.2 [74] was employed to identify barriers to gene flow within the data set (all sixteen sites). The procedure using the FST and DCE described in [75] was followed. The tested number of probable barriers was between 1 and 10. The relationships between FST, Rousset’s distance measure FST/(1- FST), RST and geographic distances among sites were analysed using Mantel test with 104 randomizations, available at Isolation By Distance Web Service v. 3.23 (IBDWS [76]). Geographic distances were calculated in Geographic Distance Matrix Generator v. 1.2.3 [77].

Effective population size was estimated using ONeSAMP [78], which relies on Approximate Bayesian Computation. The demographic history of noble crayfish was inferred through a hierarchical Bayesian model implemented in MsVar v. 1.3 [79], using the number of genetic clusters derived from the population structure analysis (software Structure v. 2.3). For each genetic cluster, five independent runs were conducted under an exponential model, with different seeds and hyperpriors (Additional file 3). For each data set, 105 steps and 9 x 104 thinning were used. The generation time was considered as per Pianka [80], using the formula t2 = (α + ω)/2, where α is the age at maturation and ω the longevity. Since the longevity of the species is debated (see Background) two approaches were used: 1) the sexual maturity of the species (mean 3.5 years) and 2) the most probable longevity of 7.5 years (generation time 5.5). MsVar analysis was performed utilizing Aristotle University of Thessaloniki, Scientific Computing Office, IT Center, Institutional Computer Cluster resources (consuming a total of 91,561 CPU hours and 38 GB RAM capacity). 50% of each chain was discarded as burn-in. The convergence among MCMC runs was assessed visually and via the Gelman & Rubin [81] and Brooks & Gelman [82] statistics in R v. 3.1.2, using the package Coda v. 0.16-1 [83]. Point estimates of less than 1.2 were considered as a good indicator of convergence [84]. Hedge’s d and mean effective size per cluster, along with their 95% CI, were calculated as described in Paz-Vinas et al. [85].

The relative probability of the demographic event was assessed using approximate Bayes Factors (BF). A generation time of 5.5 for the noble crayfish was used in the analysis. Different levels of BF were considered as in Salmona et al. [86]. BFs were computed every 50 years for a time interval between 0 and 12,000 years (Holocene). Additionally, BFs were also computed for equal length intervals (500, 100 and 50 years) covering the past 100,000, 50,000 and 7000 years, respectively. R-scripts were based on those provided by V. Sousa (personal communication) and those included in Paz-Vinas et al. [85].

Results

Mitochondrial analysis

mtDNA variation

The full dataset comprised 66 concatenated haplotypes derived from the combination of 53 unique haplotypes for COI and 17 unique haplotypes for 16S. In COI alignment no gaps were present, whereas in 16S alignment one or two gaps were observed. The aligned final dataset had a total length of 684 bp, with 334 bp for 16S and 350 bp for COI. The best fit model generated for 16S was HKY + I [87] and for COI was HKY for codon positions 1 and 2, and HKY + Γ for codon position 3. The haplotype diversity was 0.6857 for 16S and 0.985 for COI. The number of variable and parsimony informative sites, were 17 and 11 for 16S and 49 and 32 for COI, respectively. There was no evidence for heterogeneity in base frequencies for neither mtDNA genes (16S: χ2 = 1.51, df = 165, p > 0.05 and COI: χ2 = 9.26, df = 207, p > 0.05). The substitution saturation test showed no significant saturation on either gene, since the index of substitution saturation (Iss for 16S = 0.20 and Iss for COI = 0.025) was lower than the critical value (Iss.c for 16S = 0.685 and Iss.c for COI = 0.686), with T = 41.104, df = 333, p < 0.001 for 16S and T = 73.571, df = 349, p < 0.001 for COI.

Phylogenetic inference

The tree inferred from the Bayesian phylogenetic analysis revealed six distinct genetic clades (Fig. 1). Phylogenetic clades G1 and G2 represent haplotypes found in Greece. Lineage G3 consists of haplotypes originated from Croatia, Montenegro and Germany. G4 contains haplotypes from Romania, Kosovo, Czech Republic and Austria. Phylogroup G5 comprises haplotypes belonging to several European countries (Romania, Belgium, Hungary, Bulgaria, Germany, Poland, Finland, Czech Republic and Austria). Finally, phylogenetic clade G6 is composed of Croatian, Montenegrin, German and Romanian haplotypes. It should be noted that minor discrepancies (identified in the shallow branches) have been observed between the independent gene trees produced for 16S and COI (data not shown); however, the major phylogroups were recovered in both trees.

Fig. 1
figure1

Bayesian inference topology of phylogenetic relationships among 66 concatenated mtDNA haplotypes of the noble crayfish, using Beast v. 1.8. Values at nodes represent posterior probabilities. Major nodes for the noble crayfish are given (A to E). Phylogenetic clades are represented by symbols from G1 to G6, for more information see Results

All newly identified haplotypes were unique to Greece, and represent 36.4% of the total number of haplotypes (24 out of 66 haplotypes). The percentage of the phylogroup-specific haplotypes divided by the number of individuals (V) was greater than 58.06% in both phylogenetic lineages G1 and G2. In contrast, the percentage V for phylogroups G3 to G6 was between 4.87% and 16.32% (Fig. 2). Furthermore, the between genetic distances of the haplotype groups ranged from 0.1% to 1.9% for 16S and 0.6% to 4.1% for COI. The within genetic distances of the phylogroups were between 0.13% and 1.02% for 16S, and 0.37% and 3.5% for COI (Additional file 4).

Fig. 2
figure2

Graphical representation of the composition of the phylogenetic clades (G1 to G6), based on: (a) the number of haplotypes (red color) and samples (green color); (b) the percentage of the phylogroup-specific haplotypes divided by the number of individuals (V; blue color). Where: *, based on [29]; **, based on [29, 35]; *** based on [6, 29, 35]

Age estimates

Estimates of divergence times for the concatenated mtDNA dataset showed that the noble crayfish diversified during the Pleistocene. Similar time estimates were produced by the substitution rate approaches (Decapoda and Arthropoda; Table 2). Therefore, estimated time divergences using the mutation rate of Decapoda are further discussed. Node A, represents the separation of the southernmost distribution of Astacus astacus (i.e. Greece; phylogroups G1 and G2) from the remaining groups of Europe (G3 to G6) at 1.765 ± 0.0036 MYA (Table 2). The separation of the two Greek haplogroups (G1 and G2) occurred at 1.245 ± 0.0029 MYA (node B). The first group to diversify from the remaining groups of Europe is G3, at 1.284 ± 0.003 MYA (node C). Node D followed (separating G4 from groups G5 and G6) with an estimated divergence time at 0.792 ± 0.0022 MYA. Phylogenetic groups G5 and G6 shared for the last time a common ancestor at 0.556 ± 0.0016 MYA (node E). The divergence within groups (G1 to G6) was estimated at 1.765 ± 0.0036, 0.705 ± 0.0021, 0.761 ± 0.0022, 0.356 ± 0.015, 0.468 ± 0.0014 and 0.285 ± 0.0012 MYA, respectively.

Table 2 Estimates of divergence time for the concatenated mtDNA (16S and COI) data set, based on the substitution rate of Decapoda or Arthropoda

Microsatellite analysis

General summary statistics

Putative null alleles were observed in 3 microsatellite loci and 5 sites (Table 3). Omitting Aas8, Aas2498 or Aas1198 from the analysis did not have any substantial effect on the results (following the procedure described in [39]; data not shown). Therefore, all loci were included in the subsequent analysis. All microsatellite loci were polymorphic, with the number of alleles per locus ranging from 11 (Aas766) to 39 (Aas1198) (Additional file 5). A total of 111 different alleles were recognised with an average of 18.5 alleles per locus (Additional file 5). For all pairs of loci across all sites, 5 out of 240 global tests for genotypic linkage disequilibrium were statistically significant (none after Bonferroni correction) and none among locus pairs (Fisher’s method). Deviation from Hardy-Weinberg equilibrium was present at 9 out of 96 tests (none after Bonferroni correction).

Table 3 Population genetic parameters inferred from 6 microsatellite loci for 16 sites of the noble crayfish

The mean allelic richness of each site varied from 1.86 (KRN) to 5.57 (KAS) and the mean observed heterozygosity ranged from 0.19 (KRN) to 0.77 (KAS and DOX) (Table 3; see also Additional file 5). Overall inbreeding coefficient FIS for each sampled site showed heterozygote excess in several sites (AA2, LOG, NEO, KRN and KPA) due to their negative values (between −0.271 in LOG and −0.033 in AA2). Maximum value of FIS was recorded at AA1 (0.152; Table 3). The total number of private alleles was 33. Sites DOX and TSV had the highest number of private alleles with 8 counts, while in AA1, KNS, SKR, KRN and KPA no private alleles were observed.

Exact tests based on allele frequencies for all pairs of sites showed that 6 out of 240 tests were significant (none after Bonferroni correction). The overall genetic heterogeneity was high (FST = 0.4 and RST = 0.39), with the estimated values of FST generally greater than RST for most of the loci (except Aas3040 and Aas1198; Additional file 6). Pairwise estimates of FST ranged from −0.0069 (KLV-NEO) to 0.7499 (KRN-SKR) and RST from −0.0325 (KLV-NEO) to 0.8778 (KRI-KRN) (Additional file 7).

Population structure

The inferred number of clusters (K) was nine for the analysis using Structure software (Additional file 8). Sites were grouped together as follows: AA1-AA2-KPA (cluster 1), KRI-SKR (cluster 2), KLV-NEO (cluster 6), KAS-TZA (cluster 8) and LOG-KNS-KRN (cluster 9), while each one of the other sites (PLF, TSV, PRT and DOX) formed different and distinct genetic clusters (cluster 3, 5, 7 and 4, respectively; Fig. 3). The proposed number of clusters for the DAPC analysis varied between 6 and 9 (Additional file 8). It is important to notice that assignment probability for each K identified by DAPC was >0.978 (Additional file 9). In both programs, K = 9 had a similar grouping and was based on the geographic proximity of the sampling sites (Fig. 3).

Fig. 3
figure3

Graphical representation of the clustering analysis, using: (a) Structure v. 2.3 (for K = 9); (b) R package Adegenet v. 1.4-2 (for K between 6 and 9). In both analyses each bar represents an individual that is colored based on its assignment probability (colors used in both analysis for K = 9 are similar). All K are sorted based on the sampling sites. In (a) the abbreviation of the sites (Sites, below) and a number representing the cluster (Clusters 1 to 9, on top) identified are given. In (b) for each K identified by the DAPC every cluster is identified by a number (the same numbering is used in Structure)

The majority of the genetic structure in DAPC for nine clusters was captured in the first principal component (screeplot of the eigenvalues; Fig. 4). Visual inspection of the first two principal components, revealed two distinct genetic clusters (clusters 2 and 9; Fig. 4). Similar results were observed for K between 6 and 8 (Additional file 10). No evidence of isolation by distance was retrieved by Mantel tests (no statistically significant correlation between genetic and geographical distances was observed - Additional file 11).

Fig. 4
figure4

DAPC scatterplot of the first two principal components for K = 9. Clusters are represented by several distinguishable colors (same as in Structure, Fig. 4). The inset shows the discriminant analysis (DA) eigenvalues

Barriers

Barrier software revealed the occurrence of breaks in the genetic flow across continental Greece. Analysis of all microsatellite loci simultaneously using DCE or FST identified eight probable barriers. DCE approach revealed at least 93% bootstrap support for every barrier. The locus by locus analysis of FST showed that five out of six microsatellite loci (with the exception of Aas3040) agreed with the occurrence of all barriers (Additional file 12). In Aas3040 the observed differences are due to the lack of barriers between sites KRN-LOG-KNS and NEO-KLV, as well as between AA1-AA2-KPA and KAS-TZA. Since the number of data derived from the analysis was large (totally 100; Additional file 12) and the outcomes were similar, it was decided to provide the output data for eight putative barriers, 16 sites and all microsatellite loci using both FST and DCE (Fig. 5). The derived genetic barriers are closely related with the presence of mountains and river systems, corroborating the findings of structure analysis for K = 9. For instance, sites KRI and SKR are delimited by Mt. Olympus, Mts. Pieria and Aliakmon river. The names of the mountains and major rivers are provided in Additional file 13.

Fig. 5
figure5

Genetic barriers created via Delaunay triangulation (green lines) and Voronoi tessellation (blue polygons), as predicted by Barrier software for all sixteen sites. Red lines constitute the genetic barriers detected through (a) the bootstrap analysis (104 bootstraps) of DCE and (b) the analysis of the pairwise FST matrix for all microsatellite loci. Circles represent the sites (black letters). For (a) the thickness of each genetic barrier is proportional to the bootstrap support (green numbers). A background map, created using ArcGIS® and ArcMap™ by Esri, is provided in order to visualize the landscape (Copyright © 2014 Esri and its licensors. All rights reserved)

Demographic history and effective population size

Point estimates of the genetic clusters varied between 1 and 1.09 (Additional file 14), indicating a good convergence of MCMC runs. MsVar analyses indicated a decline in all genetic clusters, with no apparent overlap of the posterior distributions of past log10(N1) and present log10(N0) population sizes (Additional file 15). Analyses of the global data set (random samples from different demes were merged into one dataset) gave similar results excluding the presence of a false bottleneck signal (Additional files 14 and 15; as suggested by [56]). Since the descriptive statistics (Table 4) and the statistics describing the magnitude of demographic change (Table 5) are similar for either generation time used (3.5 and 5.5), it was decided to further discus the analysis assuming a generation time of 5.5. The ancestral population size log10(N1) among the genetic clusters ranged from 4.97 to 5.51 (Table 4). The contemporary population size log10(N0) varied between 0.5 and 1.42 (or mean N0 3.17 and 26.6, respectively) (Table 4). The current effective size estimated for each genetic cluster in ONeSAMP was small, with values varying between 22.069 and 44.191 (Table 6). Mean effective population size (Mean log10(N0/N1)) of the genetic clusters decreased by 14 to 20% (average of approximately 17%, Table 5). The magnitude of the bottleneck was strong in every genetic cluster, with negative values for Hedge’s d (between −6.02 and −9.91) (Table 5). Time since the population decline log10(T) (Additional file 16) varied between 2.27 and 2.77 among genetic clusters (Table 4). BF analysis indicated that the decline occurred during the last around 5000 years (BF > 3, Fig. 6; see also Additional file 17).

Table 4 Mean values of log10(N0), log10(N1) and log10(T) for every genetic cluster (1 to 9) with a generation time of 3.5 and 5.5
Table 5 Statistics describing the magnitude of the demographic change for each genetic cluster (1-9) and generation time (3.5 or 5.5)
Table 6 Mean estimated contemporary size (N0) for each genetic cluster of the noble crayfish, along with their 95% Confidence Intervals (CI) (OneSAMP software)
Fig. 6
figure6

BF values for equal time intervals (100 years) for the past 12,000 years (Holocene). Each genetic cluster (1 to 9) is represented by a different colored-solid line (same as population structure analysis). Horizontal dotted line represents BF values of 3

Discussion

Genetic diversity of the noble crayfish

The observed haplotypic diversity of the noble crayfish (Astacus astacus) in Europe (0.6857 for 16S and 0.985 for COI, respectively) is in accordance with the wide geographic range of the species in this area [18]. Six haplotypic lineages (G1 to G6, Fig. 1) were identified using a concatenated mitochondrial data set (16S and COI). The inclusion of samples from the southern Balkans (i.e. Greece) revealed the existence of two novel phylogenetic clades (G1 and G2), separated from the rest of the European haplogroups. The between pairwise genetic distances showed that the phylogroups associated to the southernmost Balkan Peninsula (i.e. Greece, G1 and G2) were clearly differentiated from the remaining European haplotype groups (G3 to G6; Additional file 4). Furthermore, a high number of unique haplotypes in Greece were observed (24 haplotypes from 37 samples; Fig. 2), indicating a high genetic diversity. It must be noted that the number of unique haplotypes and genetic diversity of the southern distribution of A. astacus may be biased by the distribution of the sites within each haplogroup and the number of individuals per haplogroup in the Balkan Peninsula. These findings may corroborate the taxonomical division of the noble crayfish, proposed by Karaman [88]; three subspecies are present in Europe, A. astacus astacus, A. a. colchicus and A. a. balcanicus. Based on our results, it is impelling to suggest that the Greek phylogroups may represent the subspecies A. a. balcanicus. Haplotype groups G3 to G6 (Fig. 1) are consistent with a previous study [29]. Additionally, a lower genetic diversity was observed in central and northern Europe (G4 to G6) compared to the Greek (G1 and G2) and the Croatian (G3) lineages (Fig. 2; Additional file 4), probably due to bottleneck and founder effects [16, 17, 29]. Several studies indicated that this pattern could be the result of a post-glacial range expansion from south to north regions [7, 8, 11]. Genetic distances (16S and COI; Additional file 4) were lower to those observed in Austropotamobius torrentium in Europe [35], probably reflecting different life-history traits. Phylogroup G5 indicates a pattern typical of a recent range expansion, based on the number of observed haplotypes. This pattern of noble crayfish distribution range is identified as a typical Southern glacial refugium [9]. Therefore, a “Mediterranean” molecular biogeographical pattern [11] and a “grasshopper” [7] or “chub” [8] model of expansion, can be deduced for the noble crayfish. The most likely colonization route of the noble crayfish towards European regions is through the Danube river (as suggested by Schrimpf et al. [29]) and observed, also, in temperate freshwater fishes [8]. In contrast, it has been assumed that the non-endemic Plio-Pleistocene freshwater fishes of Greece originated from further north regions, e.g. the Danube river [89]. The present analysis indicates that Greece is an older glacial refugium than the eastern Black Sea (phylogenetic clade G4) and Western Balkans (phylogenetic clade G3), with the latter two regions already identified by Schrimpf et al. [29]. Nevertheless, the scenario of multiple refugia is confirmed, as suggested for the noble crayfish [17, 29] and postulated for other species (e.g. [90, 91]).

The phylogroups of the noble crayfish diverged during Early and Middle Pleistocene. Specifically, the estimated average divergence times were between ca. 0.3 and 1.8 MYA (Table 2). In the Mediterranean, major cooling events occurred between 2.15 and 2.73 MYA ([92] and references therein), covering the upper HPD limits. During the Middle Pleistocene, two global and prominent climatic transitions seems to be the core of the induced phylogenetic divergence of the species, the “Early Middle Pleistocene Transition” (0.8 to 1.2 MYA) and the “Mid-Brunhes Event” (~0.4 MYA) [93]. Moreover, a correlation can be observed between time of lineage divergence, increased freshwater inputs (last ca. 2.3 MYA [94]) and climatic oscillations (last 1.3 MYA [95]) of the eastern Mediterranean. It has been suggested that there was a connection between the Danube basin and the Greek river basins during the late Pliocene to Pleistocene; this was based on studies related to central European and Danubian freshwater fishes [89]. A geological event that might have influenced the noble crayfish is the opening of the Danube corridor between 0.7 and 1 MYA [96]. In Greece, the delta formation of Pinios river could have, also, affected the noble crayfish genetic structure by dividing the haplotypes into two well-defined clusters (~1.1 MYA; Faugères, 1977 in Caputo [97] and Rook and Martínez-Navarro [98]). It has, also, been suggested that during the glacial-interglacial cycles, freshwater or oligosaline conditions were prevailing in the upper Aegean areas [99]. These events may have played an important role in the expansion of the noble crayfish in northern regions.

Micro-scale evolution of the noble crayfish

The noble crayfish in Greece is confined in the upper parts of river systems, with a limited, scattered, discontinuous and grouped distribution (21 crayfish sites from a total of 284 investigated sites; Additional file 12). Similar distribution pattern was observed in a former study in Greece [36], but also in other European countries [20, 100]. Anthropogenic activities (agricultural land use and settlements) are negatively affecting the species distribution in riverine systems of the Czech Republic [100]. Agriculture intensification in lowlands of the Mediterranean [101] and water quality of rivers in Greece [102] could explain the distribution pattern of the Greek noble crayfish.

The genetic heterogeneity in the southernmost distribution range of the noble crayfish (i.e. Greece) (global FST = 0.4 and between sites −0.007 < FST < 0.75; Additional file 6) was comparable to wider geographical studies (i.e., global FST = 0.264 [17] and 0.008 < FST < 0.723 [29]). Structure analysis revealed a high number of genetic clusters (Fig. 3). The higher sub-structuring of Greece compared to central and northern regions of Europe [17, 29] corroborates the findings of the mitochondrial analysis.

The identified nine clusters (K = 9) revealed a relation between them and the geographic landscape of Greece (Figs. 3 and 5). Furthermore, eight putative barriers were inferred revealing that mountains played an important role in the gene flow (Fig. 5). The majority of the inferred clusters (K = 9) also revealed a close relation with their river system, suggesting geographic isolation. These findings are in accordance with field observations on wild populations of noble crayfish (indigenous to the location). Those studies indicated a sedentary nature of the noble crayfish, with low dispersal ability [20, 103]. However, discrepancies can be observed in two clusters (i.e. 1 and 9, Fig. 3), with their components belonging to two different river systems. The most likely explanation is due to human interference. Cluster 1 (sites AA1, KPA and AA2; Fig. 3) may be the result of indirect human activities, since deviation of water by the Aoos hydroelectric power plant has been reported [102]. Cluster 9 is a typical case of a recent human translocation (sites KNS and LOG translocated from site KRN; Fig. 3 and Table 1). In fact, the Geographic Distance Matrix Generator shows that in several cases it involves long-distance translocations (e.g. 53 km in a straight line, between sites KRN and LOG). The same observation was inferred for the river Kalamas, by a former study [36]. It is important to notice, that the sedentary nature of the species depends on food availability, sex, external factors and familiarity of the noble crayfish to the habitat [20, 103]. Low pairwise FST and RST values between sites from different river catchments (Additional file 6), absence of correlation between genetic and geographical distances (Additional file 11) and close genetic relation of the majority of the clusters (Fig. 4), indicated that human translocations have shaped the genetic structure of A. astacus in Greece (and not, so much, relevant biogeographic events). Information on the status of the species (Table 1) disclosed that several identified distinct clusters are due to human translocation (e.g. cluster 8, Kalamas river).

Person(s) involved in translocation of the noble crayfish in Greece had not followed good conservation practices, such as similarities or dissimilarities of habitats, riverine systems or genetic make-up of the species (e.g. sites KNS, LOG and KRN, Table 1 and Fig. 3). Similar observations are also reported in other European regions [16, 17]. Another issue that must be addressed is the nature of the translocations. Unfortunately, translocations may occur without scientific guidance and/or governmental consensus (for instance sites LOG and KNS; Table 1). The absence of appropriate controls may be an issue, since overharvesting and violation of regulations from local and recreational fishermen is common practice. For example, 671 noble crayfish with highly variable lengths (approximately between 5 and 12 cm Total Length or TL) were collected in a 24 h period in Karya by a local recreational fisherman (personal observation). It should be noted that, in Greece, trading of A. astacus is prohibited for less than 7 cm TL (Ministerial Decision A2-3354/2007, Government Gazette 2207/B/14.11.2007), while the minimum allowed fishing length is 10 cm TL (Royal Decree 142/1971, Government Gazette 49A/1971). Similar phenomena have also been observed in Norway [104]. Prohibiting and restricting measures (e.g. restriction of fishing gear) could improve management practices.

Analysis of the demographic history in the southernmost distribution of the species (i.e. Balkan Peninsula, Greece) insinuates a severe decrease in the population size during the Holocene (last ~5000 years), resulting in small effective population sizes (Fig. 6; Tables 4, 5 and 6; Additional files 14, 15 and 16). The culminating arid conditions of the last approximately 5.5 kyr BP [105] may have played an important role in the recent evolutionary history of the noble crayfish (potential and reoccurring bottlenecks). During that period, several climatic oscillations (e.g. [106]) and variations in river flood activity (e.g. [107]) have been reported (see also Additional file 17E). Furthermore, several geological changes influencing freshwater systems of Greece have been recorded during the mid- and late Holocene [97, 108].

However, for the last ~5000 years and due to the advent of the Bronze Age (~3.1-5.3 kyr BP), it is difficult to distinguish climate-induced variations from human-mediated influence on riverine systems of Greece [109, 110]. Historical records showed that from the Bronze Age (~5.3 kyr BP) until nowadays several anthropogenic activities affected the continental Greece. Changes in settlements, hydraulic works (irrigation, e.g. aqueducts, dams, water diversions and creation of artificial lakes) and wars (e.g. [111, 112]) are some of the factors that may have directly or indirectly influenced the noble crayfish populations (potential alteration or destruction of habitats and crayfish consumption by humans). The major hydraulic/engineering project in Copais lake (Central - South Greece) during the Bronze Age by the Mynians (Mycenaean civilization [111]) illustrates the intensity of the human activity in Greece. Notably, it is known that crustaceans are consumed since the Bronze Age [113].

Conclusions

The noble crayfish revealed a high genetic diversity, with its southernmost range limit playing a crucial role. Mitochondrial analysis (16S and COI) inferred six phylogroups; two of them (corresponding to the southernmost distribution of the species) were identified for the first time. Genetic analyses revealed that the genetic diversity of A. astacus populations from Southern Balkans (Greece) is high and they form a distinct group when compared with their European counterparts. The divergence time analysis implied a potential differentiation of the species during the Early and Middle Pleistocene (between ca. 0.3 and 1.8 MYA). During that period, several Pleistocene climatic oscillations and geological processes corroborate the inferred divergence. Therefore, we consider that the area of the Southern Balkans may have played an important role in shaping the genetic make-up of A. astacus populations (old centre of expansion and/or southern glacial refugium). Microsatellite analysis focusing on the noble crayfish from its southernmost distribution (Balkan Peninsula, Greece) disclosed six to nine distinct clusters. The combination of the inferred genetic clusters with the genetic barriers and the status of each sampled site (natural or translocated) unraveled the influence of several factors. Specifically, it was deduced that the landscape of the region (mountains and river systems) and human related interventions (mostly translocations) were the main factors affecting the genetic structure of the noble crayfish. Finally, it was figured out that during mid to late Holocene, a severe decrease in the effective population size of the Greek noble crayfish took place. During that time period, several climatic, geological and anthropogenic processes occurred influencing the distribution of noble crayfish in the southernmost Balkan Peninsula (i.e. Greece). It is concluded that the noble crayfish has been primarily influenced by physical processes (climate and geology), while recently, anthropogenic implications (mainly translocations by humans and niche degradation) may have affected its genetic and geographic structure.

Abbreviations

BF:

Bayes factors

BIC:

Bayesian information criterion

CI:

Confidence interval

COI:

Cytochrome oxidase subunit I

DA:

Discriminant analysis

DAPC:

Discriminant analysis of principal components

ESS:

Effective sample size

HPD:

Highest posterior density

Iss:

Index of substitution saturation

Iss.c:

Index of substitution saturation critical value

IT:

Information technology

K:

Number of clusters

MCMC:

Markov chain Monte Carlo

mtDNA:

Mitochondrial DNA

PCR:

Polymerase chain reaction

TL:

Total length

V:

Percentage of the phylogroup-specific haplotypes divided by the number of individuals

References

  1. 1.

    Tzedakis PC, Lawson IT, Frogley MR, Hewitt GM, Preece RC. Buffered tree population changes in a quaternary refugium: evolutionary implications. Science. 2002;297:2044–7.

  2. 2.

    Médail F, Diadema K. Glacial refugia influence plant diversity patterns in the Mediterranean Basin. J Biogeogr. 2009;36:1333–45.

  3. 3.

    Deffontaine V, Libois R, Kotlík P, Sommer R, Nieberding C, Paradis E, et al. Beyond the Mediterranean peninsulas: evidence of central European glacial refugia for a temperate forest mammal species, the bank vole (Clethrionomys glareolus). Mol Ecol. 2005;14:1727–39.

  4. 4.

    Jablonski D, Jandzik D, Mikulíček P, Džukić G, Ljubisavljević K, Tzankov N, et al. Contrasting evolutionary histories of the legless lizards slow worms (Anguis) shaped by the topography of the Balkan peninsula. BMC Evol Biol. 2016;16:99.

  5. 5.

    Trontelj P, Machino Y, Sket B. Phylogenetic and phylogeographic relationships in the crayfish genus Austropotamobius inferred from mitochondrial COI gene sequences. Mol Phylogenet Evol. 2005;34:212–26.

  6. 6.

    Jelić M, Klobučar GIV, Grandjean F, Puillandre N, Franjević D, Futo M, et al. Insights into the molecular phylogeny and historical biogeography of the white-clawed crayfish (Decapoda, Astacidae). Mol Phylogenet Evol. 2016;103:26–40.

  7. 7.

    Hewitt GM. Post-glacial re-colonization of European biota. Biol J Linn Soc. 1999;68:87–112.

  8. 8.

    Hewitt GM. Genetic consequences of climatic oscillations in the quaternary. Philos T Roy Soc B. 2004;359:183–95.

  9. 9.

    Stewart JR, Lister AM, Barnes I, Dalén L. Refugia revisited: individualistic responses of species in space and time. Proc Biol Sci. 2010;277:661–71.

  10. 10.

    Hughes JM, Schmidt DJ, Finn DS. Genes in streams: using DNA to understand the movement of freshwater fauna and their riverine habitat. Bioscience. 2009;59:573–83.

  11. 11.

    Schmitt T. Molecular biogeography of Europe: Pleistocene cycles and postglacial trends. Front Zool. 2007;4:11.

  12. 12.

    Macklin MG, Lewin J, Woodward JC. Quaternary fluvial systems in the Mediterranean basin. In: Lewin J, Macklin MG, Woodward JC, editors. Mediterranean quaternary river environments. Rotterdam, Brookfield: A. A. Balkema; 1995. p. 1–25.

  13. 13.

    Alexandri P, Triantafyllidis A, Papakostas S, Chatzinikos E, Platis P, Papageorgiou N, et al. The Balkans and the colonization of Europe: the post-glacial range expansion of the wild boar, Sus scrofa. J Biogeogr. 2012;39:713–23.

  14. 14.

    Karaiskou N, Tsakogiannis A, Gkagkavouzis K, Papika S, Latsoudis P, Kavakiotis I, et al. Greece: a Balkan subrefuge for a remnant red deer (Cervus elaphus) population. J Hered. 2014;105:334–44.

  15. 15.

    Tzedakis PC. Museums and cradles of Mediterranean biodiversity. J Biogeogr. 2009;36:1033–4.

  16. 16.

    Schrimpf A, Schulz HK, Theissinger K, Pârvulescu L, Schulz R. The first large-scale genetic analysis of the vulnerable noble crayfish Astacus astacus reveals low haplotype diversity in central European populations. Knowl Manag Aquat Ecosyst. 2011;401:35.

  17. 17.

    Gross R, Palm S, Kõiv K, Prestegaard T, Jussila J, Paaver T, et al. Microsatellite markers reveal clear geographic structuring among threatened noble crayfish (Astacus astacus) populations in northern and Central Europe. Conserv Genet. 2013;14:809–21.

  18. 18.

    Kouba A, Petrusek A, Kozák P. Continental-wide distribution of crayfish species in Europe: update and maps. Knowl Manag Aquat Ecosyst. 2014;413:5.

  19. 19.

    Skurdal J, Taugbøl T. Astacus. In: Holdich DM, editor. Biology of freshwater crayfish. Oxford: Blackwell Science; 2002. p. 467–510.

  20. 20.

    Cukerzis J. La biologie de l’écrevisse (Astacus astacus L.). Paris: Institut National de la Recherche Agronomique Publications; 1984. (In French)

  21. 21.

    Skurdal J, Hessen DO, Garnås E, Vøllestad LA. Fluctuating fecundity parameters and reproductive investment in crayfish: driven by climate or chaos? Freshw Biol. 2011;56:335–41.

  22. 22.

    Lundberg U. Behavioural elements of the noble crayfish, Astacus astacus (Linnaeus, 1758). Crustaceana. 2004;77:137–62.

  23. 23.

    Edsman L, Füreder L, Gherardi F, Souty-Grosset C. Astacus astacus. IUCN Red List of Threatened Species. http://dx.doi.org/10.2305/IUCN.UK.2010- 3.RLTS.T2191A9338388 (2010). Accessed 23 Jan 2016.

  24. 24.

    Keller M. Ten years of trapping Astacus astacus for restocking. Freshwater Crayfish. 1998;12:518–28.

  25. 25.

    Pârvulescu L, Pacioglu O, Hamchevici C. The assessment of the habitat and water quality requirements of the stone crayfish (Austropotamobius torrentium) and noble crayfish (Astacus astacus) species in the rivers from the Anina Mountains (SW Romania). Knowl Manag Aquat Ecosyst. 2011;401:3.

  26. 26.

    Reynolds J, Souty-Grosset C, Richardson A. Ecological roles of crayfish in freshwater and terrestrial habitats. Freshwater Crayfish. 2013;19:197–218.

  27. 27.

    Legakis A. Invertebrates. In: Legakis A, Maragou P, editors. Red book of threatened fauna of Greece. Athens: Hellenic Zoological Society; 2009. p. 428–509. (In Greek).

  28. 28.

    Jussila J, Ojala K, Mannonen A. Noble crayfish (Astacus astacus) reintroduction project in the river Pyhäjoki, western Finland: a case study. Freshwater Crayfish. 2008;16:51–6.

  29. 29.

    Schrimpf A, Theissinger K, Dahlem J, Maguire I, Pârvulescu L, Schulz HK, et al. Phylogeography of noble crayfish (Astacus astacus) reveals multiple refugia. Freshw Biol. 2014;59:761–76.

  30. 30.

    Bláha M, Žurovcová M, Kouba A, Policar T, Kozák P. Founder event and its effect on genetic variation in translocated populations of noble crayfish (Astacus astacus). J Appl Genet. 2016;57:99–106.

  31. 31.

    Fevolden SE, Taugbøl T, Skurdal J. Allozymic variation among populations of noble crayfish, Astacus astacus L., in southern Norway: implication for management. Aquacult Fish Manage. 1994;25:927–35.

  32. 32.

    Schulz HK, Šmietana P, Schulz R. Assessment of DNA variations of the noble crayfish (Astacus astacus L.) in Germany and Poland using inter-simple sequence repeats (Issrs). Bull Fr Pêche Piscic. 2004;372–373:387–99.

  33. 33.

    Schulz R. Status of the noble crayfish Astacus astacus (L.) in Germany: monitoring protocol and the use of RAPD markers to assess the genetic structure of populations. Bull Fr Pêche Piscic. 2000;356:123–38.

  34. 34.

    Alaranta A, Henttonen P, Jussila J, Kokko H, Prestegaard T, Edsman L, et al. Genetic differences among noble crayfish (Astacus astacus) stocks in Finland, Sweden and Estonia based on the ITS1 region. Bull Fr Pêche Piscic. 2006;380–381:965–76.

  35. 35.

    Klobučar GIV, Podnar M, Jelić M, Franjević D, Faller M, Štambuk A, et al. Role of the Dinaric Karst (western Balkans) in shaping the phylogeographic structure of the threatened crayfish Austropotamobius torrentium. Freshw Biol. 2013;58:1089–105.

  36. 36.

    Koutrakis E, Perdikaris C, Machino Y, Savvidis G, Margaris N. Distribution, recent mortalities and conservation measures of crayfish in Hellenic fresh waters. Bull Fr Pêche Piscic. 2007;385:25–44.

  37. 37.

    Zhang D-X, Hewitt GM. Nuclear DNA analyses in genetic studies of populations: practice, problems and prospects. Mol Ecol. 2003;12:563–84.

  38. 38.

    Westman K, Pursiainen M, Vilkman R. A new folding trap model which prevents crayfish from escaping. Freshwater Crayfish. 1978;4:235–42.

  39. 39.

    Estoup A, Largiader CR, Perrot E, Chourrout D. Rapid one-tube DNA extraction for reliable PCR detection of fish polymorphic markers and transgenes. Mol Mar Biol Biotechnol. 1996;5:295–8.

  40. 40.

    Crandall KA, Fitzpatrick JF. Crayfish molecular systematics: using a combination of procedures to estimate phylogeny. Syst Biol. 1996;45:1–26.

  41. 41.

    Folmer O, Black M, Hoeh W, Lutz R, Vrijenhoek R. DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Mol Mar Biol Biotechnol. 1994;3:294–9.

  42. 42.

    Maniatsi S, Baxevanis AD, Abatzopoulos TJ. The intron 2 of p26 gene: a novel genetic marker for discriminating the two most commercially important Artemia franciscana subspecies. J Biol Res - Thessalon. 2009;11:73–82.

  43. 43.

    Baxevanis AD, Kappas I, Abatzopoulos TJ. Molecular phylogenetics and asexuality in the brine shrimp Artemia. Mol Phylogenet Evol. 2006;40:724–38.

  44. 44.

    Hall TA. BioEdit: a user-friendly biological sequence alignment editor and analysis program for windows 95/98/NT. Nucleic Acids. 1999;41:95–8.

  45. 45.

    Larkin MA, Blackshields G, Brown NP, Chenna R, Mcgettigan PA, McWilliam H, et al. Clustal W and Clustal X version 2.0. Bioinformatics. 2007;23:2947–8.

  46. 46.

    Buhay JE. ‘COI-like’ sequences are becoming problematic in molecular systematic and DNA barcoding studies. J Crustac Biol. 2009;29:96–110.

  47. 47.

    Xia X. DAMBE5: a comprehensive software package for data analysis in molecular biology and evolution. Mol Biol Evol. 2013;30:1720–8.

  48. 48.

    Lanfear R, Frandsen PB, Wright AM, Senfeld T, Calcott B. PartitionFinder 2: new methods for selecting partitioned models of evolution for molecular and morphological phylogenetic analyses. Mol Biol Evol. 2017;34:722–73.

  49. 49.

    Daniels SR. Reconstructing the colonisation and diversification history of the endemic freshwater crab (Seychellum alluaudi) in the granitic and volcanic Seychelles archipelago. Mol Phylogenet Evol. 2011;61:534–42.

  50. 50.

    Pedraza-Lara C, Alda F, Carranza S, Doadrio I. Mitochondrial DNA structure of the Iberian populations of the white-clawed crayfish, Austropotamobius italicus italicus (Faxon, 1914). Mol Phylogenet Evol. 2010;57:327–42.

  51. 51.

    Tonini J, Moore A, Stern D, Shcheglovitova M, Ortí G. Concatenation and species tree methods exhibit statistically indistinguishable accuracy under a range of simulated conditions. PLoS Currents Tree Of Life. 2015; doi:10.1371/currents.tol.34260cc27551a527b124ec5f6334b6be.

  52. 52.

    Meiklejohn KA, Danielson MJ, Faircloth BC, Glenn TC, Braun EL, Kimball RT. Incongruence among different mitochondrial regions: a case study using complete mitogenomes. Mol Phylogenet Evol. 2014;78:314–23.

  53. 53.

    Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol. 2012;29:1969–73.

  54. 54.

    Rambaut A, Drummond AJ. Tracer version 1.5. http://beast.bio.ed.ac.uk/Tracer (2009). Accessed 21 Mar 2015.

  55. 55.

    R Core Team. R: A language and environment for statistical computing. http://www.r-project.org/ (2014). Accessed 25 Sep 2015.

  56. 56.

    Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: molecular evolutionary genetics analysis version 6.0. Mol Biol Evol. 2013;30:2725–9.

  57. 57.

    Brower AV. Rapid morphological radiation and convergence among races of the butterfly Heliconius erato inferred from patterns of mitochondrial DNA evolution. Proc Natl Acad Sci U S A. 1994;91:6491–5.

  58. 58.

    Schubart CD, Diesel R, Hedges SB. Rapid evolution to terrestrial life in Jamaican crabs. Nature. 1998;393:363–5.

  59. 59.

    Kõiv K, Gross R, Paaver T, Kuehn R. Isolation and characterization of first microsatellite markers for the noble crayfish, Astacus astacus. Conserv Genet. 2008;9:1703–6.

  60. 60.

    Kõiv K, Gross R, Paaver T, Hurt M, Kuehn R. Isolation and characterization of 11 novel microsatellite DNA markers in the noble crayfish, Astacus astacus. Anim Genet 2009;40:124–124.

  61. 61.

    Van Oosterhout C, Hutchinson WF, Wills DPM, Shipley P. MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Mol Ecol Notes. 2004;4:535–8.

  62. 62.

    Nei M. Estimation of average heterozygosity and genetic distance from a small number of individuals. Genetics. 1978;89:583–90.

  63. 63.

    Belkhir K, Borsa P, Chikhi L, Raufaste N, Bonhomme F. GENETIX 4.05, logiciel sous Windows™ pour la génétique des populations. http://kimura.univ-montp2.fr/genetix/ (2016). Accessed 1 Jan 2016.

  64. 64.

    Rousset F. Genepop’007: a complete re-implementation of the genepop software for windows and Linux. Mol Ecol Resour. 2008;8:103–6.

  65. 65.

    Rice WR. Analyzing tables of statistical tests. Evolution. 1989;43:223–5.

  66. 66.

    Goudet J. FSTAT, a program to estimate and test gene diversities and fixation indices (version 2.9.3) http://www2.unil.ch/popgen/softwares/fstat.htm#top (2016). Accessed 23 Feb 2016.

  67. 67.

    Goodman SJ. RST calc: a collection of computer programs for calculating estimates of genetic differentiation from microsatellite data and determining their significance. Mol Ecol. 1997;6:881–5.

  68. 68.

    Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155:945–59.

  69. 69.

    Earl DA, VonHoldt BM. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv Genet Resour. 2012;4:359–61.

  70. 70.

    Jakobsson M, Rosenberg NA. CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics. 2007;23:1801–6.

  71. 71.

    Rosenberg NA. DISTRUCT: a program for the graphical display of population structure. Mol Ecol Notes. 2004;4:137–8.

  72. 72.

    Jombart T, Devillard S, Balloux F. Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genet. 2010;11:94.

  73. 73.

    Jombart T. Adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics. 2008;24:1403–5.

  74. 74.

    Manni F, Guerard E, Heyer E. Geographic patterns of (genetic, morphologic, linguistic) variation: how barriers can be detected by using Monmonier’s algorithm. Hum Biol. 2004;76:173–90.

  75. 75.

    Chambers JL, Garant D. Determinants of population genetic structure in eastern chipmunks (Tamias striatus): the role of landscape barriers and sex-biased dispersal. J Hered. 2010;101:413–22.

  76. 76.

    Jensen JL, Bohonak AJ, Kelley ST. Isolation by distance, web service. BMC Genet. 2005;6:13.

  77. 77.

    Ersts PJ. Geographic Distance Matrix Generator (version 1.2.3). http://biodiversityinformatics.amnh.org/open_source/gdmg/ (2017). Accessed 10 Mar 2017.

  78. 78.

    Tallmon DA, Koyuk A, Luikart G, Beaumont MA. ONeSAMP: a program to estimate effective population size using approximate Bayesian computation. Mol Ecol Resour. 2008;8:299–301.

  79. 79.

    Storz JF, Beaumont MA. Testing for genetic evidence of population expansion and contraction: an empirical analysis of microsatellite DNA variation using a hierarchical Bayesian model. Evolution. 2002;56:154–66.

  80. 80.

    Pianka ER. Evolutionary ecology. 2nd ed. New York, USA: Harper and Rowe Publishers; 1978.

  81. 81.

    Gelman A, Rubin DB. Inference from iterative simulation using multiple sequences. Stat Sci. 1992;7:457–511.

  82. 82.

    Brooks SPB, Gelman AG. General methods for monitoring convergence of iterative simulations. J Comput Graph Stat. 1998;7:434–55.

  83. 83.

    Plummer M, Best N, Cowles K, Vines K. CODA: convergence diagnosis and output analysis for MCMC. R News. 2006;6:7–11.

  84. 84.

    Gelman A, Hill J. Data analysis using regression and multilevel/hierarchical models. Policy anal. Cambridge: Cambridge University Press; 2007.

  85. 85.

    Paz-Vinas I, Quéméré E, Chikhi L, Loot G, Blanchet S. The demographic history of populations experiencing asymmetric gene flow: combining simulated and empirical data. Mol Ecol. 2013;22:3279–91.

  86. 86.

    Salmona J, Salamolard M, Fouillot D, Ghestemme T, Larose J, Centon J-F, et al. Signature of a pre-human population decline in the critically endangered Reunion Island endemic forest bird Coracina newtoni. PLoS One. 2012;7:e43524.

  87. 87.

    Hasegawa M, Kishino H, Yano T. Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. J Mol Evol. 1985;22:160–74.

  88. 88.

    Karaman MS. Studie der Astacidae (Crustacea, Decapoda) II. Teil. Hydrobiologia. 1963;22:111–32. (In German)

  89. 89.

    Economidis PS, Banarescu PM. The distribution and origins of freshwater fishes in the Balkan peninsula, especially in Greece. Int Revue Ges Hydrobiol. 1991;76:257–83.

  90. 90.

    Salvi D, Harris DJ, Kaliontzopoulou A, Carretero MA, Pinho C. Persistence across Pleistocene ice ages in Mediterranean and extra-Mediterranean refugia: phylogeographic insights from the common wall lizard. BMC Evol Biol. 2013;13:147.

  91. 91.

    Thanou E, Tryfonopoulos G, Chondropoulos B, Fraguedakis-Tsolis S. Comparative phylogeography of the five Greek vole species infers the existence of multiple South Balkan subrefugia. Ital J Zool. 2012;79:363–76.

  92. 92.

    Rohling EJ, Foster GL, Grant KM, Marino G, Roberts ΑP, Tamisiea ME, et al. Sea-level and deep-sea-temperature variability over the past 5.3 million years. Nature. 2014;508:477–82.

  93. 93.

    Maslin MA, Brierley CM. The role of orbital forcing in the early middle Pleistocene transition. Quatern Int. 2015;389:47–55.

  94. 94.

    Kroon D, Alexander I, Little M, Lourens LJ, Matthewson A, Robertson AHF, et al. Oxygen isotope and sapropel stratigraphy in the eastern Mediterranean during the last 3.2 million years. In: Robertson AHF, Emeis K-C, Richter C, Camerlenghi A, editors. Proceedings of the ocean drilling program, scientific results, 160. Earth and planetary science research institute publications; 1998. p. 181–9.

  95. 95.

    Tzedakis PC, Hooghiemstra H, Pälike H. The last 1.35 million years at Tenaghi Philippon: revised chronostratigraphy and long-term vegetation trends. Quaternary Sci Rev. 2006;25:3416–30.

  96. 96.

    Fitzsimmons KE, Marković SB, Hambach U. Pleistocene environmental dynamics recorded in the loess of the middle and lower Danube basin. Quaternary Sci Rev. 2012;41:104–18.

  97. 97.

    Caputo R, Bravard J-P, Helly B. The Pliocene-quaternary tecto-sedimentary evolution of the Larissa plain (eastern Thessaly, Greece). Geodin Acta. 1994;7:219–31.

  98. 98.

    Rook L, Martínez-Navarro B. Villafranchian: the long story of a Plio-Pleistocene European large mammal biochronologic unit. Quatern Int. 2010;219:134–44.

  99. 99.

    Bianco PG. Potential role of the palaeohistory of the Mediterranean and Paratethys basins on the early dispersal of euro-Mediterranean freshwater fishes. Ichthyol Explor Fres. 1990;1:167–84.

  100. 100.

    Římalová K, Douda K, Štambergová M. Species-specific pattern of crayfish distribution within a river network relates to habitat degradation: implications for conservation. Biodivers Conserv. 2014;23:3301–17.

  101. 101.

    Caraveli H. A comparative analysis on intensification and extensification in Mediterranean agriculture: dilemmas for LFAs policy. J Rural Stud. 2000;16:231–42.

  102. 102.

    Skoulikidis NT. The environmental state of rivers in the Balkans - a review within the DPSIR framework. Sci Total Environ. 2009;407:2501–16.

  103. 103.

    Bohl E. Motion of individual noble crayfish Astacus astacus in different biological situation: in-situ studies using radio telemetry. Freshwater Crayfish. 1999;12:677–87.

  104. 104.

    Skurdal J, Garnås E, Taugbøl T. Management strategies, yield and population development of the noble crayfish Astacus astacus in lake Steinsfjorden. Bull Fr Pêche Piscic. 2002;367:845–60.

  105. 105.

    Finné M, Holmgren K, Sundqvist HS, Weiberg E, Lindblom M. Climate in the eastern Mediterranean, and adjacent regions, during the past 6000 years - a review. J Archaeol Sci. 2011;38:3153–73.

  106. 106.

    Büntgen U, Myglan VS, Ljungqvist FC, McCormick M, Di Cosmo N, Sigl M, et al. Cooling and societal change during the late antique little ice age from 536 to around 660 AD. Nat Geosci. 2016;9:231–6.

  107. 107.

    Luterbacher J, García-Herrera R, Akcer-On S, Allan R, Alvarez-Castro MC, Benito G, et al. A review of 2000 years of paleoclimatic evidence in the mediterranean. In: Lionello P, editor. The Climate of the Mediterranean region: From the Past to the Future. Amsterdam: Elsevier; 2012. p. 87–185.

  108. 108.

    Vött A, Brückner H, Schriever A, Handl M, Besonen M, van der Borg K. Holocene coastal evolution around the ancient seaport of Oiniadai, Acheloos alluvial plain, NW Greece. Coastline Rep. 2004;1:43–53.

  109. 109.

    van Andel TH, Zangger E, Demitrack A. Land use and soil erosion in prehistoric and historical Greece. J Field Archaeol. 1990;17:379–96.

  110. 110.

    Lespez L. Geomorphic responses to long-term land use changes in eastern Macedonia (Greece). Catena. 2003;51:181–208.

  111. 111.

    Viollet P-L. Water engineering in ancient civilizations: 5,000 years of history. Boca Raton: CRC Press; 2007.

  112. 112.

    Bintliff J. The complete archaeology of Greece: from hunter-gatherers to the 20th century a.D. 1st ed. Wiley-Blackwell: Chichester; 2012.

  113. 113.

    Mylona D. Aquatic animal resources in prehistoric Aegean. Greece J Biol Res - Thessalon. 2014;21:2.

  114. 114.

    Bathrellos GD, Skilodimou HD, Maroukian H, Gaki-Papanastassiou K. Late quaternary evolution of the lower reaches of Ziliana stream in south Mt. Olympus (Greece). 8th IAG International Conference on Geomorphology “Geomorphology and Sustainability”. August 2013. Paris, France.

  115. 115.

    Association of Environmental Protection Vrontismenis. http://www.ekke.gr/estia/gr_pages/mko_po/organoseis/Grenorag/Hpiros/152.htm (2012). Accessed 4 Mar 2016. (In Greek).

  116. 116.

    Stavropoulou N, Papaconstantinos C, Bantaraki M, Kizilou C, Oikonomou M. The wetland of lake Tsivlou. Akrata: Center for Environmental Education Akrata; 2003. (In Greek)

Download references

Acknowledgements

The authors would like to acknowledge the support provided by Andreas Laggis throughout the progress of this study. Also, we would like to express our gratitude to Kostas Papathanasiou for the help provided during the sampling, I. Kappas for his valuable comments during the writing of this work and V. Sousa and T. Theodosiou for their help with R. Also, we would like to thank C. Perdikaris and several fishermen for the crayfish samples provided. Finally, we kindly acknowledge the effort and time invested by two anonymous reviewers for improving the manuscript.

Funding

Not applicable.

Availability of data and materials

The majority of the data generated or analysed during this study are included in this published article [and its supplementary information files]. The few not available data generated during the current study are available from the corresponding author on reasonable request.

Authors’ contributions

AL conceived the experimental design, sampling, performed genetic experiments, genetic and bioinformatic analysis, wrote the paper; ADB contribute in data analysis and writing of the manuscript; AC bioinformatic analysis and technical support; SM genetic analysis; AT contributed in cross-checking of genotypes and writing; TJA supervised AL Ph.D. thesis, conceived the experimental design and analytical work and contributed to the writing of the manuscript. All authors read and approved the final manuscript.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Ethics approval and consent to participate

Not applicable.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Author information

Correspondence to Theodore J. Abatzopoulos.

Additional files

Additional file 1:

List of sequences used in the mitochondrial analysis. A detailed table with the information of each sequence used is given. The information comprises the name of the site, country of origin, haplotype, GenBank Accession numbers and bibliographic references. (DOC 362 kb)

Additional file 2:

The modified conditions used for the PCR and amplification of each microsatellite loci are given. (DOC 32 kb)

Additional file 3:

Starting values (α, σ, β and τ) for hyperpriors (log(N0), log(N1), log(Θ) and log(T)) for each independent MCMC run. (DOC 30 kb)

Additional file 4:

Pairwise within and between genetic distances of the haplotype groups (G1 to G6) for COI and 16S. (DOC 56 kb)

Additional file 5:

Genetic diversity of the microsatellite loci and sampling sites. Table with the number of alleles, allelic richness, number of private alleles, expected heterozygosity, unbiased expected Nei’s heterozygosity, observed heterozygosity, inbreeding coefficient and exact P-value for Hardy-Weinberg equilibrium test are listed for the microsatellite loci and the sampling sites. (DOC 140 kb)

Additional file 6:

Fixation indices for all microsatellite loci. Values of the fixation indices for all microsatellite loci based on infinite allele model (FIT, FIS and FST) and stepwise mutation model (RST). (DOC 30 kb)

Additional file 7:

Estimated pairwise FST and RST values based on sampling sites. (DOC 53 kb)

Additional file 8:

Inference of number of clusters. The file contains graphical representations of a) DeltaK for each K (1 to 16) produced by Structure Harvester and b) Bayesian information criterion (BIC) for every number of clusters, using DAPC. (DOC 108 kb)

Additional file 9:

Assignment probability (K between 6 and 9) and assignment probability per cluster, as identified by DAPC analysis. (DOC 35 kb)

Additional file 10:

Graphical representations of the DAPC scatterplots of the first two principal components for K between 6 and 8. (DOC 325 kb)

Additional file 11:

Mantel test correlation between geographical distance (kilometers) and genetic distances (given as FST, FST/(1-FST) or RST) for all sites. (DOC 637 kb)

Additional file 12:

Graphical representation of genetic barriers (1 to 10) based on: a) sixteen sampling stations, all microsatellites loci and the genetic distance DCE, b) sixteen sampling stations, all microsatellites loci and the FST, c) sixteen sampling stations, each microsatellite locus (Aas8, Aas766, Aas1198, Aas2498, Aas3040 and Aas3950) and FST, d) nine genetic clusters (of the Structure software), all microsatellites loci and the genetic distance DCE, and e) nine genetic clusters (of the Structure software), all microsatellites loci and the FST. (DOC 10235 kb)

Additional file 13:

Map of the sampled noble crayfish in Greece (Copyright © 2014 Esri and its licensors. All rights reserved) where abbreviation of the sampling sites, names of major rivers, mountains and sierra, type of genetic markers used (mitochondrial and/or nuclear; mtDNA and nDNA, respectively) and elevation are given. (DOC 2355 kb)

Additional file 14:

Point estimates with their upper Confidence Interval (CI) for all parameters (log(N0), log(N1), log(Θ) and log(T)) and every genetic cluster for a generation time of 3.5 and 5.5. (DOC 48 kb)

Additional file 15:

Posteriors and priors densities plots for the past (N1) and present (N0) population sizes of each genetic cluster (cluster 1 to 9) and the global data set (samples from different demes were merged into one dataset; black color). (DOC 859 kb)

Additional file 16:

Posteriors and priors densities plots for time since the population decline (T) for each genetic cluster (cluster 1 to 9) and the global data set (samples from different demes were merged into one dataset; black color). (DOC 674 kb)

Additional file 17:

The relative probability of different hypotheses was assessed using approximate Bayes Factors (BF). A generation time of 5.5 for the noble crayfish was used in the analysis. Hypotheses were examined for the past 12,000 and 7000 years. (DOC 2096 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Keywords

  • Noble crayfish
  • mtDNA
  • 16S
  • COI
  • Microsatellites
  • Europe
  • Balkans
  • Phylogeny
  • Populations