Plio-Pleistocene sea level and temperature fluctuations in the northwestern Pacific promoted speciation in the globally-distributed flathead mullet Mugil cephalus

Background The study of speciation in the marine realm is challenging because of the apparent absence of physical barriers to dispersal, which are one of the main drivers of genetic diversity. Although phylogeographic studies using mitochondrial DNA (mtDNA) information often reveal significant genetic heterogeneity within marine species, the evolutionary significance of such diversity is difficult to interpret with these markers. In the northwestern (NW) Pacific, several studies have emphasised the potential importance of sea-level regression during the most recent glaciations as a driver of genetic diversity in marine species. These studies have failed, however, to determine whether the period of isolation was long enough for divergence to attain speciation. Among these marine species, the cosmopolitan estuarine-dependent fish Mugil cephalus represents an interesting case study. Several divergent allopatric mtDNA lineages have been described in this species worldwide, and three occur in sympatry in the NW Pacific. Results Ten nuclear microsatellites were surveyed to estimate the level of genetic isolation of these lineages and determine the role of sea-level fluctuation in the evolution of NW Pacific M. cephalus. Three cryptic species of M. cephalus were identified within this region (NWP1, 2 and 3) using an assignment test on the microsatellite data. Each species corresponds with one of the three mtDNA lineages in the COI phylogenetic tree. NWP3 is the most divergent species, with a distribution range that suggests tropical affinities, while NWP1, with a northward distribution from Taiwan to Russia, is a temperate species. NWP2 is distributed along the warm Kuroshio Current. The divergence of NWP1 from NWP2 dates back to the Pleistocene epoch and probably corresponds to the separation of the Japan and China Seas when sea levels dropped. Despite their subsequent range expansion since this period of glaciation, no gene flow was observed among these three lineages, indicating that speciation has been achieved. Conclusions This study successfully identified three cryptic species in M. cephalus inhabiting the NW Pacific, using a combination of microsatellites and mitochondrial genetic markers. The current genetic architecture of the M. cephalus species complex in the NW Pacific is the result of a complex interaction of contemporary processes and historical events. Sea level and temperature fluctuations during Plio-Pleistocene epochs probably played a major role in creating the marine species diversity of the NW Pacific that is found today.


