Background

Hirundichthys oxycephalus is an important flyingfish resource in eastern Taiwan and northwestern Japan. A substantial catch decline in Taiwan has caused serious concerns on stock status of the fish, prompting the government to impose a set of regulations on flyingfish egg fishery since 2008. However, the regulations were set in a precautionary manner, without considering the fundamental understanding of the population genetic structure. This study aims to investigate the population genetic structure of H. oxycephalus in the region based on mtDNA cytochrome oxidase I (COI) gene and to thus provide scientific information for sustainable management of the resource.

Results

Tissue samples (156) from six localities of eastern Taiwan and western Japan were collected, and 616 bp of mtDNA COI gene were sequenced. Seventy haplotypes were determined, and the haplotype diversity and nucleotide diversity were estimated as 0.93% and 0.57%, respectively. Results of various statistical analyses suggested that the genetic differentiations among the six localities were small and most variation occurred within populations, indicating a high gene flow in the region with undergoing population expansion. Although the study showed that the fishes were genetically divided into two groups, the support was low and the separation was not geologically evident.

Conclusions

The study revealed two groups of H. oxycephalus in the northwestern Pacific Ocean. However, due to high gene flow, an association of either group to a spatial distribution was not observed, and so the two groups may be considered as one population. Thus, the results favored the conclusion that H. oxycephalus from eastern Taiwan and western Japan belong to the same population and, consequently, that the management unit of the current regulations only covering eastern Taiwan does not match the spatial structure of the population. Rather, the results suggest that joint efforts from countries within the population boundary are necessary to maintain a sustainable exploitation.

The flyingfish family Exocoetidae is a group of small fish comprising 52 species and is distributed globally in the tropical and subtropical areas of three oceans (Nelson 2006; Chang et al. 2012a; Xu et al. 2013). Exocoetids are low-tropical level species (Wu et al. 2010) and are a food source for dolphinfish, tuna, and swordfish (Oxenford and Hunte 1999; Wu et al. 2006). They are known as flyingfish due to their ability to glide in the air over the water using large and elongated pectoral and asymmetric caudal fins (Davenport 1994). This uncommon flying behavior serves as an escape mechanism from predators and is uniquely present in extant fish species of the Exocoetidae group as a homoplasy of the extinct family Thoracopteridae (Davenport 1994; Kutschera 2005; Xu et al. 2013).

The flyingfish is economically important in the Caribbean and the western Pacific Ocean (Rennie 2002; Potts et al. 2003). It is the national fish and one of the national cuisines of Barbados (Cumberbatch and Hinds 2013) and provides a commercially important harvest in Japan, especially in the southwestern regions (Masuda et al. 1984). In Taiwan, it supports the stability of the Kuroshio ecosystem as an intermediary in the bioenergetics transfer, the sustainability of artisanal fisheries as a locally consumed dish and major bait for fishing, and the continuity of aboriginal Tao culture as a respectable emblem (Chang et al. 2012a).

At least 29 species of flyingfish in seven genera are reportedly found in the waters off Taiwan and adjacent areas (Chang et al. 2012a). Among these, Hirundichthys oxycephalus (bony flyingfish) is one of the dominant flyingfishes in the Kuroshio Current (Chang et al. 2012b). The fish can be found from south of Taiwan to northwestern Japan (Ichimaru 2007; Chang et al. 2012b); however, it has not been found in the Philippines according to a 3-year sampling program conducted in Batanes (SK Chang and CH Lin, unpublished data). The fish was utilized as the main source of eggs for flyingfish egg fisheries in northern Taiwan and the bait of longline and troll fisheries in eastern Taiwan (Chang et al. 2012a, b). Juveniles of the fish also constituted an important target of boat seine and setnet fisheries in northwestern Japan (Ichimaru 2007). A decline of more than 60% in catches of the fish and its eggs occurred in Taiwan during the 2006 to 2007 period (Fisheries Agency 2010). This decline has caused serious concerns regarding the fish stock status and prompted the Taiwan government to implement a set of regulations on flyingfish egg fishery since 2008, including an annual catch limit. The catch limit regulation is precautionary, however, and has no scientific support regarding the optimal fish stock abundance or, more basically, on the boundary of the fish population that is crucial for effective fishery management (Garcia 2005). Therefore, the population structure of H. oxycephalus needs to be resolved urgently in order to confirm if the population is trans-boundary and is exploited by more than one country.

