Characterization of the Complete Mitochondrial Genome Sequences of Three Croakers (Perciformes, Sciaenidae) and Novel Insights into the Phylogenetics

The three croakers (Nibea coibor, Protonibea diacanthus and Argyrosomus amoyensis, Perciformes, Sciaenidae) are important commercial species inhabiting the Eastern Indian Ocean and Western Pacific. Molecular data employed in previous research on phylogenetic reconstruction have not been adequate and complete, and systematic and comprehensive phylogenetic relationships for these fish are unresolved. We sequenced the complete mitochondrial genomes of the three croakers using next-generation sequencing for the first time. We analyzed the composition and phylogenies between 19 species in the family Sciaenidae using the mitochondrial protein coding sequences of 204 species in the Series Eupercaria. We present the characterization of the complete mitochondrial genome sequences of the three croakers. Gene arrangement and distribution of the three croakers are canonically identical and consistent with other vertebrates. We found that the family Sciaenidae is an independent branch that is isolated from the order Perciformes and does not belong to any extant classification. Therefore, this family is expected to belong to a new classification at the order level and needs further analysis. The evolution of Sciaenidae has lagged far behind the Perciformes differentiation. This study presents a novel insight into the phylogenetics of the family Sciaenidae from the order Perciformes and facilitates additional studies on the evolution and phylogeny of Series Eupercaria.

They feed on crustaceans and small fish and have high commercial value due to their nutritional and medicinal properties [4][5][6][7].
Mitochondrial DNA (mtDNA) is an important model system for studying molecular evolution, phylogeny, and genome structure [8,9]. The size of the mitochondrial genome is approximately 15-20 kb and contains 13 protein-coding genes (PCGs), two ribosomal RNAs genes (12S rRNA gene and 16S rRNA gene), 22 transfer RNA (tRNA) genes, and a major non-coding region that contains the initial sites for mtDNA replication, and RNA transcription [10]. Fish mitochondrial genomes have a similar organization to other vertebrates [11]. Because of a constant gene content, multiple copy status in a cell, and a lack of recombination and paralogues, this is of particular interest as the mitochondrial genome is a simpler system than the nuclear genome for studying the molecular dynamics and mechanisms of rearrangements that underlie variations in the genome. Other features of mtDNA are its small size, cellular abundance, maternal inheritance, compact gene arrangement, and high rate of evolution. Therefore, these make the mitochondrial genome a particularly attractive tool for studying evolutionary biology [12,13].
Thus, mtDNA is generally considered a good molecular marker for phylogenetic analyses among fish taxa. However, short mitochondrial gene fragments exhibit limitations in resolving complicated phylogenetic relationships in many fish lineages [14]. The additional informative sites from longer DNA sequences (e.g., mitochondrial genomes) allow these deeper branches and higher-level relationships to be more fully resolved [15]. Hence, the mitochondrial genomes provided in this study may help resolve the evolutionary relationships of the family Sciaenidae and the Series Eupercaria and lead to taxonomic revisions.
There are only 9 studies on N. coibor, P. diacanthus, and A. amoyensis [16][17][18][19][20][21][22][23] involving active material extraction, nosography, and mitochondrial phylogenetic and microsatellite variation. There have been some studies concerning mitochondrial genome in the genera Nibea, Protonibea, and Argyrosomu, the family Sciaenidae and order Perciformes. However, investigations of phylogenetic relationships have been limited because they are based on partial mitochondrial sequence fragments, particularly COXI and the 16S rRNA gene, since few mitochondrial genomes have been sequenced [24]. Only a portion of the mitochondrial genomes is used for the taxonomic classifications, and none have yet been comprehensively described [25][26][27][28][29][30]. Systematic and comprehensive phylogenetic relationships within the Series Eupercaria have not been resolved and the molecular data employed in previous research on phylogenetic reconstruction is incomplete. Therefore, we sequenced the complete mitochondrial genomes of three croakers using a next-generation sequencing (NGS) strategy in this study.
Firstly, our interest in this study was to compare the complete mitochondrial genome sequences of these fish which were sequenced in our study for insights into their characteristics, including their genome organization and composition, gene content, functional regions and the evolutionary dynamics and molecular mechanisms that have shaped them. Beyond this, these completely sequenced mitochondrial genomes provide a broad and useful molecular resource.
Secondly, the family Sciaenidae is a diverse and commercially important family, one of the largest within the Perciforme, comprising 68 genera and about 311 species [31]. We took a detailed view of the composition and phylogeny in 19 species including representatives from 12 of the entire16 genera in family Sciaenidae, based on the mitochondrial genomes presented in GenBank.
Thirdly, Series Eupercaria (=Percomorpharia in previous versions of the classification) possesses more than 6000 species arranged in 161 families and at least 17 orders [32], and is by far the largest series of Percomorphs. Some of the most diverse orders of fish are included in this group. Previous molecular studies obtained monophyletic groups with a combination of taxa assigned to Eupercaria, but including a far more limited sampling [33][34][35]. Although most family-level and ordinal groups within this series receive high nodal support, interrelationships among them are largely unresolved [36]. As currently circumscribed, the order Perciformes is the largest group within Eupercaria and has been traditionally regarded as a "taxonomic waste basket" [31]. Therefore, we preformed the phylogenetic analyses and divergence times based on the concatenated alignment of 13 PCGs amino acid sequences of 204 species of Series Eupercaria in this study.
To sum up, according to the above description and situation, not only does this article provide new insight into novel mitochondrial genome organization among the family Sciaenidae, but it also proposes a monophyletic definition based on robust molecular analyses and a revised taxonomic circumscription of Perciformes and Series Eupercaria. All information reported in this article may facilitate further investigation of the population genetic structures, phylogenetic studies and molecular evolution of Sciaenidae. This work also provides important new molecular resources for the species identification, fishery management, and conservation biology with respect to Sciaenidae.
The mitogenomes of A. amoyensis, N. coibor, and P. diacanthus were considerably shorter than those of Johnius belangerii and Johnius grypotus but were consistent with the average of other Sciaenidae species (Table 1). Intergenic spacers in N. coibor (13), P. diacanthus (14), and A. amoyensis (11) involving 90, 130, and 67 bp, respectively, indicated a greater diversity in the locations and intergenic nucleotides compared to the overlaps. A tRNA-Phe-12S rRNA intergenic spacer (45 bp) of P. diacanthus was the largest and distinctly different with those of A. amoyensis (0 bp) and N. coibor (0 bp), which considerably contributed to the longest complete mitogenome of P. diacanthus among the three croakers in this study ( Table 2). The 34-36 bp intergenic spacers located in a common cluster of five tRNA genes between tRNA-Asn and tRNA-Cys were secondary and almost identical between the three croakers (WANCY region, Table 2).
The overlaps of the mitochondrial genes were 20 bp for N. coibor, 29 bp for P. diacanthus, and 29 bp for A. amoyensis and occurred in 7 different locations at 1-10 bp, resulting in a compact mitochondrial genomic structure. Most overlaps were detected in tRNA-Ile-tRNA-Gln, tRNA-Gln-tRNA-Met, COXI-tRNA-Ser, ND4L-ND4, and ND5-ND6. Generally, we allowed gene overlaps between adjacent genes but seldom between PCGs and tRNAs. When a full termination codon (TAA or TAG) caused an overlap between a protein-coding gene and a tRNA, we annotated this gene using an incomplete termination codon T/TA rather than as an overlap. We annotated full termination codons, however, in overlapping PCGs (Table 2).