Background
In the marine environment, fluctuations in sea level, water temperatures and sea ice caused by glacial cycles are believed to have had a major influence on species distributions and the population connectivity of marine species [1][2][3]. For instance, lowered sea levels during glacial maxima led to the emergence of land bridges, which fragmented marine ecosystems and isolated populations of aquatic species [4][5][6]. At the same time, temperature fluctuations in the Pleistocene epoch influenced population dynamics by promoting bottlenecks and the loss of genetic heterogeneity. In fact, much of the intraspecific genetic diversity that is exhibited by contemporary species is generally assumed to have derived from such events in the Pleistocene (for a review, see [7]).
The northwestern (NW) Pacific is characterised by marginal seas, which were particularly impacted by Plio-Pleistocene glacial cycles. The lowering of sea level caused the recurrent closure of the Japan Sea, the semiclosure of the South China Sea and the partial or full exposure of the East China and Yellow Seas [8]. The closure of the Japan Sea caused fluctuations in sea temperatures, by halting the influx of the warm Tsushima Current, a branch of the Kuroshio Current (Figure 1) that supplies a large amount of heat to northern areas [9]. As a consequence, seascape dynamics of the region are believed to have profoundly influenced both intraspecific genetic diversity and species diversity of marine organisms, as exemplified by the exceptional diversity found in SE Asia [10][11][12]. This hypothesis has been confirmed by recent phylogeographic investigations in the NW Pacific area ( Table 1). The significance of the lowering of the sea level and glaciations for species diversity has generally, however, not been assessed. The overlapping distribution of multiple divergent lineages is a striking phylogeographic pattern that is common to most of the organisms that have been investigated in this region, and it raises the question of whether some of these lineages are in fact cryptic species. For example, some authors [13][14][15] considered the various divergent lineages to be different species of Pleistocene origin, while others [16,17] have interpreted these lineages as intraspecific polymorphisms, consistent with the existence of isolated populations during glacial maxima and secondary contact at the end of the Pleistocene. Finally, the comparison of phylogeographic information among different marine species inhabiting the same area can provide for conflicting interpretations of the role of decreasing sea level and Pleistocene glaciations, and highlights the difficulty of distinguishing species in the early stages of differentiation from populations experiencing secondary contact after an extended isolation.
The flathead mullet, Mugil cephalus Linnaeus, 1758, is a euryhaline teleost distributed worldwide in coastal waters, lagoons, bays and estuaries between latitudes 42°N and 42°S [18]. It is subjected to intense and unregulated inshore fisheries in the NW Pacific because of the economic value of its roe. There have been drastic declines in the landings of M. cephalus since 1980 (6 865 metric tons in 1980 versus 159 metric tons in 2008; Taiwan Fisheries Agency), which has been suggested to be the result of overfishing and the impact of global warming [19,20]. Thus, a sound management program is urgently needed, which requires precise knowledge about the taxonomy and population structure of this over-exploited species.
Previous studies have shown that up to three highly divergent mitochondrial lineages exist in sympatry in NW Pacific M. cephalus populations [14,17,21]. This high inter-lineage divergence raises questions about the taxonomic status of M. cephalus in the NW Pacific, but studies have not reached consensus on this issue. The inter-lineage divergence of 48% observed at the mtDNA control region probably indicates a species complex, because such divergence coupled with the existence of different reproductive behaviours (resident versus migratory) demonstrates that these lineages are independent [14]. By contrast, an inter-lineage divergence of 5.3 to 6.7% was estimated using cytochrome b sequences, which would argue for the existence of a single species because the genetic divergence between mugilid species exceeds this by a factor of two [17]. However, because the mtDNA genome used in both studies is maternally inherited, it is not possible to determine whether the presence of divergent mitochondrial lineages in the same sample is a result of secondary contact after an extended period of isolation and/or the presence of two sibling species. To our knowledge, only two studies on M. cephalus in this region have used biparentally inherited markers [22,23]. The first study used an allozyme locus, GPI-A (glucose-6-phosphate isomerase), which may be able to distinguish migratory from resident M. cephalus populations [22]. More recently, Liu et al. [23] found a high level of genetic structure in M. cephalus using AFLP, and four populations were identified among a set of 6 samples collected along the Chinese coast, with the southern samples (Hainan Island) being the most divergent. This marked genetic structuring in NW Pacific M. cephalus contrasts with findings from the Gulf of Mexico and the northwest Atlantic, which did not reveal any genetic heterogeneity over a similar geographic scale, and found only one mitochondrial lineage [24,25].
Therefore, given all of the above uncertainties, this study investigated the level of both historical and contemporary gene flow among M. cephalus samples collected in the NW Pacific, using a set of ten microsatellites and one mtDNA (COI) locus. The objective was to estimate the level of independence of the different M. cephalus mtDNA lineages, and reveal the existence of cryptic species. This provided an opportunity to address both the evolutionary forces and reproductive boundary uncertainties among and within the mitochondrial lineages. Finally, we aimed to better understand the importance of Pliocene and Pleistocene geologic and climatic events in the diversification of marine species in the NW Pacific.

Sampling and DNA extraction
A total of 713 sub-adult to adult M. cephalus were collected from 12 locations in the NW Pacific (Table 2, Figure 1). Additional temporal sampling (2005)(2006)(2007)(2008) was conducted at three localities around Taiwan: the coast of Keelung (KL), which is located midway along the M. cephalus spawning migration route from the Eastern China Sea; the offshore waters of Kaohsiung (KS), where spawning occurs; and the Kaoping Estuary (KP), a tropical system located in southwestern Taiwan. Sampling locations, dates of collection, sample sizes for microsatellites analyses, mtDNA lineage (highlighted either by sequencing or multiplex COI haplotype-specific PCR), mean fork length (mm) and gonadosomatic index (GSI) data for the M. cephalus from Taiwan are shown in Table 2. Muscle tissues were preserved in 95% ethanol prior to DNA extraction. Genomic DNA was extracted from muscle tissue using a DNA Purification Kit (Bioman, Taipei, Taiwan), preserved in TE buffer and then quantified and diluted to 1 ng/μl for PCR.