Genetic analysis is a useful tool for species identification and the determination of population structure which could provide a scientific basis for fishery management (Ward 2000). The technique has been successfully utilized for the sustainable management of many economically and ecologically important fish, such as albacore and bigeye tunas (Ward 2000). For flyingfish, the molecular phylogeny of selected species of the Exocoetidae group has been reconstructed based on the mitochondrial cytochrome b and the nuclear RAG2 genes (Lewallen et al. 2011). Other than this exception, research on population genetics of flyingfish is rare and previous studies all focused on Hirundichthys affinis (Gomes et al. 1998; Gomes et al. 1999, 2000). This is the first study that aims to investigate the population structure of H. oxycephalus in the northwestern Pacific Ocean inferred from mitochondrial cytochrome oxidase I (COI), and to provide a scientific basis for proper fishery management.

From 2008 to 2011, samples of H. oxycephalus were collected from four localities along the eastern coasts of Taiwan, including Keelung (KL), Yilan (YL), Hualien (HL), and Ludao (LD), as well as two localities along the northwestern coasts of Japan, including Nagasaki (NS) and Yakushima (YS). Tissue samples were obtained from 150 adult fish caught by drive-in net, gillnet, and setnet in the six localities, four juveniles from KL, and two eggs from YL using hand-held net (Figure 1 and Table 1). Tissue samples were stored frozen at −20°C before the extraction of DNA.

DNA sequencing

Crude DNA was extracted from the muscle tissue by proteinase K using a Wizard Genomic DNA Purification Kit (Promega, Madison, WI, USA). The mitochondrial COI gene was amplified with the following two primers (Ward et al. 2005): FishF1 (5′TCAACCAACCACAAAGACATTGGCAC3′) and FishR1 (5′TAGACTTCTGGGTGGCCAAAGAATCA3′). PCR was performed using a GeneAmp PCR System 2400 thermal cycler (Applied Biosystems, Carlsbad, CA, USA) with an initial denaturation at 95°C for 4 min, followed by 35 cycles of amplification (denaturing at 94°C for 0.5 min, annealing at 55°C for 0.5 min, and extension at 72°C for 50 seconds) and a final extension at 72°C for 10 min. Autosequencing was performed on ABI 3100 Avanter Genetic Analysis with an ABI PRISM BigDye™ Terminators v 2.0 Cycle Sequence Kit (Applied Biosystems, Carlsbad, CA, USA). The sequences were double checked with forward and reverse strands.

A minimum spanning network was constructed using the Arlequin 3.5.1.3 program (Excoffier and Lischer 2010). In order to measure genetic variation, haplotype diversity (h) and nucleotide diversity (π) were calculated. In order to estimate the nucleotide differences between populations, the Nm and fixation indexes (Fst) were calculated. All of the above indices were calculated by DNASP v5. To analyze if the flyingfish from Taiwan and Japan are from the same population, all samples were hierarchically clustered into Taiwan and Japan groups. Then, a nested analysis of molecular variance (AMOVA) (Excoffier et al. 1992) was performed using Arlequin to estimate population differentiation from the genetic variation of each hierarchical level. Furthermore, neutrality tests, including Tajima’s D, a frequency-based test to identify the genes under selection (Tajima 1989), and Fu’s Fs (Fu 1997), a test sensitive to population demographic expansion, were performed using Arlequin. Mismatch distribution was analyzed to examine the model of sudden population expansion (Harpending 1994). Demographic history and the divergence time for the H. oxycephalus COI sequence were further estimated using Bayesian evolutionary analysis sampling trees (BEAST) ver. 1.6.2 (Drummond and Rambaut 2007). Bayesian skyline analysis (which was calculated using BEAST) was used to infer the changes in effective population size (Ne) across time. The analysis was run with 50 million steps in a Markov chain Monte Carlo (MCMC) simulation under the HKY + I model, relaxed molecular clock model (uncorrelated log-normal), and 1.2% nucleotide substitution rate (Bermingham et al. 1997). Other operators were optimized automatically. BEAST under the Yule process model was also used to estimate the time to the most recent common ancestor (TMRCA). The other model parameters were the same with Bayesian skyline analysis. The results were viewed using Tracer v1.5 (Rambaut and Drummond 2003).

