The unequivocal identification and classification of living organisms to the species level frequently relies on genetic evidence. Specific DNA sequences act as unrepeatable signatures and therefore constitute a unique DNA barcode for each species. Hebert et al. (2003) proposed that a single gene sequence would be sufficient to differentiate all, or at least the vast majority of animal species, and proposed the use of the mitochondrial DNA gene cytochrome oxidase subunit I (COI) as a global bioidentification system for animals. Initiatives, such as The Barcode of Life Database (http://www.barcodinglife.org) including The Fish Barcode of Life (http://www.fishbol.org), use a DNA-based identification system based on a relatively small fragment of COI.

However, the approach is not without controversy (Lipscomb et al., 2003; Moritz and Cicero, 2004). For a barcoding approach to species identification to succeed, within-species DNA sequences need to be more similar to one another than to sequences in different species. Recent studies show that this is generally the case, but there are exceptions. Hybridization among species would create taxonomic uncertainty: mitochondrial DNA is maternally inherited and any hybrid or subsequent generation would have the maternal species DNA only. Thus, the idea of a multi-locus DNA barcoding approach is progressively emerging, and limitations of mtDNA underline the requirement of nuclear regions (Frezal and Leblois, 2008). It is clear that longer length DNA barcodes will provide more efficient identification labels. Barcode efficiency can be further improved by the simultaneous use of genes showing different evolutionary rates and genomic positions.

RHOopsin, a nuclear gene encoding a transmembrane Gprotein- coupled receptor for visual transduction cascade, proved to be very suitable for barcoding purposes because it has a single copy in the genome, does not contain introns in Teleostei fishes, its evolutionary rate in fish is less than 2-fold lower than, for example, cyt b gene (0.167 vs. 0.247), and RHO trees contain a high number of well supported Teleostei clades (Chen et al., 2003).

Here we examine the diversity of two genes with different phylogenetic signal: relatively variable mitochondrial COI and conservative nuclear RHOopsin (RHO) within and among 9 Mugilid species belonging to 4 genera, many of which have been examined from multiple specimens, with the goal of determining whether DNA barcoding can achieve unambiguous species recognition of Mugilid species.

Grey mullets of the family Mugilidae (Pisces, Mugiliformes) are distributed worldwide and inhabit marine, estuarine, and freshwater environments in all tropical and temperate regions. Various Mugilid species are commercially important species in fishery and aquaculture of many countries. Mugilid taxonomy still has not been finalized. The family includes from 14 to 20 genera recognized as valid according to the most recent revisions and 55 currently recognized species (Thomson, 1997; Nelson, 2006). Most of them are representatives of Liza, Chelon, and Mugil genera.

One to seven specimens were collected for each of the nine species of Mugilidae (Table 1). DNA extracts were prepared from either muscle or heart tissue or fin clips preserved in 95% ethanol alcohol according to the protocol of Sambrook et al. (1989). Amplification of the COI gene fragments was carried out using the primers: FishF1 (5′-TCAACCAACCACAAAGACATTGGCAC-3′) and FishR1 (5′-TAGACTTCTGGGTGGCCAAAGAATCA-3′) (Ward et al., 1994). Primers RHOO545 (5′-GCAAGCCCATCAGCAACTTCCG-3′) and RHOO1039 (5′-TGCTTGTTCATGCAGATGTAGA-3′) (Chen et al., 2003) were used to amplify the RHO gene fragments.

PCRs were conducted for both genes with a total volume of 25 μL consisting of approximately 50 ng of template DNA, 0.25 μM of each primer, 2.5 μL of 10× PCR reaction buffer, 2 mM of each dNTP, and 1U of Taq DNA polymerase (Sibenzyme, Russia). PCR consisted of an initial denaturation step at 94℃ for 5 min., followed by 30 cycles of denaturation at 94℃ for 30 s, annealing for 30 s at 60℃ for COI gene fragment and 55℃ for RHO gene fragments and extension at 72℃ for 1 min. The terminal extension was at 72℃ for 10 min. Aliquots (3 μL) of amplicons were examined in 1.5% gels, stained with ethidium bromide, and photographed upon transillumination. DNA fragments were cut from the agarose gels, frozen-thawed, and re-PCRed under the same conditions. The rePCRed products were precipitated with ethyl alcohol and vacuum dried. Then they were sequenced under conditions recommended by the manufacturer, using the BigDye Terminator v.3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, CA, USA). Electrophoreses of sequencing products were made in ABI Prism 3130 DNA sequencer in 50 cm capillary array with polymer POP-7.