Mitochondrial DNA analysis
Partial mtDNA COI gene sequences (627 bp) were amplified by polymerase chain reaction (PCR) using the universal primers FishF1 and FishR1 [26]. To identify M. cephalus cytochrome b lineages [17], cytochrome b gene was also sequenced for 20 individuals belonging to three COI lineages using the cytochrome b primers MCglu-1.F: GGCTTGAAAAACCACCGTTG and MCcytbR: AGTACTGTGGCAGAGCTTGG. PCR was performed in a Biometra TGradient Thermocycler with a 15-μL reaction volume containing 0.2 μm dNTPs, 1.5 μL of 10× PCR buffer (Bioman, Taipei, Taiwan), 0.5 μm each of forward and reverse primer, 0.2 U Taq DNA polymerase (Bioman, Taipei, Taiwan), and 1.0 μL of template DNA. MtDNA amplification was carried out using the following PCR conditions: 35 cycles of denaturation at 94°C for 15 s, annealing at 58°C for 15 s and extension at 72°C for 30 s after heating at 94°C for 5 min. The PCR products were electrophoresed on a 1.0% agarose gel (Bioman, Taipei, Taiwan) and stained with ethidium bromide for band characterisation via ultraviolet trans-illumination. All sequencing reactions were performed according to the manufacturer's protocol (Applied Biosystems, Foster City, CA, USA) using the same forward primer as used for PCR.
All COI sequences were automatically aligned using MAFFT version 6 [27] and manually corrected. The characterisation of the genetic variability of the COI gene sequences and the number of nucleotide substitutions, including transition/transversions, were conducted using ARLEQUIN version 3.01 [28]. Phylogenetic trees were reconstructed using the Bayesian, neighbour joining (NJ), maximum likelihood (ML) and maximum parsimony (MP) methods. The Bayesian method was applied using the program MrBayes 3.0 [29], and NJ, ML and MP were applied using the program PAUP 4.0 [30]. The transversion model with invariable sites HKY + I + G (I = 0.097, G = 0.244, Ts/Tv ratio = 5.4465) was selected for construction of the phylogenetic tree using Mugil curema as an outgroup (GenBank accession numbers: EU715464, [31]) and the Akaike information criterion (AIC) in MODELTEST 3.7 [32]. Nodes with bootstrap values ≥ 70% were considered well supported [33].
The fixation index (Φ ST ) for all pair-wise comparisons among different populations was calculated to investigate genetic diversity among M. cephalus populations, and a permutation test (10 000 permutations) was performed using ARLEQUIN version 3.01 [28]. Population genetic diversity was measured within each of the populations based on the number of distinct haplotypes, gene diversity (h) and mean nucleotide diversity (π), using ARLEQUIN version 3.01 [28]. The divergence time for the M. cephalus COI gene was estimated using Bayesian Evolutionary Analysis Sampling Trees (BEAST) ver. 1.5.2 [34] with 20 million steps in a Monte Carlo Markov Chain (MCMC) simulation (2 million step burn-in time). COI sequences of M. cephalus samples from disjunct locations in the western Atlantic (South Carolina: GenBank accession numbers: HQ149710 and Florida: GenBank accession numbers: HQ149711) and eastern Pacific (Peru: GenBank accession numbers: HQ149714 and Mexico: GenBank accession numbers: HQ149715) were used for the calibration, and assuming that the divergence of these populations corresponds to the closure of the Isthmus of Panama, 2.8 MY ago [35]. A strict molecular clock was assumed for this final run. The effective sample sizes (ESS) of parameters sampled from the MCMC were > 500 (acceptable ESS is >200). The results were viewed using TRACER 1.4 [36]. Historical demographic/spatial expansions of M. cephalus were explored using two different approaches. Tajima's D [37] and Fu's F S [38] tests were used to test population equilibrium. Deviations from the sudden population expansion model were tested using the Harpending's raggedness index RI of mismatch distributions [39].