Exocoetus volitans was designated as an out-group for tree rooting. In addition to H. oxycephalus, Cypselurus poecilopterus was included in all molecular analyses. An NJ tree was reconstructed with a Kimura 2-parameter model using MEGA 5.02 (Tamura et al. 2011). Bayesian analysis was conducted using MrBayes v3.1.2 (Ronquist et al. 2012) with the GTR + I + G model as suggested by MrModeltest (Nylander 2004) and partitions based on codon position (first, second, and third). All parameters, with the exception of topology and branch length, were allowed to vary independently with sampling for 27 million generations (two concurrent analyses, nruns = 2; three heated chains, nchains = 4; chain temperature 0.1; sample frequency 1,000; burn in = 7 million generations). Convergence was checked using the software Tracer v1.5 (Rambaut and Drummond 2003) as well as the SUMP and PLOT commands of MrBayes. Maximum parsimonious analysis was analyzed using TNT (Goloboff et al. 2008) by the heuristic algorithm with zero-length branches permitted. Reliability of NJ and MP trees were examined by bootstrapping with 1,000 replications (Felsenstein 1985) using MEGA 5.02 and TNT, respectively. Bayesian posterior probability was calculated by using the MCMC method with MrBayes v3.1.2.

A total of 616 bp of mtDNA COI gene were sequenced from the 156 H. oxycephalus tissue samples. There are 62 variable sites, including 28 parsimoniously informative sites and 34 singletons. The mean nucleotide diversity (π) is 0.57%, with values ranging from 0.52% for HL to 0.63% for YL (Table 1). Seventy haplotypes were determined, and the mean haplotype diversity for all samples is 0.93. Two of them were found in almost every locality in Taiwan and Japan. Furthermore, they were shared by 37% of the total individuals; 13 haplotypes were shared by at least two localities and 28% of the total individuals; and the remaining 55 haplotypes were found in only one locality. The frequency is 35% of the total individuals. The overall mean genetic distance of H. oxycephalus is 0.58%. The minimum spanning network of all samples revealed two distinct groups. Group I comprises most samples (71%) collected from all localities, while group II was smaller (29%) and comprised of samples from all localities except YS (Figure 2A).

All 15 Fst values for the six localities were smaller than 0.15; more than half of them (8) were 0, indicating high gene flow of H. oxycephalus in the northwestern Pacific Ocean. The AMOVA results (Table 2) did not show differentiation between Taiwan (HL, KL, LD, and YL) and Japan (NS and YS; ФCT = −0.0306, p = 0.1), and the variation was −3.0613%. Most variation was within populations (98.92%) and the ФST value was 0.0108 (p = 0.006).

In addition, the mean Tajima’s D and Fu’s Fs values over all samples and three localities (e.g., NS, YS, and LD) were negative, indicating a deviation from neutrality with population expansion (Table 1). Mismatch distribution analysis of all samples of H. oxycephalus revealed a unimodal shape. The sum of squared deviation (SSD) and Harpending’s raggedness index (RI) were not significant (SSD = 0.0076, p = 0.349; RI = 0.0221, p = 0.514), which indicates that the observed results fit with the growth expansion model (Figure 3). Both group I and group II showed similar results for the mismatch distribution analysis with the whole group of H. oxycephalus (Table 1). The Bayesian skyline plot for total H. oxycephalus revealed Pleistocene demographic expansion, which occurred between 40,000 and 140,000 years ago (Figure 4). TMRCA for total H. oxycephalus was 0.531 million years ago (HPD 95% confidence interval was 0.317 to 0.743 million years ago).

