© EDP Sciences 2017

1 Introduction

It would be rational to exploit biodiversity at different levels sustainably, including genetic diversity. Yet, this goal has not been achieved and worse, genetic diversity continues to decrease (Small, 2011). Freshwater eels (Anguillidae) form a unique group of fishes very valuable both commercially and scientifically (Van Ginneken and Maes, 2005; Arai, 2014a; Zhaoqun et al., 2015). According to the latest taxonomy of freshwater eels, the genus Anguilla includes 19 species/subspecies (Arai, 2016). Most species (11) live in the tropics but five are temperate. During the past decades, all temperate freshwater eel species have gone into decline and today some of them are threatened of extinction (Arai, 2014b). The declines of these species probably were caused by the synergistic effect of overfishing, habitat loss or fragmentation, mortality from hydropower turbines, chemical contamination, parasites, and climate warming affecting oceanic spawning sites (Castonguay and Durif, 2016). Regional atmospherically driven ocean current variations in the spawning grounds also could drastically affect eel recruitment (Bonhommeau et al., 2008; Aoyama, 2009; Baltazar-Soares et al., 2014). The first claim that the continental stock of European eel Anguilla anguilla (L.) was in decline dates back to the early 1800s (Castonguay and Durif, 2016), but serious actions were taken only in 1980s. Currently, management and conservation of the critically endangered (Stacey et al., 2014) European eel is a challenging task that requires international collaboration among managers and scientists (Kettle et al., 2011; Capoccioni et al., 2014). Most difficulties emerge due to a complex life cycle of the species (Shiao et al., 2006).

The life cycle of the European eel was a mystery until early 20th century and even nowadays scientists lack important information about the biology and genetics of this species (Aarestrup et al., 2009; Miller et al., 2015; Arai, 2016; Righton et al., 2016). It begins when silver eels arrive from Europe or North Africa to the Sargasso Sea to spawn (Avise, 2011). After spawning, leptocephali enter the main ocean currents, such as the Gulf Stream, and travel to the continental shelf (Cagnon et al., 2011; Baltazar-Soares et al., 2014), where they metamorphose into glass eels. At this stage, the eels are capable of active swimming and try to enter freshwater systems, though some remain oceanic (Harrod et al., 2005). After arriving to the continental feeding habitats the glass eels transform into elvers and then yellow eels (Arai, 2016). Finally, after growing for years (5–50) in the continental ecosystems, yellow eels become silver eels, which swim to the spawning grounds to reproduce (Dekker, 2000; Aarestrup et al., 2009; Righton et al., 2016). So far breeding eels and their eggs have not been detected in the wild and artificial breeding of this species is not profitable economically. To overcome the many limitations of the ecological studies of migrating eels in the Atlantic Ocean, especially the direct observations in the spawning grounds presumably located in the Sargasso Sea, it is possible to use indirect methods, such as otolith microchemistry and population genetics (Van Ginneken and Maes, 2005). Comprehensive investigations into genetic diversity of this species, as well as the American eel Anguilla rostrata and the Japanese eel Anguilla japonica, are clearly necessary both for practical purposes and fundamental research, as genetic data may help evaluate the number of genetic stocks, clarify the spatio-temporal stability of genetic structure, define the influences of oceanic conditions on genetic variability, and evaluate the effect of population decline on genetic variability and the overall fitness of eels (Maes and Volckaert, 2007).

Great progress in the population genetics of the European eel has been achieved over the years, but the panmixia hypothesis, which states that silver eels reproduce randomly at the spawning site independently of their growing localities in the continent, is still debated among eel researchers (Wirth and Bernatchez, 2001; Maes and Volckaert, 2002; Maes et al., 2006; Pujolar et al., 2006; Als et al., 2011; Avise, 2011; Kettle et al., 2011; Pujolar et al., 2011; Baltazar-Soares et al., 2014; Jacobsen et al., 2014; Pujolar et al., 2014; Baltazar-Soares and Eizaguirre, 2016; Baltazar-Soares et al., 2016). Wielgoss et al. (2014) presented an extensive summary of the molecular and population genetics literature for Atlantic eels (both the European eel and the American eel) until 2013. Based on the most extensive DNA microsatellites (Dannewitz et al., 2005; Palm et al., 2009; Als et al., 2011) and genomic studies (Gagnaire et al., 2012; Jacobsen et al., 2014; Pujolar et al., 2014) it could be argued that the European eel is a panmictic species experiencing within-generation selection. It is worth mentioning that Als et al. (2011) managed to generate the most comprehensive European eel dataset required for the verification of the panmixia hypothesis: larvae of European eel sampled directly in the spawning area together with glass eels were investigated using 21 microsatellite loci. However, the European eel continue to be a slippery research subject (Avise, 2011). For instance, application of in silico population genetics showed that observable genetic structure can result from both the panmixia and female philopatry scenarios, especially in years of low recruitment (Baltazar-Soares et al., 2014). In addition, the population genetic models support the existence of matriline-driven demes over a completely panmictic mode of reproduction and inferences of gene flow amongst matrilines (a mtDNA line of descent transmitted to new generations by maternal side) revealed that migration is asymmetric, inconsistent with a panmictic reproduction (Baltazar-Soares and Eizaguirre, 2016). Consequently, more investigations that aim to examine the potential maternal population structure in the European eel are necessary.