Nucleotide Composition
The overall nucleotide compositions of the H-strand of the three croakers in descending order were 31.24% C, 27.13% A, 25.28% T, and 16.35% G, with an A+T content of 52.41% (varying between 52.07-52.88%), indicating significant strand asymmetry or strand-specific bias. This is commonly observed in most vertebrates [29]. The GC-/AT-skews analysis indicated that the H-strand possessed an overrepresentation of C and A, and consequently contained a lower number of G and T bases. Strand asymmetry due to nucleotide composition was also reflected in the codon usage of genes oriented in opposite directions. Genes encoded on H-strand showed a clear preference for C in codon wobble position, whereas G or T ending codons were overrepresented in the genes on the L-strand. The underlying mechanism responsible for the strand bias has generally been interpreted as evidence of an asymmetrical directional mutation pressure. This is associated with replication processes as one strand remains transiently in a single-stranded state, making it more vulnerable to DNA damage [49,50].
The highest A+T content was detected in the control region (64.08%) of P. diacanthus consistent with previous reports on other teleosts, but this occurred in tRNA-His (65.22%) in N. coibor and tRNA-Pro (65.71%) in A. amoyensis [51,52].
We compared the mitochondrial genomic nucleotide compositions of 19 Sciaenidae species. Strand asymmetry due to nucleotide composition could be described by AT and GC skew values, and the AT skews among these 19 species near zero, while their GC skews were all negative, except for J. grypotus and J. belangerii, which were the converse ( Table 1). The A content of most species was only slightly higher than T, whereas C was considerably more prevalent than G. Such skews towards a particular nucleotide were attributed to differential mutational pressures imposed on the L-and H-strands resulting from the asymmetric replication of mtDNA. It was apparent that A and C were more prevalent across the entire mitogenomes of most of these Sciaenidae species, similar to previous observations of a bias against the use of G [47,53]. The diversity of the PCGs showed a different situation. The AT skews of 16 of these 19 species were barely less than zero and was the converse for the entire genome. The GC skews were all negative except for J. grypotus and J. belangerii, which were consistent with the entire genome. The T content of most species was only slightly higher than A, whereas C was considerably more prevalent than G in the PCGs. It was apparent that T and C were more prevalent in the PCGs of most of the Sciaenidae species. This was again similar to previous observations of a bias against the use of G [47,53]. Notably, J. grypotus and J. belangerii had the longest complete mitochondrial genomes and possessed unusual nucleotide skews ( Table 1). The reasons for this are unclear at present.
The AT skews in the 13 PCGs of the three croakers were waving near the zero, and the majority of them were negative. The GC skews were all negative except for the ND6 gene. This situation was consistent with the other 16 species, but, interestingly, the AT skews for ND6 ranged from −0.4757 to −0.333 and GC skews ranging from 0.4353 to 0.4684 ( Figure 2). These values were unusual but similar to other observations of strand asymmetry [47,53].
The diversity of the PCGs showed a different situation. The AT skews of 16 of these 19 species were barely less than zero and was the converse for the entire genome. The GC skews were all negative except for J. grypotus and J. belangerii, which were consistent with the entire genome. The T content of most species was only slightly higher than A, whereas C was considerably more prevalent than G in the PCGs. It was apparent that T and C were more prevalent in the PCGs of most of the Sciaenidae species. This was again similar to previous observations of a bias against the use of G [47,53]. Notably, J. grypotus and J. belangerii had the longest complete mitochondrial genomes and possessed unusual nucleotide skews ( Table 1). The reasons for this are unclear at present.
The AT skews in the 13 PCGs of the three croakers were waving near the zero, and the majority of them were negative. The GC skews were all negative except for the ND6 gene. This situation was consistent with the other 16 species, but, interestingly, the AT skews for ND6 ranged from −0.4757 to −0.333 and GC skews ranging from 0.4353 to 0.4684 ( Figure 2). These values were unusual but similar to other observations of strand asymmetry [47,53].