Rapid screening of the different mtDNA lineages
A multiplex COI haplotype-specific PCR (MHS-PCR) was designed following the recommendations of previous studies [40,41] to develop a rapid screening method capable of detecting different mitochondrial lineages of M. cephalus. Specific mutations of the different M. cephalus mtDNA lineages were identified using 448 COI sequences of M. cephalus produced in this study. NWP1,2F (5' GCTTTTCCCCGAATAAAT 3') was the forward primer for both lineage 1 and lineage 2, while NWP3F (5'TACTGCCCTAAGCCTACTC 3') was the forward primer for lineage 3. NWP1,3R (5' CGATCTGTTAGGAGTATGG 3') was the reverse primer for lineage 1 and lineage 3, while NWP2R (5' CTCATACGAAAAGGGGTGTT 3') was the reverse primer for lineage 2. PCR was performed using a Biometra TGradient Thermocycler with a 15-μL reaction volume containing 0.2 μm of each dNTP, 1.5 μL 10× PCR buffer (Bioman, Taipei, Taiwan), 0.5 μm forward and reverse primers, 0.2 U Taq DNA polymerase (Bioman, Taipei, Taiwan), and 1.0 μL of template DNA. MtDNA amplification was carried out using the following PCR program: 35 cycles of denaturation at 94°C for 15 s, annealing at 55°C for 15 s and extension at 72°C for 30 s after heating at 94°C for 5 min. The PCR products were electrophoresed in a 1.0% agarose gel (Bioman, Taipei, Taiwan) and stained with ethidium bromide. Lineage 1 had a PCR product of 362 bp, while lineage 2 and lineage 3 had PCR products of 283 bp and 549 bp, respectively (Additional file 1, Figure S1). The rest of the samples that were not sequenced were all rapidly screened to determine their lineages.

