Open Access

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

BMC Evolutionary BiologyBMC series – open, inclusive and trusted201717:122

https://doi.org/10.1186/s12862-017-0971-6

Received: 2 December 2016

Accepted: 17 May 2017

Published: 30 May 2017

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.

Keywords

Noble crayfish mtDNA 16S COI Microsatellites Europe Balkans Phylogeny Populations

Background

Quaternary climatic oscillations have influenced the flora (e.g. [13]) and fauna (e.g. [46]) 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 [79]. 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 [79, 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

Site

Abbreviation

Habitat type

Water Basina

Geographic coordinates

Collection method

Year

Sample size

mtDNA

nDna

Status

Arahthos

ARX

River

Arahthos

N39°15′ E21°00′

Trap (recreational fisherman)

2011

1

1

-

Unknown

Kalamas

KAS

River

Kalamas

N39°39′ E20°31′ b

Trap (from Perdikaris)

2004

20

2

20

Transfer from 1) Larissa (80’s)d and 2) unknown origins (years 1990 and 1991)e

Tzaravina

TZA

Lake

Kalamas

N39°54′ E20°30′

Hand

2011

6

2

6

Transfer from 1) Larissa (80’s)d and 2) unknown origins (years 1990 and 1991)e

Chani Kaber Aga

KPA

Stream

Arahthos

N39°43′ E20°57′

Hand

2011

6

2

6

Unknown

Aoo1

AA1

Lake

Aoos

N39°49′ E21°07′ b

Trap (fisherman)

2009

20

2

20

Unknown

Aoo2

AA2

Lake

Aoos

N39°48′ E21°06′ b

Trap (fisherman)

2011

20

2

20

Unknown

Fragkades

FRA

Stream

Aoos

N39°50′ E20°48′

Hand

2011

2

2

-

Unknown

Kalivia

KLV

Stream

Acheloos

N39°18′ E21°42′

Hand

2011, 2012

16

2

16

Unknown

Neochori

NEO

Stream

Acheloos

N39°17′ E21°43′

Hand

2011, 2012

11

2

11

Unknown

Karya

KRI

Stream

Zilianac

N39°58′ E22°25′

Hand

2011, 2012

20

2

20

Native

Skotina

SKR

Stream

Zilianac

N39°59′ E22°27′

Hand

2012

20

2

20

Native

Koniskos

KNS

Stream

Pinios

N39°48′ E21°48′

Hand

2011, 2012

17

2

17

Transfer from Krania (last 10 to 15 years)g

Loggas

LOG

Lake

Pinios

N39°49′ E21°55′

Hand

2011, 2012

17

2

17

Transfer from Krania (last 10 to 15 years)g

Begoritida/Agra

BGR

Lake

Aliakmon

N40°46′ E21°46′ b

Trap (fisherman)

2008, 2009

1

1

-

Transfer from Orhomenos and Edessa Aquaculture stationd

Palaifyto

PLF

River

Axios

N40°47′ E22°17′

Hand/Trap

2011, 2012

20

2

20

Transfer probably from Begoritida/Agrag

Tsivlo

TSV

Lake

Streams of North Peloponnese Beach

N38°04′ E22°13′

Hand

2011, 2012

20

2

20

Transfer from unknown origin f,g

Doxa

DOX

Lake

Streams of North Peloponnese Beach

N37°55′ E22°17′

Hand

2011, 2012

23

2

23

Unknown

Perivoli

PRV

Stream

Pinios

N39°03′ E22°11′

Hand

2011

1

1

-

Native

Kefalovriso

KEF

River

Pinios

N39°53′ E22°04′

Hand-land

2012

1

1

-

Transfer from unknown origing

Krania

KRN

Stream

Aliakmon

N39°51′ E21°18′

Hand

2012

20

2

20

Native

Pertouli

PRT

Stream

Acheloos

N39°32′ E21°28′

Hand

2013

22

1

22

Unknown

abased on the National Commission of Water (Government Gazette 1383/8/2-9-10 and 1572/Β/28-9-10).

bgeographic coordinates obtained from Google earth.

cfrom [114].

dfrom [36].

efrom [115].

ffrom [116].

gfrom local resident and recreational fisherman

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

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

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

Substitution rate

Decapoda

Arthropoda

Phylogroup or Node

Time

SD

95% HPD

Time

SD

95% HPD

A

1.765

0.0036

(1.1433 - 2.3995)

1.807

0.0049

(0.9810 - 2.7321)

B

1.245

