Durand et al.’s [13] mitochondrial phylogeny of the Mugilidae has uncovered a number of deeply divergent lineages within nominal species. Several of the lineages were paraphyletic with other species; other lineages represented reproductively isolated sympatric species, as demonstrated by genotypic frequencies at nuclear loci or inferred from karyotypes. Last, in some instances, deeply divergent sister-lineages characterized geographically separate populations within a species. Thus, there was substantial evidence for cryptic species in Mugilidae, based on the tree topology, on independent data from nuclear markers, and on the geographic distribution of sister lineages. We used W.N. Eschmeyer's fish database [17] as the reference for the current nomenclature. The current nomenclature was maintained for a lineage when its geographic distribution was compatible with the type locality of the species. By analogy with these cases where specific status was documented, we considered other deep lineages in Mugilidae, i.e. lineages whose distance to its nearest neighbour exceeded the gap between infra-specific and inter-specific pairwise distances (see section 2.5), to potentially represent additional cryptic species. We maintained the current nomenclature to designate those lineages that unambiguously correspond to the type material, based on the type locality, and we arbitrarily assigned capital letters to the other lineages. The other lineages were thus provisionally denominated “sp. A”, “sp. B”, etc.

We emphasize that our approach is not one of DNA barcoding, but one of molecular taxonomy, where molecular diagnoses of species and putative cryptic species are provided. We use gaps in the distributions of pairwise genetic distances as a means to distinguish deep lineages, which is where one may find analogy with barcoding. Nevertheless, the utility of CO1 barcoding for identifying species in the family Mugilidae will be the topic of a separate paper.

In this paper, genus nomenclature accords with our recent revision [18], where the following changes have been made, relative to the previous nomenclature: Moolgarda seheli and Valamugil buchanani have been placed together with Crenimugil crenilabis under Crenimugil, and Moolgarda cunnesius, Moolgarda engeli, Moolgarda perusii, and Valamugil robustus have been placed under the resurrected genus Osteomugil; likewise, Liza aurata, Liza bandialensis, Liza dumerili, Liza ramada, Liza richardsonii, Liza saliens, and Liza tricuspidens have been placed together with Chelon labrosus under Chelon; likewise, Chelon macrolepis, Chelon melinopterus, Chelon subviridis, Liza abu, Liza affinis, Liza alata, and Liza haematocheila have been placed under the resurrected genus Planiliza; C. planiceps has since then been synonymized with Liza tade [17] and placed under Planiliza; also, Sicamugil cascasia, Agonostomus monticola, Liza argentea, Rhinomugil nasutus, and Oedalechilus labiosus have been placed, respectively, under the resurrected genera Minimugil, Dajaus, Gracilimugil, Squalomugil, and Plicomugil whereas Xenomugil thoburni has been placed under Mugil; the genus names Liza, Moolgarda, Valamugil and Xenomugil have been dismissed; three new genera have been erected: Neochelon (for Liza falcipinnis), Parachelon (for Liza grandisquamis), and Pseudomyxus (for Myxus capensis).

Durand et al. [18] have also synonymized the genus Paramugil [19] with Planiliza. We must acknowledge that this was an error as explained in the following. We erroneously used as reference specimen for Planiliza parmatus individual MNHN-IC-2011-0212, numbered 118 in [18], which had been collected in south Java by S. Kleinertz. On the basis of photographs that he kindly agreed to examine, H. Senou identified this specimen as a Planiliza (“Chelon”), and not a Paramugil. This specimen was subsequently examined by J. Ghasemzadeh who also rejected our identification as Paramugil and identified it as Planiliza (“Liza”) melinoptera based on its external morphological features.

Durand et al.’s comprehensive mitochondrial phylogeny of the Mugilidae [13,18], which is based on the concatenated partial 16S rRNA, COI and cytochrome b gene sequences (3885 bp long in total) of 257 reference specimens (including 120 vouchers deposited in museum collections), was used for the present investigation. Zooms on regions of interest in this phylogeny are presented in Figs. 1 and 2.