Microsatellite analysis
Ten microsatellite loci were used to screen the genetic diversity of the M. cephalus samples [42]. Reverse primers for each locus were labelled with fluorescent dyes (6-FAM, HEX and TAMRA), and multiplex PCR was performed in a 15-μl reaction volume containing 0.1 ng DNA, 1.25 pmole each of the three reverse primers labelled with different fluorescent dyes, 1.25 pmole of each forward primer, 5 mM dNTP, 1.5 mM MgCl 2 and 0.5 U of Taq polymerase (Bioman, Taipei, Taiwan). Amplification was conducted as follows: initial denaturation at 95°C for 4 min, followed by 35 cycles at 94°C, 54-58°C and 72°C for 30 sec each. Locus polymorphisms were screened using an ABI PRISM 377 auto DNA sequencer (Applied Biosystems, Foster City, California, USA). Lengths of microsatellite alleles were determined using a TAMRA-labelled 100 bp standard (Perkin-Elmer, Waltham, Massachusetts, USA).
The calculation of the number of alleles, estimated (H E ) and observed heterozygosities (H O ) for each locality and for all samples, the genetic differentiation index (θ an estimator of f analogues of Wright's fixation indices [43]) among samples, locus-by-locus AMOVA and the exact test for deviations from Hardy-Weinberg equilibrium at every locus for each locality were performed using ARLEQUIN version 3.01 [28] and GENETIX software [44]. The significance level for multiple simultaneous comparisons was adjusted using the sequential Bonferroni technique [45]. The allelic richness of the minimum population size by a rarefaction method and Fis was estimated with FSTAT 2.9.3 [46]. The significance of the differences in allelic richness was tested by the Wilcoxon signed-rank test for paired observations (see [47] for power analysis). The program GENEPOP version 3.1 [48] was used to test for linkage disequilibrium among the ten loci analysed in this study (10 000 permutations, 1 000 dememorisation steps).
The presence of intraspecific genetic structure was tested using the model-based clustering method [49], as implemented in STRUCTURE VER. 2.1. For each value of K, which is the number of genetically distinct populations, the Markov chain Monte Carlo scheme was run with a burn-in period of 10 000 steps and a chain length of 100 000 replicates following the non-admixture model. Twenty runs were performed to evaluate the reliability of the results, with the number of populations being determined from posterior probabilities of K calculated as K = {1~6}. The K values could be incorrectly estimated if the migration rates between populations are not equal, so the values of ΔK were also calculated for each value of K [50]. Individuals were regarded as correctly assigned to a population when their q-value (i.e. the posterior probability) was at least 80% after subtracting the posterior probability assignment of another population [47].

Mitochondrial DNA
A total of 627 nucleotides of the COI gene were obtained from 448 individuals. In these sequences, 54 positions were variable (Additional file 2, Table S1), defining 36 haplotypes (Genbank accession numbers GU260664-GU260697, HQ149082-HQ149083). Phylogenetic trees reconstructed using either COI (Figure 2) or cytochrome b (Additional file 3, Figure S2) sequences [17] all supported the existence of three highly supported monophyletic lineages in NW Pacific M. cephalus. Seventeen COI haplotypes were found in lineage 1, including 15 nucleotide transitions and three nucleotide transversions. Of these, HT2 was the most frequent haplotype, being shared by the samples from Taiwan, Japan, Qingdao and Russia, and it was connected to the other 16 haplotypes (HT1-HT15, HT35-HT36) by one to two steps in a "starlike" haplotype network (Figure 3a). The other 16 COI haplotypes were minor and mainly found in Keelung and Kaohsiung. Fifteen COI haplotypes were found in lineage 2, including 11 nucleotide transitions and two nucleotide transversions. Two major haplotypes (HT17 and HT25) were shared by samples from Keelung, the Kaoping River, Okinawa, Yokosuka, the Pearl River, Shantou and the Philippines. These two main haplotypes had only two nucleotide differences and were mainly found at the Kaoping River site. Lineage 3 consisted of only four COI haplotypes, with the main haplotype (HT32) occurring at the Kaoping River and Hainan sites and in a few samples from Shantou and Keelung. The fixation index (Φ ST ) showed significant differentiation among the three lineages (Φ ST = 0.973, P < 0.001). The pair-wise Φ ST values were 0.960 (P < 0.001) between lineage 1 and lineage 2; 0.956 (P < 0.001) between lineage 2 and lineage 3 and 0.994 (P < 0.001) between lineage 1 and lineage 3. Fifteen variable sites were found within the amplified 627-bp COI fragment (2.5% of sites) between lineage 1 and lineage 2, with approximately 6.7% of the inferred changes being transversions. While 24 and 30 variable sites were found between lineages 2 and 3 (3.8%) and between lineages 1 and 3 (4.8%), respectively, with approximately 21% and 13% of the inferred changes being transversions. The calibrated mean mutation rate of the COI gene in M. cephalus was 1.6% mutations/ site/MY. Therefore, the calibrated divergence time using BEAST was 1.607 MY (HPD 95% confidence interval were 2.383-0.892 MY, ESS = 1475) between lineages 1 and 2 and 4.200 MY (HPD 95% confidence interval were 5.708-2.801 MY, ESS = 876) between lineage 3 and the ancestor of lineages 1 and 2. BEAST also suggested that lineage 3 is much older than the others, even compared to M. cephalus populations present in the eastern Pacific and western Atlantic, and does not share the most recent common ancestor of lineages 1 and 2.
The mismatch distributions of lineage 1 and 3 were unimodal, with almost no pair-wise differences between them (Figure 3b). Lineage 3 exhibited a more uneven mismatch distribution compared to lineage 1. By contrast, the mismatch distribution of lineage 2 was bimodal ( Figure 3b), with peaks at zero and two substitutions. RIs were not significant for any lineage. On the other hand, both Tajima's D and Fu's F S were negative and highly significant for lineage 1. For lineage 2, Tajima's D values were all negative, but only significant for the control region, while Fu's F S were all negative and significant for all genes. Finally, both Tajima's D and Fu's F S were negative, but not significant for the lineage 3 for any gene considered (Table 3).

Microsatellites
A total of 180 alleles with an average of 18.0 alleles per locus were observed across the ten microsatellite loci (Additional file 4, Table S2). All loci were polymorphic, except for locus Mce-8, which was monomorphic in some of the samples. Twenty-one of 180 sample-locus combinations deviated significantly from Hardy-Weinberg proportions after Bonferroni correction. However, the deviations of these loci were not significant in most locations (including Kaoping and Hainan) when individuals were sorted according to one of the three mtDNA lineages, suggesting a Wahlund effect, which is the reduction of heterozygosity in a sample due to the mixture of genetically differentiated populations. No linkage disequilibrium was detected within the M. cephalus populations (P > 0.05).
AMOVA indicated that the genetic differentiation of M. cephalus among 18 spatial and temporal samples was highly significant (θ = 0.060, P < 0.001). A locus-bylocus AMOVA indicated that all loci could detect population differentiation among the samples. However, no temporal variation was observed in allelic frequencies among the KS, KL and KP locations where the temporal samples were obtained. During three to four years of sampling, no genetic heterogeneity was observed among the temporal samples collected at the same location, except in KP, and the genetic differentiation between locations remained stable ( Table 4).
The results from the STRUCTURE assignment test supported three clusters based both on the log probability of the data [L(K)] and the statistic ΔK, as described in [49,50], respectively (Additional file 5, Table S3).
With K = 3, almost all individuals were clearly assigned to one of the three clusters ( Figure 4). The posterior probabilities corresponding to the assignment of individuals to cluster 3 (red bar, hereafter denoted as NWP3) were all unambiguous, while the assignments to one of the two other clusters (green bar for NWP1 and blue bar for NWP2) included several ambiguous individuals (7 over 645 individuals belonging either to NWP1 and NWP2). The individuals from Kaohsiung, Qingdao, the Ariake Sea and Russia were mostly assigned to NWP1; the individuals from the Philippines were all assigned to NWP2, while at least two clusters coexist in Keelung

Genetic variability in mtDNA and microsatellites
There was congruence between the two molecular markers when the results of the assignment test were compared to the COI sequence data (Figure 4). The results indicated that 385 of 390 individuals assigned to NWP1 harboured mtDNA representing lineage 1, 253 of 255 individuals assigned to NWP2 harboured mtDNA of lineage 2, and all 68 individuals assigned to NWP3 harboured lineage 3 mtDNA. This resulted in a successful assignment rate of 99.0%. Although there were seven ambiguous individuals in NWP1 and NWP2 below our assignment criteria, the mean q-values of assignment test for these two clusters were still over 0.99 (Additional file 6, Table S4).
In each location sampled, individuals belonging to the same mtDNA lineages were grouped together to further investigate the genetic structure of the three Operational Taxonomic Units (OTUs). Genetic variability analysis showed that only two of the 30 combinations, loci Mce-3 and Mce-8 of NWP3, deviated significantly from Hardy-Weinberg proportions after Bonferroni correction (Additional file 7, Table S5). The allelic richness per locus based on the minimum sample size (68 individuals of NWP3) was significantly lower at NWP3 than NWP2 (P < 0.001) and NWP1 (P < 0.001) and slightly lower at NWP2 than NWP1 (P < 0.05) for the ten microsatellite loci (H O was lower at NWP3 than the other two clusters). The genotypes at Mce-14 of NWP2 and Mce-10 of NWP3 were all homozygous.
Population genetic differentiation among these three clusters was highly significant (θ = 0.179, P < 0.001). Pair-wise θ comparison indicated that NWP3 was the most differentiated (θ NWP3-NWP1 = 0.296, P < 0.001; θ NWP3-NWP2 = 0.305, P < 0.001), while level of differentiation was lower but still highly significant between NWP1 and NWP2 (θ = 0.115, P < 0.001). In addition, no genetic heterogeneity was detected among samples belonging to any one of the three OTUs, with the , nucleotide diversity (π ± standard deviation) for each group of samples. Tajima's D and Fu's Fs, the corresponding P value, and the raggedness index (RI) of these lineages were also indicated. 1 data from [14], 2 data from [17]. *P < 0.05, ** P < 0.01, ***P < 0.001 Table 4 Temporal pair-wise θ test, using microsatellite loci between different sampling locations and years in Taiwan exception of the Ariake Sea, which had some genetic heterogeneity among the samples of NWP1 (Table 5).

Historical demography and phylogeography
Since the Plio-Pleistocene, successive periods of sea-level regression have directly impacted the species connectivity of the marginal seas of the NW Pacific. This has clearly been demonstrated by previous phylogeographic studies (Table 1), which have identified Pleistocene refugia in all marginal seas. Our results argue for the existence of three mtDNA lineages of M. cephalus in the NW Pacific but only two, NWP1 and NWP2, are of Pleistocene origin. The NWP3 lineage diverged from the NWP1 and NWP2 common ancestor during the Miocene or the Pliocene epoch, and we assume that the origin of NWP3 does not relate to geological events in the NW Pacific because this lineage is also observed in the southwestern (SW) Pacific (in Fiji and New Caledonia Islands; Durand, pers. com.). While it appears likely that the refugium of the NWP2 lineage was located in the South China Sea, the location of the NWP1 refugium is more difficult to assess. The closure of the Tsushima Straits (between Korea and Japan, Figure 1), which isolated the Japan Sea, may have been the powerful vicariant event that was at the origin of the NWP1 lineage. Reconstruction of the history of the Tsushima Current [9] demonstrated that during the interval between 3.5-1.7 MY, the current only flowed into the Sea of Japan periodically and the volume and salinity of water was lower than at present. After 1.7 MY, the Tsushima Current flowed throughout almost all interglacial periods and, more importantly, with a volume and water salinity similar to the present day. This date matches the calibrated dates of divergence estimated between NWP1 and NWP2, suggesting that the NWP1 ancestor was not able to enter in the Japan Sea before 1.7 MY. During each glacial episode, ancestors of the NWP1 lineage that were trapped in the Japan Sea probably experienced several demographic crashes due to decreases in surface water temperatures [51][52][53]. This hypothesis is consistent with the observed genetic diversity and Tajima's D and Fu's F S indices, which indicate bottleneck events in NWP1 (Table 3). Climatic oscillations during the Pleistocene epoch were probably less marked in more southerly latitudes, where the influence of the Kuroshio Current was still exerted, which would explain why NWP2 presents a higher genetic diversity than NWP1. Similar trends can be observed in the genetic diversity pattern of other NW Pacific species, such as Chelon haematocheilus and Odontamblyopus lacepedii sensu lato, suggesting that this scenario may have affected all species inhabiting the area [15,16].  Figure 4 Correspondence between mtDNA lineages and nuclear Structure clusters. First row is the distribution of the 3 mtDNA lineages (green: lineage 1, blue: lineage 2, red: lineage 3) in 18 samples of Mugil cephalus in NW pacific. Second to forth row are the results of the assignment test using STRUCTURE [49] for M. cephalus microsatellite data. Each cluster (K) is designated by a different colour with vertical bars representing individuals and the proportion of a bar assigned to a single colour representing the posterior probability that an individual is assigned to that cluster. Assignment results are shown with K = 3, 4, 5 and 6.
Lastly, the East China Sea, which was exposed during the Pleistocene glaciations, seems to be a post-glacial contact zone for M. cephalus. While the Yellow Sea was recolonised by NWP1 from a northern route of dispersion through the Tsushima Strait, the Ryukyu Islands and East Japan were recolonised by both NWP1 and NWP2 from southern routes that followed the Kuroshio Current (but see below). If the East China Sea experienced a post-glacial colonisation wave from the Japan Sea, it appears that in species such as Chelon haematocheilus and Eriocheir sensu stricto, older colonisation waves were also possible, which allowed the emergence of a specific East China Sea lineage in these species [16,54].

Cryptic M. cephalus species in the NW Pacific
It has previously been hypothesised that M. cephalus is a species complex rather than a unique panmictic unit. To date, however, this has not been proven, due to inadequate sampling schemes with samples separated by thousands of km, or the use of only maternally inherited markers [14,17,25,31,[55][56][57]. The present study demonstrates that M. cephalus is indeed composed of at least three genetically divergent species in the NW Pacific. While COI inter-lineage divergence in M. cephalus exceeds intra-lineage divergence (0.4 to 0.2%) by a ratio greater than 10, which is an indicator of cryptic species according to the 10× rule [58], the nuclear bi-parentally inherited markers provide the most definitive evidence of cryptic species. The assignment tests strongly supported the existence of three groups where 99% of M. cephalus were significantly assigned (posterior probabilities <10% and >90%). Furthermore, the complete congruence between nuclear groups and mitochondrial lineages suggests that lineage sorting of ancestral polymorphisms has been completed. Lastly, further to the argument that the M. cephalus mitochondrial lineages are true species, the gene flow estimated with microsatellite loci appeared to be more limited among sympatric individuals belonging to different mtDNA lineages than   See Table 1 for the complete name of locations. * < 0.05, ** < 0.01, *** < 0.001 in bold significant values after bonferroni correction [45].
between geographically remote individuals sharing the same mtDNA lineage. These results reveal the marked genetic isolation of the three M. cephalus species in the NW Pacific. If, as suggested earlier, the genetic differentiation of these cryptic species is related to past physical barriers to gene flow, the absence of genetic introgression suggests active reproductive isolation where they coexist sympatrically. Previous studies have reported various different reproductive behaviours in M. cephalus in the NW Pacific. While some individuals migrate over large geographic scales, others appear to be year-round residents of estuaries [59][60][61][62][63][64]. Migration was observed over the winter solstice (middle of December-early January), individuals of of 3-4 years of age (FL = 43-55 cm) migrate with the cold North China Coastal Current from the Eastern China Sea to south of the Taiwan Strait to spawn offshore of Kaohsiung, Taiwan [59][60][61]. The newly hatched larvae drift with the coastal currents and then recruit into the estuaries of the western coast of Taiwan at one to two months post-hatching [59,[62][63][64]. The presence of M. cephalus juveniles in Taiwanese estuaries earlier in the year (end of autumn) was interpreted, therefore, as evidence of the existence of a resident M. cephalus population [22].
In this study, all individuals sampled from the supposed spawning ground of migrant M. cephalus (offshore waters of Kaohsiung) in the spawning season (end of December), over three successive years, belonged to NWP1. The high GSI of these individuals suggested that they were spawning or would spawn soon. By contrast, M. cephalus collected in estuarine or coastal waters in the same period belonged to NWP1, NWP2 and NWP3, but all exhibited low GSI (Table 2). This suggests that genetic isolation among M. cephalus species is maintained by spatial and temporal differences in spawning migration patterns.

Population structure and contemporary gene flow
No nuclear genetic heterogeneity was recorded within the cryptic M. cephalus species among the locations sampled within the NW Pacific. This findings contrast with those of Liu et al. [23], who documented up to four regional groups of M. cephalus along the Chinese coast. However, because we analysed samples from the same locations (Hainan, Shantou and Qingdao, which is a station close to Lianyungang), it is evident that the population structure revealed using AFLP [23] relied mostly on the abundant difference of the three cryptic species. Whereas the Hainan, Shantou and Qingdao samples were characterised by high abundances of NWP3, NWP2 and NWP1, respectively, Liu et al. [23] described them as belonging to 3 populations: IV, II and III, respectively. If there was no experimental artefact, according to Liu et al. [23], a northern population of NWP2 is expected to be isolated in the Bohai Sea (population I), as this population is genetically close to population II.
Finally, if all cryptic M. cephalus species revealed in this study are able to disperse and migrate at sea, low genetic structuring is not surprising because large oceanographic current systems in the NW Pacific are expected to greatly facilitate larval dispersal over considerable distances. However, dynamic oceanographic systems can also profoundly restrict connectivity among groups, genetic heterogeneity is generally increased in marine zones with predictable oceanographic fronts [65][66][67].
The distributional range of the three NW Pacific M. cephalus species appears mainly to be shaped by three major oceanographic current systems, namely the South China, North China Coastal and Kuroshio Currents. The range of the NWP2 species appears to match the circulation pattern of the Kuroshio Current, which brings warm tropical water northeastward from the Luzon Straits, past both the west and east coasts of Taiwan and, finally, to Japan, where a branch of the current (the Tsushima Current) enters the Japanese Sea [68,69]. As the world's second-largest ocean current, the Kuroshio has a major effect on other current systems in the area. It, however, mainly remains offshore, whereas the other currents are located on the Taiwanese and North China shelves, and are more coastal. Thus, the South China Current, which flows northward from the South China Sea to the Taiwan Strait, brings warm water which mixes with cold water from the North China Coastal Current flowing southward from the Yellow Sea and East China Sea. The NWP3 distribution range seems to follow the warm water of the South China Current, while the NWP1 species appears to be restricted to the cold waters of the North China Coastal Current, which flows southward (with the northeastern monsoon) during the winter. Such close relationships between oceanographic current systems and fish species distributions presumably reflect temperature preferences, because differences between the Kuroshio Current and the China Coastal Current can range from 6 to 11°C ( [70], Figure 5a). While the Kuroshio Current remains warm throughout the year (23-26°C), the China Coastal Current exhibits a wide variation in temperature, from 12°C in March to 20°C in September. In this context, NWP1 may provide a particularly good indicator species for global warming. NWP1 individuals are assumed to spawn at the front (19-21°C) between the cold southward flowing North China Coastal Current and a warm branch of the northward flowing Kuroshio Current, with the front shifting according to the strength of the North China Coastal Current [20,59,60].

Conclusions
This study successfully identified three cryptic M. cephalus species in the NW Pacific Ocean, using both microsatellites and mitochondrial genetic markers. The genetic architecture and current distribution ranges of these species suggest a complex interaction of contemporary and evolutionary factors. Allopatric and adaptive speciation probably occurred during the Pliocene and Pleistocene epochs, when the Sea of Japan was periodically isolated and influenced by cold northern oceanographic currents. After the last glacial maxima, species range expansions were probably facilitated by changes in the oceanographic current structure and the species' temperature preferences. Currently, these species overlap in their distribution ranges, and it is likely that the genetic integrity of each is maintained by temporal and spatial isolation during the spawning period. This study emphasizes the importance of historical separations of marginal seas, in conjunction with fluctuating temperatures, in creating aquatic biodiversity in the NW Pacific. Isolation and demographic fluctuations are expected to constitute a powerful evolutionary force that increases speciation processes in the marine environment. These regional results also shed new light on the potential genetic diversity of M. cephalus worldwide.