0.0029

(0.7819 - 1.7700)

1.276

0.0037

(0.6481 - 1.9634)

C

1.284

0.0030

(0.7610 - 1.8444)

1.317

0.0040

(0.6953 - 2.1040)

D

0.792

0.0022

(0.4428 - 1.2068)

0.809

0.0028

(0.3496 - 1.3083)

E

0.556

0.0016

(0.3042 - 0.8348)

0.569

0.0019

(0.2598 - 0.9207)

G1

1.048

0.0025

(0.6251 - 1.4989)

1.072

0.0032

(0.5545 - 1.6842)

G2

0.705

0.0021

(0.3464 - 1.1095)

0.727

0.0027

(0.3034 - 1.2180)

G3

0.761

0.0022

(0.3957 - 1.1816)

0.781

0.0028

(0.3307 - 1.2955)

G4

0.356

0.0015

(0.1254 - 0.6256)

0.364

0.0016

(0.1157 - 0.6708)

G5

0.468

0.0014

(0.2453 - 0.7109)

0.478

0.0016

(0.2188 - 0.7847)

G6

0.285

0.0011

(0.1002 - 0.4974)

0.293

0.0013

(0.0854 - 0.5367)

For each phylogenetic group (G1 - G6) and node (A-E) (defined in Fig. 1) the mean divergence time (Time), standard deviation (SD), 95% HPD lower and upper are given

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

 

AA1

AA2

KAS

LOG

KNS

KLV

NEO

PLF

KRI

SKR

KRN

KPA

TZA

DOX

TSV

PRT

A

3.667

4.333

8.167

2.500

2.500

4.167

2.833

5.333

2.667

2.167

2.333

2.333

2.500

8.167

7.333

7.333

AR

2.844

3.295

5.567

1.969

2.069

2.760

2.281

4.082

2.021

2.056

1.856

2.333

4.833

5.333

5.019

5.080

PA

0.000

3.000

5.000

1.000

0.000

1.000

1.000

3.000

1.000

0.000

1.000

0.000

1.000

8.000

8.000

4.000

HE

0.468

0.559

0.753

0.204

0.236

0.365

0.299

0.648

0.310

0.310

0.174

0.414

0.641

0.758

0.735

0.761

HN

0.480

0.573

0.772

0.210

0.243

0.377

0.314

0.664

0.318

0.318

0.179

0.452

0.700

0.775

0.754

0.779

HO

0.408

0.592

0.767

0.265

0.216

0.333

0.356

0.592

0.308

0.300

0.192

0.472

0.667

0.768

0.717

0.727

FIS

0.152

−0.033

0.007

−0.271

0.114

0.118

−0.144

0.112

0.030

0.058

−0.074

−0.049

0.051

0.009

0.050

0.067

PH-W

0.040

0.103

0.402

0.840

0.344

0.169

0.943

0.047

0.804

0.825

0.850

0.689

0.769

0.069

0.805

0.127

L0

Aas2489

Aas8

   

Aas8

       

Aas8

 

Aas1198, Aas8

Where the number of alleles (A), allelic richness (AR), number of private alleles (PA), expected heterozygosity (HE), unbiased expected Nei’s heterozygosity (HN), observed heterozygosity (HO), inbreeding coefficient (FIS), exact P-value for Hardy-Weinberg equilibrium test (PH-W) and loci with putative null alleles (L0) are given

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

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

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

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

Genetic cluster

1

2

3

4

5

6

7

8

9

Generation time

3.5

5.5

3.5

5.5

3.5

5.5

3.5

5.5

3.5

5.5

3.5

5.5

3.5

5.5

3.5

5.5

3.5

5.5

Mean log10(N0)

0.885

1.006

0.522

0.661

0.889

0.781

1.289

1.425

1.298

1.273

0.837

0.582

1.315

1.233

1.411

1.201

0.509

0.501

S.D.

0.711

0.685

0.952

0.826

0.711

0.685

0.827

0.873

0.647

0.519

0.622

0.632

0.856

0.877

0.581

0.655

0.555

0.664

Mean log10 (N1)

5.153

5.132

5.355

5.349

4.970

4.971

5.381

5.387

5.502

5.508

5.036

5.272

5.342

5.338

5.462

5.464

5.281

5.269

S.D.

0.382

0.380

0.490

0.495

0.378

0.375

0.324

0.322

0.308

0.308

0.411

0.482

0.333

0.332

0.304

0.303

0.482

0.478

Mean log10 (T)

2.286