The first attempts of DNA sequencing of mtDNA D-loop (Lintas et al., 1998) or cyt b (Daemen et al., 2001) regions were not enough to show the actual potential of these molecular markers (Van Ginneken and Maes, 2005). New information and insights about the population genetic structure and intraspecific evolution of the European eel were obtained by two recently conducted, independent studies, which used the mtDNA D-loop region (Ragauskas et al., 2014) and the ND5 region (Baltazar-Soares et al., 2014) sequence analysis. These papers re-opened eel research using direct sequencing of different mtDNA regions and raised new important questions concerning the population genetic structure of this species (Baltazar-Soares and Eizaguirre, 2016). Daemen et al. (2001) suggested that the population genetic structure of the European eel occurs in a genetic mosaic consisting of isolated groups, but only after the studies of Baltazar-Soares et al. (2014), Ragauskas et al. (2014) and Baltazar-Soares and Eizaguirre (2016) the existence of reproductively isolated groups that are represented by phylogenetically related individuals has been re-evaluated. The existence of a genetic mosaic in the population of the European eel could be explained by a cryptic philopatric behavior of females spawning in distinct areas in the spawning site. With males maintaining gene flow among philopatric demes, apparent panmixia could still be detected. This cryptic organization may contribute to the maintenance of the adaptive potential of the species (Baltazar-Soares and Eizaguirre, 2016). Despite the fact that significant genetic differentiation was not detected using the mtDNA D-loop region (Ragauskas et al., 2014), a strong and significant genetic structure was observed using the ND5 region (Baltazar-Soares et al., 2014). Therefore, it seems that the ND5 region is more sensitive to the natural selection. In order to get more valuable information about the population genetic structure of the European eel additional mtDNA genetic markers should be used. It is reasonable to select cyt b as type I molecular marker, i.e., a marker associated with a gene of known function (Liu and Cordes, 2004), to complement the mtDNA D-loop region data, which were recently used for the evaluation of the population genetic structure of the European eel (Ragauskas et al., 2014). Thus far the analysis of mtDNA cyt b sequences has revealed no clear population genetic structure based on the distribution of different haplotypes of 107 individuals (Daemen et al., 2001) and no more attempts have been made to carry out similar studies into this species on the basis of the utilization of this genetic marker.

In this paper we present the analysis of the genetic structure of the European eel using the mtDNA cyt b molecular marker. New results were compared with previous results of the mtDNA D-loop sequences analysis. Moreover, in order to start comprehensive studies into the genetic effects of eel translocation the parameters of the genetic variability of naturally recruited and introduced eels in Lithuania and Latvia are presented.

2 Materials and methods

2.1 Sampling, DNA extraction and amplification

All European eels analysed in this study (n = 222) were yellow or silver eels with the exception of three individuals that were collected at least three decades ago and were at the glass eel stage. The particular site where those glass eels were caught is unknown, but it is known that they were fixed in formaline in Moscow (later a bottle with glass eels was transferred to Vilnius, Lithuania). Most of the studied eels were caught in Lithuania, Latvia and the Mediterranean Sea region (Fig. 1 and Table 1). Similarly to the former Ragauskas et al. (2014) research, eels caught in the Baltic Sea, the Curonian Lagoon and the Gulf of Riga represented a mixture of naturally recruited to Lithuania and Latvia eels (NRE) with estimated 2–20% of introduced eels (Shiao et al., 2006), whereas the eel samples taken from Lake Balsys, Lake Dringis and Lake Siesartis represented the introduced (INTR) eel group. Eels introduced in Lithuanian inland lakes are able to migrate to the Curonian Lagoon and the Baltic Sea via water networks due to the absence of serious migration barriers.