The sequences of COI and RHO gene fragments were analyzed both separately and together. DNA sequences were aligned using Clustal W 1.8 in MEGA 5.05. Evolution models and parameters were separately selected for COI and RHO and were then used in the following analyses.

COI sequence divergences were calculated using the Hasegawa, Kishino, and Yano (HKY85+ G) model (Hasegawa et al., 1985) as the appropriate model of sequence evolution chosen on the basis of hierarchical likelihood ratio tests (LRTs) as implemented in ModelTest 3.8 (Posada and Crandall, 1998). RHO sequence divergences were estimated using MEGA 5.05 by Kimura 2-parameter model. Saturation levels of each gene fragment and the combined data set were assessed by plotting transitions and transvertions accumulation for each pair of haplotypes against Hasegawa, Kishino, and Yano distances for COI and Kimura 2 distances for RHO in DAMBE 4.5.2 and MEGA 5.05.

Data were analyzed using three approaches, neighbour joining (NJ), maximum parsimony (MP), and Bayesian inference (BI). NJ, MP and BI trees were created to provide a graphic representation of the patterning of divergence between species computed in MEGA ver. 5.05 and MrBayes 3.2.1. MP analyses were performed using heuristic searches with 50 random stepwise additions and tree bisection-reconnection branch swapping. Bootstrap analyses were used to assess the relative robustness of branches with 1,000 replicates. BI analyses were conducted in MrBayes 3.2.1 with the selected best fit models and parameters. Each BI analysis was run over 4,000,000 generations using four Markov Chain Monte Carlo chains and every 200th tree saved. Four thousand trees were discarded as burn-in. The robustness of trees was tested using posterior probability.

Scomberomorus cavalla (Liu, 2010) and Tribolodon hakonensis were included in all analyses as outgroups. The sequences of Scomberomorus cavalla were taken from GenBank (Table 1).

All sequences are available in Genbank under the accession

numbers shown in the Table 1.

A total of 9 Mugilid species were analyzed, generating (because of multiple specimens for most species) a total of 60 sequences. Approximately 652 bp were amplified and sequenced from the 5′ region of the COI gene from mitochondrial DNA. No insertions, deletions or stop codons were observed in any sequence. One hundred ninty-eight variable sites were revealed, most of them (88.9%) were attributable to the the 3rd codon base, 10.1% to the 2nd, and 1% to the 1st codon base. Fourteen nucleotide substitutions resulted in amino acid replacement.

The GC content of the Liza-Chelon-Moolgarda species was higher of M. cephalus (48.1% vs. 45.0%) due to a higher GC content of the first (58% vs. 55%) and, especially, second (54% vs. 43%) and third (44% vs. 32%) codon base in the former.

The HKY85 distance of individuals within L. haematocheila was 0.17% compared to 3.22% within M. cephalus specimens. Mean divergence among congeneric species was 12.54%. Mean divergence among species within the family increases to 21.8% Mugil vs. Liza, and Mugil vs. Chelon, 20.5% Mugil vs. Moolgarda, 18.2% Moolgarda vs. Liza, and Moolgarda vs. Chelon, and only 10.6% Liza vs. Chelon.

RHOopsin sequences were 515 bp length. All RHOopsin sequences were species-specific. We observed no intraspecific variability, except one nucleotide substitution in 3rd codon position in L. aurata individuals from Azov and Mediterranean Seas that did not result in amino acid replacement.