Bayesian skyline plot for totalHirundichthys oxycephalusin the northwestern Pacific. The Bayesian skyline plot is derived from mtDNA COI sequences, where the x-axis is time in years and the y-axis is the product of effective population size (Ne) and generation time (T). The black line indicates the median estimate, and dashed lines indicate the 95% highest posterior density (HPD) region.

The neighbor-joining tree (Figure 2B) showed that the haplotypes of H. oxycephalus were divided into two groups with high bootstrap value but low genetic distance (0.9%). Despite varying in abundance, haplotypes from most localities were shared by the two groups, with the exception that those from YS were only present in group I. Bayesian and maximum parsimony trees showed similar topology, but group I was not recovered in these two analyses.

In this study, we used a 616-bp fragment of the mtDNA COI gene to analyze the population genetic structure of H. oxycephalus. The mean haplotype diversity of the six localities was high, but nucleotide diversity was low (Table 1). Marine fish populations can be classified into four categories according to the different combinations of values for nucleotide diversity and haplotype diversity (Grant and Bowen 1998). In H. oxycephalus, the high nucleotide diversity and low haplotype diversity seem to fit the population of the second category, i.e., with a scenario of population bottleneck followed by a rapid population growth and accumulation of mutations. Negative and significant values of Tajima’s D and Fu’s Fs tests, as well as the unimodal shape of mismatch distribution and the Bayesian skyline plot (Figure 3; Table 1), also support the above classification. During the last glacial maximum (LGM), the sea surface temperature was maximally 3°C cooler than it is today (Stott et al. 2002). The cooler temperature could have further limited the distribution of tropical and subtropical species, such as flyingfish, and thus eventually reducing the population size. The fit to an expansion model might be due to its recent colonization from warm refuge after Pleistocene cooling had occurred (Díaz-Jaimes et al. 2006). A similar episode was also observed among dolphinfish, yellowfin tuna, and whiskered velvet shrimp (Díaz-Jaimes et al. 2006; Dammannagoda et al. 2008; Chu et al. 2012).

According to the Bayesian skyline plot for H. oxycephalus, demographic expansion was estimated to occur during the Pleistocene period (Figure 4). Similar observations were noted for some marine fauna, such as mitten crab, mottled spinefoot, and ice goby in the western Pacific, whose Bayesian skyline plots also indicated that the populations had experienced Pleistocene demographic expansion (Xu et al. 2009; Ravago-Gotanco and Juinio-Meñez 2010; Kokita and Nohara 2011). In addition, the TMRCA estimates for all samples of H. oxycephalus indicated that the expansion occurred approximately in the middle Pleistocene period, which is much later than the flathead mullet in this region which occurred in the early Pleistocene period or even in the Pliocene period (Shen et al. 2011).

The low genetic differentiation among six localities in the northwestern Pacific Ocean, along with the variation mostly found within populations rather than among them (AMOVA, Table 2), implies a high gene flow of H. oxycephalus among populations. The strong and constant Kuroshio Current might facilitate the gene flow of species distributed within the northwestern Pacific Ocean and result in a weak genetic structure (Hohenlohe 2004). This may account for the low genetic differentiation of H. oxycephalus and the other species distributed in the range of the Kuroshio Current, including albacore (Wu et al. 2009), Japanese eel (Han et al. 2010), Japanese anchovy (Liu et al. 2006), crocker (Han et al. 2008a,b), mackerel (Chiou and Lee 2004), mullet (Liu et al. 2007; Liu et al. 2009; Shen et al. 2011), goby (Iida et al. 2010; Ju et al. 2013), and whiskered velvet shrimp (Chu et al. 2012).

