Genetic diversity and structure of the noble crayfish populations in the Balkan Peninsula revealed by mitochondrial and microsatellite DNA markers

Sample collection and DNA isolation

This study used more than 400 noble crayfish samples from 18 populations in the western part of the Balkans (Fig. 1, Table 1). The specimen collection was conducted in Croatia, Romania, Serbia and Slovenia with the approval of local authorities (Croatia UP/I-612-07/18-48/148; Romania 408/CJ/27.11.2018; Serbia 324-04-10/2021-04 and 03 No. 026-419/2; Slovenia 35601-1262 150/2006-6 and 35601-135/2010-9). Crayfish were collected by hand or with trapped baited traps (Westman, Pursiainen & Vilkman, 1978; Kozák et al., 2015). Crayfish samples from Albania were purchased from a fisherman on the Prespa Lake. The Prespa Lake and Ohrid Lake are approximately 10 km apart, and waters from the Prespa Lake feed Ohrid Lake through underground karst channels. Out of the Ohrid Lake springs the only outlet, the Black Drim River, which flows in a north direction into Albania and thus into the Adriatic Sea. A pereiopod from each specimen was removed and preserved in 96% ethanol. Wild caught crayfish were then released back into the water. This sampling method does not harm crayfish as these appendages regenerate after the next moulting. Total genomic DNA was extracted from the pereiopod muscle tissue with the Sigma GenElute Mammalian Genomic DNA Miniprep Kit (Sigma-Aldrich, USA) following the manufacturer’s protocol and stored at −20 °C.

Mitochondrial DNA analyses

