Pelagic larval duration, growth rate, and population genetic structure of the tidepool snake moray Uropterygius micropterus around the southern Ryukyu Islands, Taiwan, and the central Philippines

Sample collection

One hundred and ninety-five juvenile and adult U. micropterus were collected by dip net and hand-lining at rocky intertidal zones along Ishigaki Island, Japan (n = 14), six Taiwanese sites (n = 142), and Badian, the Philippines (n = 39) during 2014–2016 (Table 1; Fig. 1). Shitiping (n = 36), Jihui (n = 32), and Green Island (n = 20) are located on the east coast of Taiwan and are influenced by the strong, constant Kuroshio Current that flows northward year-round. Checheng (n = 31), Wanlitong (n = 11), and Liuqiu (n = 12) are located on the southwest coast of Taiwan, in a more complex ocean environment that is affected by numerous water masses on a seasonal basis (Shaw, 1991; Farris & Wimbush, 1996; Hu et al., 2010). Ishigaki is also affected by the Kuroshio Current, and Badian is a relatively closed environment located on Cebu in the central Philippines.

Otolith preparation and analysis

Left sagittal otoliths were extracted under stereo microscope, cleaned with deionized water, dried at 55 °C overnight, embedded in Epofix resin, and fixed on microscope slides. The prepared otoliths were grounded along the sagittal plane with 2000 and 2400 grit sandpapers until the core was revealed on the surface, and then polished until smooth with 0.05 µm alumina powder. The ground otoliths were photographed and the growth increments were counted under a compound light microscope (BX-51; Olympus, Tokyo, Japan) from the first feeding check (FFC) to the growth check (GC) (Fig. 2). The otolith growth increments during the leptocephalus stage of U. micropterus were distinguishable and the narrowest rings (approx. 1 µm) were larger than the resolution limitation of the compound light microscope. The GC was defined by Ling, Iizuka & Tzeng (2005) as the prominent check at which the growth increments change from a circular to a radiating pattern, accompanied by a lowered Sr/Ca ratio. This prominent check has been assumed to be associated with the beginning of metamorphosis in leptocephalus and has been used in several studies of marine eel species (Ling, Iizuka & Tzeng, 2005; Lee et al., 2008). Otolith growth increments from the first feeding ring to the GC were expressed as T_GC (i.e., PLD). Meanwhile, the radius of GC was measured along the longest axis and divided by T_GC to calculate the mean otolith growth rate (µm increment⁻¹, as G_GC). Some sub-increments between the wide increments before GC may be ignored, and the blurred rings before FFC were excluded from the count, which may have led to a slight underestimation of the actual number of increments. Fourteen otoliths were randomly chosen for the Sr/Ca ratio analysis to assist in the judgment of GC. Polished otoliths were coated with a layer of carbon and analyzed by an electron probe microanalyzer (EPMA, JXA-8900R; JEOL, Peabody, MA, USA). The Sr/Ca ratios were measured from the core to the edge of each otolith under electron beam conditions of 15 kV and 3 nA, beam size 5 ×4 µm, and 10 µm of spot intervals. Since the otolith growth increments of anguillid species are typically deposited daily in the early leptocephalus and glass eel stages (Sugeha et al., 2001; Shinoda et al., 2004), the otolith increments counted in this study were assumed to be deposited daily (i.e., T_GC = days and G _GC = µm day⁻¹).

The daily ages and otolith growth rates of U. micropterus at different latitudes were divided into three groups (Ishigaki, Taiwan, and Badian) for statistical analysis. The Taiwanese sampling sites were additionally divided into eastern and southwestern groups, representing different oceanic current conditions at similar latitudes, to test whether ocean currents affected the early life characteristics of U. micropterus despite the small geographic scale. The statistical differences in daily ages and mean daily otolith growth rates between groups were tested by analysis of variance (ANOVA) and the post hoc Tukey HSD test in R (R Development Core Team, 2013). The mean increment widths of sagittal otoliths for intervals of every three rings were also estimated in ImageJ (Abràmoff, Magalhães & Ram, 2004) along the longest axis. The percentages of daily ages for five-day intervals are shown in bar charts.