In general, DNA was extracted from frozen or ethanol-preserved muscle tissues following the method of Aljanabi and Martinez (1997) with slight variations. Few samples of DNA were also obtained from the parts of glass eels that were fixed in formalin. DNA fragments of mtDNA cyt b were amplified using the originally designed Ang2-F (5′-TGGCAAACCTACGAAAAACC-3′) and Ang2-R (5′-AAGTGGAATGCGAAGAATCG-3′) primer pair. The PCR volume for each sample was 25 μl and consisted of 1× PCR buffer (with 50 mM KCl), 0.2 mM dNTP, 0.2 μM of each primer, 2.5 mM MgCl₂, 0.75 U Taq DNA polymerase LC (MBI Fermentas, Vilnius, Lithuania), and 0.2 μg template DNA. Amplification started with an initial denaturation step for 5 min at 95 °C, followed by 35 cycles (denaturation for 45 s at 94 °C, annealing for 45 s at 52 °C, elongation for 1 min at 72 °C), and ended with a final elongation step for 5 min at 72 °C. The length of the amplified fragment was approximately 550 bp. PCR products were purified with the help of ExoI and FastAP (Thermo Fisher Scientific, Vilnius, Lithuania) and sequenced. Sequencing was carried out at the Sequencing Centre of the Institute of Biotechnology (Lithuania) or in the Laboratory of Molecular Ecology of the Institute of Ecology of Nature Research Centre (Lithuania) using Big-Dye^® Terminator v3.1 Cycle Sequencing Kit and 3130xl or 3500 Genetic Analyzer (Applied Biosystems, USA), respectively.

Fig. 1 The distribution of A. anguilla mtDNA cyt b haplogroups in Europe and North Africa: the eel samples collected during the earlier studies (Aoyama and Tsukamoto, 1997; Wolf et al., 2000; Aoyama et al., 2001; Daemen et al., 2001; Lin et al., 2001; Minegishi et al., 2005; Jacobsen et al., 2014) and data from GenBank (EF427617–EF427618, EU223996–EU223997, EU492326–EU492327) are indicated in squares, whereas the eel samples collected during the current work are indicated in circles; the diameter of pie charts is proportional to the number of sequences (1–52 sequences); eel samples: AR − the Arno River; BAL − Lake Balsys; BC − the Bristol Bay; BJ − the Baltic Sea; BR − the Burrishoole River; BRA − Lake of Bracciano; BU − Burrishoole; BUR − Burano; DRI − Lake Dringis; EFR − Eastern France; FR − the Coast of France; KM − the Curonian Lagoon; LE − Lough Erne; LGL − Lac Grand Lieu; MO − Moulouya; RH − Ringhals; RĮ − the Gulf of Riga; RM − Rio Minho; SE − the Estuary of Severn; SIE − Lake Siesartis; SP − the Coast of Spain; SR − the Sebou River; SW − the Coast of Sweden; SWZ − Switzerland; TUN − the Coast of Tunisia; UK − the Coast of England; VA − Valencia; VIS − the Viskan River.

2.2 MtDNA cyt b sequence analysis

A total of 394 mtDNA cyt b sequences were used for the intraspecific genetic analysis of the European eel. Sequences were truncated to 392 bp in order to investigate new data together with available Daemen et al. (2001) data. Two hundred and twenty-two out of these 394 sequences were obtained during the present study and 172 sequences were retrieved from the GenBank (Table 2). Due to the fact that catching locations of 4 eels were unknown, a set of 390 eel sequences was chosen for the statistical analysis of three groups representing the macrogeographic regions (the Atlantic Ocean Region [AOR], the North Sea Region [NSR] and the Mediterranean Sea Region [MSR]) and two groups of eels representing the eels naturally recruited to Lithuania and Latvia (NRE) or the eels introduced to the inland lakes of Lithuania (INTR). In one analysis NSR and NRE samples were combined into one sample, representing the hypothetical Northern subpopulation (NP). In order to directly compare the quantitative parameters of genetic diversity using different mtDNA molecular markers and to reveal similarities and differences between the genealogies of the A. anguilla mtDNA cyt b and D-loop region, a set of 146 eel individuals that were investigated both during the current study and previous study carried out by Ragauskas et al. (2014) was selected.