Mitochondrial 16S and COI gene fragments were amplified and sequenced with primers 16Sar/16Sbr (Simon et al., 1994) and LCO-1490/HCO-2198 (Folmer et al., 1994) allowing comparisons to be made with 16S and COI sequences from previous studies on this species (Schrimpf et al., 2014; Laggis et al., 2017; Mrugała et al., 2017) available from GenBank. Polymerase chain reactions (PCR) for COI were prepared in a total volume of 25 µl containing 10–50 ng/µl DNA template, 1.5 mM Promega Buffer, 0.04 U HotStart Polymerase, 0.15 mM of each dNTP, 0.7 mM MgCl₂ and 0.4 µM of each primer. PCR conditions for COI were as follows: initial denaturation at 94 °C for 3 min, followed by 35 cycles of denaturation at 94 °C for 45 s, annealing at 48 °C for 60 s and extension at 72 °C for 60 s, and the final extension of 10 min at 72 °C. The final reaction mix in a total volume of 10 µL for 16S gene contained 0.05 U GoTaq G2 HotStart Polymerase, 1.5 mM GoTaq FlexiBuffer, 0.2 mM of each dNTP, 0.275 µM of each primer, and 10–50 ng/µl of DNA template. PCR conditions for 16S were as follows: initial denaturation at 95 °C for 3 min, followed by 40 cycles of denaturation at 95 °C for 1 min, annealing at 52 °C for 1 min and extension at 72 °C for 1 min, and the final extension of 5 min at 72 °C. The purification of PCR products was performed with EXOAnP Mix (20 U/µl of Exonuclease I (New England Biolabs), 5 U/µl of Antarctic phosphatase (New England Biolabs). The sequencing of purified PCR products was prepared by Macrogen, Inc. (Seoul, South Korea), with the same forward primers used for the gene amplifications.

Sequences were edited in SEQUENCHER v. 5.3 (Gene Codes Corporation, Ann Arbor, USA) and aligned using MAFFT v.7.187 (Katoh & Standley, 2013). The final alignment for the COI gene fragment was 655 bp long, while for 16S was 475 bp long. New 16S and COI sequences were submitted to GenBank and BOLD data bases, and their GenBank accession numbers are MW726211 –MW726336 and MW726338 –MW726635 for 16S and COI sequences, respectively (Table S1 in Supplements). Additionally, all available sequences of 16S and COI genes of A. astacus were downloaded from GenBank (COI sequences from Schrimpf et al. (2014) and Laggis et al. (2017) were 350 bp-long, while sequences from Mrugała et al. (2017) were 635 bp long, and 16S sequences were 475 bp long). GenBank accession numbers of the haplotypes obtained in the present study, as well as the ones from Schrimpf et al. (2014), Laggis et al. (2017) and Mrugała et al. (2017) are reported in Tables S1 and S2.

In order to reconstruct phylogenetic tree that will be comparable with trees obtained in previous studies (Schrimpf et al., 2014; Laggis et al., 2017) we concatenated COI and 16S sequences, and used only those samples for which sequences of both genes were available. This made concatenated data set (used in phylogenetic reconstruction) smaller compared to the COI data set that was used for other analyses. The 350 bp long COI sequences and 475 bp 16S sequences from the same individual were concatenated and collapsed to unique haplotypes using FaBox (Villesen, 2007). The full concatenated COI/16S data set included 83 haplotype combinations (Table S2) and the final alignment was 825 bp long including a single gap-containing position observed in 16S fragment. The optimal models of nucleotide evolution for each partition of the concatenated data set were selected under the Bayesian information criterion (BIC) using the jModelTest (Darriba et al., 2012). The selected model for 16S was HKY+I, while for COI was HKY+G. Phylogenetic tree was reconstructed using the concatenated haplotypes in BEAST v.2.5.2 (Bouckaert et al., 2019). Since the null hypothesis of equal evolutionary rate throughout the tree was rejected at a 5% significance level, we used a relaxed uncorrelated lognormal clock model and the arthropod substitution rate of 2.3% pairwise sequence divergence for COI (0.0115 substitutions/s/Ma/l) (Brower, 1994) along with an estimated molecular clock for the 16S. The tree prior was set as the birth-death and independent substitution models were assigned to each mtDNA gene. The Markov Chain Monte Carlo (MCMC) analysis run comprised 300,000,000 generations, sampled every 10,000 generations. In order to determine convergence, the Effective Sample Size (ESS) values were checked in Tracer (Rambaut et al., 2018). The best fit tree was produced using the Maximum clade credibility tree option in TreeAnnotator 2.5.2 after the 20% of the sampled trees was discarded as burn-in.

Additionally, the phylogenetic relationships were estimated using Bayesian analysis (BA) in MrBayes ver.3.2. (Huelsenbeck & Ronquist, 2001; Ronquist & Huelsenbeck, 2003) with priors set according to the suggested model for each partition (16S: HKY+I, COI: HKY+G). Two separate runs with four Metropolis-coupled Monte Carlo Markov chains (MMCM) were performed for 10,000,000 generations, and trees were sampled every 1,000 generations. After effective sample size (ESS values > 200) for each parameter was confirmed with Tracer, the first 25% of sampled trees were eliminated as burn-in, and a 50% majority-rule consensus tree was constructed, with nodal values representing the posterior probabilities. Sequences of Pontastacus leptodactylus (Acc. No. KX279350) were used as an outgroup.

The median-joining network (MJ) approach (Bandelt, Forster & Röhl, 1999) was used in order to establish non-hierarchical phylogeographic and phylogenetic relationships among samples. To that end, three data sets were prepared; two including only COI sequences, and one including concatenated (16S+COI) sequences. The COI sequence data set I comprised 655 bp sequences obtained in this study, while data set II comprised 350 bp long sequences published in Schrimpf et al. (2011); Schrimpf et al. (2014)) and Laggis et al. (2017), combined with sequences obtained in this study and the study by Mrugała et al. (2017). The sequences obtained in this study and the study by Mrugała et al. (2017) were trimmed to 350 bp in order to match size of data sets from Schrimpf et al. (2011); Schrimpf et al. (2014)) and Laggis et al. (2017). This approach enabled us to associate COI haplotypes obtained in the present study to the COI haplotypes obtained in their research and indirectly to the lineages sensu Schrimpf et al. (2014) and groups sensu Laggis et al. (2017). Since one of our aims was to position our new samples into previously constructed phylogenies (Schrimpf et al., 2014; Laggis et al., 2017), we adopted their naming of lineages or/and groups. The MJ networks were generated using the program PopART v.1.7 (Leigh & Bryant, 2015) with all parameters set to default values.

Genetic diversity indices (number of segregating sites, number of haplotypes, haplotype diversity, nucleotide diversity, average number of nucleotide differences) were calculated in program DNASP v.6 (Rozas et al., 2017).

