Assessing Temporal Patterns and Species Composition of Glass Eel ( Anguilla spp.) Cohorts in Sumatra and Java Using DNA Barcodes

: Anguillid eels are widely acknowledged for their ecological and socio-economic value in many countries. Yet, knowledge regarding their biodiversity, distribution and abundance remains superﬁcial—particularly in tropical countries such as Indonesia, where demand for anguillid eels is steadily increasing along with the threat imposed by river infrastructure developments. We investigated the diversity of anguillid eels on the western Indonesian islands of Sumatra and Java using automated molecular classiﬁcation and genetic species delimitation methods to explore temporal patterns of glass eel cohorts entering inland waters. A total of 278 glass eels were collected from monthly samplings along the west coast of Sumatra and the south coast of Java between March 2017 and February 2018. An automated, DNA-based glass eel identiﬁcation was performed using a DNA barcode reference library consisting of 64 newly generated DNA barcodes and 117 DNA barcodes retrieved from BOLD for all nine Anguilla species known to occur in Indonesia. Species delimitation methods converged in delineating eight Molecular Operational Taxonomic Units (MOTUs), with A. nebolusa and A. bengalensis being undistinguishable by DNA barcodes. A total of four MOTUs were detected within the glass eel samples, corresponding to Anguilla bicolor , A. interioris, A. marmorata, and A. nebulosa/A. bengalensis . Monthly captures indicated that glass eel recruitment peaks in June, during the onset of the dry season, and that A. bicolor is the most prevalent species. Comparing indices of mitochondrial genetic diversity between yellow/silver eels, originating from several sites across the species range distribution, and glass eels, collected in West Sumatra and Java, indicated a marked difference. Glass eels displayed a much lower diversity than yellow/silver eels. Implications for the management of glass eel ﬁsheries and species conservation are discussed.


