Variation in Population Structures of Tropical Freshwater Eels, Anguilla Marmorata, A. Bicolor Bicolor and A. Bengalensis Bengalensis, in the Indo-Pacic

Freshwater eels, genus Anguilla, have a distinctive catadromous life history, which could be associated with certain oceanic current systems and offshore spawning sites. Thus, migration and dispersion patterns are believed to be important factors inuencing the population structures of each species. Temperate eel species are well studied, while little research has been conducted on the tropical counterparts that comprise two-thirds of all eel species. The population structures of three tropical species, A. marmorata, A. bicolor bicolor and A. bengalensis bengalensis, which are distributed widely in the Indo-Pacic region, were explored by means of DNA sequence analysis of mitochondrial cytochrome c oxidase subunit I (COI). This study found that A. bicolor bicolor had two genetically distinct populations, and these different populations co-occur geographically in the Indo-Pacic region, while A. marmorata and A. bengalensis bengalensis showed a panmictic-population structure in this region. The populations of A. bicolor bicolor were also found to have lower genetic variability than the populations of A. marmorata and A. bengalensis bengalensis. This study is the rst to explore the population genetic structure of A. bengalensis bengalensis. The present results also suggest plausible dispersion and migration of these tropical species into their continental habitats.


Introduction
The freshwater eels of the genus Anguilla Schrank, 1798 have a unique catadromous life history, consisting of 19 species and subspecies [1][2][3] . These freshwater eels have various ranges of oceanic migration from less than one hundred to more than ve thousand km 4 . Their distribution ranges are almost worldwide except for the continental habitats near to the eastern Paci c and South Atlantic oceans 1 . The growth habitats and/or spawning sites of each eel species either partly or completely overlap with at least one congener. Freshwater eels are categorised into either temperate or tropical species based on their geographical habitats 5 . The temperate and tropical eels consist of 6 and 13 species/subspecies, respectively. The fundamental biology and ecology of temperate eels are generally well understood, but those of tropical eels, which make up two-thirds of all freshwater eel species are very rudimentary 4 .
Each species of freshwater eel has been hypothesized to have its own speci c migration loop, and subsequently, any temporal and/or spatial changes in the migration loop could lead to speciation 6 . This hypothesis could in principle explain the divergence between congener eel species such as Atlantic eels.
The spawning sites of two Atlantic eels, A. anguilla and A. rostrata overlap in the Sargasso Sea 7,8 . However, the differences in their life history characteristics, such as lengths of larval duration, recruitment periods and peaks of spawning seasons, diverge and maintain their migration loops [9][10][11] . Similarly, the spawning sites of a temperate eel, A. japonica, and a tropical eel, A. marmorata, also overlap in the western North Paci c Ocean. However, these eels are allopatrically distributed in Southeast and East Asian countries due to the different optimal temperatures required upon recruitment to continental habitats 12 . These ndings indicate that both biotic and abiotic factors in uence migration loops and divergence in freshwater eels. However, the divergence mechanisms in most of the eel species are still unclear.
The Indo-Paci c region is a suitable location to study the divergence of freshwater eels because there are multiple species, including A. marmorata, A. bicolor and A. bengalensis, which are sympatrically distributed in this region 1 . A. marmorata is the most widespread species 1,13 and is distributed from the east coast of Africa to the Indo-Paci c Ocean 5 . A. marmorata has also been discovered in the Caroline Islands 14 , the Galapagos Islands 15 and the Palmyra Atoll 16 , indicating that its distribution range is much farther east than previously thought. A. bicolor is the second most widespread species and is distributed from the eastern coast of Africa to the Paci c Ocean 1 . The distribution range of A. bengalensis is from the east coast of southern Africa to the west coast of northwest Sumatra of Indonesia and Peninsular Malaysia 1, 3, 17 and is believed to have the third widest geographical distribution. Previous research on the life history characteristics has found that these three species have similar features, such as comparatively shorter larval stages, faster larval growth and longer spawning seasons than temperate species [17][18][19][20][21][22][23] .
A multiple-population structure was previously observed in A. marmorata from the Indian and Paci c oceans [24][25][26] . Furthermore, there were two genetically distinct populations in the Indian Ocean, one on the western side (southwestern Indian Ocean) and another on the eastern side (Sumatra Island of Indonesia) 26,27 . A. bicolor has been differentiated into two subspecies: A. bicolor bicolor in the Indian Ocean and A. bicolor paci ca in the western Paci c Ocean 1 . There was no substantial difference observed in A. bicolor bicolor populations between the western and eastern sides of the Indian Ocean by means of morphological and DNA analyses [28][29][30] . However, two mitochondrial sublineages of A. bicolor bicolor were instead found sympatrically in the Indian Ocean 29,30 . In contrast to A. marmorata and A. bicolor, however, the population genetic structure of A. bengalensis has never been examined, although a study based on morphological observation (total number of vertebrae) found that A. bengalensis has two subspecies, A. bengalensis labiata and A. bengalensis bengalensis, in the western and eastern Indian Oceans, respectively 1, 31 . The genetic population structures of these three species could provide clues to the ecological and environmental factors that might affect population divergence and gene ow in freshwater eels and could provide a better understanding of the distribution and biogeography of eels, whether physical migration barriers in the ocean exist or not. Furthermore, research on the population structure of freshwater eels is necessary for establishing their protection, conservation and management, as well as for understanding the taxonomy and evolution of diadromous shes.
The aim of the present study was to explore the population structures of A. marmorata, A. bicolor bicolor and A. bengalensis bengalensis in the Indo-Paci c region by means of mitochondrial DNA sequence analyses. We also discuss the potential dispersion and migration mechanisms and the importance of their population structures in the region.