Analysis of molecular variance (AMOVA) (Excoffier, Smouse & Quattro, 1992) was performed in Arlequin v.3.5 (Excoffier & Lischer, 2010) in order to estimate hierarchical distribution of genetic diversity of A. astacus. Populations were grouped on the basis of major tributaries of the Danube River (Drava, Sava, Vuka, Velika Morava, Timok, Tisa) and Prespa Lake. Standard AMOVA computations were performed with three hierarchical levels: among groups (river basins/tributaries), among populations within groups, and within populations. The variance components were tested statistically by non-parametric randomisation tests using 10,000 permutations. Genetic differentiation among populations and river catchments was analysed through estimation of pairwise values of Φ_ST.

Microsatellite DNA analyses

A total of 19 species-specific tetranucleotide repeat microsatellite loci were amplified following the protocols and procedures described by Gross et al. (2017) with the following modifications; initial screening of a few Balkan A. astacus populations revealed that they possess much higher variability and wider allele size range at many loci than the eastern European (Czech Republic and Estonia) populations used in the study of Gross et al. (2017). Therefore, we split the single 19-plex microsatellite panel into two multiplex (10-plex and 9-plex) panels (Table S3). Only 15 loci were used for data analysis, as it later became apparent that at four loci, the allele size ranges still overlapped in some populations (Aast4_26, Aast4_47, Aast4_10 and Aast4_30) (Table S3). The PCR reaction (10 μl) contained 1x Type-it Multiplex PCR Master Mix (QIAGEN, Germany), 200 to 400 nM of each primer (Table S3 ), and ca 5 ng/ μl of DNA template. Touchdown program was used for PCR amplification: initial activation of 5 min at 95 °C, followed by 20 cycles of 30 s at 95 °C, 90 s at 60 °C, 30 s at 72 °C, with the annealing temperature decreasing 0.5 °C per cycle, followed by 10 cycles of 30 s at 95 °C, 90 s at 50 °C, 30 s at 72 °C and a final extension for 30 min at 60 °C. The multiplex panels of microsatellite loci were genotyped on Applied Biosystems 3500 Genetic Analyser (Life Technologies, USA) using internal GeneScan 600 LIZ Size Standard v2.0 (Life Technologies, USA) and microsatellite genotypes were scored using GeneMapper v.5 software (Life Technologies, USA). Microsatellite genotypes were initially scored automatically and were double-checked manually by two experts. Subsample of 50 individuals (12% of total 413 analysed individuals) was genotyped twice to check the genotyping consistency.

MICRO-CHECKER 2.2.3 (Van Oosterhout et al., 2004) was used to assess the potential presence of genotyping errors due to scoring of stutter peaks, large allele dropouts and null alleles. FSTAT v. 2.9.3.2 programme package (Goudet, 2001) was used for calculating allele frequencies, F_IS and pair-wise F_ST values (as estimated by Weir and Cockerham’s θ), for estimating expected (unbiased genetic diversity) and observed heterozygosities (H_E, H_O) and rarefied allelic richness (A_R), and for testing the significance of differences in average values of A_R, H_E and H_O among groups of populations (1,000 permutations, two-sided tests). The private allelic richness (A_PrRar) was estimated by rarefaction method using HP Rare v.Feb-2-2009 (Kalinowski, 2005) and multiplied by number of loci in order to make it comparable to the number of observed private alleles (A_Pr). GENEPOP v. 3.3 (Raymond & Rousset, 1995a) was used to test genotypic distributions for conformance to Hardy-Weinberg (HW) expectations and to test the loci for genotypic disequilibria. All probability tests were based on the Markov chain method (Guo & Thompson, 1992; Raymond & Rousset, 1995b) using 1,000 de-memorization steps 100 batches and 1,000 iterations per batch. Sequential Bonferroni adjustments (Rice, 1989) were applied to correct for the effect of multiple tests. Analysis of molecular variance (AMOVA) incorporated in ARLEQUIN v. 3.5.1.2 (Excoffier & Lischer, 2010) was used to partition genetic variance hierarchically between population groups, between populations within groups and among individuals within the populations. Populations were grouped based on major tributaries of the Danube River (Drava, Sava, Tisa, Velika Morava, Timok and Vuka) and the Prespa Lake.