Protein-Coding Genes
The 13 PCGs in the mitogenomes of the three croakers were similar to those of other vertebrates. These included 3 cytochrome c oxidase complexes (COX I-III), 7 NADH ubiquinone oxidoreductase complex (ND1-6, ND4L), a cytochrome b oxidoreductase complex (Cyt b), and 2 ATP synthases (ATP6 and ATP8). The coding regions ranged in size from 168 (ATP8) to 1839 bp (ND5) comprising 11,404-11,435 bp accounting for approximately 69% of the entire mitogenomes. The majority of these genes had similar lengths except for ATP6 and ND5 that possessed 27 and 15 bp gaps, respectively. All protein-coding genes were located on the H strand except for ND6.
Initiation and termination codons were determined based on alignments with the corresponding genes and proteins of other Sciaenidae fish. The mitogenomes of the three croakers exhibited a canonical genetic code shared by all vertebrates. An orthodox initiation codon methionine (ATG) was used for most of the PCGs, while GTG was utilized in ND1 of N. coibor and in ATP6 of P. diacanthus, ATA in ATP6 of N. coibor, and ATC in ND5 of P. diacanthus. These were also canonical mitochondrial initiation codons for vertebrate mitogenomes [54].
The termination codons AGA, T, TA, TAA, and TAG were present in all three species, and T and TAA were the most prevalent. Incomplete termination codons using either TA or T would be completed as a TAA via post-transcriptional polyadenylation. A diverse pattern of codon usage within stop codons seems to be a common tendency in fish mitogenomes. Since the three croakers were AT-rich organisms, we expected that A or T ending codons would predominate. The codon usage patterns in the 13 PCGs of the three croaker mitogenomes was similar to bony fish (Osteichthyes) [55].
In all 13 PCGs of the three croakers, the average non-synonymous/synonymous mutation ratio (Ka/Ks) was 0.0671 and varied from 0.0100 (COXI between AA-PD) to 0.2714 (ND5 between the NC-PD) (Table S2, Figure 3). Two non-adaptive forces, random genetic drift, and mutation pressure define the fundamental features of genome evolutional though functional constraints imposed burdens on mutation events [56]. Therefore, the mutation-associated disadvantages were difficult to establish under purifying selection. The selection processes maintained long-term stability of the biological structure. Our Ka/Ks ratio indicated that the functional genes evolved under strong purifying selection, which meant natural selection against deleterious mutations with negative selective coefficients [57]. The selection pressures differed among the genes and were likely to evolve in different ways [56]. Additionally, the ND5 genes had the highest Ka/Ks values, indicating that the selection pressures were strand-independent.