Page 4/23
A total of 511 specimens from 23 sites in the Indo-Paci c region were used in the present study (Table 1-3; Fig. 1). From these 511 specimens, 35 specimens were ampli ed and sequenced in this study, the BLAST search showed that 23 and 12 specimens were identi ed as A. marmorata and A. bengalensis bengalensis, respectively, with high identity matches of 99 to 100%, and 476 specimens were obtained from the GenBank database (Table 1-3; Fig. 1). Therefore, by molecular analysis, this study con rmed the occurrence of A. marmorata in Lombok and Sulawesi of Indonesia, and A. bengalensis bengalensis in Penang and Perak of Peninsular Malaysia. Haplotype analysis showed the presence of 23, 23 and 18 different haplotypes in A. marmorata, A. bicolor and A. bengalensis bengalensis, respectively (Table 1-3).   The population structures of the three tropical eels were analysed by using a phylogenetic tree and haplotype network. For A. marmorata, both the MP and NJ trees seemed to show two clades or groups, but they were not supported by the bootstrap probabilities (Fig. 2). The ML tree also did not show the occurrence of these two groups. This suggests that A. marmorata from the Indo-Paci c region formed a panmictic population. Similarly, the haplotype network (Fig. 2) for A. marmorata also suggests a panmictic population structure, as no distinct grouping could be observed, some haplotypes, especially H4, were shared between several localities. It is interesting to note that for H2, this haplotype was shared by Indonesia (Aceh and Bengkulu on Sumatra Island, and Java Island), French Polynesia and Vietnam, which were of considerably distant from each other geographically.
For A. bicolor, both phylogenetic trees and haplotype networks supported the divergence of this species into two subspecies, A. bicolor bicolor and A. bicolor paci ca (Fig. 3). This is supported by the high bootstrap probabilities observed in the MP, NJ and ML trees. For A. bicolor bicolor, two haplogroups or populations could be clearly observed in the MP and NJ trees, with higher bootstrap probabilities shown by the MP tree. However, the ML tree did not support the occurrence of these two haplogroups. Two haplogroups of A. bicolor bicolor were also clearly observed in the haplotype network (Fig. 3). Both groups could be found together in the same localities -Aceh and Bengkulu on Sumatra Island, Java Island and Peninsular Malaysia.
For A. bengalensis bengalensis, no distinct grouping could be observed in the MP tree (Fig. 4). Similarly, neither NJ nor MP showed any well-supported grouping. This suggests that this species formed a panmictic population in the Indo-Paci c region. Similarly, the haplotype network for A. bengalensis bengalensis also did not show any distinct grouping. Furthermore, two haplotypes, H4 and H6, were shared by several localities. Interestingly, three haplotypes, H1, H3 and H7, were more divergent than the other haplotypes.
As shown in Table 4, the panmitic populations of A. marmorata and A. bengalensis bengalensis were found to have higher haplotype diversity and nucleotide diversity than the two populations of A. bicolor bicolor. Haplogroup 2 of A. bicolor bicolor had slightly higher genetic variability than the haplogroup 1. Tajima's test showed negative Tajima's D values in all populations, which suggests recent population expansion, although the D statistic was not signi cant for one of the A. bicolor bicolor populations.