Population genetic analysis

DNA was extracted from muscle tissue using a Qiagen DNA extraction kit (Qiagen, Hilden, Germany) following the manufacturer’s protocols. Polymerase chain reactions (PCRs) were run in a total volume of 50 µL, including 6 µL of TaKaRa (http://www.clontech.com) 10 × buffer, 4 µL of 2.5 mM dNTPs, 4 µL of 10 µm of each primer, 0.25 µL of TaKaRa Ex Taq DNA polymerase, 6 µL of template DNA at 50 ng/ μL, and 25.75 µL of deionized water. The fragments of cytochrome b (cyt b) 680-bp and cytochrome oxidase subunit I (COI) 656-bp were respectively amplified using the primers cyt b: L14725 (5′-GTG ACT TGA AAA ACC ACC GTT G-3′) (Song, Near & Page, 1998) and H15573 (5′-AAT AGG AAG TAT CAT TCG GGT TTG ATG-3′) (Taberlet, Meyer & Bouvet, 1992); and COI: FishF2 (5′-TCG ACT AAT CAT AAA GAT ATC GGC AC-3′) and FishR2 (5′-ACT TCA GGG TGA CCG AAG AAT CAG AA-3′) (Ward et al., 2005). The annealing temperatures of cyt b and COI were 47 °C and 50 °C, respectively. The thermal profiles of PCR were 94 °C for 5 min, followed by 37 cycles of 94 °C for 1 min, annealing temperature for 45 s, and 72 °C for 1 min, with a final extension at 72 °C for 10 min. The quality of PCR products was checked by electrophoresis with 1.5% agarose gel and then purified using a Macherey-Nagel purification kit (http://www.mn-net.com) according to the manufacturer’s protocols. DNA sequences were generated by an ABI 3730 automated sequencer at the Center for Biotechnology, National Taiwan University. Sequences were assembled and edited manually and aligned using MEGA version 6.0 (Tamura et al., 2013).

Sequences of cyt b and COI were concatenated as a single genetic marker and the analyses that followed were based on this data set. The genetic diversity indexes of haplotype diversity (h) and nucleotide diversity (π) were calculated in DnaSP version 5.0 (Librado & Rozas, 2009) according to Nei (1987). Pairwise Φ_ST comparisons among sampling sites and among groups with different classes of T_GC were estimated in Arlequin version 3.5 (Excoffier & Lischer, 2010), and 10,000 permutations were used to estimate the departure from the null hypothesis of genetic homogeneity. The statistical significance of pairwise Φ_ST values was adjusted with Bonferroni correction (Rice, 1989) for multiple comparisons. The hierarchical levels of genetic diversity were tested through analysis of molecular variance (AMOVA; Excoffier, Smouse & Quattro, 1992), and the proportions of variations among groups (Φ_CT), among populations within groups (Φ_SC), and within populations (Φ_ST) were calculated in Arlequin; 10,000 permutations were used to estimate statistical significance. Two hypothetical grouping treatments were used for AMOVA: (1) based on the three different latitudinal groups of Ishigaki, Taiwan, and Badian; and (2) based on the sampling sites that were associated with three different ocean current conditions, namely the Kuroshio Current system (Ishigaki, Shitiping, Jihui, and Green Island), mixed effect by numerous water masses (Checheng, Wanlitong, and Liuqiu), and the interior current systems of the Philippine archipelago (Badian). The minimum spanning network (MSN) of haplotypes was built using Arlequin version 3.5 and HapStar version 0.7 (Teacher & Griffiths, 2011) to connect haplotypes based on the minimum differences between sequences.

Otolith microstructure and microchemistry

The otolith microstructure of U. micropterus was similar to those of other marine eels (Correia et al., 2004; Ling, Iizuka & Tzeng, 2005; Lee et al., 2008). After polishing, the otolith core became a hole surrounded by a thick ring, referred to as a hatch check (HC; Fig. 2). The first feeding check (FFC) was assumed to form when yolks were absorbed completely and the leptocephali began to ingest external food. There were three to five blurry increments between HC and FFC in some individuals. The increments beyond FFC were circular, and the increment widths gradually increased to a peak of 1.1–2.0 µm at approximately the 10th to 20th increments, followed by a gradual decrease to a minimum of 0.5–0.9 µm. Then, the growth increment width abruptly increased to 1.5–4.0 µm by three to 28 increments and formed a profound growth check (GC). The growth increments after GC were wider (5–15 µm), diffused, and radiative. The Sr/Ca ratios of 14 U. micropterus fluctuated between 3 and 16 ×10⁻³ from the core to the GC, with no apparent pattern. The Sr/Ca ratios then dropped rapidly, accompanied by the appearance of the GC in all but one individual (Fig. 3).

Genetic data

One hundred and thirty-two haplotypes from 1,336 bp concatenated mtDNA sequences from 195 U. micropterus were identified (GenBank accession number MF190188–MF190364). In total, 179 polymorphic sites, 112 parsimony informative sites, and 67 singleton variable sites were found. Haplotype diversity (h) and nucleotide diversity (π) ranged from 0.9921 to 1 (average = 0.9938) and 0.0060 to 0.0073 (average = 0.0065), respectively (Table 1). The minimum spanning network showed many unique haplotypes, with only 28 of the 132 haplotypes shared by more than one individual. These unique haplotypes were connected to the center haplotype that occurred in Ishigaki, Shitiping, Checheng, and Badian. The most common haplotype consisted of seven individuals from all locations except Liuqiu and Wanlitong. Closely related haplotypes consisted of individuals from distinct regions, revealing no obvious geographic pattern (Fig. 4).

The pairwise Φ_ST values among sampling sites ranged from −0.025 to 0.121. Liuqiu revealed low but significant genetic variations with all sampling sites except Green Island. Only Liuqiu versus Checheng was statistically significant after the most conservative Bonferroni correction (Table 2). All pairwise Φ_ST values among groups with different classes of T_GC remained low and insignificant (Table S1). The AMOVA results showed that over 99% of variations occurred at the population level under both groupings. Only the Φ_SC of different ocean current conditions revealed low but statistically significant structure (Φ_SC = 0.017, P < 0.05) (Table 3). The results of genetic analysis support the conclusion that U. micropterus in the study areas should be considered genetically homogeneous with weak genetic structure.

Pelagic larval duration and growth rate

The pelagic larval duration represented by the T _GC ranged from 33 to 98 days for the eels examined (Table 4). The maximum T_GC of 98 days occurred in Liuqiu and the minimum of 33 days occurred in Badian. Liuqiu specimens had the longest and most variable T_GC. Ishigaki and Taiwanese specimens had significantly longer T_GC than those from Badian (Tukey HSD, P < 0.01). The otolith growth rate from the first feeding ring to the GC represented by G_GC ranged from 0.91 to 2.40 µm day ⁻¹. The maximum G_GC of 2.40 µm day⁻¹ occurred in Badian and the minimum value of 0.91 µm day⁻¹ occurred in Jihui. Significantly lower G_GC were observed in Ishigaki and Taiwanese specimens compared with Badian (Tukey HSD, P < 0.01).

The eels were divided into three latitudinal groups according to sampling site, namely Ishigaki (n = 14), Taiwan (Shitiping, Jihui, Green Island, Checheng, Wanlitong, and Liuqiu; n = 141) and Badian (n = 39). The respective mean ± SD daily age and otolith growth rate were 54.6 ± 7.7 days with 1.28 ±0.16 µm day⁻¹, 54.6 ± 9.3 days with 1.35 ± 0.18 µm day⁻¹, and 43.9 ± 4.9 days with 1.60 ±0.19 µm day⁻¹. Ishigaki specimens had the lowest growth rate, whereas the Badian specimens had the highest growth rate and shortest mean T_GC (Table 5, Figs. 5A & 6). There were no significant differences in T_GC and G _GC between the Ishigaki and Taiwanese specimens (Tukey HSD, P >0.05), and both were significantly different from the Badian specimens (Tukey HSD, P < 0.001). The percentage of T_GC peaked at 56–60 days (29%) in the Ishigaki specimens, 46–50 days (24%) in the Taiwanese specimens, and 41–45 days (44%) in the Badian specimens (Fig. 5A). U. micropterus tended to have longer T_GC and lower G_GC at higher latitudes.

The Taiwanese specimens were further divided into eastern (Shitiping, Jihui, and Green Island; n = 87) and southwestern groups (Checheng, Wanlitong, and Liuqiu; n = 54), which had mean ± SD days and otolith growth rates of 52.6 ± 8.0 days (1.36 ± 0.16 µm day⁻¹) and 57.9 ± 10.5 days (1.33 ± 0.19 µm day⁻¹), respectively. The specimens collected in southwestern Taiwan had significantly longer T_GC than those collected in eastern Taiwan (Tukey HSD, P < 0.01), whereas there was no significant difference in G_GC (Tukey HSD, P > 0.05) (Table 6). The percentage of T_GC peaked at 46–50 days (30%) for the eastern Taiwan group and 51–55 days (22%) for the southwestern Taiwan group (Fig. 5B). U. micropterus from southwestern Taiwan had longer T_GC.

Early life history traits

The diameters of fertilized eggs and total lengths of newly hatched preleptocephalus larvae of U. micropterus were 3.0–3.5 and 10.0 mm (CH Chen, 1994, unpublished data). The total lengths of U. micropterus leptocephali are approximately 60.0 mm at metamorphosis, based on observations of two Uropterygius spp. in the early metamorphosis stage with total lengths of 51.1 and 56.5 mm (Okiyama, 2014). Therefore, the larval growth rates of U. micropterus with PLD of 33–98 days are likely to be between 0.5 and 1.5 mm day⁻¹. This fast larval growth rate greater than 1 mm day⁻¹ was also observed in four other eel species in the eastern Gulf of Mexico (Gymnothorax saxicola, Ophichthus gomesii, Ariosoma balearicum, and Paraconger caudilimbatus) (Bishop, Torres & Crabtree, 2000).

In the present study, the commencement of leptocephalus metamorphosis was defined by the GC where the Sr/Ca ratio drastically decreased and the otolith growth increment widths abruptly increased, similar to findings in other muraenids and eel species (Marui et al., 2001; Correia et al., 2004; Ling, Iizuka & Tzeng, 2005). Although the biological meaning of GC remains unclear, the first growth increment to this check should contain most of the leptocphalus stage and sufficiently represent the PLD of U. micropterus among sampling locations.

Insignificant pairwise Φ_ST values among groups with different PLDs (Table S1) indicate that the plasticity of the PLD is likely due to individual acclimatization rather than different genotypes. Although PLD variation can be explained by the different growth conditions or birthplaces of the larvae, clearly characteristic PLDs were still observed among the groups in this study, indicating that regional environmental factors may influence tendencies in early life history traits (Searcy & Sponaugle, 2000; Sponaugle & Pinkard, 2004; Bay et al., 2006). Leis et al. (2013) suggested three general classes of factors that might lead to differences in larval dispersal, including biological (e.g., spawning mode and PLD), physical (water movement and habitat fragmentation) and biophysical differences (principally temperature). In this study, U. micropterus at relatively higher latitudes tend to have lower otolith growth rates and longer PLDs compared with specimens at lower latitudes, potentially owing to water temperature. The sea surface temperatures (SSTs) during the spawning season and pelagic leptocephalus stage of U. micropterus (Chen, 1997) were derived from Global Sea Temperature (http://www.seatemperature.org/) (Fig. S1). The SST of Itoman (Japan) is 5 °C lower than that of Guihulngan (the Philippines) in May and December, indicating that the leptocephali in the different study areas experienced different SSTs regardless of whether the eels were self-recruited or transported from other habitats. The effects of temperature on larval development were discussed in rearing experiments that demonstrated positive correlations in Anguilla japonica glass eels and elvers (Fukuda et al., 2009). Leptocephali at lower latitudes with higher SSTs may grow faster and require less time to reach the minimum size for metamorphosis (Reveillac et al., 2008). Because the latitude of Ishigaki is similar to that of eastern Taiwan, neither the mean PLDs nor the otolith growth rates were expected to be different.

Since the effect of SST on larval growth can be excluded from the sampling sites in Taiwanese waters, the relatively longer PLD in the southwestern specimens might be attributable to different current systems or recruitment routes. Larval dispersal route is considered to be strongly influenced by ocean currents (Kim et al., 2007). For instance, newly hatched larvae of C. myriaster can be retained in mesoscale eddies for several months, resulting in longer PLD than expected (Kurogi et al., 2012). The strong Kuroshio Current flows northward along the east coast of Taiwan year-round (Rudnick et al., 2011). When leptocephali competent to metamorphose drift to this area, they can instantly metamorphose and settle. However, the current system in southwestern Taiwan is affected by different water masses seasonally, including the intrusion of the Kuroshio branch, the substitution of South China Sea surface water mass into the Taiwan Strait, and other transient oceanographic events (Shaw, 1991; Farris & Wimbush, 1996; Hu et al., 2010). In addition, a weak anticyclonic eddy with a diameter of 100–200 km was identified in this area (Fig. 1, Fig. S2). The leptocephali might entrain in the anticyclonic eddy or complicated current system and therefore require a longer time to settle along the southwest coast of Taiwan.

Population structure and larval dispersal

Previous studies have shown that the dispersal or retention of larvae may greatly influence gene flow for demersal fishes (Taylor & Hellberg, 2003). During the field collection stage, we found that U. micropterus inhabited shallow water close to the upper tidal zone, which was usually shallower than 1 m and was exposed to air during low tide. These habitats are usually fragmented, but each has an extremely high density of individuals, suggesting that some larvae hatched in the upper tidal zone may remain in nearshore areas and contribute to local habitats. However, the fast-flowing Kuroshio Current is frequently used by the larvae of many eel species for transportation to a wide range of areas (Miller et al., 2002), which may facilitate U. micropterus gene flow. The sampling site at Badian is relatively isolated from the main oceanic current systems around the Philippine archipelago. Compared with the Taiwanese sites in this study, there may be fewer U. micropterus leptocephali drifting away from Badian per generation. It is likely that a handful of migrants have sufficiently contributed to the genetic homogenization among the study regions over timescales of tens to hundreds of thousands of years (Reece et al., 2010).

Nevertheless, it seems that limited larval exchange over time is inadequate to explain the significant genetic divergence of Liuqiu from all locations but Green Island with Φ _ST before Bonferroni correction. Unknown environmental mechanisms that restrict larval dispersal and enhance massive self-recruitment may be a possible explanation. Previous studies have suggested that the self-recruitment of larvae is a common phenomenon for marine fishes and invertebrates, disregarding larval dispersal capabilities (Taylor & Hellberg, 2003; Teske et al., 2015). Self-recruitment of leptocephali is also surmised more often than dispersal in some areas and is attributed to local current systems, semi-enclosed ocean environments, or the swimming ability of the larvae (Miller et al., 2011; Miller et al., 2016). Furthermore, the PLD of Liuqiu specimens was the longest and most variable (46 to 98 days) of all sampling sites, indicating that U. micropterus larvae had more complicated composition and different transportation routes to Liuqiu. The longer PLD implies that some larvae might recruit from places with cooler water or remote populations not included in this study. Further studies with larger-scale sampling among a range of habitats may provide more details on the population genetic structure of U. micropterus.