Historical and ecological factors may influence the demography of marine fishes and result in significant population structure (Riginos and Nachman 2001; Yu et al. 2002; Lourie and Vincent 2004; Liu et al. 2007). On the contrary, species with high dispersal potential and longer pelagic larval duration, such as bigeye tuna (Chiang et al. 2008), albacore (Wu et al. 2009), Atlantic mackerel (Nesbø et al. 2000), and billfish (Graves and McDowell 1995), tend to have high level of gene flow and a lack of population genetic structure (Palumbi 1994, 2003; Shulman 1998; Bohonak 1999; Planes et al. 2001; Purcell et al. 2006, 2009). Despite their distant dispersal and pelagic stage, H. affinis of the central western Atlantic can be divided into three lineages with distinct geographic distribution. Furthermore, its population genetic structure is considered to be associated with spawning behavior (Gomes et al. 1999). In the present study, phylogenetic trees and minimum spanning network of H. oxycephalus revealed two distinct clades with minor genetic distance. The subdivision implied that the specimens of H. oxycephalus were composed of two populations. However, the interrelationship of either stock to any sampling locality, the Kuroshio, or other minor currents was not detected. Similar results have also been observed in other marine fishes in the northwestern Pacific (Chiou and Lee 2004; Liu et al. 2007; Han et al. 2008a,b; Liu et al. 2009; Wu et al. 2009; Han et al. 2010; Iida et al. 2010; Shen et al. 2011; Chu et al. 2012; Ju et al. 2013), in which major ocean current patterns are not good predictors of gene flow. On the other hand, historical factors such as Pleistocene glaciations might account for the genetic pattern of H. oxycephalus. Further study that includes more comprehensive samplings throughout its distribution range would facilitate classification of the genetic pattern.

There are two major definitions for a biological population: the ecological definition that reflects mainly the co-occurrence and demographic interaction of individuals in time and space, and the evolutionary definition that is based on genetic structure due to reproductive interactions among individuals (Waples and Gaggiotti 2006; Reiss et al. 2009). If a significant and reproducible genetic differentiation can be detected, the populations could be considered as separate evolutionarily defined populations with demographical independence (Bentzen 1998). Accordingly, the present study suggested that genetically speaking, there were two populations of H. oxycephalus that existed in the studied region; however, the genetic differentiation was not ‘significant’ enough according to the results of the genetic analysis (Figure 2B).

On the other hand, the present study showed that the haplotypes from most localities were shared by the two groups with the exception of those from YS, which were only present in group I (Figure 2). There were only eight specimens available from YS for the present study, and the limited sample size might have inevitably introduced some uncertainties to the results of the analysis. YL, a locality with only seven specimens, also showed similar pattern, i.e., all haplotypes except one were placed in group I (Figure 2). Therefore, if temporarily ignoring samples with possible sampling bias, the results tended to suggest that the H. oxycephalus from eastern Taiwan and from western Japan was a single population. A separate 3-year study on the spatio-temporal distributions of fish density, body length, and gonadosomatic index (GSI) data of the fish from the same region also suggested one population in the region with two possible migration routes (SK Chang, CH Lin, S Jan, T Ichimaru, and CE Chou, unpublished data). These observations suggested that H. oxycephalus from the northwestern Pacific Ocean belonged to one ecologically defined population.

Fish populations are the fundamental biological unit. Ignoring the congruence of spatial scales between the population structure of fish species and management units can result in reduced productivity and local reduction of population (Kenchington et al. 2003; Waples and Gaggiotti 2006; Worm et al. 2006). Currently, the only known spawning ground of H. oxycephalus is north of Taiwan (Chang et al. 2012b). The spawning activities of the fish take place from May to July, which is the major spawning season based on the back-calculation performed in a daily growth increments study (Chang et al. 2012a). A reproductive biology study on H. oxycephalus suggested that the fish was a batch spawner with asynchronous oocyte development patterns (CH Lin, YT Tsuei, and SK Chang, unpublished data). Based on the assumption of one ecological population, there might be another unknown distinct spawning ground for the same population during the second spawning season in winter time. Currently, there are no implemented management measures on this species in Japan, and the measures implemented by Taiwan are only for flyingfish egg fishery in eastern Taiwan. Although the hypothesis of either one or two populations has not been confirmed, it is evident that the current management scheme covering only the eastern waters of Taiwan mismatches with the spatial structure of the biological population. To properly and effectively manage the resource of H. oxycephalus, negotiations and joint efforts between the Taiwanese and Japanese governments are needed.

