Complete Mitochondrial Genome of Three Species (Perciformes, Amblyopinae): Genome Description and Phylogenetic Relationships

The vertebrate mitochondrial DNA is a double-stranded circular molecule of 16-18kb in length consisting of genes for 22 transfer RNA genes (tRNA), 2 ribosomal RNA genes (rRNA) and 13 protein-coding genes as well as two non-coding regions: origin of light strand replication (OL) and control region (Bibb et al., 1981). The usefulness of mitogonomes was becoming increasingly important in solving the long-term controversial phylogenetic relationship and confusing taxon (Cheng et al., 2012). In recent years, longer mitochondrial DNA sequences were used to reconstruct higher level phylogenetic relationships (Boore et al., 2005) and those studies were necessary to resolve these controversial problems.

The suborder Gobioidei consisted of six divided families (Thacker, 2009), and Gobiidae was the largest family including five subfamilies (Gobiidae, Gobionellinae, Sicydiinae, Oxudercinae and Amblyopinae) (Hoese, 1984). Moreover, studies on Gobiidae taxonomy have been controversial and attract extensive attention for a long time. The subfamily Amblyopinae are usually inshore, mud-dwelling fishes which are generally given the vernacular name of “eel gobies” or “worm gobies”. The genus Odontamblyopus within Amblyopinae was firstly proposed by Bleeker (1874). Since then, Norman (1966) was the only recent author who clarified Odontamblyopus into more than two species. Murdy and Shibukawa (2001) demonstrated that Odontamblyopus comprised 4 species. In 2003, O. rebecca was found in Vietnam by Murdy and Shibukawa. Up to now, the genus O. consists of 5 species including O. rubicundus, O. roseus, O. lacepedii, O. tenuis and O. rebecca. Due to the similar morphological characters, some taxon problems still exist.

In this paper, we sequenced the complete mitogenome sequences and examined the genomic structure of three species (O. lacepedii, O. rebecca and Odontamblyopus sp.). We also analyzed the main features in terms of the genome organization, gene arrangement and codon usage. Furthermore, combined with 12 protein-coding genes of other Gobiidae fishes, the phylogenetic relationships were conducted using Bayesian and ML analyses. This study aims to solve confusing taxon problems using longer sequences and comparison between gobies.

The three species (O. lacepedii, O. rebecca and Odontamblyopus sp.) were collected from Ariake Bay in Japan and from Zhujiang and Zhoushan in China, respectively (Table I). All samples were morphologically discriminated to the species level based on Wu and Zhong (2008) and Murdy and Shibukawa (2001, 2003). Muscles excised were used for DNA extraction and preserved in 95% ethanol until use.

Whole genomic DNA was extracted from muscle tissue by proteinase-K digestion followed by the standard phenol/chloroform method and used as a template for subsequent PCR reactions to determine the complete mitogenome sequences of O. lacepedii, O. rebecca and Odontamblyopus sp. Six sets of primers were designed for long-PCR amplification based on the aligned mitogenome sequences of Trypauchen vagina (NC_016693) as previously determined. New pairs of primers were designed for the subsequent amplification by primer walking method. Thirty-two normal PCR primer sets were used to accomplish the entire mitogenome. Essentially every contiguous segment overlapped by at least 50 bp to explicit the accuracy of the sequences.

All PCRs were performed in an Eppendorf thermal cycler. TaKaRa Ex-Taq and LA-Taq Kits (Takara Biomedical) were used for normal and long-PCR reactions, respectively. Long-PCR reactions were carried out in 25 μl reaction mixture containing 15.25 μl of sterile distilled H2O, 2.5 μl of 10×Buffer, 4 μl of dNTP, 1μl of each primer (5μM), 0.25 μl of LA Taq polymerase (1 unit/μl, Takara), and 1 μl of DNA template. The long-PCR reactions consisted of an initial denaturing step at 94°C for 2 min, followed by 30 cycles of denaturing at 94°C for 30 s, annealing at 58°C for 3 min and a final extension at 72°C for 10 min. The normal PCR was performed following the standard procedure (Liu et al., 2007). Negative controls were included in all PCR amplifications to conform the absence of contaminants. PCR product was purified with a Gel Extraction Mini Kit (Watson BioTechnologies Inc., China). The purified product was then sequenced on ABI Prism 3730 (Applied Biosystems) from both strands with the same primers as those used for PCRs.

