DNA barcoding of a tropical anguillid eel, Anguilla bicolor (Actinopterygii: Anguilliformes), in Indo-Pacific region and notes on its population structure

The tropical anguillid eel, Anguilla bicolor McCelland, 1844, includes two subspecies, Anguilla bicolor bicolor McCelland, 1844 and Anguilla bicolor pacifica Schmidt, 1928, and is distributed across the Indo-Pacific region. Although A. bicolor is widely distributed and recognized as an important fish resource in the Indo-Pacific region, few studies have been conducted on its genetic variation and population structure. DNA barcoding of A. bicolor specimens collected in the Indo-Pacific region was carried out in this study using mitochondrial cytochrome c oxidase subunit I. Anguilla bicolor was found to diverge genetically, which supported its classification into two different subspecies. In addition, our study showed that A. bicolor bicolor had two genetically distinct populations/groups, and these different populations co-occur geographically in Indonesia and Malaysia in the eastern Indian Ocean. Our findings suggest that the eel larvae might be transported from at least two geographically different spawning grounds in the Indian Ocean, and then recruited to and settled in the same habitats in Indonesian and Malaysian waters. The molecular evidence calls for further research on the life history, stock assessment and protection of the populations of A. bicolor bicolor in Indonesia and Malaysia.

Due to the wide geographical distribution, offshore spawning migration, passive larval transportation and the absence of significant physical barriers in the Indo-Pacific Ocean, a weak or an absence of population structure was expected for A. bicolor (Minegishi et al. 2012, Fahmi et al. 2015, as observed in temperate species such as Anguilla anguilla Linnaeus, 1758 (Palm et al. 2009) and Anguilla japonica Temminck & Schlegel, 1846 (Han et al. 2008). However, based on its allopatric distribution and slight differences in morphological characteristics, A. bicolor has been subdivided into A. bicolor bicolor (Indian Ocean) and A. bicolor pacifica (western Pacific Ocean) (Ege 1939, Watanabe et al. 2014.
Anguilla bicolor was previously sampled across its entire geographic range and analyzed using microsatellites and mitochondrial control region sequences (Minegishi et al. 2012). They found that A. bicolor was divided into (i) a homogeneous eastern and western side of the Indian Ocean (A. bicolor bicolor), and (ii) a western side of the Pacific Ocean (A. bicolor pacifica). Analysis of mitochondrial cytochrome b sequences and microsatellite markers further confirmed the differentiation into two subspecies with the microsatellite loci, suggesting a moderate differentiation between the subspecies (Fahmi et al. 2015). Although A. bicolor bicolor in the Indian Ocean was found to have two mitochondrial sublineages with random occurrence among and within sampling locations (Minegishi et al. 2012, Fahmi et al. 2015, little work has been done on the population structure of A. bicolor bicolor in the Indian Ocean.
Knowledge on the life history of A. bicolor has been gradually accumulating (e.g., Arai et al. 1999a, 1999b, 2020, Robinet et al. 2003, 2019, Arai and Abdul Kadir 2017a, 2017b, Chai and Arai 2018. Hatching and recruitment are estimated to be almost year-round and at a constant age, respectively . Tropical glass eels are recruited to a river mouth throughout the year (Arai et al. 1999c, Sugeha et al. 2001, which corresponds to a year-round spawning Abdul Kadir 2017a). In contrast, temperate eels do not show year-round spawning and recruitment. The broad ranges of otolith strontium/calcium (Sr:Ca) ratios found in A. bicolor bicolor and A. bicolor pacifica suggest that these subspecies have flexible migratory behaviors in ambient waters (Chino and Arai 2010a, 2010b, 2020, Arai and Chino 2019. The aim of the present study was to investigate the genetic variation through DNA barcoding of A. bicolor in Indonesian, Malaysian and Philippine waters. All eel specimens were examined for their morphological characteristics and cytochrome c oxidase subunit I (COI) gene sequence. We also discuss the potential dispersion mechanism and importance of the population structure of A. bicolor in these regions.
A total of 18 anguillid eels were collected from Southeast Asia: nine specimens in Indonesia (Cilacap and Yogyakarta, central Java), four specimens in Malaysia (Penang, northwestern Peninsular Malaysia) and five specimens in the Philippines (Mindanao) through traps, hooks and lines from 2009 to 2019 ( Fig. 1, Table 1). The external morphometric characteristics were determined ( Table 1). The fin difference index (FDI), which is the relative ratio of the ano-dorsal (the origin of the dorsal fin to the anus) length (Z) and the total length (L T ) (Ege 1939), was calculated as follows: FDI = 100 Z L T -1 . DNeasy Blood & Tissue Kit (QIAGEN, Germany) was used to extract genomic DNA from a dorsal fin clip. Approximately 600 bp region of the mitochondrial COI gene was amplified using different combinations of universal primers: FishF1 (5'TCA ACC AAC CAC AAA GAC ATT GGC AC3'), FishF2 (5'TCG ACT AAT CAT AAA GAT ATC GGC AC3'), FishR1 (5'TAG ACT TCT GGG TGG CCA AAG AAT CA3'), and FishR2 (5'ACT TCA GGG TGA CCG AAG AAT CAG AA3') (Ward et al. 2005). The samples that could not be amplified by these primers were subsequently amplified with the primers, AngF (5'TCA CCC GTT GAT TCT TTT CT3') and AngR (5'CCG ATA GCC ATT ATT GCT CAG3'), which were designed in this study to amplify 835 bp of COI (starting from position 8 until 842) from Anguilla eels. The primers were designed by targeting the conserved sites of the COI genes from different species of sequenced Anguilla eels (Minegishi et al. 2005). Each PCR reaction contained 2 µl of DNA sample, 2.5 µl of each 10 µM primer, 25 µl of 2x Taq PCR Master Mix (QIAGEN, Germany) and 18 µl of distilled water. The PCR conditions were initially 95°C for 2 mins, then 35 cycles of 94°C for 30 s, 50°C for 30 s and 72°C for 60 s, and finally 72°C for 10 minutes. PCR amplicons were purified using QIAquick Gel Extraction Kit (QIAGEN, Germany), and sequenced bi-directionally with the same primers.
MEGA version X (Kumar et al. 2018) was used to edit and assemble DNA sequence trace files. To validate the species identity, the contig sequences were compared for percent similarity with the reference sequences in the GenBank database (Benson et al. 2013) by using BLAST search (Altschul et al. 1990). The contig sequences of eels from Indonesia, Malaysia and Philippines were submitted to the GenBank database with accession numbers MT107677 -MT107685, MT107673 -MT107676 and MT107686 -MT107690, respectively. In addition to these eighteen sequences, the COI sequences of thirty A. bicolor bicolor specimens from other localities that were found in the GenBank were also included in the analysis: Banda Aceh in northern Sumatra of Indonesia (22 specimens; KY618768 -KY618795), Sukabumi in western Java of Indonesia (one specimen; KU692247), Andhra Pradesh in southeast India (two specimens; MG675613, KP979655), Langkawi (one specimen; KF182304) and Penang islands of northwest Peninsular Malaysia (3 specimens; KM875503 -KM875505), and Yangon of Myanmar (one specimen; AP007236). The COI sequences of A. bicolor pacifica from Cebu in the Philippines (one specimen; AP007237) and Anguilla megastoma Kaup, 1856 (as outgroup; AP007243) were also included ( Fig. 1).
GUIDANCE version 2 (Sela et al. 2015) was used to create multiple sequence alignment via its inbuilt Clustal W version 2 (Larkin et al. 2007), which also assigned the alignment with a robust score of 1. The multiple sequence alignment was trimmed  Phillipps, 1925(Ege 1939, Watanabe et al. 2004, Arai 2016. To further determine which species from the list, diagnostic geographical distribution (Ege 1939) was used. Thirteen specimens from Indonesia and Malaysia and five specimens from the Philippines were identified as A. bicolor bicolor and A. bicolor pacifica, respectively (Table 1).
DNA barcoding using the COI gene confirmed that the 13 specimens from Indonesia and Malaysia and five specimens from the Philippines were A. bicolor bicolor and A. bicolor pacifica, respectively (Table 1), with 99-100 % maximum identity matches with the reference sequences in the GenBank database. Apart from the universal primers used, this study showed that DNA barcoding could also be carried out using novel primers that were designed for Anguilla eels.
From our phylogenetic analyses (Fig. 2), A. bicolor bicolor samples were observed to form a clade (bootstrap proportion = >60 %) that was distinct from the clade (bootstrap proportion ≥ 90 %) formed by A. bicolor pacifica samples. This confirms the differentiation of A. bicolor into two subspecies. The MP tree strongly supports the presence of two subgroups within the A. bicolor bicolor clade (bootstrap proportion: subgroup 1 = 80%, subgroup 2 = 66%). These two subgroups are also found in the NJ tree, though weakly supported by bootstrap proportion (subgroup 1 = 50%, subgroup 2 = 58%). However, these two subgroups were not clearly observed in the ML tree.
Haplotype analysis revealed a total of 18 different haplotypes in A. bicolor bicolor (H1 -H18 in Fig. 3) and 2 different haplotypes in A. bicolor pacifica (H19-H20 in Fig. 3). The constructed haplotype network (Fig. 3) further confirmed the presence of two mitochondrial sublineages in A. bicolor bicolor, as the haloptypes could be grouped into two haplogroups or populations.
Anguilla bicolor from the Indian Ocean (A. bicolor bicolor) and Pacific Ocean (A. bicolor pacifica) were found to diverge from each other, and the outcomes of our COI gene sequence were similar to those previous studies using microsatellites, mitochondrial control region and cytochrome b sequences (Minegishi et al. 2012, Fahmi et al. 2015. The previous studies suggest that A. bicolor has diverged into two subspecies in the two oceans. In the eastern and western sides of the Indian Ocean, the genetic profiles of A. bicolor bicolor were found to be similar. However, the genetic profile of either the eastern or western side of the Indian Ocean was observed to have two divergent mitochondrial sublineages (Minegishi et al. 2012, Fahmi et al. 2015. These two sublineages coexist in the Indian Ocean, occurring randomly among and within sampling locations. The present results also supported the occurrence of two mitochondrial sublineages, as there were two distinctive groups in the eastern Indian ocean (Figs 2, 3).
Based on the COI sequence analysis, the results suggest that different A. bicolor bicolor populations could co-occur geographically in Indonesian and Malaysian waters. Spawning areas of A. bicolor bicolor are thought to be located in the waters off west Sumatra (Jespersen 1942), and off east Madagascar (Jubb 1961, Robinet et al. 2003. These two geographically different spawning areas might explain the occurrence of two divergent populations of A. bicolor bicolor in the Indian Ocean. Thus, the present study supports the hypothetical occurrence of the two spawning grounds of A. bicolor bicolor in the Indian Ocean. Furthermore, the ages of recruitment of this species between the east (148-202 days in Java of Indonesia) and west coasts (68-96 days in Mascarene Islands) of the Indian Ocean definitely differed (Arai et al. 1999a, Robinet et al. 2003. This suggests that A. bicolor bicolor could be reasonably assumed to have different populations in the east and west sides of the Indian Ocean. Anguilla bicolor bicolor in the eastern Indian Ocean has a longer leptocephalus period (119-171 days) (Arai et al. 1999a) and is spawning throughout the year Abdul Kadir 2017a). These life history strategies might enable A. bicolor bicolor to disperse across the Indian Ocean. Seasonal changes of the South Equatorial Current and the long larval migration period in the Indian Ocean, which was estimated to be more than one year in anguillid eels (Pous et al. 2010), might be capable of mixing A. bicolor bicolor larvae from different spawning areas. Further studies are needed to determine the exact spawning grounds and larval transportation mechanism and routes in A. bicolor bicolor. Nevertheless, these particular life history characteristics of A. bicolor bicolor and the oceanographic features of the Indian Ocean could explain the sympatric occurrence of different populations. Figure 3. Haplotype network constructed with Anguilla bicolor COI gene sequences. Each colour represents a sample site. The size of the circle is proportional to the number of samples that belong to each haplotype. White circle refers to median vector which is a hypothesized or missing haplotype. Each dash, which appears on the line that connects two haplotypes together, symbolizes one mutational step.
The specific choice of spawning area might be imprinted in anguillid eels including A. bicolor bicolor, although there is no physical barrier to migration to the different spawning areas. The occurrence of two populations in the eastern Indian Ocean suggests that population divergence in anguillid eels could still occur even in the absence of physical barriers during the spawning migration. Minegishi et al. (2012) presented a different hypothesis, as this previous study found that A. bicolor bicolor diverged into two sublineages based on the mitochondrial control region but not when using the microsatellite markers. The study hypothesizes that there could be an incomplete lineage sorting or a possible secondary contact between the two groups after an initial population splitting.
This study focused only in the eastern side but did not include the western side of Indian Ocean. The sample sizes were also small in all localities. Therefore, more comprehensive research is needed to understand the population structure of A. bicolor bicolor across the Indian Ocean. This study has revealed a unique population structure in A. bicolor bicolor, revealing a clear genetic divergence in the eastern Indian Ocean. The molecular evidence calls for further research on the life history, stock assessment and protection of the populations of A. bicolor bicolor in Indonesia and Malaysia.