Protein-Coding Genes
The 13 PCGs in the mitogenomes of the three croakers were similar to those of other vertebrates. These included 3 cytochrome c oxidase complexes (COX I-III), 7 NADH ubiquinone oxidoreductase complex (ND1-6, ND4L), a cytochrome b oxidoreductase complex (Cyt b), and 2 ATP synthases (ATP6 and ATP8). The coding regions ranged in size from 168 (ATP8) to 1839 bp (ND5) comprising 11,404-11,435 bp accounting for approximately 69% of the entire mitogenomes. The majority of these genes had similar lengths except for ATP6 and ND5 that possessed 27 and 15 bp gaps, respectively. All protein-coding genes were located on the H strand except for ND6.
Initiation and termination codons were determined based on alignments with the corresponding genes and proteins of other Sciaenidae fish. The mitogenomes of the three croakers exhibited a canonical genetic code shared by all vertebrates. An orthodox initiation codon methionine (ATG) was used for most of the PCGs, while GTG was utilized in ND1 of N. coibor and in ATP6 of P. diacanthus, ATA in ATP6 of N. coibor, and ATC in ND5 of P. diacanthus. These were also canonical mitochondrial initiation codons for vertebrate mitogenomes [54].
The termination codons AGA, T, TA, TAA, and TAG were present in all three species, and T and TAA were the most prevalent. Incomplete termination codons using either TA or T would be completed as a TAA via post-transcriptional polyadenylation. A diverse pattern of codon usage within stop codons seems to be a common tendency in fish mitogenomes. Since the three croakers were AT-rich organisms, we expected that A or T ending codons would predominate. The codon usage patterns in the 13 PCGs of the three croaker mitogenomes was similar to bony fish (Osteichthyes) [55].
In all 13 PCGs of the three croakers, the average non-synonymous/synonymous mutation ratio (Ka/Ks) was 0.0671 and varied from 0.0100 (COXI between AA-PD) to 0.2714 (ND5 between the NC-PD) (Table S2, Figure 3). Two non-adaptive forces, random genetic drift, and mutation pressure define the fundamental features of genome evolutional though functional constraints imposed burdens on mutation events [56]. Therefore, the mutation-associated disadvantages were difficult to establish under purifying selection. The selection processes maintained long-term stability of the biological structure. Our Ka/Ks ratio indicated that the functional genes evolved under strong purifying selection, which meant natural selection against deleterious mutations with negative selective coefficients [57]. The selection pressures differed among the genes and were likely to evolve in different ways [56]. Additionally, the ND5 genes had the highest Ka/Ks values, indicating that the selection pressures were strand-independent. We evaluated the conservation of mtDNA genes based on the overall p-genetic distance among 19 Sciaenidae species (Figure 4) [58]. Of the 13 PCGs, the COXIII gene has the lowest overall p-genetic distance (0.036) and the ATP8 and ND5 genes had the highest value (0.133) based on data of the first and second nucleotides of codons. A full-length sequence comparison of each gene resulted in the lowest value for COXIII (0.147) and the highest for ND5 (0.229). These data were consistent with the values based on data of the first and second nucleotides of codons. According to these results, ND5 most likely had the fastest evolutionary rate among Sciaenidae species, while COXIII had the lowest. For the third or wobble nucleotide, all genes had a high overall p-genetic distance value, ranging from 0.372 to 0.460. As was the case for the other fishes, most of the differences in the mtDNA PCGs occurred at the wobble position [48]. We evaluated the conservation of mtDNA genes based on the overall p-genetic distance among 19 Sciaenidae species (Figure 4) [58]. Of the 13 PCGs, the COXIII gene has the lowest overall p-genetic distance (0.036) and the ATP8 and ND5 genes had the highest value (0.133) based on data of the first and second nucleotides of codons. A full-length sequence comparison of each gene resulted in the lowest value for COXIII (0.147) and the highest for ND5 (0.229). These data were consistent with the values based on data of the first and second nucleotides of codons. According to these results, ND5 most likely had the fastest evolutionary rate among Sciaenidae species, while COXIII had the lowest. For the third or wobble nucleotide, all genes had a high overall p-genetic distance value, ranging from 0.372 to 0.460. As was the case for the other fishes, most of the differences in the mtDNA PCGs occurred at the wobble position [48].