Sequences of overlapping fragments were assembled manually and aligned initially against the complete mitochondrial genome sequences of T. vagina using SEQMAN software and DNASTAR software with default parameters and further adjusted manually. The boundaries for rRNAs and protein-coding genes were determined by DOGMA comparing with published mitochondrial sequences. The codon usage of the 13 protein-coding genes and the base composition of 37 genes for three species were analyzed using the program MEGA 4.0. Translation initiation and translation termination codons were identified using genetic codon table for mitochondrion in MEGA 4.0. The putative OL and control region were determined by sequence homology and proposed secondary structures. The secondary structure of the putative OL was analyzed with the software RNA structure 5.4. The complete mitochondrial sequences of three species (O. lacepedii, O. rebecca and Odontamblyopus sp.) have been deposited in GenBank with accession number KR815520, KT633953 and KT633954.

Phylogenetic relationship using homologous sequences downloaded from GenBank (Table II) were performed. Only 12 protein-coding regions were used for the subsequent phylogenetic analyses with the exception of ND6 (Fig. 1) because of its heterogeneous base composition and consistently poor performance (Miya et al., 2003). Multiple alignments of sequences were performed using Clustal X version 1.8. Alignment of all protein coding genes was triplet and no internal stop codons were found in any fragments. To explore the degree of saturation present in the datasets, we plotted sequence divergence (GTR-based distance) vs. number of transition and transversion substitutions for all pairwise comparisons among taxa for each codon position in DAMBE. If the codon position sites were saturated, we would expect to see a plateau in such a plot, where little or no additional substitution is detectable with increased p distance.

Phylogenetic analyses of the concatenated genes were conducted under Bayesian Inferences (BI) and maximum likelihood (ML) with the programs MrBayes and PAUP 4.0, respectively. Appropriate evolution substitution models were screened by the ‘‘decision-theoretic performance-based’’ approach (DT; Minin et al., 2003) selection strategy in jModeltest 2.1.1 with partition strategies of all nucleotide sites and each codon. The best-fit model of all nucleotide sites was used for reconstructed phylogenetic trees in BI and ML analyses. BI was also conducted for the combined dataset using the optimal models with each codon position. Micropercops swinhonis and Perccottus glenii were chosen as outgroups. A total of 100 ML bootstrap replicates (MLBS) were performed using PAUP. Bayesian inferences using Markov chain Monte Carlo (MCMC) sampling were performed with four chains which ran simultaneously for 2,000,000 generations with tree sampled every 100 generations. The average standard deviation of split support (ASDOSS) and effective sample size (ESS) were used to assess convergence of MCMC. When value of ASDOSS was less than 0.01, the runs were stopped; ESSs detected by Tracer v.1.5 were more than 200, achieving convergence. The burn-in trees sampled prior to convergence (25% of trees from each run; 5000 sampled trees) were discarded, and the subsequent trees sampled independently from the posterior probability distribution were combined to produce phylogram and 50% majority-rule consensus trees.

The complete mitogenome sequences were 17245 bp, 17009 bp and 17004 bp in O. lacepedii, O. rebecca and Odontamblyopus sp., respectively (Supplementary Table I; Fig. 1). The length and gene order of protein-coding genes were highly in agreement with previous reports (Cui et al., 2009; Kim et al., 2004; Wang et al., 2008). The length variation of mitochondrial DNA was mainly due to the different length of control region for most vertebrates (Randi et al., 1998; Ketmaier and Bernardini, 2005).

The overall base compositions of three species were listed in Supplementary Table II. The G content was 15.3% (O. lacepedii), 15.3% (O. rebecca) and 15.4% (Odontamblyopus sp.) showing an obvious bias against G, supporting by previous studies (Miya et al., 2003; Mabuchi et al., 2007; Wang et al., 2008). The A+T content exhibited higher values than G+C content indicating that the codon usage with A and T nucleotides preferred to C and G nucleotides at the third codon position (Supplementary Table II).