Principal Coordinates Analysis (PCoA), implemented in PAST v. 4.05 (Hammer, Harper & Ryan, 2001), was used to explore and to visualize D_A genetic distance (Nei, Tajima & Tateno, 1983) matrix between populations in multidimensional space. The D_A distances were calculated using POPULATIONS v1.2.3.1 software (Langella, 1999).

Long-term effective population size (N_e) was assessed from the microsatellite data using the coalescent-based approach implemented in Lamarc v.2.1 (Kuhner, 2006), considering a mutation rate of 5 × 10⁻⁴ and a mixed mutational model with 30% KAM as in Gouin et al. (2011). We ran three replicates of two initial chains retaining 10,000 genealogies sampled every 200 and discarding the first 1,000 as burn-in, followed by one final chain retaining 30,000 genealogies sampled every 200 and discarding the first 5,000 as burn-in.

Combined mtDNA and microsatellite analyses

In order to visualise the distribution of A. astacus genetic diversity across geographic space and connect genetic diversity patterns to geographic features we applied Alleles in Space (AIS) package (Miller, 2005), with software implementation of Mantel test and interpolation of genetic landscape shape (GLS). The Mantel test was used to test for correlations between genetic and geographic distances using three datasets: (a) COI data set including all sequences (350 bp long), (b) COI data set including only sequences obtained in this study (655 bp long), (c) genotypes of 15 microsatellite loci. Mantel tests were performed at the individual level using an analogue of Nei’s genetic distances (Nei, Tajima & Tateno, 1983) between pair of individuals. For the Mantel test and other analyses run in AIS package, significance was tested using 10,000 permutations.

GLS interpolation was performed using two datasets: (a) full COI data set including all sequences (350 bp long) aiming to get insight into diversity on the European level, and (b) genotypes from 15 microsatellite loci in order to get insight into local diversity on the Balkan Peninsula. First, sampling sites were connected through a network based on the Delaunay triangulation, in which the simple mismatch molecular distances between connected sampling sites were calculated based on the molecular data obtained from all individuals. Surface calculation was based on midpoints of edges derived from Delaunay triangulation. Surface heights were calculated based on residual genetic distances.