Mitochondrial Gene Codon Usage
The amino acids Leu and Ser were utilized by six different codons, while all other amino acids were encoded by either two or four. The frequencies in the PCGs of Leu (CTC, CTA) and Ala (GCC) were the most frequently used and Ser (TCG), Arg (CGT), and Cys (TGT) were least frequently used. However, relative synonymous codon usage (RSCU) analysis revealed that the codons coding Ala (GCC), Pro (CCC), Ser (TCC), and Arg (CGA) were the most frequently present, while those coding Pro (CCG), Leu (TTG), Thr (ACG), and Ala (GCG) were rare (Table S1 and Figure 5). Therefore, the codon frequency and RSCU were totally different indexes to value the mitochondrial gene codon usage. Only when the species had a similar codon usage preference would codon distribution and amino acid content be consistent among them. Moreover, differences in A+T content and the AT and GC skew of the PCGs were reflected in the codon usage (Table 1; Figure 2). They were closely related and mutually consistent [47,48,53]. The RSCU values indicated that codons with A or C in the third position were more used than T or G, so the codons NNA and NNC were in majority, while the synonymous codons NNT and NNG were in the minority.

Mitochondrial Gene Codon Usage
The amino acids Leu and Ser were utilized by six different codons, while all other amino acids were encoded by either two or four. The frequencies in the PCGs of Leu (CTC, CTA) and Ala (GCC) were the most frequently used and Ser (TCG), Arg (CGT), and Cys (TGT) were least frequently used. However, relative synonymous codon usage (RSCU) analysis revealed that the codons coding Ala (GCC), Pro (CCC), Ser (TCC), and Arg (CGA) were the most frequently present, while those coding Pro (CCG), Leu (TTG), Thr (ACG), and Ala (GCG) were rare (Table S1 and Figure 5). Therefore, the codon frequency and RSCU were totally different indexes to value the mitochondrial gene codon usage. Only when the species had a similar codon usage preference would codon distribution and amino acid content be consistent among them. Moreover, differences in A+T content and the AT and GC skew of the PCGs were reflected in the codon usage (Table 1; Figure 2). They were closely related and mutually consistent [47,48,53]. The RSCU values indicated that codons with A or C in the third position were more used than T or G, so the codons NNA and NNC were in majority, while the synonymous codons NNT and NNG were in the minority.

Transfer and Ribosomal RNA Genes
The mitogenomes of the three croakers contained the 22 tRNA genes that are typically found in metazoan mitogenomes. They varied in size from 66 bp (tRNA-Cys) to 75 bp (tRNA-Lys) and the anticodons were identical to those reported in other vertebrates. Serine were determined by two types of anticodon (TGA, GCT) and leucine by TAG and TAA. With the exception of tRNA-Ser (GCT) of N. coibor, tRNA-Ser (GCT), and tRNA-Asn (GTT) of P. diacanthus, all of the 22 tRNA genes were predicted to be folded into canonical cloverleaf secondary structures, although there were numerous non-complementary and T-G base pairs in the stem regions. Stem mismatches seemed to be a common phenomenon for mitochondrial tRNA genes and are probably repaired via a post-transcriptional editing process [59]. The tRNA-Ser (GCT) found in the N. coibor and P. diacanthus mitogenomes had no recognizable DHU stem, which was the case for almost all vertebrate mitogenomes [60,61]. This tRNA-Ser has been previously described as a pseudogene to explain its unusual structure. That suggested that it could also perform the function by adjusting its structural conformation to fit the ribosome in a similar way to that of usual tRNAs in the ribosom [62].
As in all other mitogenomes described so far, two rRNA genes were present, a small (12S) and a large (16S) subunit. The 12S rRNA gene ranged in size from 907 to 953 bp and the 16S rRNA gene from 1696 to 1720 bp. Both genes were located on the H-strand between tRNA-Phe (GAA) and tRNA-Leu (TAA), separated by tRNA-Val as is the case in most vertebrates [54].
The AT content of the 12S rRNA gene was between 51.42% and 51.95% with an AT-skew of 0.1846-0.2 and a GC-skew of −0.0842 to −0.0502. The AT content of the 16S rRNA gene was between 52.41% to 53.26% with an AT-skew of 0.2174-0.2251 and a GC-skew of −0.1692 to −0.1339. The A and C content was more prevalent in the 12S rRNA and 16S rRNA genes of the three croakers, as has been reported in other bony fish [63]. That was similar to the entire mitogenome, and different with the PCGs. This implied that different kinds of mitogenome genes preferred to different bias.