The total length of 13 protein-coding genes was 11415 bp, 11412 bp and 11415 bp in O. lacepedii, O. rebecca and Odontamblyopus sp., respectively. As with the common features of other bony fishes (Miya et al., 2003; Cheng et al., 2012), three reading-frames overlap were noted on the same strand: ATPase8 and ATPase6 overlapped by 4 nucleotides, ATPase6 and COIII overlapped by 1 nucleotide, and ND4L and ND4 overlapped by 7 nucleotides (Supplementary Table I). The pair of genes ND5-ND6 encoding on the different strands overlapped by 4 nucleotides.

All the initiation and stop codons were examined according to the corresponding genes and proteins of other goby fishes (Cheng et al., 2012). Most protein-coding genes began with ATG codon (Supplementary Table I) except COI gene and ATP6 gene used the initiation codon GTG and ATA. Differing from initiation codons (Catanese et al., 2010; Oh et al., 2007), complete stop codons and incomplete codons were identified. ATG seemed to be the most common initiation codon, though there were exceptions (Miya et al., 2003). The open reading frames for ND6 genes ended with TAG in O. lacepedii and O. rebecca, however, Odontamblyopus sp. ended with TAA. All other reading frames of three species ended with TAA (ND1, ND2, ND4L and ND5), AGG (COI) and TAG (ATPase8). Six protein-coding genes ended with incomplete stop codons TA (ATPase6), as well as in complete stop codons T (COII, COIII, ND3, ND4 and Cytb). Based on Ojala et al. (1981), the presence of such noncanonical stop codons can be converted into a fully functional TAA stop codon via post-transcriptional polyadenylation.

Comparative base composition analysis at each codon position (Supplementary Table II) reflected that the proportion of G at third codon position showed a relatively low value (Miya et al., 2003; Mabuchi et al., 2007; Oh et al., 2007) and a clearly excess of pyrimidines over purines probably as a result of the hydrophobic character of the proteins (Naylor et al., 1995). The patterns of codon usage were exhibited in Figures 2, 3 and 4. Codons for Leucine accounted for the highest percent value, however codon usage for Cysteine were the least (Fig. 2). The overall codons ending with A or C were all used more frequently than those ending with T and G (Fig. 3). Moreover, evidence also can be observed codon usage of genes oriented in opposite directions. For amino acids with fourfold degenerate codons, third position ending with A were the most frequently employed, followed by codon families ending with C and U. Otherwise, among twofold degenerate codons, C seemed to be more often used than T in pyrimidine codon families; purine codons encoded mostly with A (Fig. 4). The codon usage pattern of 13 protein-coding genes of the three species was in common with other fishes (Cheng et al., 2010).

The mitogenomes of three species examined included 22 typital tRNA genes dispersed between rRNA and protein-coding genes (Table II). The tRNA genes varied from 64 bp to 76 bp and gene inversion of tRNA genes has been observed in some fishes (Lin et al., 2006; Zhu and Yue, 2008). Among these tRNA genes, two forms of tRNA-Leu (UUR and CUN) and tRNA-Ser (UCN and AGY) were determined (Table II). The three tRNA clusters-IQM (isoleucine, glutamine and methionine), WANCY (tryptophan, alanine, asparagines, cysteine and tyrosine) and HSL (histidine, serine and leucine) were well conserved. The A:T:C:G base composition of the 22 tRNA genes was exhibited in Supplementary Table II, among which A+T content was the highest in O. rebecca (56.6%).

Similar to other Gobiidae fishes (Jin et al., 2012, 2015), small encoding subunit (12S rRNA) and large subunit (16S rRNA) located between the tRNA-Phe and tRNA-Leu (UUR) genes on the H-strand and separated by the tRNA-Val genes (Supplementary Table I; Fig. 1) were identified among three species. The 12S rRNA genes were 950 bp in O. rebecca and 949bp in O. lacepedii and Odontamblyopus sp. However, the 16S rRNA genes were in different length (1693 bp in O. lacepedii, 1695 bp in O. rebecca and 1692 bp in Odontamblyopus sp., respectively). The rRNA genes showed a slightly lower A+T contents comparing with the protein-coding genes and tRNA genes (Supplementary Table II), however, slightly richer than other bony fishes (Kim et al., 2004).