The sequence alignment was performed using the ClustalW algorithm (Thompson et al., 1994) implemented in the MEGA6 (Tamura et al., 2013) program. Quantitative parameters of genetic diversity, neutrality tests (Tajima, 1989; Fu and Li, 1993) and genetic differentiation (Hudson et al., 1992a,b) were estimated using the DNASP 4.50.3 program (Rozas et al., 2003). Genetic differentiation was estimated with the statistic Φ_ST between pairs of samples. The significance of the null hypothesis Φ_ST = 0 was tested by a permutation-based procedure, using Monte Carlo simulations to estimate significance levels, and is particularly suited for data on nucleotide variation in samples from two or more localities (Hudson et al., 1992a). The haplotype network was constructed using the Median Joining (MJ) (Bandelt et al., 1999) algorithm, implemented in NETWORK 4.6.1.3. software. Different haplogroups were distinguished and haplotypes were attributed to one or another haplogroup on the basis of the following criteria: sequences representing different rare haplotypes separated from the most frequent haplotype by minimum mutational steps were considered belonging to the same haplogroup. The frequency of haplogroups across regions was tested by Pearson's χ² test using R v. 3.3.3. (https://www.r-project.org/) program. Unless stated otherwise, none of the deletions was considered to be the 5th position. The Mantel test was used to test the concordance of the haplotype group structure of cyt b and D-loop region. In order to do it a diagonal 146 × 146 matrix with all the possible pairs of individual eels were made. Then each term of the matrix was considered and we gave value 1 if the individuals in the pair shared the same cyt b haplogroup and 0 if they do not. The same procedure was applied for D-loop region data. Finally, correlation between the two different matrices was tested by Mantel test/Two-tailed test using Microsoft XLSTAT 2017 software. The P value has been calculated using the distribution of r(AB) estimated from 10 000 permutations.

3 Results

3.1 Data on mtDNA sequence variation and haplotype network

On the basis of the 392 bp fragment of the mtDNA cyt b gene, a total of 62 different haplotypes out of 394 aligned sequences was ascertained (Fig. 2). Transversions were detected in 8 different positions in the analyzed DNA fragments, but indels were not found. Singletons accounted for 11.17%. The links between haplotypes in the haplotype network were without median vectors and were suitable for grouping these 62 haplotypes into 9 haplogroups, designated with capital letters (Figs. 1 and 2). The greatest genetic diversity was detected in the haplogroup A (represented by 19 different haplotypes). The prevailing mtDNA cyt b haplotypes of the European eel were haplotypes A and B, detected in 40.1% and 20.3% of scored eels, respectively. Only these most frequent haplotypes were detected in the glass eels caught several decades ago (UNK; n = 3). Despite the fact that all common haplotypes of the haplogroups A, B, C, D, E, F and G were found in the samples collected while carrying out this work, 19 rare haplotypes that were previously detected in Europe and North Africa were not observed this time. In addition, the following 32 new haplotypes, distributed among 9 haplogroups, were determined: A3, A4, A5, A6, A7, A8, A10, A11, A12, A13, A18, B3, B4, B5, B6, B7, B8, B9, B10, B11, C3, C4, C5, D2, E4, F2, F3, F4, G2, H1, H2 and I (Table 2).

Fig. 2 MJ haplotype network of 394 A. anguilla mtDNA cyt b sequences: mtDNA cyt b sequences of eels caught in other European locations, as well as their frequency, are marked in white; the numbers between the nodes indicate particular mutations (red − transitions; blue − transversions) between the most closely related sequences; the radius of the circles is proportional to the frequency of haplotypes; BJ − the Baltic Sea; KM − the Curonian Lagoon; RĮ − the Gulf of Riga; UK − the Coast of England; BAL − Lake Balsys; BRA − Lake of Bracciano; BUR − Burano; DRI − Lake Dringis; SIE − Lake Siesartis; TUN − the Coast of Tunisia; UNK − glass eels fixed in formaline in Moscow.

3.2 Distribution patterns of mtDNA cyt b haplogroups of the European eel in Europe and North Africa

In general, the results of the study revealed similar patterns of the distribution of mtDNA cyt b haplogroups of the European eel in different sampling sites of Europe and North Africa (Fig. 1). The tendencies of composition of haplogroups in macrogeographic and microgeographic regions were similar. However, the representatives of haplogroup B in the North Sea Region were rare.

The statistical significance of the frequency of haplogroups across regions using Pearson's χ² test presented in Table 3. Significant P values were revealed by pairwise comparisons between regions for AOR/NSR (P = 0.002587; χ² = 21.9544; df = 7), NSR/NRE (P = 0.02148; χ² = 14.8477; df = 6) and NSR/INTR (P = 0.007638; χ² = 20.821; df = 8).

3.3 Quantitative parameters of genetic diversity