CG View Comparison Tool (CCT) Map
We compared the complete mitochondrial genomes of N. albiflora (Sciaenidae) with other available mitogenomes of the family Sciaenidae. Both the nucleotide ( Figure 6) and coding DNA sequences (Figure 7) of the mitochondrial genomes illustrated the gene similarity between the three croakers sequenced in our study and the other Sciaenidae species, with the N. albiflora randomly selected as the reference sequence.
There was less identity in the mitogenome nucleotides than the coding DNA sequences (compare Figures 6 and 7). In the 37 genes and the control region (D-loop) of the mitogenome, ND5 was the least conservative in all species. This was appropriate for the phylogenetic and evolutional analyses within the family Sciaenidae. On the contrary, Cytb was the most conservative, which is appropriate for the analyses above the family Sciaenidae (Figure 7).

Phylogenetic Analyses
Based on the concatenated alignment of 13 PCGs amino acid sequences of 204 species in the Series Eupercaria analyzed in this study (Table S3), Bayesian and ML analyses produced almost identical topologies with similar branch lengths and strong bootstraps (ML analysis) and posterior probabilities (Bayesian inference) values (Figure 8). We preferred the phylogenetic analyses based on the 13 amino acid sequences of PCGs rather than the mitochondrial genome nucleotides because the 204 species of the Series Eupercaria belonged to different orders and possessed wide-ranging taxonomies. By multiple sequence alignment of their amino acid and nucleotide sequences, we figured out that nucleotide sequences were more conservative and suitable for the phylogenetic analyses of the intensive taxonomy [64][65][66]. To summarize, whether the 13 amino acid sequences or nucleotides of the 13 PCGs of mitochondrial genome is appropriate for the phylogenetics dependent upon the specific species adopted.

Phylogenetic Analyses
Based on the concatenated alignment of 13 PCGs amino acid sequences of 204 species in the Series Eupercaria analyzed in this study (Table S3), Bayesian and ML analyses produced almost identical topologies with similar branch lengths and strong bootstraps (ML analysis) and posterior probabilities (Bayesian inference) values (Figure 8). We preferred the phylogenetic analyses based on the 13 amino acid sequences of PCGs rather than the mitochondrial genome nucleotides because the 204 species of the Series Eupercaria belonged to different orders and possessed wide-ranging taxonomies. By multiple sequence alignment of their amino acid and nucleotide sequences, we figured out that nucleotide sequences were more conservative and suitable for the phylogenetic analyses of the intensive taxonomy [64][65][66]. To summarize, whether the 13 amino acid sequences or nucleotides of the 13 PCGs of mitochondrial genome is appropriate for the phylogenetics dependent upon the specific species adopted. The phylogenetic tree analysis indicated that all the species of the Sciaenidae clustered on the same branch, and all the posterior probabilities were 1 except for Miichthys (0.8041) and Argyrosomus (0.9884). Furthermore, the intergeneric and interspecific taxonomic positions were explicit and clear.
Unexpectedly, the overall taxonomy of the whole Sciaenidae branch based on the mitochondrial genome classification was not consistent with the traditional Chinese fish taxonomy references [7,67]. The family Sciaenidae is traditionally classified as Perciformes consistent with Fishbase. However, the NCBI taxonomy software based on the DNA sequences revealed that the Sciaenidae position in the Eupercaria is uncertain. Algorithms used in the DeepFin program based on 20 nuclear genes and a mitochondrial gene show that the Sciaenidae have an uncertain placement in the Percomorpharia [32,36,68]. In addition, the Percomorpharia are at a higher taxonomic level than the Eupercaria. Therefore, although there are great differences between traditional DNA molecular marker taxonomies, the analyses using the coding gene sequences from the mitochondrial genome in our manuscript were closer to the ones in NCBI and DeepFin. The phylogenetic tree analysis indicated that all the species of the Sciaenidae clustered on the same branch, and all the posterior probabilities were 1 except for Miichthys (0.8041) and Argyrosomus (0.9884). Furthermore, the intergeneric and interspecific taxonomic positions were explicit and clear.
Unexpectedly, the overall taxonomy of the whole Sciaenidae branch based on the mitochondrial genome classification was not consistent with the traditional Chinese fish taxonomy references [7,67]. The family Sciaenidae is traditionally classified as Perciformes consistent with Fishbase. However, the NCBI taxonomy software based on the DNA sequences revealed that the Sciaenidae position in the Eupercaria is uncertain. Algorithms used in the DeepFin program based on 20 nuclear genes and a mitochondrial gene show that the Sciaenidae have an uncertain placement in the Percomorpharia [32,36,68]. In addition, the Percomorpharia are at a higher taxonomic level than the Eupercaria. Therefore, although there are great differences between traditional DNA molecular marker taxonomies, the analyses using the coding gene sequences from the mitochondrial genome in our manuscript were closer to the ones in NCBI and DeepFin.
The Sciaenidae have evolved phyletically and were closer to the orders Lobotiformes and Ephippiformes and at the same classification level (Figure 8). Tracing back at the branch nodes of the evolutionary tree from Sciaenidae to Perciformes, they have experienced six branch nodes with node bootstrap values of 1, 1, 0.9982, 0.9791, 0.9974, and 1. Sciaenidae are the most distant from the order Perciformes and the high posterior probabilities of each branch support that the Sciaenidae are an independent branch, isolated from the order Perciformes. Sciaenidae have been included in the orders Acanthuriformes, Acanthuridae, Emmelichthyidae, Luvaridae, and Zanclidae [31]. However, there are exceptions to this scheme where the placement of Emmelichthyidae and Sciaenidae in the Acanthuriformes was not supported, consistent with our viewpoint [32].
In summary, synthetically considering the taxonomy of Sciaenidae and the relationship between the Sciaenidae and the Perciformes, we predict that the Sciaenidae belongs to a new unknown order rather than the extant order Perciformes.