Introduction
The freshwater eel family Anguillidae consists of 20 species and two genera [1], all well known for their catadromous life-cycles. Adults spawn in marine environments, and hatched larvae, known as leptocephalus, migrate to inland waters. Upon approaching continental shelfs, larvae become competent (glass eels) before entering inland waters, where they acquire pigmentation (elvers). Eels further grow up in freshwaters (yellow eels) until they become sexually mature (silver eels) and return to marine spawning grounds [2][3][4][5][6].
Anguillid eels are a significant food resource around the world, and contribute greatly to many national economies [2,[7][8][9][10]. For example, approximately 150,000 tons of eels were consumed each year both in Japan (2000)(2001)(2002) and China (2012China ( -2013 [11]. In Indonesia, the annual eel trade is estimated to be worth 100 million USD [12]. Demand for anguillid species has steadily increased in recent years; however, eel farming solely relies on wildcaught juvenile eels (elvers, silver eels), as breeding in captivity is not yet commercially viable [11,13,14]. The growing demand, in combination with habitat destruction and fragmentation caused by dams and other infrastructures, has led to a precipitous decline in populations across the globe, and some species have become endangered [2,[15][16][17][18][19]. Unfortunately, scientific knowledge regarding eel diversity, distribution and abundance is sparse [15,19]. In particular, little is known about the ecology of anguillid eels in the tropics, even though such regions typically support relatively high abundances and species richness' of anguillid eels [7,[20][21][22]. Indeed, two-thirds of the known species of the genus, Anguilla, occur in the tropical Indian and Pacific Ocean, while the remaining species occur in temperate parts of both the Pacific and Atlantic Ocean [1]. For Indonesia, a total of seven to nine species has been reported, as the biological status of two species pairs (A. bengalensis vs. A. nebulosa, A. borneensis vs. A. malgumora) is still under debate [7,[22][23][24]. This diversity is among the highest in the Indo-Pacific Ocean; however, few attempts have been made to document diversity and distribution of Anguilla spp. in Indonesia. Previous studies have mostly focused on either specific life stages or selected species [7,22,25,26]. Therefore, there is a clear and pressing need to assess spatial and temporal patterns of all life stages of tropical anguillid eels to gain a comprehensive understanding of their ecology.
Studies of biodiversity and distribution of tropical anguillid eels were long impeded by difficulties in identifying species using morphological characters, particularly for early life stages [27]. The development of standardized molecular approaches to automated specimen identification, such as DNA barcoding [28,29], opened new perspectives for monitoring and the exploration of species boundaries. Several studies have already successfully applied DNA barcoding to characterize the Indonesian ichthyofauna in various contexts, ranging from island inventories [30] over lineage-specific re-appraisals [31][32][33][34][35][36] to larval identification [37,38].
This study investigated the diversity of anguillid eels on the western Indonesian islands of Sumatra and Java using DNA barcoding. Anguillid glass eels are currently being intensively harvested in the western parts of Indonesia, particularly A. bicolor and A. marmorata on Java [22]. In addition, the rapid development of dams and other river infrastructures is compromising upstream glass eel migration. Therefore, this study aimed to: (1) establish a DNA barcode reference library for all known Anguilla species occurring in Indonesia based on both available and newly generated Cytochrome Oxidase I (COI) sequences, (2) apply this library to automated glass eel identification to examine temporal patterns of recruitment and species composition of glass eel cohorts in western Indonesia, where the fishing pressure is currently high.

Sampling
Sampling consisted of two distinct campaigns in Indonesia targeting freshwater life stages (yellow/silver eels) and early life stages (glass eels/elvers) entering freshwaters. Yellow/silver eels were sampled across several sites on Sumatra, Java, Bali, Lombok, Sulawesi and Ambon Islands between November 2012 and July 2019 using electrofishing and fishing rods (DS-ANGUILLA; doi:xx). Specimens were photographed, individually labeled and voucher specimens were preserved in a 5% formalin solution. Prior to fixation, a fin clip or a muscle biopsy was taken and fixed separately in a 96% ethanol solution for further Diversity 2021, 13,193 3 of 15 genetic analyses. Both tissues and voucher specimens were deposited in the National Collections at the Museum Zoologicum Bogoriense (MZB), Research Center for Biology (RCB), Indonesian Institute of Sciences (LIPI). Specimens were further identified using field guides [39][40][41]. Glass eels and elvers were sampled using a systematic and standardized procedure, developed for the purpose of this study, on Sumatra and Java between March 2017 and February 2018 (Figure 1). A trap device (dimensions 3.5 m × 2 m × 0.8 m) was used for 6-h periods after dark once every week ( Figure 2). Glass eel and elver samples were collected and preserved in 70% ethanol solution for further genetic analysis.
Sampling consisted of two distinct campaigns in Indonesia targeting freshwater life stages (yellow/silver eels) and early life stages (glass eels/elvers) entering freshwaters. Yellow/silver eels were sampled across several sites on Sumatra, Java, Bali, Lombok, Sulawesi and Ambon Islands between November 2012 and July 2019 using electrofishing and fishing rods (DS-ANGUILLA; doi:xx). Specimens were photographed, individually labeled and voucher specimens were preserved in a 5% formalin solution. Prior to fixation, a fin clip or a muscle biopsy was taken and fixed separately in a 96% ethanol solution for further genetic analyses. Both tissues and voucher specimens were deposited in the National Collections at the Museum Zoologicum Bogoriense (MZB), Research Center for Biology (RCB), Indonesian Institute of Sciences (LIPI). Specimens were further identified using field guides [39][40][41]. Glass eels and elvers were sampled using a systematic and standardized procedure, developed for the purpose of this study, on Sumatra and Java between March 2017 and February 2018 ( Figure 1). A trap device (dimensions 3.5 m × 2 m × 0.8 m) was used for 6-h periods after dark once every week ( Figure 2). Glass eel and elver samples were collected and preserved in 70% ethanol solution for further genetic analysis.

DNA Extraction, Amplification and Sequencing
Total genomic DNA was extracted from muscle tissue using the Geneaid extraction procedure following manufacturer's specifications (www.geneaid.com, accessed on 12 March 2020). A partial fragment of 652 bp of the mitochondrial COI coding gene, corresponding to the standard DNA Barcoding fragment [42], was amplified using the universal primers, Fish-COI-F and COI-Fish-R [43]. The primer sequences were as follows: Fish The PCR products were visualized on a 1% agarose gel and purified (Thermo Scientific PCR purification kit). Yellow/silver eel amplicons were bi-directionally sequenced and larval amplicons were sequenced with the reverse primer (COI-Fish-R). Sequencing was done using the EZ-Seq service (Macrogen) for glass eels, and the Centre for Biodiversity Genomics (University of Guelph) for yellow/silver eels. Sequences were deposited in BOLD in the dataset DS_ANGUILLA (doi:XX) and GenBank

DNA Extraction, Amplification and Sequencing
Total genomic DNA was extracted from muscle tissue using the Geneaid extraction procedure following manufacturer's specifications (www.geneaid.com, accessed on 12 March 2020). A partial fragment of 652 bp of the mitochondrial COI coding gene, corresponding to the standard DNA Barcoding fragment [42], was amplified using the universal primers, Fish-COI-F and COI-Fish-R [43]. The primer sequences were as follows: The PCR products were visualized on a 1% agarose gel and purified (Thermo Scientific PCR purification kit). Yellow/silver eel amplicons were bi-directionally sequenced and larval amplicons were sequenced with the reverse primer (COI-Fish-R). Sequencing was done using the EZ-Seq service (Macrogen) for glass eels,

Data Analysis
As one main objective was to implement automated identification of Anguilla glass eels and elvers through DNA barcodes, we first established a reference library. The reference library was based both on sequences produced from yellow/silver eel specimens, identified using morphological characters, and on DNA barcode records retrieved from BOLD. To validate each reference, we examined the match of species boundaries, delineated either by morphological characters, or Molecular Operational Taxonomic Units (MOTU), representing diagnosable molecular lineages [45][46][47]. It is well accepted that a combination of different delimitation approaches is capable of overcoming potential pitfalls arising from uneven sampling [33,[48][49][50][51]. Therefore, we used a scheme based on a 50 percent consensus between four algorithms: (1) Refined Single Linkage (RESL) as implemented in BOLD and used to produce Barcode Index Numbers (BIN) [52], (2)  (mPTP) rates version as implemented in the stand-alone software mptp_0.2.3 [54,55], and (4) General Mixed Yule-Coalescent (GMYC) in its single (sGMYC) and multiple (mGMYC) rates version as implemented in the R package Splits 1.0-19 [56].
As the mPTP algorithm uses a phylogenetic tree as an input file, a maximum likelihood (ML) tree was first reconstructed using RAxML [57] based on a GTR+I+Γ substitution model. An ultrametric, fully resolved tree was reconstructed using the Bayesian approach implemented in BEAST 2.6.2 [58], to be later applied with the GMYC algorithm. Duplicated sequences were pruned prior to reconstructing the ultrametric tree based on a strict-clock of 1.2% per million year [59]. Two Markov chains, each with a length of 10 million, were run independently using Yule pure birth and GTR+I+Γ substitution models. Trees were sampled every 5000 states, after an initial burnin period of 1 million. Both runs were combined using LogCombiner 2.6.2 and the maximum credibility tree was constructed using TreeAnnotator 2.6.2 [58].
A final COI gene tree was reconstructed using the SpeciesTreeUCLN algorithm of the StarBEAST2 package [60]. This approach implements a mixed-model, including a coalescent component within species and a diversification component between species, that allows accounting for variations of substitution rates within and between species [61]. StarBEAST2 jointly reconstructs gene trees and species trees, and as such requires species designations, which were determined using the consensus of our species delimitation analyses. The StarBEAST2 analysis was performed using the same parameters as the BEAST analysis described above.
Finally, unknown sequences from early life stages were classified to the species level using the R package BarcodingR [62]. A single sequence was randomly selected for each of the MOTUs for classification using the following algorithms: (1) BP [62], (2) fuzzy-set [63], and (3) Bayesian [64]. Similarly, the final assignment of an unknown to a known sequence was established based on a 50% consensus.

Results
A total of 64 and 278 COI sequences were generated for all yellow/silver and glass eels sampled, respectively. The yellow/silver eel COI sequences were combined with 117 sequences retrieved from BOLD into a DNA barcode reference library comprising 181 COI mostly full-length sequences (652 bp) for the 9 species of Anguilla occurring in Indonesia [1,23].
DNA-based species delimitation analyses resulted in congruent delimitation schemes with 7 MOTUS for mPTP, 8 MOTUs for RESL, ABGD and mGMYC, 10 MOTUs for sPTP and 11 MOTUs for sGMYC; with a consensus scheme consisting of 8 MOTUs (Figure 3; Table S1). All MOTUs matched the species boundaries defined by morphology-based identifications, with the exception of Anguilla bengalensis and A. nebulosa, which displayed tightly intertwined mitochondrial lineages (Figure 3). The coalescent tree of this species pair was the youngest of the Anguilla mitochondrial gene trees, with a Most Recent Common Ancestor (MRCA) dated around 500,000 years. By contrast the mitochondrial MRCA of the Anguilla analyzed in our study dated back 5.3 million years (Myr). Species showed a distinct barcoding gap, which is defined as the lack of overlap between the distributions of maximum intraspecific and minimum interspecific genetic distance. Maximum intraspecific distances were ten-fold higher on average than minimum interspecific distances (Table 1).
Common Ancestor (MRCA) dated around 500,000 years. By contrast the mitochondri MRCA of the Anguilla analyzed in our study dated back 5.3 million years (Myr). Specie showed a distinct barcoding gap, which is defined as the lack of overlap between the di tributions of maximum intraspecific and minimum interspecific genetic distance. Max mum intraspecific distances were ten-fold higher on average than minimum interspecif distances (Table 1).

Figure 3.
Mitochondrial gene tree of Anguilla spp. inferred with SpeciesTreeUCLN, including 95% HPD interval for node age estimates, species boundaries according to morphology-based identifications, and genetic species delimitation results for six methods and their 50% consensus. Assignment of unknown glass eels/elvers sequences to known yellow/silver eels sequences resulted in fully congruent classifications across the three methods (Table 1 and  Table S2). Fuzzy average yielded the lowest support, with assignment probabilities ranging from 0.753 to 0.989 per species on average; while Bayesian assignments were all supported by probabilities of 1 (Table 1). Among the 278 larvae sequences, 15 were assigned to A. interioris, 26 to A. marmorata, 24 to the species pair A. bengalensis/A. nebulosa and 213 to A. bicolor (Table S2). During the sampling period from February 2017 to March 2018, glass eels were collected only between May and November 2017 (Figure 4). A peak of abundance was observed in May and June, a period corresponding to the onset of the dry season and tightly following the first peak of annual temperatures in May (Figure 4). Multi-species glass eel assemblages were observed throughout the dry season between May and August, and later in November, with a peak of species richness in June when 4 species were collected. eels were collected only between May and November 2017 (Figure 4). A peak dance was observed in May and June, a period corresponding to the onset of the son and tightly following the first peak of annual temperatures in May (Figure 4 species glass eel assemblages were observed throughout the dry season between August, and later in November, with a peak of species richness in June when were collected. All estimates of mitochondrial genetic diversity indicated that glass eels/e play a much lower diversity than their species gene pools ( Table 2). The numbe lotypes (h) and the haplotypic diversity (Hd) were four-fold and three-fold lowe tively, in glass eels/elvers than in yellow/silver eels. A similar trend was observ nucleotide level, with nucleotide diversity (π) and Theta (θw) being two-and-aand ten-fold lower, respectively, in glass eels/elvers than in yellow/silver eels. T was also reflected by noticeable differences in the Tajima D test, which shifted positive D value in larvae indicating a scarcity of rare alleles. However, none of th D tests performed were significant, except for the A. bengalensis/A. nebulosa spe among yellow/silver eels with a significant negative D value. All estimates of mitochondrial genetic diversity indicated that glass eels/elvers display a much lower diversity than their species gene pools ( Table 2). The number of haplotypes (h) and the haplotypic diversity (Hd) were four-fold and three-fold lower, respectively, in glass eels/elvers than in yellow/silver eels. A similar trend was observed at the nucleotide level, with nucleotide diversity (π) and Theta (θw) being two-and-a-half-fold and ten-fold lower, respectively, in glass eels/elvers than in yellow/silver eels. This result was also reflected by noticeable differences in the Tajima D test, which shifted toward a positive D value in larvae indicating a scarcity of rare alleles. However, none of the Tajima D tests performed were significant, except for the A. bengalensis/A. nebulosa species pair among yellow/silver eels with a significant negative D value. Table 2. Summary statistics of genetic diversity per species, including yellow/silver and glass eels/elevers, for the four species detected among glass eel cohorts. N, number of individuals; h, number of haplotypes; Hd, haplotypic diversity; π, nucleotide diversity; θw, theta; Tajima's D test including D value and significance of the test (* significant). Regional refers to estimates at the scale of species distribution ranges. Local refers to estimates at the scale of glass eel sampling sites.

Discussion
The present study further highlights the usefulness of DNA barcoding in capturing boundaries between Anguilla species [7,22,[73][74][75][76]. Among the nine species included, seven displayed well-differentiated mitochondrial lineages, with minimum genetic distances to the nearest neighbor largely exceeding maximum intraspecific genetic distance. DNA barcodes failed to differentiate between the species-pair Anguilla bengalensis and A. nebulosa. The young age of their mitochondrial MRCA may be indicative of a recent divergence. However, the persistence of ancestral polymorphisms between diversifying lineages [77] tend to inflate genetic diversity in diverging pairs [78]. Here, genetic diversity indices show no differences between this species pair and the other species, thus questioning their distinctiveness. Although A. nebulosa and A. bengalensis are currently both considered as valid species, they have a complex taxonomic history and were synonymized in the past [23]. The lack of a declared type locality in the original description of A. bengalensis is further confusing its taxonomic status. From a mitochondrial perspective, this species pair behaves as a single cohesive lineage, which suggests that further studies are needed across their distribution range to establish a more comprehensive assessment of their genetic diversity.
Tropical river systems support assemblages of multiple anguillid species [26,79,80]. The co-occurrence of multiple species in the tropics, particularly Southeast Asia, is thought to be the result of shorter migration routes, easing dispersal [24], and of higher productivity leading to better growth rates and increased carrying capacities of aquatic habitats [81]. The detection of multiple species in glass eel cohorts in Sumatra and Java is consistent with the co-occurrence of multiple species of anguillid eels. Although the four species detected have been previously reported from Western Indonesia [1,23,82], little information is available on diversity patterns and the abundance of yellow/silver eels upstream. In particular, our results question how species abundance in glass eel cohorts translates into yellow/silver eel communities. Ecological dynamics during growth may translate into important shifts in species composition from larval to later stage assemblages, a trend that was observed in previous DNA barcoding studies of early life stages of marine and freshwater fishes [38,[83][84][85]. An alternative explanation might be that anguillid eels are typical r-strategists, producing a large number of offspring at the cost of low survival and recruitment rates [5,21,86]. Such species are subject to important stochastic fluctuations of their demographic parameters depending on temporal heterogeneity and environmental disturbances [87]. Recent new records of amphidromous fishes (i.e., species for which the adult breeds in freshwater and the larva drifts into the sea) in Java, particularly in the family Gobiidae, based on a single or a few individuals, seem to corroborate the stochastic nature of migratory fish recruitment in Western Indonesia [30].
A. bicolor, for instance, was extremely abundant during our glass eel sampling campaign. Yet, yellow/silver eels of A. bicolor are known to be uncommon throughout the species' range distribution [88]. A. bicolor may have several spawning sites, one of which is supposed to be located off the Southwest cost of Sumatra [88]. The abundance of A. bicolor in glass eel samples may corroborate this hypothesis; however, it questions the supposed scarcity of yellow/silver eels. A similar trend has been observed for A. interioris, in which most yellow/silver eel specimen records have been reported from New Guinea, while glass eels have been observed in several western islands of Indonesia [25,81,88]. Due to this knowledge gap on distribution and biology, A. interioris is currently listed as data deficient in the IUCN Red List [6].
Comparison of mitochondrial genetic diversity between yellow/silver eels originating from their known species range and glass eels from local streams in Sumatra and Java indicated large differences. Several species, such as A. marmorata, A. bicolor and A. bengalensis/A. nebulosa, display wide range distributions throughout the Indian and Pacific oceans. These differences might be explained by the existence of several spawning sites and associated genetic population structuring [21,79]. Likewise, mechanisms inducing glass eel migration choices may account for these large differences between glass and yellow/silver eel mitochondrial genetic diversity. Natal homing, that is migration directed towards parental habitats, is currently observed for A. anguilla and A. americana, which seemingly share spawning grounds in the Atlantic Ocean and hybridize but whose yellow/silver eels do not co-occur in Europe and North America [89][90][91]. Alternatively, inter-annual variability in reproductive success and larval survival might be responsible for these large differences between glass and yellow/silver eel genetic diversity in A. marmorata, A. bicolor and A. bengalensis/A. nebulosa. This result is more intriguing for A. interioris, as its yellow/silver eels have been mainly observed in the northern rivers of New Guinea and it is supposed to have a single spawning site offshore of northern New Guinea. The presence of glass eels of A. interioris at low frequency in western Indonesia suggests random dispersal of glass eels for A. interioris and very low recruitment outside of New Guinea. These results suggest interspecific differences in life history traits, and phenotypic plasticity seems to be an important driver of species-specific life history traits and adaptation in Anguilla [92].
The general trend towards lower mitochondrial genetic diversity among glass eels in western Indonesia than among yellow/silver eels over species range distributions suggests small effective population sizes as a consequence of either fragmentation and the existence of several spawning sites or stochastic fluctuations in breeding and recruitment outcomes, as well as potential migration preferences for some species. This pattern suggests that Anguilla populations in South Sumatra and Java should probably be considered small populations. If commercial harvesting of glass eels were to be intensified in the future, this could bring them potentially close to an extinction vortex, i.e., a situation in which reduced population size and increased demographic variance either induce spatial fragmentation or a reduction in population adaptive potential [93,94]. Furthermore, the rapid development of infrastructure during the last two decades in this part of Indonesia has already resulted in fragmentation and reduced habitat availability, as well as water pollution associated with industrial waste and agricultural run-off [95][96][97][98].

Conclusions
Given the intense harvesting of glass eels in Java, the present study calls for a broader assessment of their recruitment dynamics to inform the sustainable management of these fisheries. The presence of only a small portion of the species' mitochondrial genetic diversity in Sumatra and Java glass eels sampled suggests that they have small effective population sizes and are potentially endangered in some regions. It also suggests important shifts of species abundance between glass and yellow/silver eels, pointing to major gaps in our understanding of tropical Anguilla life cycle ecology, and highlighting inadequacies in knowledge about glass eel recruitment per se for the sustainable management of Anguilla fisheries. The yellow and silver eel life stages of these species should receive high priority for conservation, given that they play an inherently important role in recruitment, particularly in the context of fragmented and disconnected aquatic landscapes.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10 .3390/d13050193/s1, Table S1: Results of the genetic species delimitation analyses. Table S2: Results of the assignment analyses of unknown sequences to known species.