The values of the quantitative parameters of genetic diversity of different samples of eels representing three macrogeographic regions and NRE and INTR groups of eels were quite similar (Table 4). The calculated parameters of genetic diversity h, K, S and π ranged from 0.741 ± 0.048 (MSR) to 0.827 ± 0.032 (NRE), from 1.2769 (MSR) to 1.8588 (INTR), from 11 (NSR) to 27 (INTR) and from 0.00326 ± 0.00039 (MSR) to 0.00474 ± 0.00048 (INTR), respectively. With the exception of the h value, all the remaining estimates of the quantitative parameters of genetic diversity were the highest in the INTR sample as compared to those of other studied samples. In general, samples of the eels collected in Lithuanian lakes differed considerably from one another according to the distribution of different haplotypes determined in each lake. For instance, 8, 7 and 8 unique haplotypes were detected exclusively in Lake Balsys, Lake Dringis and Lake Siesartis, respectively (Fig. 2). Haplotypes A, B, D and E were found in all the studied samples of the eels collected from Lithuanian lakes.

3.4 Neutrality tests

Both Tajima's D and Fu's Fs statistics were negative (−2.2594 and −90.068) and significant using all available sequence data (Table 5). This indicate that population of the European eel have an excess of rare nucleotide site variants compared to the expectation under a neutral model of evolution and most likely underwent a population expansion.

3.5 Genetic differentiation

Significant genetic differentiation was revealed by pairwise comparisons between regions for AOR/NSR (n = 135; Φ_ST = 0.02768; P = 0.0071; χ² = 45.595; df = 25) and NSR/NRE (n = 109; Φ_ST = 0.01144; P = 0.0440; χ² = 35.723; df = 23) (Table 6). After Bonferroni correction the result of determined genetic differentiation between the Atlantic Ocean and the North Sea remained statistically significant (0.0071 < 0.0083). No genetic differentiation was detected for the AOR/NSR/MSR/NRE/INTR comparison (n = 390; Φ_ST = 0.00729; P>0.05) or among stocked Lithuanian lakes (n = 106; Φ_ST = 0.00754; P > 0.05).

3.6 Quantitative parameters of genetic diversity using different mtDNA molecular markers and genealogies of the A. anguilla mtDNA cyt b and D-loop region

The comparison revealed that the quantitative parameters of genetic variability in the mtDNA D-loop region is several times greater as compared to that found in mtDNA cyt b sequences (Table 7). The analysis of sequences of both mitochondrial regions also revealed the existence of high haplotype diversity but relatively low nucleotide diversity.

The comparison presented in Figure 3 show that several genetic lineages do exist in the population of the European eel because certain cyt b haplotypes tend to be related with respective haplogroups of mtDNA D-loop region that is represented in the haplotype–haplogroup network. For instance, in the mtDNA D-loop region haplotype–haplogroup network haplogroup H is separated by many mutational steps from the haplogroups C7 and C8, as well as from haplogoups K, L and L1. In the same haplotype–haplogroup network E and E4 mtDNA cyt b haplotypes were observed exceptionally among the eel individuals attributed to the haplogroup H of mtDNA D-loop region. Similarly, F and F2 mtDNA cyt b haplotypes were observed exceptionally among the eel individuals attributed to the haplogroups C7 and C8 of mtDNA D-loop region. Finally, certain cyt b haplotype relation with respective haplogroup of mtDNA D-loop region was also observed in the case of haplotype D of mtDNA cyt b haplotype. The presented examples of relations between different mtDNA regions represent the existence of three distinct genetic lineages. However, mentioned tendency of direct relation between different mtDNA regions was not proved in some other cases. The distribution of the most common mtDNA cyt b haplotypes (A, B and G) in the mtDNA D-loop region haplotype–haplogroup network seems to be less predictable. In agreement with these variable results, the Mantel test for the concordance of the haplogroup structure of the two markers was significant but gave a low value for the correlation coefficient (r(AB) = 0.323; P = 0.0001).

Fig. 3 The distribution of mtDNA cyt b haplotypes plotted on the maximum parsimony (MP) haplotype–haplogroup network of A. anguilla mtDNA D-loop region (Ragauskas et al., 2014): numbers between the nodes are mutational steps between the most closely related sequences; the radius of the circles is proportional to the frequency of mtDNA D-loop region haplotypes; mv represent hypothetical sequences that were not found during the study but are important to the construction of the haplotype network; gray shapes represent different haplogroups of mtDNA D-loop region; different colors represent different mtDNA cyt b haplotypes; white circles represent mtDNA D-loop region haplotypes detected in the eels that were not examined using mtDNA cyt b.

4 Discussion