Divergence Times of the Series Eupercaria
We estimated divergence times based on phylogenetic trees constructed with the Bayesian method, and the species differentiation time estimated using the ML method of RelTime MEGA 7.0 [69] ( Figure 9). In the ProtTest 3 model preferences [70], the model MtMam+I+G+F was optimal with an AIC value of 650,366.12. The model JTT+I+G+F was suboptimal with an AIC value of 653,040.01, similar to the one for model MtMam+I+G+F. Our software did not support the MtMam+I+G+F model so we used the JTT+I+G+F model for the molecular clock. The Sciaenidae have evolved phyletically and were closer to the orders Lobotiformes and Ephippiformes and at the same classification level (Figure 8). Tracing back at the branch nodes of the evolutionary tree from Sciaenidae to Perciformes, they have experienced six branch nodes with node bootstrap values of 1, 1, 0.9982, 0.9791, 0.9974, and 1. Sciaenidae are the most distant from the order Perciformes and the high posterior probabilities of each branch support that the Sciaenidae are an independent branch, isolated from the order Perciformes. Sciaenidae have been included in the orders Acanthuriformes, Acanthuridae, Emmelichthyidae, Luvaridae, and Zanclidae [31]. However, there are exceptions to this scheme where the placement of Emmelichthyidae and Sciaenidae in the Acanthuriformes was not supported, consistent with our viewpoint [32].
In summary, synthetically considering the taxonomy of Sciaenidae and the relationship between the Sciaenidae and the Perciformes, we predict that the Sciaenidae belongs to a new unknown order rather than the extant order Perciformes.

Divergence Times of the Series Eupercaria
We estimated divergence times based on phylogenetic trees constructed with the Bayesian method, and the species differentiation time estimated using the ML method of RelTime MEGA 7.0 [69] (Figure 9). In the ProtTest 3 model preferences [70], the model MtMam+I+G+F was optimal with an AIC value of 650,366.12. The model JTT+I+G+F was suboptimal with an AIC value of 653,040.01, similar to the one for model MtMam+I+G+F. Our software did not support the MtMam+I+G+F model so we used the JTT+I+G+F model for the molecular clock.  The divergence time calculations indicated that the Perciformes began to differentiate from other species about 73.49 million years ago and then subsequently evolved as an independent branch. The Sciaenidae differentiation from the Ephippiformes and Lobotiformes lagged far behind Perciformes and began around 51.52 million years ago. Therefore, there was a large gap of 21.97 million years in the differentiation initiation times of the Perciformes and Sciaenidae. The Sciaenidae initiated their differentiation after 4 important differentiation nodes at 65.76, 63.38, 60.15, and 55.92 million years ago and all lagged behind Perciformes differentiation. In summary, the Sciaenidae was an independent branch of fish evolution that was entirely independent with, and lagged far behind, Perciformes differentiation.