Forty-four variable sites were revealed, 28 of them were found to be parsimony-informative, most of them (67%) were attributable to the 3rd codon base, 5% to the 2nd, and 28% to the 1st codon base. Eight nucleotide substitutions resulted in amino acid replacement. GC content was about the same in all the species studied-53-55%. The level of divergence due to conservativeness of RHO gene was much lower compared to COI ranging from only 0.4% between L. saliens and C. labrosus and 5% between L. aurata and M. cephalus.

Although DNA barcoding aims to develop species identification systems, some phylogenetic signal was apparent in the data. No saturation was detected for each dataset of COI and RHO sequences. The plots of transversion and transitions against genetic distances were almost linear. Therefore, all substitutions in COI and RHO gene fragments were used for phylogenetic reconstructions.

In all cases, no significant differences among the NJ, MP, and BI tree topologies based on COI and RHO gene sequences were found. The consensus tree obtained from Bayesian analysis of COI dataset was completely identical in topology to the one produced using NJ and MP analysis. The topologies of NJ, MP, and BI trees based on RHO gene analyses were very similar to each other. Thus only the Bayesian tree based on combined COI and RHO gene fragments was presented together with major nodal support values for NJ and MP (Fig. 1).

Haplotypes of the same species were always clustered together in both phylogenetic reconstructions based on COI, RHO and combined COI and RHO gene fragments dataset. All haplotypes of M. cephalus sampled from the seas of Azov, Mediterranean and Japan and near the water of Taiwan formed the monophyletic cluster on the phylogenetic tree based on COI dataset as well as RHO sequences with high bootstrap and posterior probabilities values (BP, PP). The haplotypes of M. cephalus appeared basal as sister group to all other Mugilidae on both phylogenetic trees (Fig. 1). L. haematocheila haplotypes sampled in the East Sea and Okhotsk and near Taiwan were clustered together with high nodal support. The closest relative of L. haematocheila was C. subviridis (Fig. 1). Chelon macrolepis was intermediate between Mediterranean mullets and L. haematocheila-C. subviridis cluster. Chelon labrosus grouped with L. aurata, L. ramada, L. saliens to form a monophyletic clade (Fig. 1), which turned out to exclusively comprise species distributed in Mediterranean- Azov waters. Liza aurata haplotypes sampled in Mediterranean and Azov seas clustered together (Fig. 1). The other Chelon species C. subviridis together with L. haematochalia sampled, all from the Pacific as well as the specimens acclimatized on Azov Sea from the East Sea, formed a distinct clade.

Although the phylogenetic trees based on COI and RHO sequences were significantly similar to each other, there were differences between them. Haplotypes of M. cephalus sampled in Azov-Mediterranean, the East Sea and Taiwan resolved into three clusters on the phylogenetic tree based on COI sequences with high BP and PP support. Haplotypes of L. haematocheila consisted of three lineages, one of which was acclimatized individuals sampled in Azov Sea, and two were from the East Sea and near Taiwan. No variation on RHO sequences was revealed within species M. cephalus and L. haematoheila. Topology of C. macrolepis and M. cunnesius branches was unstable in the phylogenetic tree based on RHO sequences, BP and PP support were low.

Molecular markers are used for explanation and refining the taxonomic and phylogenetic relationships in the groups of the species. Furthermore, molecular markers are used for the accurate and unambiguous identification of fish and fish product from eggs to adult, and it is important in many areas.

Although phylogenetic analyses of Mugilidae species were carried out based on 16S and 12S rRNA mtDNA sequences (Liu et al., 2010), Phe 12SrRNA, cyt b and COI (Heras et al., 2009), three mtDNA loci 16SrRNA, COI and cyt b (Durand et al., 2012), the author supposed that employing slowerevolved and independent nuclear genes might prove helpful to provide molecular systematic of Mugilidae.

For the barcoding we proposed to use not only variable mitochondrial cytochrome oxidase subunit I, but also the highly conserved nuclear RHOopsin gene, which has not previously been used for such purposes. The use of the nuclear genes for the barcoding is especially important in the case of hybridization among species. The mitochondrial COI gene is the most efficient gene marker. However, in some cases 5′-COI “barcoding” region is not enough and even be misleading for purposes of identification due interspecific COI recombination reflecting putative historical hybridization events between species (Balakirev et al., 2012).