The OL region located between tRNA-Asn and tRNA-Cys comprised 35 bp in length within the WANCY region (Supplementary Table I) and predictably had the capacity to fold into a stable stem-loop hairpin structure consisting of 13 bp in the stem and 11 bp in the loop. While, both O. lacepedii and Odontamblyopus sp. had single base difference with O. rebecca in the loop (Supplementary Fig. S1). In addition, the OL region had a strong asymmetry in the codon usage in the stem with an obvious over-representation of pyrimidines in the 5’ side in the sequence (Supplementary Fig. S1). The conserved motif (5’-GCCGG-3’) at the base of the stem within tRNA-Cys among three species seemed to be in connection with the transition from RNA synthesis to DNA synthesis (Hixson et al., 1986). T-rich or C-rich loop was highly in agreement with previous reports (Zardoya et al., 1995; Cheng et al., 2010). It was suggested that tRNA genes can be functional as origins of replication as they can form OL-like structure without loss of other main functions (Desjardins and Morais, 1990).

Tandemly repeated sequences were a common feature for vertebrate mitochondrial DNA control regions (Bentzen et al., 1998; White and Martin, 2009). Three species (1571 bp, 1336 bp and 1332 bp in O. lacepedii, O. rebecca and Odontamblyopus sp., respectively) were A+T-rich (Supplementary Table II) and identified to be variably long with different repeat number and length variation near the 3’ end of control regions (Supplementary Table I; Table I, Fig. 5A). Two types of repeats were found: perfect repeats and imperfect repeats (Table I). Such variations implied a rapid evolution of the structure in the control region, which may provide information for dissecting the structure/function relationships of the control region.

In this paper, we discussed the phylogenetic relationship based on the 12 protein-coding genes and the analysis demonstrated that clade Amblyopinae included three species (O. lacepedii, O. rebecca and Odontamblyopus sp.). In addition, O. lacepedii has been proved to be a much closer affinity to Odontamblyopus sp. than O. rebecca. The result was in agreement with previous records that Tang et al. (2010) found the cryptic species Odontamblyopus sp. which was similar to O. lacepedii based on collected molecular and morphological data; The two species were sister species and cannot easily be differentiated with naked eye. Agorreta et al. (2013) analyzed sequences data of five molecular makers (two mitochondrial DNA and three nuclear DNA) averaging 222 species of gobioids, while the phylogenetic relationship included only Odontamblyopus sp. within Odontamblyopus and the relationship among three species hasn’t been described with such numberous species studied. Briefly, the phylogenetic relationship among the Godiidae species still remains poorly understood, and our study is confirmed to be a validly supplementary proof and expected to provide significant information for future taxon studies.

The mtDNA genome of three species (O. lacepedii, O. rebecca and Odontamblyopus sp.) were all similar with other bony fishes including genome organization, gene arrangement and codon usage. An unexpected structure in the control region among the three species was discovered that only termination-associated sequence domain (TAS) and conserved sequence blocks domain (CSB-1, CSB-2 and CSB-3) were found. Phylogenetic relationships based on the concatenated 12 protein-coding genes indicated the tree topologies represented two groups, and also demonstrated clade Amblyopinae included all three goby species and that O. lacepedii has been proved to be a much closer affinity to Odontamblyopus sp. than O. rebecca. Three species are distinguished authoritatively by analysis of the complete mitochondrial DNA.

We are grateful to Mr. Xiaozhe Pan for sample collection and Mr. Wei Zhou and Ms. Nan Zhang for offering technical assistance during the paper writing. We also would like to acknowledge Ms. Lu Liu and Mr. Yuman Ju for proofing the manuscript. This study is supported by Public Science and Technology Research Funds Projects of Ocean (201305043; 201405010).

There is supplementary material associated with this article. Access the material online at: http://dx.doi.org/10.17582/journal.pjz/2018.50.6.2173.2183

Authors have declared no conflict of interest.