The data obtained in the current study combined with the available data of the earlier studies (Aoyama and Tsukamoto, 1997; Wolf et al., 2000; Aoyama et al., 2001; Daemen et al., 2001; Lin et al., 2001; Minegishi et al., 2005; Jacobsen et al., 2014 and unpublished data available from the GenBank under Accession numbers EF427617–EF427618, EU223996–EU223997, EU492326–EU492327) provide a dataset consisting of 394 sequences, which represent 62 haplotypes that belong to 9 haplogroups (Fig. 2). These data set indicates that cyt b gene in the European eel is characterized by an average haplotype diversity of 0.790 ± 0.017 and average nucleotide diversity of 0.00396 ± 0.000224. Quantitative parameters of genetic diversity (K, S and π) of European eel calculated for different mtDNA regions (Table 7) indicate that the genetic variability in the mtDNA D-loop region is several times greater in comparison to that found in mtDNA cyt b sequences. This result is expected and in agreement with the general knowledge of genetic diversity of these mtDNA regions in fishes (Meyer, 1993). Grant and Bowen (1998) reviewed the mitochondrial DNA variability in several species of fishes and found a range of 0.11 to 1.00 for h, and 0.07% to 3.20% for π. They also classified fish species according to different combinations of small and large values of h and π in four categories. Based on this classification, the European eel should be in the second category (high h and low π). High haplotype diversity but relatively low nucleotide diversity in comparison to many fish species indicates that the European eel most likely underwent a population expansion after a long-term period of low effective population size (Grant and Bowen, 1998). This conclusion is strengthened by negative and significant values of Tajima's D and Fu's Fs statistics (Table 5) and in agreement with Jacobsen et al. (2014) and Baltazar-Soares et al. (2016) findings. Alternatively, an excess of low frequency mutations in the European eel could indicate directional selection, but positive selection in the mitogenome of this species was clearly detected only in the ATP6 gene (Gagnaire et al., 2012). It would be interesting to know whether the patterns of diversity found at the mtDNA are also seen at the nuclear genome. Baltazar-Soares et al. (2016) studied mtDNA and two types of nuclear markers (microsatellites and MHC). Despite a clear discrepancy between signatures obtained for each class of markers, all of them revealed relatively high genetic diversity in the European eel. Recently conducted genomic investigations also determined relatively high genetic diversity in this species (Gagnaire et al., 2012; Jacobsen et al., 2014; Pujolar et al., 2014). High genetic diversity would provide the species with better chances to survive and reproduce.

The data of the current study (Table 4) indicate that genetic diversity of the eels sampled in the Mediterranean Sea region is slightly lower than that of the eels found in other parts of the species range, in agreement with previous studies (Daemen et al., 2001). These differences are not statistically significant due to large errors of the h and π estimates. However, and despite the fact that different samples of the eels from the Mediterranean Sea region were used in different studies, comparable tendencies were observed after the same parameters of genetic diversity obtained from mtDNA D-loop region sequences were evaluated (Ragauskas et al., 2014). Lower diversity values in the Mediterranean are unexpected for neutral markers given that all European eels are supposed to form a single mating group. In fact, catadromous silver eels from the western and eastern parts of the Mediterranean basin should likely reach the reproduction grounds in sequence thus creating many different “batches” of leptocephali (Maes et al., 2009). A possible explanation of the reduced diversity of the Mediterranean eels is offered by Capoccioni et al. (2014), on the basis of the (late) time of capture of some maturing eels from the deep bottoms of the Sicilian Strait: the eastern eels are lost for reproductive purposes, and only the western “stock” may therefore be able to mate and contribute (together with the African Atlantic individuals) to the replenishment of the eel production. In this way, only a few “batches” of haplotypes are generated from Mediterranean parents and present in southern Europe and North Africa. Also variations in the oceanic larval dispersal, described by Kettle et al. (2011) and Baltazar-Soares et al. (2014), may as well have resulted in lower diversity for the Mediterranean eels.

Our samples consisted of various phases (glass eels, yellow eels and silver eels) of eels that were caught in different time (indeed representing different generations of European eels that lived in the last decades). Due to restricted availability to investigate genetic differentiation of different eel samples in connection to time issues we focused on genetic differentiation between macrogeographical regions and NRE and INTR samples. During the current study statistically significant differentiation in haplotype frequencies was detected between the AOR and NSR regions (Φ_ST = 0.02768; P = 0.0071) and between NSR and NRE (Φ_ST = 0.01144; P = 0.0440) (Table 5). Taking into account that eel samples were collected in the continental ecosystems, it is challenging to evaluate the real biological meaning behind the observed differentiation. One possible explanation is that natural selection is acting on the mitogenome of European eel leading to a specific population structure. Baltazar-Soares et al. (2014) also detected strong and significant genetic structure across European eel populations using the ND5 region of mtDNA (F_ST values ranged between 0.058 and 0.16). The facts that the largest and most significant differentiation was detected between different pairwise comparisons containing two Irish samples and no significant genetic differentiation was detected among samples taken along the remaining European coast from Finland to Portugal suggested that positive selection could be acting over population structure. Evidence of positive selection was also observed by Jacobsen et al. (2014) in the ND5 gene in the American eel. An alternative explanation is that there is some cryptic structure in the spawning ground and that eels have some sort of philopatric behaviour (different groups of females mate in different parts of the reproductive grounds) (Wirth and Bernatchez, 2001; Baltazar-Soares and Eizaguirre, 2016). Under the hypothesis of female-structured spawning grounds, it would be expected to detect groups of related individuals in the continent. Als et al. (2011) identified two full-sib families of leptocephalus larvae at two westernmost transects. They also detected some structure in the Sargasso Sea, which showed up as a significant correlation between genetic and geographical distances. Baltazar-Soares and Eizaguirre (2016) used a combination of ocean modelling and genetics based on mitochondrial and nuclear loci to indirectly evaluate the possible existence of cryptic demes. They observed significant genetic differentiation among populations which was consistent with the existence of female-driven reproductive units combined with a low abundance of spawners due to current variations in the Sargasso Sea area.