Discussion
The present study explored the population structures of three tropical freshwater eels, A. marmorata, A. bicolor bicolor and A. bengalensis bengalensis, and only a few studies have previously reported these structures. The study showed that there were two genetic populations of A. bicolor bicolor present in the Indo-Paci c region (Fig. 3), whereas A. marmorata and A. bengalensis bengalensis each showed a panmictic population structure (Figs. 2, 4). Given that A. bicolor bicolor can be found on both sides of the Indian Ocean, it could be hypothesized that its population structure should be similar to that of A. marmorata. However, this was not the case. The panmictic population of A. marmorata is also in disagreement with the results from previous studies, whereby A. marmorata was found to have a multiple-population structure [24][25][26] . Although there were inconsistencies between these previous studies, Ishikawa et al. 24 and Minegishi et al. 25 proposed ve populations -North Paci c Ocean, South Paci c Ocean, western Indian Ocean, eastern Indian Ocean (Sumatra Island) and Guam populations -while Gagnaire et al. 26 supported the existence of three divergent populations only -North Paci c Ocean, South Paci c Ocean and western Indian Ocean populations. The present study used mitochondrial COI sequences, but the previous studies examined different DNA markers, which might explain the different results observed. In addition, the present study also did not analyse any A. marmorata specimens from the western Indian Ocean. In this study, A. marmorata from Aceh and Bengkulu on Sumatra Island did not form a distinct eastern Indian Ocean population, and A. marmorata from French Polynesia and Hawaii did not form a South Paci c population, as in the previous studies [24][25][26] . However, this study found that A. marmorata from these localities formed a panmictic population with the specimens from other localities, especially Java, Sulawesi and Vietnam (Fig. 2).
As A. marmorata was previously found to have two genetic populations in the Indian Ocean (western and eastern Indian Ocean populations), these populations were suggested to originate from different spawning areas 24,26 . Two plausible spawning sites of A. marmorata could exist in the Indian Ocean, i.e., in the southwestern Indian Ocean and off Sumatra Island 27 . If the latter is true, A. marmorata should be dominantly recruited to nearby regions. However, on the Java and Sumatra islands, A. bicolor bicolor was found to be the dominant species, accounting for 85% of eels found in these regions 32 ; thus, these regions are not part of the major distribution zone for A. marmorata. Similarly, in northwestern Peninsular Malaysia, which is adjacent to Sumatra Island, A. bicolor bicolor was found to be the most dominant species, accounting for 88.1% of the eels found in this region, whereas A. bengalensis bengalensis was the second-most dominant species at 11.7%, while A. marmorata was present at only 0.2% 33 . These ndings suggest that the transportation of the leptocephali of A. marmorata to the Sumatra and Java islands was likely not from this spawning site in the Indian Ocean but was possibly from the Paci c Ocean. A. marmorata is also distributed in Brunei Darussalam and Sabah of Malaysia on Borneo Island, and on Bali and Lombok islands of Indonesia, which are near the Sumatra and Java islands (Fig. 1).
Oceanic currents ow from the North Paci c and South Paci c oceans through off eastern Borneo Island (Makassar Strait) and off northern Australia to the Sumatra and Java islands 34 (Fig. 1). Therefore, A. marmorata on the Sumatra and Java islands might be transported from spawning site(s) in the Paci c Ocean.
One spawning site of A. marmorata in the Paci c Ocean is located in the western North Paci c 35 (Fig. 1). Brunei Darussalam, Sabah (Malaysia), Thailand, and Vietnam are located in the South China Sea (Fig. 1).
The South China Sea is subjected to seasonal monsoons; during the northwest monsoon, the surface circulation is directed southward, while during the southeast monsoon, the circulation is in the opposite direction 34 . The larvae of A. marmorata might be transported (i) from the northern South China Sea (between the northern Philippines and southern Taiwan) by the North Equatorial Current (NEC) and Kuroshio Current, and/or (ii) from the southern South China Sea through the Celebes (Sulawesi) and Sulu Seas by the NEC and the Mindanao Current (Fig. 1). A. marmorata from Ambon, Bali and Lombok could be further transported by the Indonesian Through ow via the Makassar Strait (Fig. 1). This one spawning site would ensure a homogenous population of A. marmorata in the Indo-Paci c region. However, A. marmorata from different potential spawning sites in the Paci c Ocean could still form a panmictic population through secondary contact.
Anguilla bicolor from the Indian Ocean (A. bicolor bicolor) and the Paci c Ocean (A. bicolor paci ca) were found to diverge from each other (Fig. 3), and the outcomes of our COI gene sequence analysis were similar to those of previous studies using microsatellites, mitochondrial control regions and cytochrome b sequences 29,30 . These previous ndings also showed that A. bicolor in the two oceans has diverged into two subspecies. Furthermore, they also found that on the western and eastern sides of the Indian Ocean, A. bicolor bicolor was observed to have two divergent mitochondrial sublineages 29,30 . These two populations co-occur geographically in the Indian Ocean. The present results also supported the occurrence of two populations of A. bicolor bicolor in the eastern Indian Ocean. These different A. bicolor bicolor populations could occur together within the same habitats on Java Island, Sumatra Island (Aceh and Bengkulu) and Peninsular Malaysia (Fig. 3). Spawning sites of A. bicolor bicolor are thought to be located off west Sumatra Island in the eastern Indian Ocean 32 and off east Madagascar Island in the western Indian Ocean 36,37 . However, no strong evidence is yet available to support these proposed spawning sites of A. bicolor bicolor, unlike the spawning site of A. marmorata in the North Paci c Ocean where its eggs, smaller leptocephali and fully matured adults have been found. These two geographically different spawning areas might explain the occurrence of the two divergent populations of A. bicolor bicolor in the Indian Ocean. The present study supports the hypothetical occurrence of these two spawning sites. Furthermore, the ages of recruitment of this species between the eastern (148-202 days in Java, Indonesia) and western coasts (68-96 days in Mascarene Islands) of the Indian Ocean were substantially different from each other 18,37 . Thus, A. bicolor bicolor could be reasonably assumed to have different populations on the east and west sides of the Indian Ocean.
Anguilla bicolor bicolor in the eastern Indian Ocean has a long duration of the leptocephalus stage (119-171 days) 18 , and the spawning season is year-round 17,20 . These life history strategies might enable A. bicolor bicolor to disperse widely across the Indian Ocean. Seasonal changes in 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 38 , might be capable of mixing A. bicolor bicolor larvae from different spawning sites. Two distinct populations of A. marmorata co-occur geographically, which suggests that they return to their different spawning grounds (Minegishi et al., 2008). Although there is no physical barrier that can prevent migration to these different spawning sites, the speci c choice of spawning site could be imprinted in these eels. Similarly, migration to the different spawning sites in the Indian Ocean might also be imprinted in A. bicolor bicolor. The occurrence of two populations in the eastern Indian Ocean suggests that population divergence in anguillid eels is possible even in the absence of physical barriers.
In contrast to A. bicolor bicolor, A. bengalensis bengalensis had a panmictic population in the Indian Ocean (Fig. 4). This could be because the species has a narrower distribution range than that of A. bicolor bicolor in the Indian Ocean. The distribution range of A. bengalensis bengalensis includes India, Sri Lanka, Myanmar, Andamans, the west coasts of North Sumatra of Indonesia and Peninsular Malaysia 1,31,33,39 , and the species is restricted to the eastern half of the Indian Ocean only, unlike A. bicolor bicolor, which is distributed widely throughout the ocean. Although the present study suggests a panmictic population in A. bengalensis bengalensis, several haplotypes (H1, H3 and H7) were observed to show divergence from the other haplotypes (Fig. 4). Therefore, the occurrence of sub-populations in the Indian Ocean might be possible. Information on the life history of A. bengalensis bengalensis is scarce compared with those of A. marmorata and A. bicolor bicolor 4,17,33 , and only a few leptocephali have been collected from off Sumatra Island and in the eastern Indian Ocean 32,40 . Therefore, further intensive studies on the spawning site, life history and molecular genetic analyses are needed to understand the transportation and dispersion mechanisms of this eel. This study has revealed the unique and different population structures in A. marmorata, A. bicolor bicolor and A. bengalensis bengalensis, which are sympatrically distributed in the Indo-Paci c region. A. bicolor bicolor and A. bengalensis bengalensis were proposed to be transported from spawning site(s) in the Indian Ocean. However, A. marmorata from the eastern Indian Ocean (Sumatra and Java islands) was likely not transported from the Indian Ocean but from spawning site(s) in the Paci c Ocean. The transportation and recruitment mechanisms are likely complicated and different among these three species in the Indo-Paci c region. The molecular evidence from this study calls for further research on the life history, stock assessment and protection of the populations of tropical freshwater eels in the region. To understand the dispersal and transportation mechanisms, further study is needed to investigate the spawning sites of all tropical species together with their population structure, phylogeny and life history analyses in the Indo-Paci c region.