2.584

2.454

2.770

2.167

2.267

2.358

2.680

2.441

2.611

2.218

2.732

2.430

2.545

2.430

2.432

2.482

2.657

S.D.

0.654

0.631

0.893

0.784

0.755

0.797

0.597

0.489

0.574

0.587

0.786

0.822

0.539

0.604

0.514

0.619

0.826

0.842

Standard Deviation (SD) for every parameter is given

Table 5

Statistics describing the magnitude of the demographic change for each genetic cluster (1-9) and generation time (3.5 or 5.5)

Genetic cluster

1

2

3

4

5

6

7

8

9

Generation time

3.5

5.5

3.5

5.5

3.5

5.5

3.5

5.5

3.5

5.5

3.5

5.5

3.5

5.5

3.5

5.5

3.5

5.5

Mean effective size log10(N0/N1)

-1.809

−1.734

−1.976

−1.96

−1.718

−1.764

−1.515

−1.414

−1.576

−1.555

−1.761

−1.993

−1.523

−1.576

−1.449

−1.643

−2.037

-2.030

95% CI lower

−1.809

−1.734

−1.976

−1.96

−1.718

−1.764

−1.515

−1.414

−1.576

−1.555

−1.761

−1.993

−1.523

−1.576

−1.449

−1.643

−2.037

-2.030

95% CI upper

−1.809

−1.734

−1.976

−1.96

−1.718

−1.764

−1.515

−1.414

−1.576

−1.555

−1.761

−1.993

−1.523

−1.576

−1.449

−1.643

−2.037

-2.030

S.D.

0.826

0.804

0.916

0.894

0.838

0.888

0.623

0.468

0.662

0.596

0.863

0.902

0.587

0.686

0.513

0.714

0.930

0.920

Variance

0.682

0.646

0.840

0.800

0.702

0.788

0.388

0.219

0.438

0.356

0.745

0.814

0.344

0.471

0.263

0.51

0.865

0.847

Hedges’ d

−7.479

−7.447

−6.386

−6.885

−7.164

−7.584

−6.516

−6.024

−8.301

−9.915

−7.965

−8.346

−6.199

−6.19

−8.734

−8.347

−9.183

-8.241

95% CI lower

−9.439

−9.407

−8.346

−8.846

−9.124

−9.544

−8.476

−7.984

−10.261

−11.875

−9.925

−10.306

−8.159

−8.15

−10.694

−10.307

−11.143

-10.201

95% CI upper

−5.519

−5.487

−4.426

−4.925

−5.204

−5.623

−4.556

−4.063

−6.341

−7.955

−6.005

−6.386

−4.239

−4.23

−6.774

−6.387

−7.222

−6.281

Mean effective size log10(N0/N1), Standard Deviation (SD), Variance, along with Hedges’ d and their 95% Confidence Intervals (CI) are given. Negative values indicate significant bottlenecks

Table 6

Mean estimated contemporary size (N0) for each genetic cluster of the noble crayfish, along with their 95% Confidence Intervals (CI) (OneSAMP software)

Cluster

Mean N0

95% CI lower

95% CI upper

1

29.134

20.632

44.899

2

22.069

14.988

37.702

3

23.765

18.996

33.244

4

34.107

27.110

48.675

5

38.318

27.868

61.107

6

34.361

24.070

59.253

7

27.029

22.345

35.878

8

44.191

33.663

65.395

9

23.051

13.907

39.118

Fig. 6

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

Declarations

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.

Open AccessThis 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.

Authors’ Affiliations

(1)
Department of Genetics, Development and Molecular Biology, School of Biology, Aristotle University of Thessaloniki
(2)
Scientific Computing Office, Information Technology (IT) Center, School of Sciences, Aristotle University of Thessaloniki

