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  <cit> .
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  <cit> .
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 .
the northwestern  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 semi-closure of the south china sea and the partial or full exposure of the east china and yellow seas  <cit> .
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  that supplies a large amount of heat to northern areas  <cit> .
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  <cit> .
this hypothesis has been confirmed by recent phylogeographic investigations in the nw pacific area .
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  <cit>  considered the various divergent lineages to be different species of pleistocene origin, while others  <cit>  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.
n: sample size, nl: number of lineage, % div.
: percentage of divergence.
the flathead mullet, mugil cephalus linnaeus,  <dig>  is a euryhaline teleost distributed worldwide in coastal waters, lagoons, bays and estuaries between latitudes 42° n and 42° s  <cit> .
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  <dig> , which has been suggested to be the result of overfishing and the impact of global warming  <cit> .
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  <cit> .
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  demonstrates that these lineages are independent  <cit> .
by contrast, an inter-lineage divergence of  <dig>  to  <dig> % 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  <cit> .
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  <cit> .
the first study used an allozyme locus, gpi-a , which may be able to distinguish migratory from resident m. cephalus populations  <cit> .
more recently, liu et al.
<cit>  found a high level of genetic structure in m. cephalus using aflp, and four populations were identified among a set of  <dig> samples collected along the chinese coast, with the southern samples  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  <cit> .
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  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.
methods
sampling and dna extraction
a total of  <dig> sub-adult to adult m. cephalus were collected from  <dig> locations in the nw pacific .
additional temporal sampling  was conducted at three localities around taiwan: the coast of keelung , which is located midway along the m. cephalus spawning migration route from the eastern china sea; the offshore waters of kaohsiung , where spawning occurs; and the kaoping estuary , a tropical system located in southwestern taiwan.
sampling locations, dates of collection, sample sizes for microsatellites analyses, mtdna lineage , mean fork length  and gonadosomatic index  data for the m. cephalus from taiwan are shown in table  <dig>  muscle tissues were preserved in 95% ethanol prior to dna extraction.
genomic dna was extracted from muscle tissue using a dna purification kit , preserved in te buffer and then quantified and diluted to  <dig> ng/μl for pcr.
rs: rapid pcr screening of different lineages.
the ratio of different lineages in each sample, their fork length and gonadosomatic index  are also shown.
mitochondrial dna analysis
partial mtdna coi gene sequences  were amplified by polymerase chain reaction  using the universal primers fishf <dig> and fishr <dig>  <cit> .
to identify m. cephalus cytochrome b lineages  <cit> , cytochrome b gene was also sequenced for  <dig> individuals belonging to three coi lineages using the cytochrome b primers mcglu- <dig> f: ggcttgaaaaaccaccgttg and mccytbr: agtactgtggcagagcttgg.
pcr was performed in a biometra tgradient thermocycler with a 15-μl reaction volume containing  <dig>  μm dntps,  <dig>  μl of 10× pcr buffer ,  <dig>  μm each of forward and reverse primer,  <dig>  u taq dna polymerase , and  <dig>  μl of template dna.
mtdna amplification was carried out using the following pcr conditions:  <dig> cycles of denaturation at 94°c for  <dig> s, annealing at 58°c for  <dig> s and extension at 72°c for  <dig> s after heating at 94°c for  <dig> min.
the pcr products were electrophoresed on a  <dig> % agarose gel  and stained with ethidium bromide for band characterisation via ultraviolet trans-illumination.
all sequencing reactions were performed according to the manufacturer's protocol  using the same forward primer as used for pcr.
all coi sequences were automatically aligned using mafft version  <dig>  <cit>  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  <dig>   <cit> .
phylogenetic trees were reconstructed using the bayesian, neighbour joining , maximum likelihood  and maximum parsimony  methods.
the bayesian method was applied using the program mrbayes  <dig>   <cit> , and nj, ml and mp were applied using the program paup  <dig>   <cit> .
the transversion model with invariable sites hky + i + g  was selected for construction of the phylogenetic tree using mugil curema as an outgroup  and the akaike information criterion  in modeltest  <dig>   <cit> .
nodes with bootstrap values ≥ 70% were considered well supported  <cit> .
the fixation index  for all pair-wise comparisons among different populations was calculated to investigate genetic diversity among m. cephalus populations, and a permutation test  was performed using arlequin version  <dig>   <cit> .
population genetic diversity was measured within each of the populations based on the number of distinct haplotypes, gene diversity  and mean nucleotide diversity , using arlequin version  <dig>   <cit> .
the divergence time for the m. cephalus coi gene was estimated using bayesian evolutionary analysis sampling trees  ver.
<dig> .
<dig>  <cit>  with  <dig> million steps in a monte carlo markov chain  simulation .
coi sequences of m. cephalus samples from disjunct locations in the western atlantic  and eastern pacific  were used for the calibration, and assuming that the divergence of these populations corresponds to the closure of the isthmus of panama,  <dig>  my ago  <cit> .
a strict molecular clock was assumed for this final run.
the effective sample sizes  of parameters sampled from the mcmc were >  <dig> .
the results were viewed using tracer  <dig>   <cit> .
historical demographic/spatial expansions of m. cephalus were explored using two different approaches.
tajima's d  <cit>  and fu's fs  <cit>  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  <cit> .
rapid screening of the different mtdna lineages
a multiplex coi haplotype-specific pcr  was designed following the recommendations of previous studies  <cit>  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  <dig> coi sequences of m. cephalus produced in this study.
nwp <dig> f  was the forward primer for both lineage  <dig> and lineage  <dig>  while nwp3f  was the forward primer for lineage  <dig>  nwp <dig> r  was the reverse primer for lineage  <dig> and lineage  <dig>  while nwp2r  was the reverse primer for lineage  <dig>  pcr was performed using a biometra tgradient thermocycler with a 15-μl reaction volume containing  <dig>  μm of each dntp,  <dig>  μl 10× pcr buffer ,  <dig>  μm forward and reverse primers,  <dig>  u taq dna polymerase , and  <dig>  μl of template dna.
mtdna amplification was carried out using the following pcr program:  <dig> cycles of denaturation at 94°c for  <dig> s, annealing at 55°c for  <dig> s and extension at 72°c for  <dig> s after heating at 94°c for  <dig> min.
the pcr products were electrophoresed in a  <dig> % agarose gel  and stained with ethidium bromide.
lineage  <dig> had a pcr product of  <dig> bp, while lineage  <dig> and lineage  <dig> had pcr products of  <dig> bp and  <dig> bp, respectively .
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  <cit> .
reverse primers for each locus were labelled with fluorescent dyes , and multiplex pcr was performed in a 15-μl reaction volume containing  <dig>  ng dna,  <dig>  pmole each of the three reverse primers labelled with different fluorescent dyes,  <dig>  pmole of each forward primer,  <dig> mm dntp,  <dig>  mm mgcl <dig> and  <dig>  u of taq polymerase .
amplification was conducted as follows: initial denaturation at 95°c for  <dig> min, followed by  <dig> cycles at 94°c, 54-58°c and 72°c for  <dig> sec each.
locus polymorphisms were screened using an abi prism  <dig> auto dna sequencer .
lengths of microsatellite alleles were determined using a tamra-labelled  <dig> bp standard .
the calculation of the number of alleles, estimated  and observed heterozygosities  for each locality and for all samples, the genetic differentiation index  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  <dig>   <cit>  and genetix software  <cit> .
the significance level for multiple simultaneous comparisons was adjusted using the sequential bonferroni technique  <cit> .
the allelic richness of the minimum population size by a rarefaction method and fis was estimated with fstat  <dig> .
<dig>  <cit> .
the significance of the differences in allelic richness was tested by the wilcoxon signed-rank test for paired observations .
the program genepop version  <dig>   <cit>  was used to test for linkage disequilibrium among the ten loci analysed in this study .
the presence of intraspecific genetic structure was tested using the model-based clustering method  <cit> , as implemented in structure ver.
<dig> .
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  <dig>  <dig> steps and a chain length of  <dig>  <dig> 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 = { <dig> ~ 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  <cit> .
individuals were regarded as correctly assigned to a population when their q-value  was at least 80% after subtracting the posterior probability assignment of another population  <cit> .
RESULTS
mitochondrial dna
a total of  <dig> nucleotides of the coi gene were obtained from  <dig> individuals.
in these sequences,  <dig> positions were variable , defining  <dig> haplotypes .
phylogenetic trees reconstructed using either coi  or cytochrome b  sequences  <cit>  all supported the existence of three highly supported monophyletic lineages in nw pacific m. cephalus.
seventeen coi haplotypes were found in lineage  <dig>  including  <dig> nucleotide transitions and three nucleotide transversions.
of these, ht <dig> was the most frequent haplotype, being shared by the samples from taiwan, japan, qingdao and russia, and it was connected to the other  <dig> haplotypes  by one to two steps in a "star-like" haplotype network .
the other  <dig> coi haplotypes were minor and mainly found in keelung and kaohsiung.
fifteen coi haplotypes were found in lineage  <dig>  including  <dig> nucleotide transitions and two nucleotide transversions.
two major haplotypes  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  <dig> consisted of only four coi haplotypes, with the main haplotype  occurring at the kaoping river and hainan sites and in a few samples from shantou and keelung.
the fixation index  showed significant differentiation among the three lineages .
the pair-wise Φst values were  <dig>   between lineage  <dig> and lineage 2;  <dig>   between lineage  <dig> and lineage  <dig> and  <dig>   between lineage  <dig> and lineage  <dig>
fifteen variable sites were found within the amplified 627-bp coi fragment  between lineage  <dig> and lineage  <dig>  with approximately  <dig> % of the inferred changes being transversions.
while  <dig> and  <dig> variable sites were found between lineages  <dig> and  <dig>  and between lineages  <dig> and  <dig> , 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  <dig> % mutations/site/my.
therefore, the calibrated divergence time using beast was  <dig>  my  between lineages  <dig> and  <dig> and  <dig>  my  between lineage  <dig> and the ancestor of lineages  <dig> and  <dig>  beast also suggested that lineage  <dig> 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  <dig> and  <dig>
the mismatch distributions of lineage  <dig> and  <dig> were unimodal, with almost no pair-wise differences between them .
lineage  <dig> exhibited a more uneven mismatch distribution compared to lineage  <dig>  by contrast, the mismatch distribution of lineage  <dig> was bimodal , 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 fs were negative and highly significant for lineage  <dig>  for lineage  <dig>  tajima's d values were all negative, but only significant for the control region, while fu's fs were all negative and significant for all genes.
finally, both tajima's d and fu's fs were negative, but not significant for the lineage  <dig> for any gene considered .
number of individuals , number of haplotype  or alleles , haplotype diversity  or observed heterozygosity , nucleotide diversity  for each group of samples.
tajima's d and fu's fs, the corresponding p value, and the raggedness index  of these lineages were also indicated.
<dig> data from  <cit> ,  <dig> data from  <cit> .
*p <  <dig> , ** p <  <dig> , ***p <  <dig>
microsatellites
a total of  <dig> alleles with an average of  <dig>  alleles per locus were observed across the ten microsatellite loci .
all loci were polymorphic, except for locus mce- <dig>  which was monomorphic in some of the samples.
twenty-one of  <dig> sample-locus combinations deviated significantly from hardy-weinberg proportions after bonferroni correction.
however, the deviations of these loci were not significant in most locations  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 .
amova indicated that the genetic differentiation of m. cephalus among  <dig> spatial and temporal samples was highly significant .
a locus-by-locus 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 .
kl:keelung; ks:kaohsiung; kp:kaoping river.
* <  <dig> , ** <  <dig> , *** <  <dig>  in bold significant values after bonferroni correction  <cit>
the results from the structure assignment test supported three clusters based both on the log probability of the data  and the statistic Δk, as described in  <cit> , respectively .
with k =  <dig>  almost all individuals were clearly assigned to one of the three clusters .
the posterior probabilities corresponding to the assignment of individuals to cluster  <dig>  were all unambiguous, while the assignments to one of the two other clusters  included several ambiguous individuals .
the individuals from kaohsiung, qingdao, the ariake sea and russia were mostly assigned to nwp1; the individuals from the philippines were all assigned to nwp <dig>  while at least two clusters coexist in keelung and the kaoping river of taiwan , hainan, the pearl river and shantou of china  and okinawa and yokosuka of japan  .
according to their occurrence areas, nwp <dig> is mainly observed in the north of our sampling area , whereas the nwp <dig> is mostly found at southern stations .
nwp <dig> has the widest range; it is present from the north to the south of the sampling area .
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 .
the results indicated that  <dig> of  <dig> individuals assigned to nwp <dig> harboured mtdna representing lineage  <dig>   <dig> of  <dig> individuals assigned to nwp <dig> harboured mtdna of lineage  <dig>  and all  <dig> individuals assigned to nwp <dig> harboured lineage  <dig> mtdna.
this resulted in a successful assignment rate of  <dig> %.
although there were seven ambiguous individuals in nwp <dig> and nwp <dig> below our assignment criteria, the mean q-values of assignment test for these two clusters were still over  <dig>  .
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 .
genetic variability analysis showed that only two of the  <dig> combinations, loci mce- <dig> and mce- <dig> of nwp <dig>  deviated significantly from hardy-weinberg proportions after bonferroni correction .
the allelic richness per locus based on the minimum sample size  was significantly lower at nwp <dig> than nwp <dig>  and nwp <dig>  and slightly lower at nwp <dig> than nwp <dig>  for the ten microsatellite loci .
the genotypes at mce- <dig> of nwp <dig> and mce- <dig> of nwp <dig> were all homozygous.
population genetic differentiation among these three clusters was highly significant .
pair-wise θ comparison indicated that nwp <dig> was the most differentiated , while level of differentiation was lower but still highly significant between nwp <dig> and nwp <dig> .
in addition, no genetic heterogeneity was detected among samples belonging to any one of the three otus, with the exception of the ariake sea, which had some genetic heterogeneity among the samples of nwp <dig> .
see table  <dig> for the complete name of locations.
* <  <dig> , ** <  <dig> , *** <  <dig>  in bold significant values after bonferroni correction  <cit>
discussion
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 , 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, nwp <dig> and nwp <dig>  are of pleistocene origin.
the nwp <dig> lineage diverged from the nwp <dig> and nwp <dig> common ancestor during the miocene or the pliocene epoch, and we assume that the origin of nwp <dig> does not relate to geological events in the nw pacific because this lineage is also observed in the southwestern  pacific .
while it appears likely that the refugium of the nwp <dig> lineage was located in the south china sea, the location of the nwp <dig> refugium is more difficult to assess.
the closure of the tsushima straits , which isolated the japan sea, may have been the powerful vicariant event that was at the origin of the nwp <dig> lineage.
reconstruction of the history of the tsushima current  <cit>  demonstrated that during the interval between  <dig> - <dig>  my, the current only flowed into the sea of japan periodically and the volume and salinity of water was lower than at present.
after  <dig>  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 nwp <dig> and nwp <dig>  suggesting that the nwp <dig> ancestor was not able to enter in the japan sea before  <dig>  my.
during each glacial episode, ancestors of the nwp <dig> lineage that were trapped in the japan sea probably experienced several demographic crashes due to decreases in surface water temperatures  <cit> .
this hypothesis is consistent with the observed genetic diversity and tajima's d and fu's fs indices, which indicate bottleneck events in nwp <dig> .
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 nwp <dig> presents a higher genetic diversity than nwp <dig>  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  <cit> .
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 nwp <dig> from a northern route of dispersion through the tsushima strait, the ryukyu islands and east japan were recolonised by both nwp <dig> and nwp <dig> from southern routes that followed the kuroshio current .
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  <cit> .
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  <cit> .
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  by a ratio greater than  <dig>  which is an indicator of cryptic species according to the 10× rule  <cit> , 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 .
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 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  <cit> .
migration was observed over the winter solstice , individuals of of 3- <dig> years of age  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  <cit> .
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  <cit> .
the presence of m. cephalus juveniles in taiwanese estuaries earlier in the year  was interpreted, therefore, as evidence of the existence of a resident m. cephalus population  <cit> .
in this study, all individuals sampled from the supposed spawning ground of migrant m. cephalus  in the spawning season , over three successive years, belonged to nwp <dig>  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 nwp <dig>  nwp <dig> and nwp <dig>  but all exhibited low gsi .
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.
<cit> , who documented up to four regional groups of m. cephalus along the chinese coast.
however, because we analysed samples from the same locations , it is evident that the population structure revealed using aflp  <cit>  relied mostly on the abundant difference of the three cryptic species.
whereas the hainan, shantou and qingdao samples were characterised by high abundances of nwp <dig>  nwp <dig> and nwp <dig>  respectively, liu et al.
<cit>  described them as belonging to  <dig> populations: iv, ii and iii, respectively.
if there was no experimental artefact, according to liu et al.
<cit> , a northern population of nwp <dig> is expected to be isolated in the bohai sea , 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  <cit> .
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 nwp <dig> 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  enters the japanese sea  <cit> .
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 nwp <dig> distribution range seems to follow the warm water of the south china current, while the nwp <dig> species appears to be restricted to the cold waters of the north china coastal current, which flows southward  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  <dig> to 11°c .
while the kuroshio current remains warm throughout the year , the china coastal current exhibits a wide variation in temperature, from 12°c in march to 20°c in september.
in this context, nwp <dig> may provide a particularly good indicator species for global warming.
nwp <dig> individuals are assumed to spawn at the front  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  <cit> .
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.
authors' contributions
kns carried out the molecular genetic studies, performed the statistical analysis and drafted the manuscript.
bwj and cch participated in the collection of samples and biological data.
wnt participated in the design of the study and obtained funding.
jdd conceived the study, participated in its design and coordination, and helped to draft the manuscript.
all authors have read and approved the final version of the manuscript.
supplementary material
additional file 1
table s <dig>  variable positions in the  <dig> bp mitochondrial coi gene segment of mugil cephalus from  <dig> locations in the northwestern pacific.
dots represent identical nucleotides relative to haplotype  <dig>  frequency of each haplotypes for each lineages  are also shown.
different color means the locations of the lineage specific nucleotides.
click here for file
additional file 2
table s <dig>  genetic variability at ten microsatellite loci of mugil cephalus collected in the northwestern pacific for spatial and temporal genetic structure test.
table-wide significance levels were applied using the sequential bonferroni technique  <cit> .
click here for file
additional file 3
table s <dig>  log probability and Δk  <cit>  for each number of clusters in the bayesian assignment test as implemented in structure  <cit> .
click here for file
additional file 4
table s <dig>  the mean q-values and standard deviations  for assignment test of  <dig> mugil cephalus cryptic species  as implemented in structure  <cit> .
click here for file
additional file 5
table s <dig>  genetic variability at ten microsatellite loci of mugil cephalus among  <dig> cryptic species.
table-wide significance levels were applied using the sequential bonferroni technique  <cit> .
click here for file
additional file 6
figure s <dig>  ethidium bromide stained 2% agarose gel showing the multiplex coi haplotype-specific pcr  for the rapid screening of three mugil cephalus cryptic species in the nw pacific.
1-13: unidentified individuals belonging either to lineage nwp <dig> , nwp <dig>  or nwp <dig> .
m: 100-bp dna ladder.
click here for file
additional file 7
figure s <dig>  phylogenetic relationships within mugil cephalus recovered from  <dig> bp of the cytochrome b sequences according to the neighbour-joining tree using kimura  <dig> parameter distance.
leaves of the tree correspond to haplotypes of m. cephalus observed by ke et al.
<cit>   and this study .
the values above the branches are bootstrap support .
bootstrap supports higher than 50% are displayed.
click here for file
acknowledgements
this study could not have been carried out without the generosity in collecting specimens of those listed below: j.x.
liu, x.p.
nie, x.c.
liu, h.p.
chen and x. qi from china; s. chow and h. imai from japan; m.p.
garcia from the philippines; y. kartavtsev from russia; k.t.
shao and c.w.
chang from taiwan.
this study was financially supported by the national science council of the republic of china .
the authors would like to thank z.c.
lin for microsatellite genotyping and dna sequencing, o. sadio for statistical advice and a.k.
whitfield, f. lecomte, n. leander and d. mckenzie for their comments on the manuscript.

 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. 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. 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. the northwestern 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 semi-closure of the south china sea and the partial or full exposure of the east china and yellow seas. 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 that supplies a large amount of heat to northern areas. 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. this hypothesis has been confirmed by recent phylogeographic investigations in the nw pacific area. 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 considered the various divergent lineages to be different species of pleistocene origin, while others 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. n: sample size, nl: number of lineage, % div.: percentage of divergence. the flathead mullet, mugil cephalus linnaeus, is a euryhaline teleost distributed worldwide in coastal waters, lagoons, bays and estuaries between latitudes 42° n and 42° s. 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, which has been suggested to be the result of overfishing and the impact of global warming. 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. 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 demonstrates that these lineages are independent. by contrast, an inter-lineage divergence of to % 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. 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. the first study used an allozyme locus, gpi-a, which may be able to distinguish migratory from resident m. cephalus populations. more recently, liu et al. found a high level of genetic structure in m. cephalus using aflp, and four populations were identified among a set of samples collected along the chinese coast, with the southern samples 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. 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 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. methods sampling and dna extraction a total of sub-adult to adult m. cephalus were collected from locations in the nw pacific. additional temporal sampling was conducted at three localities around taiwan: the coast of keelung, which is located midway along the m. cephalus spawning migration route from the eastern china sea; the offshore waters of kaohsiung, where spawning occurs; and the kaoping estuary, a tropical system located in southwestern taiwan. sampling locations, dates of collection, sample sizes for microsatellites analyses, mtdna lineage, mean fork length and gonadosomatic index data for the m. cephalus from taiwan are shown in table muscle tissues were preserved in 95% ethanol prior to dna extraction. genomic dna was extracted from muscle tissue using a dna purification kit, preserved in te buffer and then quantified and diluted to ng/μl for pcr. rs: rapid pcr screening of different lineages. the ratio of different lineages in each sample, their fork length and gonadosomatic index are also shown. mitochondrial dna analysis partial mtdna coi gene sequences were amplified by polymerase chain reaction using the universal primers fishf and fishr. to identify m. cephalus cytochrome b lineages, cytochrome b gene was also sequenced for individuals belonging to three coi lineages using the cytochrome b primers mcglu- f: ggcttgaaaaaccaccgttg and mccytbr: agtactgtggcagagcttgg. pcr was performed in a biometra tgradient thermocycler with a 15-μl reaction volume containing μm dntps, μl of 10× pcr buffer, μm each of forward and reverse primer, u taq dna polymerase, and μl of template dna. mtdna amplification was carried out using the following pcr conditions: cycles of denaturation at 94°c for s, annealing at 58°c for s and extension at 72°c for s after heating at 94°c for min. the pcr products were electrophoresed on a % agarose gel and stained with ethidium bromide for band characterisation via ultraviolet trans-illumination. all sequencing reactions were performed according to the manufacturer's protocol using the same forward primer as used for pcr. all coi sequences were automatically aligned using mafft version 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. phylogenetic trees were reconstructed using the bayesian, neighbour joining, maximum likelihood and maximum parsimony methods. the bayesian method was applied using the program mrbayes, and nj, ml and mp were applied using the program paup. the transversion model with invariable sites hky + i + g was selected for construction of the phylogenetic tree using mugil curema as an outgroup and the akaike information criterion in modeltest. nodes with bootstrap values ≥ 70% were considered well supported. the fixation index for all pair-wise comparisons among different populations was calculated to investigate genetic diversity among m. cephalus populations, and a permutation test was performed using arlequin version. population genetic diversity was measured within each of the populations based on the number of distinct haplotypes, gene diversity and mean nucleotide diversity, using arlequin version. the divergence time for the m. cephalus coi gene was estimated using bayesian evolutionary analysis sampling trees ver.. with million steps in a monte carlo markov chain simulation. coi sequences of m. cephalus samples from disjunct locations in the western atlantic and eastern pacific were used for the calibration, and assuming that the divergence of these populations corresponds to the closure of the isthmus of panama, my ago. a strict molecular clock was assumed for this final run. the effective sample sizes of parameters sampled from the mcmc were >. the results were viewed using tracer. historical demographic/spatial expansions of m. cephalus were explored using two different approaches. tajima's d and fu's fs 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. rapid screening of the different mtdna lineages a multiplex coi haplotype-specific pcr was designed following the recommendations of previous studies 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 coi sequences of m. cephalus produced in this study. nwp f was the forward primer for both lineage and lineage while nwp3f was the forward primer for lineage nwp r was the reverse primer for lineage and lineage while nwp2r was the reverse primer for lineage pcr was performed using a biometra tgradient thermocycler with a 15-μl reaction volume containing μm of each dntp, μl 10× pcr buffer, μm forward and reverse primers, u taq dna polymerase, and μl of template dna. mtdna amplification was carried out using the following pcr program: cycles of denaturation at 94°c for s, annealing at 55°c for s and extension at 72°c for s after heating at 94°c for min. the pcr products were electrophoresed in a % agarose gel and stained with ethidium bromide. lineage had a pcr product of bp, while lineage and lineage had pcr products of bp and bp, respectively. 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. reverse primers for each locus were labelled with fluorescent dyes, and multiplex pcr was performed in a 15-μl reaction volume containing ng dna, pmole each of the three reverse primers labelled with different fluorescent dyes, pmole of each forward primer, mm dntp, mm mgcl and u of taq polymerase. amplification was conducted as follows: initial denaturation at 95°c for min, followed by cycles at 94°c, 54-58°c and 72°c for sec each. locus polymorphisms were screened using an abi prism auto dna sequencer. lengths of microsatellite alleles were determined using a tamra-labelled bp standard. the calculation of the number of alleles, estimated and observed heterozygosities for each locality and for all samples, the genetic differentiation index 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 and genetix software. the significance level for multiple simultaneous comparisons was adjusted using the sequential bonferroni technique. the allelic richness of the minimum population size by a rarefaction method and fis was estimated with fstat.. the significance of the differences in allelic richness was tested by the wilcoxon signed-rank test for paired observations. the program genepop version was used to test for linkage disequilibrium among the ten loci analysed in this study. the presence of intraspecific genetic structure was tested using the model-based clustering method, as implemented in structure ver.. 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 steps and a chain length of 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 = { ~ 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. individuals were regarded as correctly assigned to a population when their q-value was at least 80% after subtracting the posterior probability assignment of another population.  mitochondrial dna a total of nucleotides of the coi gene were obtained from individuals. in these sequences, positions were variable, defining haplotypes. phylogenetic trees reconstructed using either coi or cytochrome b sequences all supported the existence of three highly supported monophyletic lineages in nw pacific m. cephalus. seventeen coi haplotypes were found in lineage including nucleotide transitions and three nucleotide transversions. of these, ht was the most frequent haplotype, being shared by the samples from taiwan, japan, qingdao and russia, and it was connected to the other haplotypes by one to two steps in a "star-like" haplotype network. the other coi haplotypes were minor and mainly found in keelung and kaohsiung. fifteen coi haplotypes were found in lineage including nucleotide transitions and two nucleotide transversions. two major haplotypes 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 consisted of only four coi haplotypes, with the main haplotype occurring at the kaoping river and hainan sites and in a few samples from shantou and keelung. the fixation index showed significant differentiation among the three lineages. the pair-wise Φst values were between lineage and lineage 2; between lineage and lineage and between lineage and lineage fifteen variable sites were found within the amplified 627-bp coi fragment between lineage and lineage with approximately % of the inferred changes being transversions. while and variable sites were found between lineages and and between lineages and, 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 % mutations/site/my. therefore, the calibrated divergence time using beast was my between lineages and and my between lineage and the ancestor of lineages and beast also suggested that lineage 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 and the mismatch distributions of lineage and were unimodal, with almost no pair-wise differences between them. lineage exhibited a more uneven mismatch distribution compared to lineage by contrast, the mismatch distribution of lineage was bimodal, 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 fs were negative and highly significant for lineage for lineage tajima's d values were all negative, but only significant for the control region, while fu's fs were all negative and significant for all genes. finally, both tajima's d and fu's fs were negative, but not significant for the lineage for any gene considered. number of individuals, number of haplotype or alleles, haplotype diversity or observed heterozygosity, nucleotide diversity for each group of samples. tajima's d and fu's fs, the corresponding p value, and the raggedness index of these lineages were also indicated. data from, data from. *p <, ** p <, ***p < microsatellites a total of alleles with an average of alleles per locus were observed across the ten microsatellite loci. all loci were polymorphic, except for locus mce- which was monomorphic in some of the samples. twenty-one of sample-locus combinations deviated significantly from hardy-weinberg proportions after bonferroni correction. however, the deviations of these loci were not significant in most locations 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. amova indicated that the genetic differentiation of m. cephalus among spatial and temporal samples was highly significant. a locus-by-locus 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. kl:keelung; ks:kaohsiung; kp:kaoping river. * <, ** <, *** < in bold significant values after bonferroni correction the results from the structure assignment test supported three clusters based both on the log probability of the data and the statistic Δk, as described in, respectively. with k = almost all individuals were clearly assigned to one of the three clusters. the posterior probabilities corresponding to the assignment of individuals to cluster were all unambiguous, while the assignments to one of the two other clusters included several ambiguous individuals. the individuals from kaohsiung, qingdao, the ariake sea and russia were mostly assigned to nwp1; the individuals from the philippines were all assigned to nwp while at least two clusters coexist in keelung and the kaoping river of taiwan, hainan, the pearl river and shantou of china and okinawa and yokosuka of japan. according to their occurrence areas, nwp is mainly observed in the north of our sampling area, whereas the nwp is mostly found at southern stations. nwp has the widest range; it is present from the north to the south of the sampling area. 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. the results indicated that of individuals assigned to nwp harboured mtdna representing lineage of individuals assigned to nwp harboured mtdna of lineage and all individuals assigned to nwp harboured lineage mtdna. this resulted in a successful assignment rate of %. although there were seven ambiguous individuals in nwp and nwp below our assignment criteria, the mean q-values of assignment test for these two clusters were still over. 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. genetic variability analysis showed that only two of the combinations, loci mce- and mce- of nwp deviated significantly from hardy-weinberg proportions after bonferroni correction. the allelic richness per locus based on the minimum sample size was significantly lower at nwp than nwp and nwp and slightly lower at nwp than nwp for the ten microsatellite loci. the genotypes at mce- of nwp and mce- of nwp were all homozygous. population genetic differentiation among these three clusters was highly significant. pair-wise θ comparison indicated that nwp was the most differentiated, while level of differentiation was lower but still highly significant between nwp and nwp. in addition, no genetic heterogeneity was detected among samples belonging to any one of the three otus, with the exception of the ariake sea, which had some genetic heterogeneity among the samples of nwp. see table for the complete name of locations. * <, ** <, *** < in bold significant values after bonferroni correction discussion 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, 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, nwp and nwp are of pleistocene origin. the nwp lineage diverged from the nwp and nwp common ancestor during the miocene or the pliocene epoch, and we assume that the origin of nwp does not relate to geological events in the nw pacific because this lineage is also observed in the southwestern pacific. while it appears likely that the refugium of the nwp lineage was located in the south china sea, the location of the nwp refugium is more difficult to assess. the closure of the tsushima straits, which isolated the japan sea, may have been the powerful vicariant event that was at the origin of the nwp lineage. reconstruction of the history of the tsushima current demonstrated that during the interval between - my, the current only flowed into the sea of japan periodically and the volume and salinity of water was lower than at present. after 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 nwp and nwp suggesting that the nwp ancestor was not able to enter in the japan sea before my. during each glacial episode, ancestors of the nwp lineage that were trapped in the japan sea probably experienced several demographic crashes due to decreases in surface water temperatures. this hypothesis is consistent with the observed genetic diversity and tajima's d and fu's fs indices, which indicate bottleneck events in nwp. 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 nwp presents a higher genetic diversity than nwp 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. 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 nwp from a northern route of dispersion through the tsushima strait, the ryukyu islands and east japan were recolonised by both nwp and nwp from southern routes that followed the kuroshio current. 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. 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. 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 by a ratio greater than which is an indicator of cryptic species according to the 10× rule, 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. 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 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. migration was observed over the winter solstice, individuals of of 3- years of age 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. 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. the presence of m. cephalus juveniles in taiwanese estuaries earlier in the year was interpreted, therefore, as evidence of the existence of a resident m. cephalus population. in this study, all individuals sampled from the supposed spawning ground of migrant m. cephalus in the spawning season, over three successive years, belonged to nwp 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 nwp nwp and nwp but all exhibited low gsi. 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., who documented up to four regional groups of m. cephalus along the chinese coast. however, because we analysed samples from the same locations, it is evident that the population structure revealed using aflp relied mostly on the abundant difference of the three cryptic species. whereas the hainan, shantou and qingdao samples were characterised by high abundances of nwp nwp and nwp respectively, liu et al. described them as belonging to populations: iv, ii and iii, respectively. if there was no experimental artefact, according to liu et al., a northern population of nwp is expected to be isolated in the bohai sea, 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. 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 nwp 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 enters the japanese sea. 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 nwp distribution range seems to follow the warm water of the south china current, while the nwp species appears to be restricted to the cold waters of the north china coastal current, which flows southward 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 to 11°c. while the kuroshio current remains warm throughout the year, the china coastal current exhibits a wide variation in temperature, from 12°c in march to 20°c in september. in this context, nwp may provide a particularly good indicator species for global warming. nwp individuals are assumed to spawn at the front 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.  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. authors' contributions kns carried out the molecular genetic studies, performed the statistical analysis and drafted the manuscript. bwj and cch participated in the collection of samples and biological data. wnt participated in the design of the study and obtained funding. jdd conceived the study, participated in its design and coordination, and helped to draft the manuscript. all authors have read and approved the final version of the manuscript. supplementary material additional file 1 table s variable positions in the bp mitochondrial coi gene segment of mugil cephalus from locations in the northwestern pacific. dots represent identical nucleotides relative to haplotype frequency of each haplotypes for each lineages are also shown. different color means the locations of the lineage specific nucleotides. click here for file additional file 2 table s genetic variability at ten microsatellite loci of mugil cephalus collected in the northwestern pacific for spatial and temporal genetic structure test. table-wide significance levels were applied using the sequential bonferroni technique. click here for file additional file 3 table s log probability and Δk for each number of clusters in the bayesian assignment test as implemented in structure. click here for file additional file 4 table s the mean q-values and standard deviations for assignment test of mugil cephalus cryptic species as implemented in structure. click here for file additional file 5 table s genetic variability at ten microsatellite loci of mugil cephalus among cryptic species. table-wide significance levels were applied using the sequential bonferroni technique. click here for file additional file 6 figure s ethidium bromide stained 2% agarose gel showing the multiplex coi haplotype-specific pcr for the rapid screening of three mugil cephalus cryptic species in the nw pacific. 1-13: unidentified individuals belonging either to lineage nwp, nwp or nwp. m: 100-bp dna ladder. click here for file additional file 7 figure s phylogenetic relationships within mugil cephalus recovered from bp of the cytochrome b sequences according to the neighbour-joining tree using kimura parameter distance. leaves of the tree correspond to haplotypes of m. cephalus observed by ke et al. and this study. the values above the branches are bootstrap support. bootstrap supports higher than 50% are displayed. click here for file acknowledgements this study could not have been carried out without the generosity in collecting specimens of those listed below: j.x. liu, x.p. nie, x.c. liu, h.p. chen and x. qi from china; s. chow and h. imai from japan; m.p. garcia from the philippines; y. kartavtsev from russia; k.t. shao and c.w. chang from taiwan. this study was financially supported by the national science council of the republic of china. the authors would like to thank z.c. lin for microsatellite genotyping and dna sequencing, o. sadio for statistical advice and a.k. whitfield, f. lecomte, n. leander and d. mckenzie for their comments on the manuscript.
in the northwestern 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. 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. 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. the divergence of nwp from nwp dates back to the pleistocene epoch and probably corresponds to the separation of the japan and china seas when sea levels dropped. three cryptic species of m. cephalus were identified within this region using an assignment test on the microsatellite data. this study successfully identified three cryptic species in m. cephalus inhabiting the nw pacific, using a combination of microsatellites and mitochondrial genetic markers. several divergent allopatric mtdna lineages have been described in this species worldwide, and three occur in sympatry in the nw pacific. 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. 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. each species corresponds with one of the three mtdna lineages in the coi phylogenetic tree. nwp is the most divergent species, with a distribution range that suggests tropical affinities, while nwp with a northward distribution from taiwan to russia, is a temperate species. although phylogeographic studies using mitochondrial dna information often reveal significant genetic heterogeneity within marine species, the evolutionary significance of such diversity is difficult to interpret with these markers. among these marine species, the cosmopolitan estuarine-dependent fish mugil cephalus represents an interesting case study.
 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 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 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.  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 using an assignment test on the microsatellite data. each species corresponds with one of the three mtdna lineages in the coi phylogenetic tree. nwp is the most divergent species, with a distribution range that suggests tropical affinities, while nwp with a northward distribution from taiwan to russia, is a temperate species. nwp is distributed along the warm kuroshio current. the divergence of nwp from nwp 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.  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. 