The values of molecular distance were set in the mid-points of each connection in the network using the Alleles in Space (AIS) software (Miller, 2005). The raw molecular distances were interpolated. The matrix of the ‘elevation’ values was imported into QGIS 3.10 software (available at: https://qgis.org/en/site/) to generate molecular divergence surface image using the inverse distance weighted (IDW) algorithm, plotted over a map of Europe.

Phylogenetic relationships among populations using mitochondrial COI and 16S

The phylogenetic tree inferred by BEAST using concatenated data set indicated six mostly unsupported previously described genetic lineages sensu Schrimpf et al. (2014) and groups sensu Laggis et al. (2017) (Fig. 2A), with all newly obtained concatenated haplotypes (Hap49, Hap50, Hap51, Hap52, Hap53, Hap54) nested within them. Most of the lineages diversified during Pleistocene, within the period between 1.7 and 0.5 mya (Fig. 2A). Phylogenetic reconstruction using BA in MrBayes revealed unresolved relationships among lineages/groups and weak nodal support with numerous polytomies (Fig. S1).

New haplotypes from Serbian populations (Hap49, Hap52 and Hap53) grouped with haplotypes from Kosovo, Montenegro and Germany (Fig. 2A). Haplotype 51 from Slovenia and haplotype 54 from Montenegro recovered within the same unsupported clade together with haplotypes from Croatia, Austria and Czechia (Fig. 2A). The haplotype representing specimens from the Lake Prespa in Albania (Hap50) was positioned within the clade encompassing haplotypes from Greece (Fig. 2A).

The median-joining (MJ) networks for concatenated data set (Fig. 2B), as well as for two COI data sets (Figs. 3 and 4), were mostly congruent and they depicted haplotype relatedness and distribution within A. astacus. The MJ network based on long COI sequences (655 bp) revealed the existence of 11 unique haplotypes, separated by different number of mutational steps (Fig. 3). The analysis of short-sequences data set (350 bp) produced MJ network containing 60 haplotypes, with 12 established in the present study, six of them (ssh1, ssh4, ssh5, ssh6, ssh10, ssh11) obtained for the first time (Fig. 4). The MJ network based on the concatenated data set encompassed 83 haplotypes, six of them (Hap49-Hap54) obtained for the first time (Fig. 2B).

In the long-sequences MJ network, haplotype Lshm9 (detected in the Slovenian populations BLO and KOC (Sava River, tributary of the Danube River)) and haplotype Lshm10 (discovered in Croatia, in the TOT and MAK populations (Drava and Sava, tributaries of the Danube River, respectively)) were separated by five mutational steps, whereas 10 bp changes were observed between Lshm9 and Lshm3, and 11 between Lshm9 and Lshm6 (Fig. 3). In the short-sequences MJ those two haplotypes (ssh10 and ssh11, respectively) were separated by one and three bp changes, respectively, from the closest haplotype that was recorded in Austria and Czech Republic (Fig. 4). The same was observed from the concatenated MJ; Hap51 (COI haplotype ssh10) and Hap25 (COI haplotype Aas18) were separated by one base pair change (Fig. 2B). Moreover, when BEAST tree was reconstructed with concatenated sequences the same grouping was observed; Hap51 and Hap25 form well supported subclade within unsupported lineage 4 sensu Schrimpf et al. (2014) (Fig. 2A), while this subclade was less supported in the Bayesian phylogram (Fig. S1). The Lsmh2/ssh2 COI haplotype was recorded in Croatian populations from the Drava River system (MOT and JAN) as well as in the Romanian populations (BEZ and SOM) from the Tisa River system where also Lshm4/ssh5 was recorded (Figs. 3 and 4). The former (ssh2 haplotype) is widely distributed in Europe (Fig. 4). Interestingly, in the populations MOT and JAN another COI haplotype was recorded (Lshm6/ssh7, Figs. 3 and 4) indicating presence of crayfish belonging to two distinct lineages/groups in the Drava River system. In concatenated data set, in both MJ and BEAST tree, Hap42 (COI haplotype ssh7) was closely related to Hap41 (COI haplotype Aas26) and Hap54 (COI haplotype ssh6) distributed also in Germany and Montenegro, respectively (Fig. 2). Some of the A. astacus specimens from Serbian waterbodies had ssh4 (in concatenated data set Hap53) haplotype, while most of them possessed the Lshm3/ssh3 haplotype (Figs. 3 and 4). Both haplotypes were geographically restricted to the Danube River basin tributaries Velika Morava and Timok. In addition, in concatenated data set the later formed Hap26 that is positioned close to haplotypes distributed also in Romania and Kosovo (Fig. 2). The remaining specimens from Serbia were collected from the Sava River system (KAC) and they carried haplotype Lshm1/ssh1 (Figs. 3 and 4) that in concatenated data set forms Hap49, positioned closely to haplotype Hap28 from Croatia (Fig. 2). The Lshm7/ssh8 haplotype was established in the Vuka population, and in concatenated data set (Hap43) analyses is positioned closely to haplotypes recovered from the Save River in Croatia, but also haplotypes found in Germany (Figs. 2, 3 and 4). Finally, haplotype Lshm11/ssh12 was found in Croatian populations, and in concatenated data set it formed Hap28 that was positioned central to haplotypes from Danube’s tributaries in Croatia, Serbia and Romania, but also from the Rhine River system (Figs. 2, 3 and 4).

Summarised results of genetic diversity indices of 655 bp long COI sequences revealed that populations from the Danube’s tributaries Sava, Drava and Tisa rivers possessed higher number of haplotypes and nucleotide diversity compared to other river systems that were characterised by only one haplotype (Table 2).

Genetic differentiation of populations based on pairwise comparison of COI sequences revealed very high Φ_ST-values indicating genetically isolated populations with limited gene flow (Table S4). However, some of the population pairs’ Φ_ST-values were extremely low and not statistically significant. If populations are grouped according to the Danube River tributaries (Sava, Drava, Velika Morava, Vuka, Timok, Tisa) and the Prespa Lake, values of genetic differentiation varied from 0.193 (Drava-Sava) to 1.00 (Prespa Lake vs Vuka, Velika Morava and Timok), again demonstrating genetically different groups (data not shown).

Results of AMOVA conducted on 655 bp long COI sequences grouped according to specimen’s affiliation to major tributaries of the Danube River (Drava, Sava, Vuka, Velika Morava, Timok, Tisa) and Prespa Lake, revealed that most of the variance is contained among population within groups (61.05%, P < 0.001; Table 3), followed by variance among groups (31.1%), and a small amount of the variance was found within populations (8.56%).

Microsatellites analyses

Genetic differentiation and population structure

The overall level of genetic differentiation between all studied samples was high (global F_ST = 0.501) with pairwise estimates of F_ST ranging from 0.211 (between JAR and MOT populations in Croatia) to 0.838 (between Serbian GRL and Romanian BEZ populations) (Table S4). The average level of differentiation among populations from different major tributaries of the Danube (Sava, Drava, Velika Morava, Tisa) was relatively high (average pair-wise F_ST ranging from 0.360 between Sava and Drava to 0.623 between Velika Morava and Tisa) but lower than their differentiation from the PRE population of the Adriatic Sea basin (average pair-wise F_ST from 0.578 between Drava and PRE to 0.736 between Tisa and PRE). Within the major tributaries of the Danube, populations were more differentiated in Sava, Velika Morava and Tisa (average pair-wise F_ST 0.425, 0.473 and 0.502, respectively) than in Drava (average pair-wise F_ST = 0.281) (data not shown).

PCoA analysis provided good resolution of spatial population relationships reflecting their affiliation to the river systems (Fig. 5). The Adriatic Sea basin PRE population was clearly separated from the Black Sea basin populations on both PCo1 and PCo2. Further, PCo1 clearly separated populations belonging to the Tisa, Timok and Velika Morava tributaries from the populations belonging to the Sava and Vuka rivers (Fig. 5). PCo2 separated populations from the Tisa River tributaries and populations from Serbia (both Timok and Velika Morava tributaries) (Fig. 5).

Effective population size

Effective population size estimates varied strongly between populations and within river basins, ranging from 9 to 4920 (Table 5). The highest values were found in the Croatian populations of MOT, TOT and JAN located in the north-western part of the study area in the Drava River basin, with N_e values of 4920, 4864 and 4742, respectively. Globally, the Croatian populations displayed the highest effective population sizes (mean N_e = 2951), while lower estimates appeared mainly in the eastern populations of Serbia (mean N_e = 1113; the value drops to 539 if the north westernmost population of KAC from the Sava River basin is not considered), Romania (mean N_e = 272) and Albania (N_e = 491). The lowest N_e values were found in the Romanian population of BEZ (N_e = 11), and the Serbian populations of GAZ (N_e = 18) and GRL (N _e = 9).

Spatial analysis

Mantel’s tests showed significant correlations between genetic and geographical distances (r = 0.317 (P_Mantel < 0.001) for 655 bp long COI sequences dataset; r = 0.096 (P_Mantel < 0.001) for 350 bp long COI sequences dataset; r = 0.493 (P_Mantel < 0.001) for microsatellite dataset). Alleles in Space (AIS) analysis was used to visualise the general genetic divergence on the European level for the COI data set (Fig. 6). Also, AIS was used to delineate genetic divergence in the area of interest (Balkan Peninsula) for the microsatellite data set (Fig. 7). The GLS interpolation on the COI data set showed that the area of the highest genetic divergences for A. astacus was located in the southern part of the Balkan Peninsula (Greece, Albania and North Macedonia) (Fig. 6). Moderate level of genetic divergence was observed among populations in Croatia and Slovenia, whereas the lowest divergences were indicated among populations from southern Serbia, Romania, Germany and Poland. The GLS interpolation on the microsatellite dataset indicated the highest genetic divergences was across a large area in the western part of Balkans (covering Sava and Drava River basins in Croatia and Slovenia) and the area in the southern Balkans (approximately corresponding to North Macedonia) (Fig. 7). Areas with the lowest genetic divergence were from parts of Romania, Serbia and Hungary.