References

  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.View ArticlePubMedGoogle Scholar
  2. Médail F, Diadema K. Glacial refugia influence plant diversity patterns in the Mediterranean Basin. J Biogeogr. 2009;36:1333–45.View ArticleGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  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.View ArticlePubMedPubMed CentralGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  7. Hewitt GM. Post-glacial re-colonization of European biota. Biol J Linn Soc. 1999;68:87–112.View ArticleGoogle Scholar
  8. Hewitt GM. Genetic consequences of climatic oscillations in the quaternary. Philos T Roy Soc B. 2004;359:183–95.View ArticleGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  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.View ArticleGoogle Scholar
  11. Schmitt T. Molecular biogeography of Europe: Pleistocene cycles and postglacial trends. Front Zool. 2007;4:11.View ArticlePubMedPubMed CentralGoogle Scholar
  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.Google Scholar
  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.View ArticleGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  15. Tzedakis PC. Museums and cradles of Mediterranean biodiversity. J Biogeogr. 2009;36:1033–4.View ArticleGoogle Scholar
  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.View ArticleGoogle Scholar
  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.View ArticleGoogle Scholar
  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.View ArticleGoogle Scholar
  19. Skurdal J, Taugbøl T. Astacus. In: Holdich DM, editor. Biology of freshwater crayfish. Oxford: Blackwell Science; 2002. p. 467–510.Google Scholar
  20. Cukerzis J. La biologie de l’écrevisse (Astacus astacus L.). Paris: Institut National de la Recherche Agronomique Publications; 1984. (In French)Google Scholar
  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.View ArticleGoogle Scholar
  22. Lundberg U. Behavioural elements of the noble crayfish, Astacus astacus (Linnaeus, 1758). Crustaceana. 2004;77:137–62.View ArticleGoogle Scholar
  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. Keller M. Ten years of trapping Astacus astacus for restocking. Freshwater Crayfish. 1998;12:518–28.Google Scholar
  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.View ArticleGoogle Scholar
  26. Reynolds J, Souty-Grosset C, Richardson A. Ecological roles of crayfish in freshwater and terrestrial habitats. Freshwater Crayfish. 2013;19:197–218.Google Scholar
  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).Google Scholar
  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.Google Scholar
  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.View ArticleGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  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.Google Scholar
  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.View ArticleGoogle Scholar
  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.View ArticleGoogle Scholar
  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.View ArticleGoogle Scholar
  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.View ArticleGoogle Scholar
  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.View ArticleGoogle Scholar
  37. Zhang D-X, Hewitt GM. Nuclear DNA analyses in genetic studies of populations: practice, problems and prospects. Mol Ecol. 2003;12:563–84.View ArticlePubMedGoogle Scholar
  38. Westman K, Pursiainen M, Vilkman R. A new folding trap model which prevents crayfish from escaping. Freshwater Crayfish. 1978;4:235–42.Google Scholar
  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.Google Scholar
  40. Crandall KA, Fitzpatrick JF. Crayfish molecular systematics: using a combination of procedures to estimate phylogeny. Syst Biol. 1996;45:1–26.View ArticleGoogle Scholar
  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.PubMedGoogle Scholar
  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.Google Scholar
  43. Baxevanis AD, Kappas I, Abatzopoulos TJ. Molecular phylogenetics and asexuality in the brine shrimp Artemia. Mol Phylogenet Evol. 2006;40:724–38.View ArticlePubMedGoogle Scholar
  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.Google Scholar
  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.View ArticlePubMedGoogle Scholar
  46. Buhay JE. ‘COI-like’ sequences are becoming problematic in molecular systematic and DNA barcoding studies. J Crustac Biol. 2009;29:96–110.View ArticleGoogle Scholar
  47. Xia X. DAMBE5: a comprehensive software package for data analysis in molecular biology and evolution. Mol Biol Evol. 2013;30:1720–8.View ArticlePubMedPubMed CentralGoogle Scholar
  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.Google Scholar
  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.View ArticlePubMedGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  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. 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.View ArticlePubMedGoogle Scholar
  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.View ArticlePubMedPubMed CentralGoogle Scholar
  54. Rambaut A, Drummond AJ. Tracer version 1.5. http://beast.bio.ed.ac.uk/Tracer (2009). Accessed 21 Mar 2015.
  55. R Core Team. R: A language and environment for statistical computing. http://www.r-project.org/ (2014). Accessed 25 Sep 2015.
  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.View ArticlePubMedPubMed CentralGoogle Scholar
  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.View ArticlePubMedPubMed CentralGoogle Scholar
  58. Schubart CD, Diesel R, Hedges SB. Rapid evolution to terrestrial life in Jamaican crabs. Nature. 1998;393:363–5.View ArticleGoogle Scholar
  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.View ArticleGoogle Scholar
  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.Google Scholar
  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.View ArticleGoogle Scholar
  62. Nei M. Estimation of average heterozygosity and genetic distance from a small number of individuals. Genetics. 1978;89:583–90.PubMedPubMed CentralGoogle Scholar
  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. Rousset F. Genepop’007: a complete re-implementation of the genepop software for windows and Linux. Mol Ecol Resour. 2008;8:103–6.View ArticlePubMedGoogle Scholar
  65. Rice WR. Analyzing tables of statistical tests. Evolution. 1989;43:223–5.View ArticleGoogle Scholar
  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. 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.View ArticleGoogle Scholar
  68. Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155:945–59.PubMedPubMed CentralGoogle Scholar
  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.View ArticleGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  71. Rosenberg NA. DISTRUCT: a program for the graphical display of population structure. Mol Ecol Notes. 2004;4:137–8.View ArticleGoogle Scholar
  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.View ArticlePubMedPubMed CentralGoogle Scholar
  73. Jombart T. Adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics. 2008;24:1403–5.View ArticlePubMedGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  76. Jensen JL, Bohonak AJ, Kelley ST. Isolation by distance, web service. BMC Genet. 2005;6:13.View ArticlePubMedPubMed CentralGoogle Scholar
  77. Ersts PJ. Geographic Distance Matrix Generator (version 1.2.3). http://biodiversityinformatics.amnh.org/open_source/gdmg/ (2017). Accessed 10 Mar 2017.
  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.View ArticlePubMedGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  80. Pianka ER. Evolutionary ecology. 2nd ed. New York, USA: Harper and Rowe Publishers; 1978.Google Scholar
  81. Gelman A, Rubin DB. Inference from iterative simulation using multiple sequences. Stat Sci. 1992;7:457–511.View ArticleGoogle Scholar
  82. Brooks SPB, Gelman AG. General methods for monitoring convergence of iterative simulations. J Comput Graph Stat. 1998;7:434–55.Google Scholar
  83. Plummer M, Best N, Cowles K, Vines K. CODA: convergence diagnosis and output analysis for MCMC. R News. 2006;6:7–11.Google Scholar
  84. Gelman A, Hill J. Data analysis using regression and multilevel/hierarchical models. Policy anal. Cambridge: Cambridge University Press; 2007.Google Scholar
  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.View ArticlePubMedGoogle Scholar
  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.View ArticlePubMedPubMed CentralGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  88. Karaman MS. Studie der Astacidae (Crustacea, Decapoda) II. Teil. Hydrobiologia. 1963;22:111–32. (In German)View ArticleGoogle Scholar
  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.View ArticleGoogle Scholar
  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.View ArticlePubMedPubMed CentralGoogle Scholar
  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.View ArticleGoogle Scholar
  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.View ArticlePubMedGoogle Scholar
  93. Maslin MA, Brierley CM. The role of orbital forcing in the early middle Pleistocene transition. Quatern Int. 2015;389:47–55.View ArticleGoogle Scholar
  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.Google Scholar
  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.View ArticleGoogle Scholar
  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.View ArticleGoogle Scholar
  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.View ArticleGoogle Scholar
  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.View ArticleGoogle Scholar
  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.Google Scholar
  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.View ArticleGoogle Scholar
  101. Caraveli H. A comparative analysis on intensification and extensification in Mediterranean agriculture: dilemmas for LFAs policy. J Rural Stud. 2000;16:231–42.View ArticleGoogle Scholar
  102. Skoulikidis NT. The environmental state of rivers in the Balkans - a review within the DPSIR framework. Sci Total Environ. 2009;407:2501–16.View ArticlePubMedGoogle Scholar
  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.Google Scholar
  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.View ArticleGoogle Scholar
  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.View ArticleGoogle Scholar
  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.View ArticleGoogle Scholar
  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.Google Scholar
  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.Google Scholar
  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.Google Scholar
  110. Lespez L. Geomorphic responses to long-term land use changes in eastern Macedonia (Greece). Catena. 2003;51:181–208.View ArticleGoogle Scholar
  111. Viollet P-L. Water engineering in ancient civilizations: 5,000 years of history. Boca Raton: CRC Press; 2007.Google Scholar
  112. Bintliff J. The complete archaeology of Greece: from hunter-gatherers to the 20th century a.D. 1st ed. Wiley-Blackwell: Chichester; 2012.View ArticleGoogle Scholar
  113. Mylona D. Aquatic animal resources in prehistoric Aegean. Greece J Biol Res - Thessalon. 2014;21:2.View ArticlePubMedGoogle Scholar
  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.Google Scholar
  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. Stavropoulou N, Papaconstantinos C, Bantaraki M, Kizilou C, Oikonomou M. The wetland of lake Tsivlou. Akrata: Center for Environmental Education Akrata; 2003. (In Greek)Google Scholar

Copyright

© The Author(s). 2017

Advertisement