All 9 species can be differentiated by both COI and RHO. The lack of stop codons is consistent with all amplified sequences being functional mitochondrial COI sequences, and that, together with the fact that all amplified sequences were about 628-652 bp in length, suggests that NUMTs (nuclear DNA sequences originating from mitochondrial DNA sequences) were not sequenced (vertebrate NUMTS are typically smaller than 600 bp (Zhang and Hewitt, 1996).

The data obtained showed that information based on COI sequences was diagnostic not only for species level identification but also for recognition of intraspecific units, e.g., allopatric populations of circumtropical M. cephalus, whose pronounced population genetic structure is well-known, or even native and acclimatized specimens of L. haematocheila.

COI sequences of each specimen of each sample proved to be identical. We found no difference within each of M. cephalus and L. aurata samples from Azov and Mediterranean Seas. Sequences of geographically remote L. haematocheila (Azov Sea, East Sea, and Taiwan) and M. cephalus (Azov-Mediterranean, East Sea, and Taiwan) had their own haplotypes, 1-2 nucleotide substitutions being observed between the redlip mullet specimens causing 1 amino acid substitution, while 17-31 nucleotide substitutions being revealed between the flathead mullet specimens, which resulted in 4 amino acid substitutions. COI short DNA sequence provides sufficient identification labels in terms of nucleotide positions to discriminate even between congeneric fish species (Hebert et al., 2003) despite the intraspecific variation that was especially pronounced between geographically distant M. cephalus specimens.

All RHOopsin sequences were species-specific. Even M. cephalus specimens that were geographically remote and deeply diverged by mtDNA appeared to be identical. Thus, RHO sequences can be considered useful for barcoding purposes owing to low or zero intraspecific variation and the absence of problems connected with introgressive hybridization as may emerge in the case of mtDNA.

The clades revealed after bootstrapping generally corresponded well with expectations. Topologies of NJ, MP and BI trees conformed to our previous data (Semina et al., 2007) based on PCR-restiction fragment length polymorphismanalysis of extended mtDNA segments.

The data obtained indicated that M. cephalus was the most genetically distant from the rest of the species, which conformed to the results of other authors based on mitochondrial DNA, allozyme, gemoglobin, and karyologic study (Cataudella et al., 1974; Rizotti, 1993; Caldara et al., 1996; Rossi et al., 1998, 2004; Gornung et al., 2001, 2007; Papasotiropoulos et al., 2001, 2002, 2007; Turan et al., 2005; Imsiridou et al., 2007). Three mitochondrial lineages of M. cephalus based on COI sequences had the regional distribution (Fig. 1). One corresponded the Mediterranean lineage included M. cephalus from Morocco, Tunisia and France and two conformed to the lineages revealed in the Pacific region (Durand et al., 2012). We found that all Mediterranean-Azov Chelon and Liza species clustered together as well as all Pacific Chelon and Liza. The results conformed the close genetic relation of Chelon and Liza representatives. Many other researchers also stated the close relationships of these genera. Morphologic (pharyngobranchial organ), cytogenetic, genetic (mitochondrial DNA), allozyme analyses performed by other authors, revealed no substantial difference between Mediterranean Liza species and C. labrosus (Cataudella et al., 1974; Caldara et al., 1996; Gornung et al., 2001, 2007; Murgia et al., 2002, Papasotiropoulos et al., 2002, 2007; Rossi et al., 2004; Turan et al., 2005; Imsiridou et al., 2007; Durand et al., 2012). We consider that molecular phylogeny offered no evidence for distinguishing Liza from Chelon, leading us to challenge their validity. Therefore, taking into account the scientific issue data, our results point out the necessity for taxonomic revision of Chelon and Liza genera.

Based on the data obtained, we conclude that COI, as well as RHO sequencing, can be used to unambiguously identify fish species. Topologies of phylogeny based on COI and RHO sequences coincided with each other, while together they had a good phylogenetic signal.