Over the past century, stock identification, which is crucial for fisheries stock assessments, has relied heavily on intraspecies variations to differentiate management units. In recent years, however, DNA approaches have shed additional light on some aspects of the natural history and ecology of species and stocks. The Crescent Grunter Terapon jarbua is widely distributed across the Indo-Pacific region. It inhabits coastal waters with sandy substrates and tends to aggregate in estuaries. In the present study, we collected 96 individuals of T. jarbua from 10 locations around the shoreline of Taiwan as well as near Guangdong, China. A concatenated sequence (1,753 bp) of mtDNA (cytochrome c oxidase I and cytochrome b) was obtained from all individuals. We discovered two genetically distinct clades (lineages A and B) with different historical demographies occurring sympatrically except at the Yilan collection site, which was dominated by lineage A haplotypes. Connectivity within this region is high according to FST and AMOVA tests. The genetic variation between the two clades is far below the interspecies threshold for the genus Terapon (0.007 versus 0.3156). Therefore, we suggest that they be considered different genetic stocks from a fisheries management point of view and that future stock reassessment should be conducted based on the genetic information provided in this study. Further large-scale sampling is needed to understand the mechanisms that drive genetic partitioning on regional (Northwest Pacific) and global (Indo-Pacific) scales.
Many marine ecosystems are being degraded by ongoing catastrophic impacts caused by human activities such as overfishing and pollution. The former is currently considered the main factor accounting for the accelerated shrinkage of populations and extinction of species (Jackson et al. 2001; Worm et al. 2006). According to Worm et al. (2006), the fisheries of all currently harvested taxa will collapse by the year 2048 if current exploitation rates persist. Drastic declines in population size due to overfishing could result in a loss of genetic diversity, which may in turn impair species' ability to survive in changing environments. This could produce a vicious cycle leading to failure to recover, even smaller population sizes, and the eventual collapse of marine ecosystems (Worm et al. 2006; Cano et al. 2008).
Stock identification is an interdisciplinary science for studying fisheries stocks, through which appropriate scales of harvesting and monitoring are developed to maintain sustainable fisheries and protect endangered species (Begg et al. 1999). The definition of a fish stock has varied with the approach to identifying and classifying stocks and evolved with management needs (Begg and Waldman 1999). Early definitions relied on distinct internal dynamics and fishery characteristics (Russell 1931). Over the past 100 years, they have strongly depended on different methodological approaches designed to evaluate the homogeneity of populations, e.g., measurements of meristic and reproductive traits as well as geographic distributions. However, these methods can overlook stock structure when dealing with cryptic species/ linages that cannot be distinguished by biological and morphological traits. Without proper methodologies for evaluating stock structure, multiple subpopulations of any given species have generally been managed as a single homogeneous unit. Additionally, due to the enormous size of the ocean and the lack of apparent physical barriers, local populations were considered open systems until a few decades ago (review in Swearer et al. 2002). However, recent evidence has shown that environmental and landscape features can greatly affect and shape population structure at relatively small scales in terms of the life history diversification and local adaptation of aquatic organisms (Storfer et al. 2007). Molecular genetics has made remarkable contributions to our understanding of populations, providing novel insight into previously inaccessible aspects of biological disciplines, e.g., new species descriptions (Allen and Erdmann 2012; Liu et al. 2013) as well as species and management unit identification (Liu et al. 2013; Lim et al., in press). This is especially true for specimens lacking morphological and biological differences among populations (Viñas et al. 2014) or having no morphological information, e.g., shark fins and fish filets (Cutarelli et al. 2014). Genetic markers have also been successfully applied to discriminating the stock structure of fishes (Wu et al. 2014; Singh et al., in press). There is an urgent need to conduct research of this type since species might be genetically heterogeneous over their ranges.
The Crescent Grunter Terapon jarbua is widely distributed across the Indo-Pacific region, from the Red Sea to Fiji and Japan (Vari 2001). It inhabits coastal waters with sandy substrates and is euryhaline, with juveniles commonly being found in brackish and fresh waters. It feeds mainly on small fish and invertebrates. Its planktonic larval duration is relatively short and stable, averaging 25.1 d (Lavergne 2012). The Terapon jarbua occurs along the coast of Taiwan, especially in estuaries, and is a popular game fish caught mainly for food. Larvae and juveniles enter estuaries in great abundance during May and November for feeding (Miu et al. 1990). Once fully grown, they return to deeper waters away from estuaries for spawning. Although the commercial fishery is minor in terms of total amount, T. jarbua is a high-priced species in Taiwan, with adults being caught by long lines and trawlers at depths shallower than 50 m. Its reproduction, body deformation caused by warm water discharged from nuclear power plants, and linear alkylbenzene sulfonate toxicity have been investigated for Taiwanese waters (Miu et al. 1990; Wang and Huang 1995; Chang et al. 2010). However, studies focusing on genetic stock structure have not yet been conducted.
In the present study, the cytochrome b (cyt b) and cytochrome c oxidase I (COI) regions of mtDNA were used as molecular markers to examine the level of genetic differentiation and the historical demography of T. jarbua in Taiwanese waters for future management and conservation purposes.
Sample collection and DNA extraction.—A total of 96 T. jarbua were collected from 10 locations around the shoreline of Taiwan, including New Taipei City (TP), Yilan County (YL), Hua-lian City (HL), Hua-lian County (HFB), Tai-Tung County (TT), Kaohsiung County (KH), Chiayi County (CI), Taichung County (TC), Hsinchu County (HC), Penghu County (PH), and a site near Guangdong, China (RP) (Figure 1). Sample sizes from each locality ranged from 5 to 10 individuals. Genomic DNA was extracted from muscle tissue or fin clips using commercial DNA extraction kits (Genomics BioSci and Tech, Taiwan). DNA extracts were diluted in tris-EDTA buffer and stored at -20°C until amplification by polymerase chain reaction (PCR).
Genetic diversity and phylogenetic analyses.—Complete cyt b and partial COI genes were amplified with the following primers: COI: FishF1 = 5′-TCAACCAACCACAAA GACATTGGCAC-3′, FishRl = 5′-ACTTCAGGGTGACCG AAGAATCAGAA-3′ (Ward et al. 2005); cyt b: Thr33 = 5′-TAACCTTCAGCCTCCTGCTTACA-3′, Glu31 = 5′-GTGA CTTGAAAAACCACCGTT-3′ (Davis et al. 2012). Polymerase chain reactions were run in 30-µL amounts containing 10–40 ng template DNA, 3 µL 10× buffer, 0.2 mM dNTPs, 1.5 mM MgCl2, 10 mM of each primer, and 0.2 units of Taq polymerase (MDbio, Taipei). The thermocycling profile consisted of initial denaturation at 94° C for 2 min followed by 35 cycles of denaturation at 94°C for 30 s, annealing at 48°C for 30 s, extension at 70°C for 40 s, and a final extension at 72°C for 2 min. The purification of PCR products and sequencing reactions was done by Genomics BioSci and Tech. Nucleotide sequences of forward and reverse strands were determined using an AB I 3730XL automated sequencer (Applied Biosystems, Carlsbad, California). Sequences were assembled and edited using Sequencher version 4.2 software (Gene Codes, Ann Arbor, Michigan). The sequences obtained in this study were deposited in GenBank under accession numbers KP204162-KP204259 and KP152133-KP152230 for COI and cyt b, respectively. Sequences were aligned using CLUSTAL W (Thompson et al. 1994), followed by manual editing with Sequencher 4.2 (Gene Codes). We concatenated the sequences of two loci for all individuals, and the analyses performed in this study were based on this data set. Arlequin 3.5 (Excoffier and Lischer 2010) was used to analyze genetic diversity indexes and index of population differentiation (Φst). Unique haplotypes were quantified and the genetic diversity of populations calculated. Analysis of molecular variance (AMOVA; Excoffier et al. 1992) was performed under two hypothetical groupings and the proportions of variations among regions (Φct), among populations within regions (Φsc), and within populations (Φst) were estimated. One of the groupings was based on the difference in coastal habitats between Taiwan's west coast, which is primarily characterized by muddy or sandy bottoms (RP, KH, PH, HC, TC, CI and TP), and its east coast, which is mainly rocky or reef habitat (TT, HFB, HL and YL). Another grouping included these two areas but also considered a possible genetic break between China and Taiwan, caused by either the currents of the Taiwan Strait or the emergence of the continental shelf under the Taiwan Strait during the Pleistocene Epoch (Wang 1999). Therefore, in this latter grouping RP was separated from the west coast group. Median-joining parsimony analysis was done with NETWORK (Bandelt et al. 1999) on the nucleotide sequence matrix of the 96 concatenated sequences using the program's default settings.
The general time reversible model plus the gamma rate heterogeneity model (GTR+G) was selected by MEGA 6 (Tamura et al. 2013) as the best-fitting substitution model based on the Bayesian information criterion. Maximum likelihood tree searching was conducted in RAxML-HPC BlackBox under CIPRES Science Portal version 3.1. Sequences of closely related species (T. theraps and T. puta) were downloaded from GenBank for estimating interspecies Kimura two-parameter (K2P) genetic distance with MEGA 6 (Tamura et al. 2013).
Additionally, Bayesian phylogenetic reconstructions were created by Mr. Bayes 3.12 (Ronquist and Huelsenbeck 2003) with the K2P model and default priors. Subsequently, two duplicate runs of three heated and one cold Markov chain-Monte Carlo (MCMC) chains were established, and each was initiated from a random tree and run for 1,000,000 generations. The convergence diagnostic was applied and the stop probability was set to 0.01. Trees were then sampled every 100 generations and a consensus tree built for all trees with the exclusion of the first 25% of sampled trees to allow for sufficient burn-in.
Historical demography and molecular dating.—The significance of the F-statistics for population comparisons was assessed using 1,000 permutations. In addition, a pairwise mismatch distribution comprised of the pairwise differences among all haplotypes was performed for the historical demographic test of samples of the two major clades based on the deep divergence found in the phylogenetic and minimum-spanning trees. The distribution is usually multimodal when data comply with a demographic equilibrium. In contrast, a unimodal distribution may indicate a recent demographic expansion (Slatkin and Hudson 1991).
The software BEAST version 1.7.5 (Drummond and Rambaut 2007) was used to estimate the time to the most recent common ancestor (TMRCA) using an MCMC Bayesian approach. Divergence rates of approximately 2.0% and 1.2% per million years for cyt b and COI, respectively, were used to estimate the absolute TMRCA values (Horne et al. 2012; Taillebois et al. 2013). All analyses were performed using the GTR+G model of nucleotide substitution. Two independent MCMC chains were run for 100 million generations, with sampling every 10,000 generations and 10% burn-in of the posterior samples. The effective sampling size parameter was found to exceed 200, implying that there was acceptable mixing and sufficient sampling. The analysis was run five times to test the stability and convergence of the MCMC chains in plots of posterior log likelihoods in Tracer version 1.6 (Rambaut et al. 2014).
A 1,753-bp concatenated sequence of mtDNA COI (612 bp) and cyt b (1,141 bp) was analyzed for 96 individuals obtained from 11 locations (Figure 1) in Taiwanese waters. The nucleotide composition was as follows: 24% adenine, 28.7% thymine, 31.6% cytosine, and 15.7% guanine. In total, 106 variable sites and 45 parsimoniously informative sites were found. Among the 96 individuals sequenced, 63 unique haplotypes were identified. Nucleotide diversity (π) ranged from 0.0013 ± 0.0009 (mean ± SD) to 0.0079 ± 0.0045, and haplotype diversity (h) ranged from 0.8571 ± 0.1371 to 1.0000 ± 0.0625 among locations (Table 1).
The phylogenetic analysis revealed two clades among the 96 individuals, herein named lineage A and lineage B, both having strongly supporting values on the nodes (99 maximum likelihood, 0.99 Bayesian support; Figure 2). Both lineages were found at all locations except YL, where we observed only the lineage A haplotype. The general topology of the median-joining tree corresponded with the maximum likelihood tree, and the two lineages were separated by nine mutation steps. Lineage A showed star-like branches (Figure 3). Pairwise θ estimates based on concatenated sequences ranged from —0.09184 to 0.17288, and only two comparisons (TC with HFB and TC with RP) showed marginally significant genetic structure (Table 2). The results of the AMOVA comparisons of genetic variation revealed that the variance was mainly from within the population level and that there was no significant genetic structure under these two groupings (P > 0.05; Table 3). Mismatch distribution tests on lineage A and lineage B failed to reject the hypothesis of the sudden-expansion model (Figure 3). However, the topologies of the distributions differ, with lineage A being unimodal and lineage B being bimodal. The TMRCA of lineage A and lineage B analyzed using BEAST dated to 0.4305 million years ago (effective sampling size, 5,335.203; 95% credible interval, 0.3049–0.5687 million years ago). The average K2P genetic distance (only COI) among the three Terapon species was 0.3156, and the distance between lineage A and lineage B within T. jarbua was 0.007 (Table 4).
Sampling locations and diversity indices based on concatenated sequences (1,753 bp) in 11 populations of Terapon jarbua from Taiwanese waters. Abbreviations are as follows: N = sample size, h = haplotype diversity, and π = nucleotide diversity.
High haplotype diversity (h = 0.86-1) and low nucleotide diversity (π = 0.0013-1) were found in the 11 populations of T. jarbua. This phenomenon is commonly observed in marine fishes and has been interpreted as a sign of historical population expansion (Grant and Bowen 1998). The results of the mismatch distribution tests for lineages A and B also support the expansion model. Such results have been inferred to be closely related to population expansion after glaciation in several marine fishes (Tzeng 2007; Horne et al. 2008; Reece et al. 2010). In addition, lineage A showed a typical unimodal distribution while lineage B showed a bimodal distribution, which can be interpreted as the populations' having gone through expanding and diminishing sizes and structured sizes, respectively (Excoffier and Schneider 1999). Historical population expansion is also supported by the median-joining network, which shows a bush-like pattern of haplotypes for lineage A.
Population Genetic Structure
During the Pleistocene Epoch, several land bridges formed that enclosed the Sea of Japan, partially enclosed the South China Sea, and partially or fully exposed the areas now occupied by the East China and Yellow Seas as a result of lower sea levels that had a great influence on surface currents (Ujiié et al. 1991; Wang 1999). In addition, current hydrodynamic regimes (such as deep-sea gaps along Taiwan and the Japanese Archipelago) may have had profound effects on the genetic divergence of marine populations (Kojima et al. 2006; Kawane et al. 2008). These changes in environmental and geographical regimes are suggested as the mechanism that influenced the population genetic structures of the fishes in the northwest Pacific (Liu et al. 2007; Hu et al. 2011; Shen et al. 2011). Owing to the influence of historical events, the population genetic structure of T. jarbua would be expected to exhibit distinct genetic lineages derived from sympatric individuals. Owing to the influence of hydrodynamic barriers, it should be closely associated with the distribution of populations. The present study shows that in Taiwanese waters T. jarbua is composed of two genetically divergent lineages. Based on haplotype distribution, we found no geographical gradient or separation across the 11 locations around Taiwan, including one from the coast of China. This result is also supported by AMO VA, which reveals no geographic structure among the 11 locations. Because the distributions of the two haplotype groups overlap, no geographic structure was observed. The genetic partition that we discovered in T. jarbua is unlikely to have been caused by current hydrodynamic regimes such as deep-sea gaps, more likely being due to historical events. A similar result was found in a genetic barcoding study of marine fishes in the South China Sea, with two genetically divergent clades (Zhang and Hanner 2012).
Pairwise θ estimates among Terapon jarbua populations based on concatenated mtDNA haplotype frequencies. The location codes are given in Table 1. Values in bold italics indicate P < 0.05.
Analysis of molecular variance (Schneider et al. 2000) based on concatenated mitochondrial DNA haplotype sequences according to two distinct models for the geographic partitioning of populations.
The population genetic structure of T. jarbua fits the historical scenario. Therefore, we suggest that these two distinct lineages formed in the marginal seas of the northwest Pacific (the Sea of Japan and the South China Sea) during periods of low sea level (Voris 2000), with the estimated TMRCA of lineages A and B being dated to 0.4305 million years ago. After the end of glaciation, the two lineages may have recolonized the same waters around Taiwan and now coexist. However, given the broad distribution of this species, the historical events may have occurred elsewhere. Further studies with large-scale sampling of T. jarbua may provide more information regarding the history of genetic divergence.
Terapon jarbua K2P interspecies genetic distance matrix.
Species Boundary and Connectivity
The existence of cryptic lineages is common in marine environments (Knowlton 1993, 2000; Rocha and Bowen 2008). Such lineages often represent previously undiscernible species, and several new marine fish species have been described based on cryptic lineages discovered in phylogenetic studies (Drew et al. 2010; Liu et al. 2013). Intraspecies divergence is commonly observed across different taxa in the northwest Pacific, including fishes, crustaceans, and mollusks (review in Shen et al. 2011). However, the threshold for determining cryptic species varies depending on the taxon and genetic marker used. Although morphological characters are the major criteria for defining species, genetic distance matrices are commonly used to infer the level of differentiation (Baldwin et al. 2011; Liu et al. 2013). The genus Terapon is comprised of three valid species: T. jarbua, T. puta, and T. theraps. The genetic distance between lineages A and B did not exceed the average interspecies genetic distance for the genus. Therefore, we suggest that these two lineages diverged relatively recently and can be considered two genetic stocks of the same species. Recently, Lavergne et al. (2014) revealed the complex genetic structure of T. jarbua in the wider Gulf of Aden (in the range of tens of kilometers even without a geographic barrier between adjacent locations). In contrast to what we found in our study, connectivity is low there. The high genetic connectivity that we found around Taiwan is possibly due to the shared haplotype between localities as well as the mixed haplotypes of the two lineages at most locations. As it happens, two samples collected from the Gulf of Aden shared a unique COI haplotype that can also be found in the South China Sea (Lavergne et al. 2014). This suggests that the range of population expansion after glacial retreat was not restricted to the northwest Pacific but also extended into the Indian Ocean. Broad-scale sampling is needed to test this hypothesis and to understand the historical and contemporary population connectivity within the Indo-Pacific region.
During the past few decades, mtDNA has become the most effective molecule to use for assessing intraspecific genetic variation and genealogy. A well-known case is the European Anchovy Engraulis encrasicolus, in which extensive genetic variation was discovered among Mediterranean populations by means of mtDNA analysis (Magoulas et al. 1996, 2006). Based on these studies, the Mediterranean Sea, Bay of Biscay, and Senegal coastal waters (southeastern Atlantic) are characterized by the coexistence of two highly divergent mtDNA clades. This mixed-stock phenomenon is common in fishes across their distribution, which creates some difficulties in fisheries management (Cao et al. 2014). Genetic approaches allow the contributions of different populations to mixed-population fisheries to be inferred (Ruzzante et al. 2006). Even though these approaches have been heavily applied in the past decade to resolve intraspecies issues such as the genetic structure and phylogeography of marine organisms, genetic information has only been used to reassess life history traits between genetic stocks in a few cases (von der Heyden et al. 2014).
Recent studies have indicated that late in the year the population size of California coastal Chinook Salmon Oncorhynchus tshawytscha is larger than that of the Klamath River population (Satterthwaite et al. 2014). This indicates that there are two genetic stocks that act differently during their life span. The Terapon jarbua populations in Taiwanese waters have previously been considered to be a single stock, and several biological studies, including those of reproductive biology (Miu et al. 1990; Chang 2006), feeding behavior (Lin and Lee 1991), and growth (Wei 1995), have been conducted under this assumption. Therefore, the biological parameters of T. jarbua (such as growth rate, age at maturity, and fecundity) that have been based on the single-stock assumption should now be reexamined based on genetic stock information to eliminate the effects of mixing stocks.
To summarize the present study, we found (1) that two genetically distinct stocks of T. jarbua (based on a concatenated sequence of COI and cyt b) occur sympatrically in Taiwanese waters, (2) that the average genetic distance between the two lineages is much lower than the average interspecies genetic distance for the genus Terapon (which indicates that the divergence between lineages was relatively recent), and (3) that reassessment of biological traits should be conducted based on the genetic information provided in this study. Further genetic studies with rangewide sampling using both mitochondrial and nuclear DNA markers may provide more information on the population genetic structure of this species.
We are grateful to the anonymous referees for their constructive comments. We thank Ming-Tai Chou, Zen-Wei Lin, Tzu-Wei Liu, and Chih-Hao Cheng for their assistance in sampling and Chien-Yao Lu for providing the specimens from Penghu, and we specially thank Buford Pruitt Jr. for English editing. This study was supported by grants to T.Y.L. from the Ministry of Science and Technology (NSC 102-2621-B-110-001) and Asia-Pacific Ocean Research Center, National Sun Yat-Sen University. Shang-Yin Vanson Liu and I-Hsiang Huang contributed equally to this article.
 This is an Open Access article distributed under the terms of the Creative Commons Attribution License ( http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. The moral rights of the named author(s) have been asserted.