Sampling and Materials
The wild specimens of the three croakers (N. coibor, P. diacanthus, and A. amoyensis) were collected from the South China Sea near Guangdong province (N23 • 19 , E117 • 09 ), China, on 3-18 July 2014 (identification code: NC_005072014, PD_006072014, AA_007072014). The dorsal muscle was preserved in 95% ethanol and stored at −70 • C until they were used for DNA extraction. Voucher specimens were deposited in the Germplasm Resources Lab, College of Marine Sciences, South China Agricultural University, Guangzhou, China (accession numbers: 005072014, 006072014, 007072014). Genomic DNA was isolated from muscle tissue as previously described [71]. Total DNA was eluted in sterile deionized water and was stored at −20 • C. All animal experiments were conducted in accordance with the guidelines and approval of the Animal Research and Ethics Committees of South China Agricultural University.

Library Preparation for Sequencing
DNA sample quality and quantity were characterized by gel electrophoresis and the Nano-Drop 2000 spectrometer (Thermo Scientific, Waltham, MA, USA). The high-quality genomic DNA were used to prepare DNA library, with insert sizes of 500 bp for paired-end sequencing. Paired-end reads of 100 bp were generated on an Illumina HiSeq2500 (Illumina, Inc., San Diego, CA, USA) using sequencing protocols provided by the manufacturer.
Codon usage was determined for all PCGs. The Relative Synonymous Codon Usage (RSCU) was obtained using MEGA 7 software [69]. Statistical analyses of the distributions and visualization of codon usage in the form of heat maps were conducted using R language with RSCU values, a measure of non-uniform usage of synonymous codons in a coding sequence [78]. The RSCU value was the number of times a particular codon was observed relative to the number of times that the codon would be observed for a uniform synonymous codon usage (i.e., all codons for a given amino acid exhibiting similar probabilities). The RSCU value in the absence of any codon usage bias is 1.00. A codon used less frequently than expected will have RSCU values <1.00, whereas codons used more frequently than expected will have RSCU values >1.00.

Phylogenetic Analyses
We conducted phylogenetic analyses based on 204 mitochondrial genomes sequences of the Series Eupercaria (Table S3). Cyprinus carpio and Danio rerio (Cypriniformes, Cyprinidae) were used as out-groups. All sequences were deposited in GenBank. We aligned amino acid sequences of PCGs in the 204 species using the MUSCLE program in MEGA 7. Our alignments of individual genes excluded the stop codon and the third codon.
Phylogenetic analysis was performed using the Maximum Likelihood (ML) and Bayesian Inference (BI) methods. The best model MtMam+I+G+F based on the amino acid sequences in this study was selected using AIC (Akaike Information Criteria). Minimum and best values were fitted with ProtTest 3 with optimized parameters [70]. ML analysis was conducted using the RAxML 8.1.5 with 1000 bootstrap replicates based on the best-scoring protein substitution model determined automatically by the software [79]. Bayesian phylogenetic analysis with 10,000,000 generations was carried out with MrBayes 3.2.6 software [80]. Four independent Markov chains were used at the same time with sampling every 100 generations. The BI Tree was reliable since the standard deviation of split frequencies was below 0.01. The phylogenetic trees were drawn with the Evolview (http://www.evolgenius.info/evolview/) [81,82]. Fishbase (http://fishbase.org), the NCBI taxonomy server (http://www.ncbi.nlm.nih.gov/taxonomy/), and DeepFin (http://deepfin.org/) were used in the additional analyses.
The divergence times of the Series Eupercaria was estimated using MEGA 7.0 [83] with the RelTime-ML method and JTT+F+I+G modeling. Divergence times were presented as in the Time Tree database (http://www.timetree.org/) [69]. Bayesian phylogenetic data were input to MEGA 7.0 to generate divergence times that ensured consistency between the phylogenetic tree and divergence times.
Sequence differences between the three croakers and those from other Sciaenidae species using the Nibea albiflora (NC_015205) as the reference sequence. Sequences were aligned using BLAST and annular genetic similarity mapping was visualized using the CGView Comparison Tool (http://stothard.afns.ualberta.ca/downloads/CCT).

Conclusions
We have here presented the characterization of the complete mitochondrial genome sequences of three croakers in the order Perciformes, family Sciaenidae. Gene arrangement and distribution of the three croakers are canonically identical and consistent with other vertebrates. The Sciaenidae is an independent branch because it is isolated from the order Perciformes and does not belong to any currently known order. The differentiation of the independent Sciaenidae is lagging far behind that of Perciformes. This study presents a novel insight into the phylogenetic analysis of the family Sciaenidae from the order Perciformes and facilitates additional studies on the evolution and phylogeny of the Series Eupercaria.

Conflicts of Interest:
The authors report no conflicts of interest. The authors themselves are responsible for the content and writing of the paper.