The proposal of on genetic subdivision on the eel reproductive grounds resulting in genetic differentiation across continental localities is based on the existence of several mitochondrial matrilines. In this study, several most common haplotypes (most of them separated by 1–2 base pairs) were found in the haplotype network consisting of 394 sequences of mtDNA cyt b (Fig. 2). A similar tendency was observed in the haplotype network consisting of 403 sequences of mtDNA ND5 region studied by (Baltazar-Soares et al., 2014). It deviates from the star-like shape emerging from a single central haplotype pattern, which is characteristic of panmictic populations (Komiya et al., 2011). Even so, the unique biology of the European eel, stochasticity in its environment and randomness in coalescent process suggest that we should not expect to find a simple star-like shape in the haplotype network. The mtDNA ND5 region and cyt b haplotype networks are less complicated as compared to the haplotype network consisting of 229 sequences of the mtDNA D-loop region (Ragauskas et al., 2014). While the identification of major haplogroups both in the current study and in the Ragauskas et al. (2014) work requires additional examination, our results clearly revealed that some cyt b haplotypes tend to be related with some haplogroups of mtDNA D-loop region that is represented in the haplotype–haplogroup network (Fig. 3). In addition, Mantel test/Two-tailed test detected significant relationship (r(AB) = 0.323; P = 0.0001) between the haplogroups at the two gene regions.

To sum up, the potential for maternal population structure in the species was verified using mtDNA cyt b maker and recent knowledge of the European eel. The results of the current study are in agreement with Ragauskas et al. (2014) findings indicating that the population genetic structure of the European eel should be characterized as a genetic mosaic (Daemen et al., 2001) formed by the existence of reproductively isolated groups. Signals of existence of reproductively isolated groups presuppose that the European eel was not a panmictic species in the past and/or ongoing cryptic philopatric behavior of females spawning in fragmented spawning site. However, their contribution to the maintenance of the adaptive potential of the species (Baltazar-Soares and Eizaguirre, 2016) requires additional investigations. It could be suggested that this species should be managed just as a single unit (Als et al., 2011), but genetic diversity must be investigated and preserved in eel contingents associated with the different ecotypes, as it was recommended for the American eel (Pavey et al., 2015). It is also necessary to evaluate issues related to eels stocking. For instance, in order to maintain commercial European eel fisheries, intense stocking programmes have been undertaken in the Baltic Region during the past 50 years (Westin, 2003). Today, similar activities are still pursued though to a lesser extend. However, most stocking programmes were carried out without evaluating their potential success (Shiao et al., 2006). Some scientists, focusing their attention primarily on the positive effects of restocking on the fishing yields, have assumed a neutral or hopefully positive effect to eel population (Dekker, 2000). For instance, restocking of glass eels or young yellow eels in a particular water body produces a positive effect on the yield several years later. Meanwhile ICES, focusing on the protection of the spawning stock, recommended not to use restocking of eels as a stock-rebuilding measure. Recently Stacey et al. (2014) investigated the conservation stocking of the American eel and suggested that this activity should be applied with caution, as stocked eels appear to be following life-history patterns comparable with conspecifics in the geographical range of the donor streams where they were collected. Despite the fact that stocked European eels from Lithuanian inland lakes can migrate to the Curonian Lagoon thus far it is unclear what percentage of eels reach the Baltic Sea. In addition, it is not clear whether introduced eels that are capable reach Lithuanian coastal waters can find their way to the spawning grounds. In case the restocked eels do not contribute to the spawning stock, the effect of the restocking activity on conservation is negligible. Alternatively, if they do contribute, then this activity affects the formation of the population genetic structure of the European eel. The results of the current study indicate that there are no statistically significant genetic differences between the groups of the eels naturally recruited to Lithuania and Latvia and those introduced to Lithuanian lakes (Tables 3 and 6). In order to reveal clearer genetic differences between the naturally recruited and restocked eels it is crucial to investigate their migratory history with the help of otolith microchemistry. However, this has not been done yet. It also could be suggested that the modeling of Baltazar-Soares et al. (2014) should be expanded and obtained conclusions verified, as their experiment was repeated annually for the period between 1960 and 2005 while the first claim that the continental stock of the European eel was in decline dates back to the early 1800s (Castonguay and Durif, 2016). The extended research into the similarities and differences between the genealogies of the A. anguilla mtDNA cyt b, D-loop region and ND5 region is also strongly advisable.

Acknowledgements

We thank Dr. M. Baltazar-Soares and three anonymous reviewers for valued input to the manuscript submitted. This research was supported by the Open Access to research infrastructure of the Nature Research Centre under Lithuanian open access network initiative. The authors are grateful to the Lithuanian-Latvian-Taiwan (Republic of China) Mutual Fund for financial support for joint research projects between the countries. Sincere gratitude goes to the Lithuanian State Studies Foundation and Research Council of Lithuania for providing PhD scholarships in support of this study. Special thanks go to Dr. M.S. Romdhane, Dr. L. Ložys, Dr. V. Kesminas, G. Riauba and fishermen who generously provided eel samples and S. Stropaitytė for support and technical assistance. Also, a word of thanks to Dr. P. Prakas and A. Pažusytė for their support and help in improving the quality of the manuscript prior to its submission. We are also indebted to Dr. D. Grauda and N. Krasnevska for assistance in statistical analysis.

References

All Tables

All Figures

Fig. 1 The distribution of A. anguilla mtDNA cyt b haplogroups in Europe and North Africa: the eel samples collected during the earlier studies (Aoyama and Tsukamoto, 1997; Wolf et al., 2000; Aoyama et al., 2001; Daemen et al., 2001; Lin et al., 2001; Minegishi et al., 2005; Jacobsen et al., 2014) and data from GenBank (EF427617–EF427618, EU223996–EU223997, EU492326–EU492327) are indicated in squares, whereas the eel samples collected during the current work are indicated in circles; the diameter of pie charts is proportional to the number of sequences (1–52 sequences); eel samples: AR − the Arno River; BAL − Lake Balsys; BC − the Bristol Bay; BJ − the Baltic Sea; BR − the Burrishoole River; BRA − Lake of Bracciano; BU − Burrishoole; BUR − Burano; DRI − Lake Dringis; EFR − Eastern France; FR − the Coast of France; KM − the Curonian Lagoon; LE − Lough Erne; LGL − Lac Grand Lieu; MO − Moulouya; RH − Ringhals; RĮ − the Gulf of Riga; RM − Rio Minho; SE − the Estuary of Severn; SIE − Lake Siesartis; SP − the Coast of Spain; SR − the Sebou River; SW − the Coast of Sweden; SWZ − Switzerland; TUN − the Coast of Tunisia; UK − the Coast of England; VA − Valencia; VIS − the Viskan River.

In the text

Fig. 2 MJ haplotype network of 394 A. anguilla mtDNA cyt b sequences: mtDNA cyt b sequences of eels caught in other European locations, as well as their frequency, are marked in white; the numbers between the nodes indicate particular mutations (red − transitions; blue − transversions) between the most closely related sequences; the radius of the circles is proportional to the frequency of haplotypes; BJ − the Baltic Sea; KM − the Curonian Lagoon; RĮ − the Gulf of Riga; UK − the Coast of England; BAL − Lake Balsys; BRA − Lake of Bracciano; BUR − Burano; DRI − Lake Dringis; SIE − Lake Siesartis; TUN − the Coast of Tunisia; UNK − glass eels fixed in formaline in Moscow.

In the text

Fig. 3 The distribution of mtDNA cyt b haplotypes plotted on the maximum parsimony (MP) haplotype–haplogroup network of A. anguilla mtDNA D-loop region (Ragauskas et al., 2014): numbers between the nodes are mutational steps between the most closely related sequences; the radius of the circles is proportional to the frequency of mtDNA D-loop region haplotypes; mv represent hypothetical sequences that were not found during the study but are important to the construction of the haplotype network; gray shapes represent different haplogroups of mtDNA D-loop region; different colors represent different mtDNA cyt b haplotypes; white circles represent mtDNA D-loop region haplotypes detected in the eels that were not examined using mtDNA cyt b.

In the text