Methods
Eel samples.
A total of 511 specimens from 23 sites in the Indo-Paci c region from India to French Polynesia were used in the present study. From these 511 specimens, 35 specimens were collected from February 2015 to December 2019, and 476 specimens were obtained from the GenBank database (Table 1- (Table 1, Fig. 1). There were 262 specimens of A. bicolor bicolor originating from 3 sites in Indonesia, 2 sites in Malaysia, one site each from India and Myanmar, and 5 specimens of A. bicolor paci ca from 2 sites in the Philippines ( Table 2, Fig. 1). For A. bengalensis bengalensis, a total of 80 specimens originated from 6 sites in India, 3 sites each in Indonesia and Malaysia, and one site each in Bangladesh and China (Table 3, Fig. 1). Our protocols followed the ethical guidelines for the use of animals of Universiti Brunei Darussalam (UBD) and were approved by the animal ethics committee at UBD.

Mitochondrial DNA analysis.
Genomic DNA samples were extracted using the QIAGEN DNeasy Blood & Tissue Kit according to the manufacturer's instructions. Mitochondrial cytochrome c oxidase subunit I (COI) was ampli ed by polymerase chain reaction (PCR) using QIAGEN 2x Taq PCR Master Mix according to the manufacturer's instructions. The PCR primers used were either universal primers 41 or primers that we designed for Anguilla eels, AngF (5'TCA CCC GTT GAT TCT TTT CT3') and AngR (5'CCG ATA GCC ATT ATT GCT CAG3'). The PCR products were puri ed using a QIAGEN QIAquick Gel Extraction Kit. The puri ed COI fragments were sent to a sequencing service provider, where they were sequenced bi-directionally with the same primers. MEGA X 42 was used to align, inspect and edit forward and reverse sequences. The resulting contig sequences were uploaded to the GenBank database with accession numbers MW375913-MW375947. To identify the species, a BLAST search in the GenBank database was used to match the contig sequences with the GenBank reference sequences. MEGA X was used to carry out multiple sequence alignment via ClustalW. The multiple alignment was trimmed at both ends to remove columns with missing data, resulting in lengths of 506 bp for A. marmorata, 445 bp for A. bicolor, and 529 bp for A. bengalensis bengalensis.
MEGA X was used to construct phylogenetic trees via neighbour joining (NJ), maximum parsimony (MP) and maximum likelihood (ML) algorithms. In this study, only the representative MP trees are shown. The NJ tree was constructed using the Kimura 2-parameter model. The MP tree was obtained using the Tree-Bisection-Reconnection method with search level 1, in which the initial tree was obtained by 1000 random addition replicates. The ML tree was constructed with the best-t nucleotide substitution model, the Tamura 3-parameter (T92) for A. marmorata and A. bengalensis bengalensis, and T92 with a discrete gamma distribution (+ G; 5 categories). A heuristic search starting with the initial tree (based on the NJ and BioNJ algorithms) was conducted using the nearest-neighbour-interchange method. For all trees, all codon positions were included, and a bootstrap test with 1000 replicates was also carried out.
Quanti cation of DNA polymorphisms and statistical tests were performed using DnaSP 6 43 . A haplotype network was constructed via the median joining method using Network 10 (www. uxusengineering.com).

Declarations
Map of locations of tropical freshwater eels, A. marmorata, A. bicolor bicolor and A. bengalensis bengalensis in the Indo-Paci c, which were used in our DNA analyses. The location of an offshore spawning area of A. marmorata in the North Paci c Ocean (green ellipse with M) and the oceanic currents (grey lines) in the Indo-Paci c are shown. The base map was downloaded from the OpenStreetMap (open access) at https://www.openstreetmap.org. ITF: Indonesian Through ow, NMC: Northeast Monsoon Current Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors. Each colour represents a sample site. The size of each circle is proportional to the number of samples that belong to each haplotype. Each dash, which appears on the line that connects two haplotypes together, symbolizes one mutational step.

Figure 3
Phylogenetic tree and haplotype network in Anguilla bicolor. (a) Maximum parsimony tree using COI gene sequences of A. bicolor from the Indo-Paci c. Anguilla megastoma was used as outgroup. The consistency index is 0.70, the retention index is 0.81, and the composite index is 0.71 for all sites. The percentages of the bootstrap test for maximum parsimony/neighbor joining/maximum likelihood are shown next to the branches. (b) Haplotype network constructed with A. 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. Each dash, which appears on the line that connects two haplotypes together, symbolizes one mutational step. Figure 4