Study of population structures is a significant step toward realizing the goal of management of fishery and conservation of the fish resources in their natural population (Abdul-Muneer 2014). This approach has been previously applied to H. affinis in the central western Atlantic (Gomes et al. 1999). To our knowledge, this is the first study to investigate the population genetics of flyingfish in the northwestern Pacific Ocean, particularly with a focus on an economically valuable species, H. oxycephalus, distributed among East Asian countries (Shao 2009). Neutrality tests and mismatch distribution analysis based on the fragment of the mtDNA COI gene of H. oxycephalus collected from the northwestern Pacific Ocean suggested that H. oxycephalus was undergoing population expansion, probably due to recent colonization from warm refuge after Pleistocene cooling (Díaz-Jaimes et al. 2006). AMOVA showed that most variation was within populations, implying that high gene flow of H. oxycephalus occurred among populations, which was probably facilitated by the Kuroshio Current. Phylogenetic analyses and minimum spanning network showed that the fish were genetically divided into two groups. However, the genetic divergence between the two groups was low, and the separation was not geologically evident. On the other hand, based on the genetic data of the present study, which is supported by an independent ecological study (SK Chang, CH Lin, S Jan, T Ichimaru, and CE Chou, unpublished data), H. oxycephalus from eastern Taiwan and from northwestern Japan might be considered as one population. Effective fishery management of H. oxycephalus requires joint efforts from countries within its population boundary in order to maintain sustainable exploitation.

Acknowledgements

The authors are thankful for the special efforts from Mr. Chung-Hui Lin in fish sampling and Mr. Tsung-Wei Lin in molecular analysis. Constructive comments from Dr. Tzong-Der Tzeng of the Shu-Te University on the early work and from the two anonymous reviewers on the manuscript are much appreciated. This study was financially supported by the Ministry of Science and Technology (NSC99-2611-M-110-017, NSC100-2611-M-110-016, and NSC101-2611-M-110-013) and partially by the Asia-Pacific Ocean Research Centre, National Sun Yat-sen University.

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

CEC conducted the analyses and prepared the draft. TYL provided the technical text and discussion and made contributions equal to those of CEC. HWC commented on and edited the draft. SKC provided the idea for the study and text on fisheries, management and conclusions and edited the manuscript. SKC is also the primary corresponding author of this manuscript. All authors read and approved the final manuscript.

Ichimaru T (2007) The life cycle of three species of flyingfish in the north waters of Kyusyu and the recruitment of young flyingfish to the fishing ground. Bull Nagasaki Prefectural Inst Fisheries 33:7–110Google Scholar

Rennie J (2002) Review of the social and economic status of the flyingfish fishery in Grenada. In: FAO/WECAFC (ed) Report on the second meeting of the WECAFC Ad Hoc Flyingfish Working Group of the Eastern Caribbean. Bridgetown, Barbados, 8–12 January 2001. FAO Fisheries Report No. 670, 156 pp., FAO, Rome, pp 103–157Google Scholar

Shulman MJ (1998) What can population genetics tell us about dispersal and biogeographic history of coral-reef fishes? Aust J Ecol 23(3):216–225, doi:10.1111/j.1442-9993.1998.tb00723.xView ArticleGoogle Scholar

Waples RS, Gaggiotti O (2006) What is a population? An empirical evaluation of some genetic methods for identifying the number of gene pools and their degree of connectivity. Mol Ecol 15(6):1419–1439, doi:10.1111/j.1365-294X.2006.02890.xPubMedView ArticleGoogle Scholar

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited.