The short-finned eel Anguilla bicolor is known to be subdivided in two distinct subspecies (i.e. A. bicolor bicolor and A. bicolor pacifica), each subspecies being geographically distributed in allopatry. The present survey intends to describe genetic differentiation, population structure, molecular variance and phylogeny of both subspecies of A. bicolor in Indonesian waters. The genotypes of seven microsatellite locus and sequences of the entire cytochrome b were analyzed on 180 specimens collected in 10 representative locations, where one of the two subspecies spend their freshwater life. The results showed high heterozygosity (0.767 < He < 0.891). Significant deviation from Hardy-Weinberg equilibrium were essentially detected on AjTR04 and Aro63 loci. No diagnostic microsatellite loci was observed between the subspecies which shared most of their alleles. Genetic Reynolds distances computed for each population ranged between 0.029 and 0.073 among A. bicolor pacifica populations, between 0.045 and 0.149 among A. bicolor bicolor populations and between 0.042 and 0.114 among populations of different subspecies. Both the mitochondrial and the microsatellite markers confirm the subdivision into two subspecies while microsatellite loci suggest a moderate differentiation between subspecies.
Introduction
Freshwater eels of the genus Anguilla are catadromous fish with long spawning migrations expanding from a few hundred for tropical eels to thousands of kilometers for temperate ones (Aoyama 2009). Beside distance migration, biological differences between tropical and temperate eel species occur in pelagic larval duration, distribution range and spawning season (Arai et al. 1999, 2000). Geographic distribution of tropical anguillid species partly or completely overlaps with at least one other congener in spawning area or in growing habitat (Aoyama 2009), questioning the species integrity process in a shared spawning area.
Most of freshwater eels are considered to form randomly mating panmictic populations. The mechanism producing panmixia is that individuals from different freshwater origins spawn in common area in the ocean (Avise et al. 1986). Few species have been evidenced as completely panmictic, each comprising a single and randomly mating population (A. anguilla and A. rostrata for temperate species; A. borneensis, for tropical species). However, recent evidences from studies of microsatellite variation of European eel has challenged this model (Wirth & Bernatchez 2001, Daemen et al. 2001, Maes & Volckaert 2002). Wirth & Bernatchez (2001) found that populations of European eel show slight but significant genetic differentiation and genetic distance between populations. Nevertheless, this observation was rapidly contradicted by the fact that observed disequilibrium should be due to temporal differences between samples and stochastic events (Dannewitz et al. 2005, Maes et al. 2006, Pujolar et al. 2006, 2014, Côté et al. 2013). This question is still debated.
The short-finned eel, Anguilla bicolor, has a relatively wide geographic distribution compared to the 19 species and subspecies of the genus Anguilla. The species is distributed longitudinally from the eastern coasts of Africa through the seas around Indonesia to New Guinea in Pacific Ocean (Ege 1939). The wide distribution of A. bicolor is one of the characteristics of marine fishes including large population sizes and pelagic eggs and larvae (Palumbi 1994, Ward et al. 1994). These characteristics, combined with the lack of physical barrier in the oceans, lead to the expectation of weak or no population structure (Waples 1998).
However, Ege (1939) was the first to recognize two distinct subspecies within this species: A. bicolor bicolor is restricted to the Indian Ocean and A. bicolor pacifica occurs in the western Pacific Ocean and all the seas around the northern part of Indonesia. Allopatric distribution has led to differences in the morphology of both subspecies and Ege (1939) noticed slight differences in their means and ranges of total number of vertebrae.
Indonesia is at an important biogeographic crossroads for marine ichthyofauna with its western part connected to Indian Ocean and its eastern part connected to Pacific Ocean. Concerning the widespread distribution of A. bicolor, this exceptional geographic distribution and structure made this species an important model for eel biology understanding, together with A. marmorata showing a similar range (Gagnaire et al. 2011). A. bicolor has been investigated at its whole distribution scale based on microsatellites and CR sequences (Minegishi et al. 2012). While the authors considers that more samples are necessary for a complete description, especially in its eastern distribution, they depicted, through microsatellite data, a species split into (i) an homogeneous western part in Indian Ocean (A. b. bicolor), as already observed on morphology and mitochondrial DNA (mtDNA) by Watanabe et al. (2005) and (ii) an eastern one only sampled in Philippines (A. b. pacifica).
In fact, despite the commercial importance of this widespread and abundant species, few is known in Indonesian waters, such as its exact distribution and genetic differentiation. The publication of Sriati & Sulistiono (2001) try to describe the elver abundances in Indonesia but obviously suffers of approximate taxonomy. The investigations for spawning area localization by Aoyama et al. (2007) on A. b. bicolor in front of Sumatra and of Wouthuyzen et al. (2009) on A. b. pacifica around Sulawesi were mainly inconclusive. More recently, Fahmi et al. (2012) gave a complete description of the two subspecies range in Indonesian waters.
Since a decade, hyper-variable markers such as microsatellites are widely used for studying genetic mapping, population structures and evolutionary genetics (Wirth & Bernatchez 2001, Gagnaire et al. 2009, Minegishi et al. 2012). Indonesian eels have been rarely involved. Revealing the genetic connectivity and population structure of the short-finned A. bicolor in Indonesian waters by using microsatellites is therefore useful. In comparison to temperate monopopulation species, the large distribution and the taxonomy of this species question the number of its spawning areas, its panmixia within each subspecies and the differentiation between them.
Material and Methods
Sampling
A total of 180 specimens belonging to A. bicolor bicolor and A. b. pacifica were collected from ten Indonesian localities in Indian and Pacific Oceans. Their subspecies were previously determined (Fahmi et al. 2012) using the mtDNA semi-multiplex method (Fahmi et al. 2013). Concerning the subspecies A. b. bicolor, the following samples were involved: Aceh (24 individuals), Padang (24), Bengkulu (23), Pelabuhan Ratu (24), Pangandaran (21) and Bali (24). Localities sampled for A. b. pacifica are: Bengalon (11), Donggala (18), Obi (7) and Poigar (4). All specimens were collected between 2008 and 2011 (see Fig. 1 and Table 1 for more details).
Microsatellite genotyping
Total genomic DNA was extracted using gSYNC Mini Kit Tissue (Geneaid). The protocol, according to the manufacturer's recommendations, is as follow: approximately 25 mg of eel tissue were incubated with lysis buffer and Proteinase-K at 60 °C for one hour or until lysis was completed, and then incubated with a second lysis buffer at 70 °C. The following steps included the DNA precipitation with ethanol and its separation from undesirable material by using GD column. The last step was the DNA washing successively with two distinct wash buffers. The genomic DNA was eluted by Nuclease Free-Water. The quality of DNA was checked by electrophoresis in agarose gel.
Eight microsatellite loci, which have been published for other anguillid eels species, were crossprimed (Table 2). Polymerase chain reactions (PCR) were performed on a Mastercycler (Eppendorf, Hamburg, Germany) in a 10 μl reaction volume containing 1.5 mM MgCl2, 0.2 mM of each dNTP, 1μl 10× PCR buffer (Promega, Fitchurg WI, USA), 0.5 μM of forward and reverse fluorescent-labeled primers, 0.2 units of Taq polymerase, and 2.0 μl of template DNA. Amplification program was 35 cycles of denaturation at 95 °C for 30 s, annealing for 30 s (see Table 2 for annealing temperatures) and extension at 72 °C for 30 s after heating at 95 °C for 30 s. Labeled PCR product were electrophoresed with a DNA ladder (Promega) on an 8.0 % acrylamide gel, which was later scanned on a FMBIO II Multi-View (Hitachi, Tokyo, Japan) scanner. Fragment size were determined using the software FMBIO Analysis 8.0 (Hitachi). Each molecule size produced by the software is then assimilated to an expected allele after checking the gel by eye.
Fig. 1.
Sampling location in Indian Ocean for Anguilla bicolor bicolor (square shape) and in Pacific Ocean for Anguilla bicolor pacifica (circle shape). Main regional currents: ECC = Equatorial Counter Current, IMC = Indian Monsoon Current, ITF = Indonesian Through Flow, SEC = South Equatorial Current, SJC = South Java Current.
Table 1.
List of specimen used for microsatellite analyses in the present study.
Cytochrome b sequencing
The cytochrome b was amplified through polymerase chain reaction (PCR). The PCR primers were designed especially for this study, based on sequences of all species of the genus Anguilla published in GenBank with number access AP007233-AP007249. The whole sequences of cytochrome b were aligned with Clustal W [in Molecular Evolution Genetic Analysis (MEGA) ver. 5.1 computer software package (Kumar et al. 2001)] and two invariable zones detected for primers determination. A pair of primers, FEEL-Cytb: 5′-CCA CCG TTG TAA TTC AAC-3′ and REELCytb: 5′-AAG CTA CTA GGC TTA TC-3′ allowed the amplification of a 1184 bp fragment. PCR was carried out in a total volume 25 μl containing 13 μl FastStart PCR Master mix, 4 μl of sterile distilled water, 2 μl of each primer and 5 μl of total DNA. Amplification of the fragment was performed in a Bio-Rad Thermal Cycler, programmed for a denaturation step at 94 °C for 5 min, followed by 35 cycles with one minute of denaturation at 94 °C, 1 minute of annealing at 55 °C and one minute of extension at 72 °C. The final extension was carried out at 72 °C for 10 minutes. The PCR products were checked on 1.5 % agarose gel electrophoresis and stained within 1 % Cyber Safe for 35 minutes. The migrated DNA was visualized using blue light under digital camera.
Fig. 2.
Phylogenetic relationship of A. bicolor based on cytochrome b gene (ML construction). Only bootstrap value over 50 % are indicated.
Sequencing procedure, amounts of reagents and buffers followed the manufacture's protocols (ABI, Applied Biosystem Inc.) and dye labeled fragments were analyzed on a capillary sequencer (ABI, Applied Biosystem Inc.). After checking by eye, the sequence were then stored into a MEGA ver.5.1 file.
The molecular phylogenetic tree was constructed with the MEGA ver. 5.1 computer software package by using the maximum likelihood (ML) method under the best model of DNA evolution chosen from 24 possible models and following the procedure given in MEGA 5.1. The reliability of each branch was assessed by bootstrap with the Felsenstein's (1985) bootstrap resampling technique (100 replicates) in PHYLIP 3.6 software.
Table 2.
List of loci used in this study.
Population genetic statistical analyses
Most statistical analyses were performed using the program GENETIX 4.05 (Belkhir et al. 1998): (i) calculation of the mean number of alleles by locus (A), allele frequencies, observed heterozygosity (Ho), expected heterozygosity (He) for each locus in each locality for all samples, (ii) estimation of Fis following Weir & Cockerham (1984) and its statistical significance through 5000 permutations and (iii) linkage disequilibrium among loci (significance estimated with 5000 permutations test). Multiple test significances were subjected to the Bonferroni correction (Rice 1989).
The genetic differentiations between localities and oceans were estimated using the estimator theta (Weir & Cockerham 1984) of the Fst as available in GENETIX 4.05, together with the genetic distance (Reynolds et al. 1983) as implemented in PHYLIP 3.6 (Felsenstein 2005).
The presence of intraspecific genetic structure was tested using the assignment model-based clustering method of Pritchard et al. (2000) implemented in STRUCTURE 2.1 for each value of K, which is the number of genetically distinct population-like subgroups. The Markov Chain Monte Carlo scheme was run with a burn-in period of 100000 steps and a chain length of 200000 replicates following the admixture model. Each run between K = 1 and K = 6 was repeated six times. The DeltaK method (Evanno et al. 2005) was applied using Structure Harvester online software (Earl & von Holdt 2012) in order to detect the most informative level of K.
Results and Discussion
Mitochondrial DNA phylogeny
Phylogenetic relationships of A. bicolor haplotypes were estimated with the maximum-likelihood method (ML) based on the best substitution pattern (Fig. 2). The results of the analysis with the “model selection — MEGA5″ showed the model HKY + G as the best substitution pattern with the smallest Bayesian Information Criterion value (BIC = 5347.657). The Hasegawa-Kishino-Yano model allows unequal base frequencies, distinguishes between the rate of transitions and transversions, and assumes a non-uniformity of evolutionary rates among sites modeled by a discrete Gamma distribution (+G). The estimated substitution rates were, 0.426 (A <-> G) and 0.526 (C <-> T) for transitions and 0.014 (A <-> T), 0.010 (C <-> G), 0.010 (G <-> T), 0.014 (A <-> C). The estimated transition/transversion ratio R TS/TV = 19.18.
Fig. 3.
Clustering analyses (K = 2) for the A. bicolor data set (180 individual, 10 population and 7 loci) performed using STRUCTURE program.
Fig. 4.
Neighbor Joining construction based on Reynolds distances, showing the Pacific samples clustering together (underlined location codes).
This tree is very similar to that published by Minegishi et al. (2012). The two subspecies are clearly discriminated into an Indian and a Pacific clades. The bicolor subspecies, in Indian Ocean, is composed of two distinct families of haplotypes with no geographic specificity, corresponding to haplotype groups of relict ancestral populations (Minegishi et al. 2012).
Microsatellites description
Eight pairs of primer were tested and seven of them gave easily interpretable and highly polymorphic microsatellite patterns. The genetic variability of each population is shown in Table 3. A total of 172 alleles were recorded in the seven microsatellite loci with a number by locus ranging from three in AjTR12 to 22 in Aro63, all populations included. The population parameters at the seven microsatellite loci are presented in Table 3. Multilocus expected (He) and observed (Ho) heterozygosities range from 0.767 to 0.891 and from 0.594 to 0.777 respectively. The two loci AjTR04 and Aro63 were the most polymorphic with average alleles numbers of 13.7 and 13.6 respectively.
Panmixia
Using Fis estimations, significant deviation from Hardy-Weinberg equilibrium were detected in 26 out of 70 potential cases (7 loci × 10 populations) and essentially on the loci AjTR12, Aro63 and AjTR04. According to multilocus Fis, given in Table 3, all samples show a significant deficit in heterozygotes except Poigar, probably because of its small sample size. Considering the eels vital cycle, panmixia is expected in each location together with no linkage disequilibrium nor differentiation between locations of the same subspecies.
Nuclear linkage disequilibrium
The linkage disequilibrium (LD) allow to detect non-random associations of alleles at two or more loci. LD was calculated only on samples of more than 15 specimens (Obi, Poigar and Bengalon were not included in the analysis). There is no significant linkage disequilibrium between genotypes for the seven microsatellite loci and for analyzed populations. All loci are independent at p < 0.05.
Inter-samples nuclear differentiation
Mean genetic differentiation between populations as estimated through Fst is 0.018. An Fst range value between 0.000 to 0.053 corresponds to a moderate genetic differentiation (Table 4).
Table 3.
Genetic variability of population parameters at seven microsatellite loci of tropical eel in Indonesian waters.
Table 4.
Fst estimators (Theta of Weir & Cockerham, 1984) are given above the diagonal (* and **: significant departure from zero after Bonferroni correction) and Reynolds genetic distance among population (below the diagonal). Pacific located stations are designated by underlined codes.
The genetic distances between each pair of populations were calculated according to Reynolds distance (Reynolds et al. 1983), and are given in Table 4. The stemming Neighbor Joining (NJ) tree (using PHYLIP Ver.3.6) is shown on Fig. 4 with bootstrap values based on 1000 replications. Genetic distances between populations ranged from 0.029 (between Pelabuhan Ratu and Benkulu) to 0.149 (between Obi and Poigar). Genetic distances between populations from Indian Ocean group (i.e. A. b. bicolor) ranged from 0.045 to 0.149, between populations from Pacific Ocean group (A. b. pacifica) they ranged from 0.029 to 0.073 and between populations located in Pacific and Indian Oceans from 0.042 to 0.114. According to the phenetic tree and the bootstrap values, A. bicolor is divided into three groups: the first group consists of six A. b. bicolor populations from Indian Ocean (Aceh, Padang, Pelabuhan Ratu, Pangandaran, Bali and Bengkulu, bootstrap value: 71 %); the second group consists of three populations from Pacific (Donggala, Bengalon and Obi, bootstrap value: 61 %) and finally a third cluster constituted by only one Pacific sample, Poigar (bootstrap value: 100 %). The results confirm the topology fitting with geography (Indian versus Pacific) or according to the subspecies.
Assignment tests analyzed under Structure program according to microsatellite genotypes followed by the DeltaK method (Evanno et al. 2005), indicated that the best K is three (Fig. 3). No clear inter-samples nor intersubspecies difference of the three subgroups is observed among the 10 populations of A. bicolor in Indonesia waters: for each value of K, Structure attributes various fractions of each cluster to all individuals.
Indonesian A. bicolor eels structure
The assignment analysis (STRUCTURE software) was unable to detect any geographic organization within each subspecies, which was expected, but also between both subspecies, which is quite surprising (see mtDNA differentiation). The assignment cut inside each individual, not between them. Two hypotheses can be proposed: (i) the present A. bicolor species is composed of three original clusters now perfectly mixed and (ii) this is a statistical artefact. The first hypothesis adds nothing to the question of geographic structure and the second hypothesis is far more parsimonious because more fitting the mtDNA tree and the eels biology predicting large panmictic populations.
Distances and Fst estimations gave a more expected result: the NJ tree (Fig. 4) clearly separated the samples from Pacific Ocean (Donggala, Bengalon, Poigar, Obi) to the other samples from the Indian Ocean (Aceh, Padang, Pangandaran, Bali). The populations from Indian Ocean which form this genetic cluster are divided in two subcluster, but with no bootstrap support and no geographic coherence.
According to Fis estimations, each sample within each subspecies appears in disequilibrium. This result is similar to that of Wirth & Bernatchez (2001) on European eels. It was evidenced, after several other surveys (Dannewitz et al. 2005, Côté et al. 2013) that the disequilibrium should be due to differences in sampling periods and in development stage of sampled eels. In the present A. bicolor study, sampling dates spanned along three years and the eels stages glass eels and subadults of different ages. Diachronic sampling and differences in development stages of sampled eels can be a good explanation for the observed unexpected differentiation within each Indian and Pacific supposed homogeneous populations.
The slight differentiation observed between subspecies (with genetic distances but not with assignment analysis), while mtDNA easily confirms this taxonomy, indicating that this differentiation is moderate. Microsatellites have been successful when used to solve the long controversial issue of panmixia and population genetic structure of European eels (A. anguilla) and American eels (A. rostrata) (Wirth & Bernatchez 2001, Mank & Avise 2003, Dannewitz et al. 2005) because both species have slight genetic and morphological differences. In the tropical area, using of microsatellite also have been involved to solve the population structure of wide geographical distribution of tropical eel A. marmorata (Ishikawa et al. 2004, Minegishi et al. 2008, Gagnaire et al. 2009).
The first study of population structure of A. bicolor using microsatellite marker has been done by Minegishi et al. (2012). This study confirmed the genetic differentiation among both subspecies and also revealed a genetic homogeneity within A. bicolor bicolor in the Indian Ocean. Compared to the present survey, the number of loci is similar but only four of them are common. The efficiency of the Minegishi et al. (2012) study to differentiate the two subspecies using microsatellites can be due to this difference in loci choice. It can also be due to the fact that Minegishi et al. (2012) used a Philippine sample as unique representative of the pacifica subspecies, highlighting a Philipine particularity rather than a A. b. pacifica one.
Conclusion
The present survey, when investigated on the mtDNA, confirmed the clear structure of A. bicolor into two subspecies. Using microsatellites nuclear DNA, a limited genetic divergence (Fst, genetic distances, NJ tree) was observed in A. bicolor in Indonesian water, mainly between Pacific and Indian Oceans, which is consistent with the classical subspecies designation. Other tests (assignment) confirmed the moderate differentiation between the two subspecies. Unexpected disequilibria (Fis) and differentiation within each subspecies (Fst, genetic distances) do not fit with the mono-population subspecies pattern of each subspecies. Diachronic sampling only partially explain these observations. Further sampling together with more and different nuclear markers are necessary to understand the two subspecies populations genetics like it was necessary to close the A. anguilla/A. rostrata debate (Côté et al. 2013).