Pairwise nucleotide distances between haplotypes sampled within each of seven mugilid genera (Chelon, Crenimugil, Dajaus, Ellochelon, Mugil, Osteomugil, Planiliza) were estimated under Mega5 [20] from the concatenated haplotype sequences at loci 16S rRNA, COI and cytb, which have been published previously [13,18]. Nucleotide distance was estimated according to the model of molecular evolution that, among the list of models proposed by Mega5, ranked as the most likely after the GTR-related model used to construct the phylogeny of [13], because the GTR model is not proposed by Mega5 for estimating nucleotide distances. The model thus chosen was the Tamura-Nei (TN93; [21]) with gamma distribution and invariable sites (+G + I) model. Nucleotide distances between lineages estimated according to the Kimura-2 parameter (K2P; [22]) model of molecular evolution were also presented. For each of the seven genera or species complexes we focused on, the resulting phylogenetic tree was examined together with the matrix of pairwise nucleotide distances between haplotypes. Our objective was first to determine the threshold below which distances all were infra-specific and above which they were all inter-specific. We then used this value as a yardstick to determine deep lineages that may represent cryptic species.

Further, alternative analysis of the dataset was done using the automatic gap determination algorithm proposed by N. Puillandre and co-authors to detect putative species from barcode datasets (ABGD; http://wwwabi.snv.jussieu.fr/public/abgd/; [23]). The analysis was run on each of the seven sequence datasets representing genera or species complexes, using the default settings of the program. This algorithm detects the gap in the distribution of pairwise nucleotide distances as the first significant gap beyond infra-specific distances and uses it to partition the dataset. Inference of the limit and gap detection are then recursively applied to previously obtained groups to get finer partitions until there is no further partitioning [23].

Pairwise distributions of nucleotide distances among individuals within each of the seven genera focused on in the present paper are presented Fig. 3. Detailed examination of the distribution in the genus Mugil (Fig. 3E) revealed a gap after 1%: it is therefore sensible to consider the values ≤1% separately from the rest of the distribution and to ascribe them to genetic variation at the infra-specific level. This 1% threshold value also precisely coincided with the right boundary of the first mode of pairwise nucleotide distances in the genus Chelon (Fig. 3A), and it encompassed the homologous first modes in Dajaus (Fig. 3C) and Ellochelon (Fig. 3D). Detailed examination of intraspecific distances within C. crenilabis and its morphologically distinct sister-species Crenimugil sp. B, the two most closely related lineages in the genus Crenimugil (Fig. 1D), showed no infra-specific distance greater than 1.5%. Within the other Crenimugil lineages (Fig. 1D) the highest pairwise distance was 2.2% while the lowest inter-lineage distance was 3.4%. Similarly, in the genus Osteomugil a gap in pairwise distances occurred between 2.1% and 4.3%. In Planiliza, a similar, although narrower gap was observed between 2.1%, the highest distance found within P. subviridis, and 2.6%. Thus, placing a threshold at 1% allowed the delineation of the first mode of the distribution of pairwise nucleotide distances in four (Chelon, Dajaus, Ellochelon, and Mugil) (Fig. 3A, C–E) of the seven mugilid genera tested. In the three remaining genera (Crenimugil, Osteomugil, and Planiliza) (Fig. 3B, F, G), the threshold should be placed at 2.5% based on the gap in the distribution of pairwise nucleotide distances.

In this section, we review all cases where lineages separated from the closest neighbouring lineage by a nucleotide distance larger than the threshold defined in the preceding section correspond to distinct species. In the mitochondrial phylogeny of the Mugilidae published by [13], Dajaus monticola consisted of three deeply rooted lineages (Fig. 1B), the two most recently diverged of which were geminate lineages distributed on either side of Central America, separated by 7.9% net nucleotide divergence under the K2P model and 7.0% under the TN93 + G + I model, all three markers combined. These two lineages have likely been geographically isolated from one the other for millions of years, hence are likely to represent separate species. The third lineage, from the eastern Pacific, branches externally to the two latter and is likely to represent another species. The type-locality of D. monticola is Jamaica [17]. Therefore, we here maintain the epithet monticola for the Atlantic lineage and provisionally designate the two other lineages, both from the eastern Pacific, as Dajaus sp. A and sp. B. The geographic distribution of the different lineages within D. monticola is presented in Fig. 4B. A subsequent study [15] estimated the divergence between the two geminate lineages D. monticola and Dajaus sp. A to be 31.8–11.8 million years old; it also evoked morphological differences between the two lineages from the eastern Pacific, confirming their status as separate species. The same study reported a fourth lineage currently under D. monticola from the Mexican rivers of the Gulf of Mexico, that is, geographically separated from what we here consider to be the true D. monticola [15].

The mitochondrial phylogeny of the M. cephalus/M. liza species complex (Fig. 2A) revealed 15 separate lineages, each with deep (>1%) rooting and shallow (<1%) within-lineage diversity. Three of these lineages, which occur in sympatry in Taiwan, belong to genetically distinct forms reproductively isolated from one another as demonstrated by their distinct composition at nuclear loci and by the quasi-absence of hybrids [14]. These lineages, coined NWP1-3 by [14], show 3.2%–4.8% nucleotide divergence under the K2P model, all three markers combined, and 3.3%–5.1% divergence under the TN93 + G + I model [20]. The three lineages are here assigned arbitrary species names sp. C, sp. I and sp. L., respectively (Fig. 2A). Mugil curema similarly consisted of a complex of species, where deeply rooted lineages were paraphyletic with another species, Mugil thoburni [18] (Fig. 2B). The type locality of M. curema is Bahia, Brazil [24] where only one lineage, also characterized by a chromosome complement number of 2n = 28, is present [13]. This is the “Type 2” (T2) karyological form of [25]. The topology of the tree (Fig. 2B) shows lineage T2 to root externally to the sub-clade consisting of M. thoburni and the other M. curema lineages. We here keep the name M. curema exclusively for lineage T2 and we designate the other lineages as Mugil sp. M to sp. O. The case of the three Mugil sp. M to sp. O lineages will be discussed in the following section.

A number of deeply divergent lineages potentially represent additional cryptic species in the Mugilidae. These cases are examined genus by genus in the following, where each lineage either was assigned a capital letter, or conserved its current name.

In the genus Chelon (Fig. 1A), Durand et al. [13] had sampled an unidentified Chelon sp. lineage from southeastern Africa (their specimen no. 161), which is here provisionally designated as Chelon sp. A. In the same genus, the haplogroup corresponding to Chelon dumerili actually comprised two distinct lineages separated by a net nucleotide distance, all three markers combined, of 7.5% (under the K2P model) or 8.1% (under the TN93 + G + I model). One lineage was exclusively sampled in western Africa, including Saint-Louis, Senegal [13], which is the type locality of the species [26]. We maintain epithet dumerili for this lineage. The other lineage was sampled exclusively in southeastern Africa and is here provisionally referred to as Chelon sp. B. We consider Chelon sp. B to be putatively a species distinct from C. dumerili based on the disjunction in geographic distributions and the level of nucleotide distance between the two lineages.

In the genus Crenimugil (Fig. 1D), three distinct lineages were observed within C. seheli under its current definition. These three lineages, which occur sympatrically in the Indo-West Pacific, are separated by a net nucleotide divergence, all three markers combined, of 4.5%–7.8% under the K2P model and 4.8%–8.6% under the TN93 + G + I model, whereas intraspecific nucleotide diversity under both models was ≤2.2%. The three lineages were paraphyletic with C. crenilabis and with an undescribed Crenimugil species sampled from Taiwan and Fiji, represented by individuals Nos. 238, 239 and 241 of [13]. Therefore, we consider them to characterize putative, distinct species, here designated as Crenimugil spp. A-C. The undescribed Crenimugil sp. species from Taiwan and Fiji is here designated as Crenimugil sp. D. Fig. 4A presents the geographic distribution of all four deep lineages within the C. seheli/C. crenilabis species complex.

In the genus Ellochelon (Fig. 1C), two separate lineages were observed, which diverged by 4.8% net nucleotide distance under the K2P model and 5.1% under the TN93 + G + I model, all three markers combined. One lineage included specimens from Waigeo, the type-locality [27] and French Polynesia, and another lineage was represented by a specimen from an unknown location in Australia. Epithet vaigiensis is here provisionally retained for the Ellochelon lineage sampled in Waigeo while the other lineage is provisionally assigned putative species name Ellochelon sp. A.

In the genus Mugil, the 13 distinct lineages originally uncovered within M. cephalus belonged to the same sub-clade as M. liza [13]. The average ± SD net nucleotide distance between lineages, from which M. liza was excluded was, all three markers combined, 3.6% ± 1% under the K2P model and 3.8% ± 1% under the TN93 + G + I model. Subsequently, a fourteenth lineage comprising haplotypes sampled from the Galapagos Islands was reported [18] (Fig. 2A). Three of the lineages currently within M. cephalus (i.e., sp. C, sp. I and sp. L), which occur in sympatry, belong to genetically distinct forms reproductively isolated from one another (see preceding section). Hence, basing our analogy on similarity in ratios of inter- to intra-lineage nucleotide distance, and also taking into account the current taxonomic standards that designate M. liza as a species separate from M. cephalus, we consider all other 11 lineages within M. cephalus under its current definition to be putative, distinct species. One notes that these lineages apparently have allopatric or parapatric distributions (Fig. 4C). The original description of M. cephalus [28] geographically refers to a species which “habitat Oceano Europeo”. Accordingly, we can designate without ambiguity the only lineage sampled in the Mediterranean Sea [13,29] as characterizing the actual M. cephalus. The 10 remaining lineages are here provisionally assigned putative species names Mugil sp. A, sp. B, spp. D-H, sp. J, sp. K and sp. Q. Fig. 4C presents the geographic distribution of the 15 deeply rooted lineages within the M. cephalus species complex (i.e., M. cephalus, M. liza, cryptic species Mugil spp. C, I, L, and putative cryptic species Mugil spp. A, B, D-H, J, K, Q). Three other Mugil spp. lineages within the species initially designated as M. curema were uncovered (see preceding section). Mugil sp. N and Mugil sp. O differ from M. curema by their karyotypes (respectively, 2n = 24 and 48) indicating that they are likely reproductively isolated from the latter as well as from one the other [13]. Mugil sp. M, from the Pacific Ocean, is paraphyletic with Mugil sp. N and Mugil sp. O, both from the Atlantic Ocean (Fig. 2B). These lineages differ by 3.2%–5.4% net nucleotide distance under the K2P model and 3.4–5.8% under the T93 + G + I model, all three markers combined. In comparison, nucleotide diversities within a lineage ranged from 0.2% to 0.5% (under both models). Hence, we consider the three lineages to represent distinct species. Still in Mugil, a lineage originally assigned to Mugil hospes by [13] was represented exclusively by haplotypes sampled from the Gulf of Mexico. Because the type-locality of M. hospes is Mazatlan in the eastern Pacific [30], it is sensible to provisionally assign the lineage from the Gulf of Mexico to a yet undetermined Mugil species, here designated as Mugil sp. R. Last, Mugil rubrioculus comprises two distinct lineages, one sampled from Venezuela, the type locality of the species [31], and the other one from the eastern Pacific [13]. The lineage from Venezuela retains epithet rubrioculus while the eastern Pacific lineage is here provisionally designated Mugil sp. P. This lineage may represent the same species as M. aff. rubrioculus previously mentioned from the eastern Pacific [31]. Fig. 4D presents the geographic distribution of the 7 deep lineages (i.e., M. curema, Mugil incilis, M. thoburni, and putative cryptic species Mugil spp. M-P) within the M. curema species complex.

In the genus Osteomugil (Fig. 1D), three distinct lineages representing Osteomugil cunnesius under its current definition have been found to be paraphyletic with Osteomugil perusii [13]. The type locality of O. cunnesius is the Moluccas [24]. Hence, the lineage sampled in Taiwan and in Vietnam by [13], geographically closest to the Moluccas, is provisionally retained as the actual O. cunnesius while the two other lineages, one from eastern Western Australia, the other one from South Africa, are here assigned provisional names Osteomugil sp. A and Osteomugil sp. B, respectively.

In the genus Planiliza (Fig. 1A), both Planiliza melinoptera and Planiliza tade were found to be polyphyletic ([13]; present work). Individual 118 of [13], which was initially, erroneously identified as Planiliza parmata is now recognized as a cryptic lineage of Planiliza melinoptera (see section 2.2). This lineage was separated from the P. melinoptera sampled in Fiji (individuals 109 and 111 of [13]) by 10.3% nucleotide distance (under the K2P model) and 11.4% (under the TN93 + G + I model), all three markers combined. We kept the Fiji specimens under P. melinoptera, because of the geographic proximity of Fiji with Vanikoro, the type locality [24] and we assigned provisional species name Planiliza sp. B to the lineage sampled in South Java. One of the two P. tade lineages concerned specimens sampled in Myanmar; the other lineage was sampled in northern Australia. These two regions being remote from the Red Sea, the type-locality of P. tade [32], we here followed a cautious line by designating as Planiliza sp. F the lineage from northern Australia, and Planiliza sp. I the lineage from Myanmar. Two sister-lineages were observed within P. macrolepis, separated by 3.5% (under the K2P model) or 3.7% (under the T93 + G + I model) net nucleotide divergence, all three markers combined: one lineage was exclusive to the western Indian Ocean west of the Seychelles Islands, including the Seychelles Islands and including South Africa, the type-locality [33] while the other lineage had a wide geographic distribution, as it consisted of all haplotypes sampled east of the Seychelles Islands, from the Maldives Islands to Fiji [13]. The latter is here provisionally designated as Planiliza sp. H, while its sister lineage retains the epithet macrolepis. Fig. 4E presents the geographic distribution of the two deep lineages within P. macrolepis (i.e., the actual P. macrolepis, and Planiliza sp. H). A lineage represented by a single individual from Taiwan (no. 108, “Planiliza sp.” of [18]) separated from its sister-lineage, P. melinoptera, by 7.2% nucleotide distance (under the K2P model) and 7.8% (under the TN93 + G + I model), all three markers combined, is here provisionally designated as Planiliza sp. G. Five other undetermined Planiliza species were here assigned provisional species names Planiliza sp. A (including individuals nos. 113 and 114 of [13]), Planiliza sp. C (nos. 103, 107), Planiliza sp. D (nos. 062, 062b), Planiliza sp. E (no. 057b), and Planiliza sp. J (no. 063). Our distinguishing Planiliza sp. J from its sister-lineage Planiliza sp. D is justified by the distance between the two lineages (2.6–2.8%, above the 2.5% threshold set for the genus).

Through automatic gap determination using the ABGD algorithm of [23], the present sequence dataset was found to conceal 10 separate lineages in genus Chelon, 6 in Crenimugil, 3 in Dajaus, 15 in the M. cephalus species complex, 6 in the M. curema species complex, 6 in Osteomugil, and 17 in Planiliza. The lineages designated by ABGD were all identical to those reported on Figs. 1 and 2, except for one lineage in Planiliza (Planiliza sp. J) that escaped detection using the default settings of ABGD.

The present results lead us to propose molecular diagnoses for a number of mugilid species currently considered as valid [17]. We aligned the partial sequences of the specimens characterized at the three loci (i.e., 16S, CO1 and cytb) used as phylogenetic markers [13,18] in each of 7 cases treated in the present study (i.e., Chelon spp., Crenimugil spp., Dajaus spp., species in the Mugil cephalus species complex, species in the M. curema species complex, Osteomugil spp. and Planiliza spp.). Ellochelon spp. was excluded because it consisted of two main lineages only, one of which was represented by a single individual in our dataset, thus insufficient for a comparison of inter-lineage vs intra-lineage variation. Variable nucleotide sites in each alignment were highlighted (Supplementary material, Tables S1–S21). Nucleotide sites diagnostic of species were determined. This information is summarized in Table 1. For example, M. cephalus L. is here diagnosed relative to the other species in the M. cephalus species complex by triplets (T₂₃₄, T₄₃₅, C₆₉₃) at locus CO1 and (T₁₈₃, T₄₈₃, G₅₁₀) at locus cytb, where nucleotide sites are numbered from the start of the gene, using the mitochondrial DNA sequence of Mugil sp. C (GenBank no. AP002930) as reference. Anonymous lineages designated by alphabetical letters were similarly diagnosed (Table 1). No molecular diagnosis was proposed for those species for which a single specimen was available: Chelon saliens, C. tricuspidens, Osteomugil cunnesius, and Planiliza abu. Similarly, no diagnosis was proposed for lineages Chelon sp. A, Dajaus sp. B, Mugil spp. C, K, Osteomugil sp. A, and Planiliza spp. E-G.

We are aware that future additional samples may lead to restricting the number of diagnostic sites for any given species relative to the other species in a genus. This is most likely to occur if additional cryptic species are sampled. However, the information in Table 1 may still provide the basis to future identification keys.