Samples of A. barbatulus were collected from Poyang Lake in northern Jiangxi Province, China. The experimental material for this study was provided by the Institute of Hydrobiology, Chinese Academy of Sciences. Specimens were stored in absolute ethanol at the College of Fisheries Henan Normal University. In the sampling process, the specific sampling location was not clearly marked, which led to the failure to obtain the latitude and longitude information of the sampling point. Samples from A. barbatulus were extracted using the phenol-chloroform protocol (Sambrook and Russell 2001), dissolved by the addition of 40 μl double-distilled water and stored at −20°C. The 30× genome of A. barbatulus was re-sequenced by Personalbio (Nanjing) using high-throughput sequencing technology to obtain the whole genome. The whole genome sequence of A. barbatulus was compared with the mitochondrial genome sequence of Acheilognathus tonkinensis to determine the star and termination sites of the gene. Then, the mitochondrial genome was extracted from the whole genome data. MitoZ (https://github.com/linzhi2013/MitoZ) was used to obstain the complete mitochondrial genome (GenBank format) (Meng et al. 2019). We constructed libraries with 400 inserts and sequenced them using Next-generation sequencing and paired-end (PE) sequencing, based on Illumina NovaSeq sequencing platform. The bwa (0.7.12-r1039) mem programme was used to compare the filtered high-quality data with the reference genome and the parameters were compared, based on the default parameters of bwamem. Picard 1.107 software (http://www.psc.edu/index.php/user-resources/software/picard) was used to sort and convert the sam files to bam files. We used the “FixMateInformation” command to ensure consistency between all PE reads information. The total number of reads in this specie was 149005815 and the number of reads on the reference genome accounted for 92.01% of the total number of reads.

The complete mitochondrial genome sequence was obtained by comparing it to the published mitochondrial genome of subfamily Acheilognathinae in NCBI. The position of start transfer RNA (tRNA^Phe) was determined using the MITOS website (http://mitos2.bioinf.uni-leipzig.de/index.py) (Bernt et al. 2013). Preprocessing of each gene on the mitochondrial genome was based on the mitochondrial codon of bony fish and by checking whether the mitochondrial genome was ringed. This process yielded a gene map of the A. barbatulus mitochondrial genome. Mitofish (http://mitofish.aori.u-tokyo.ac.jp/annotation/input.html) was used to determine the positions of 13 PCGs, 22 tRNAs, two ribosomal RNA (rRNAs), control regions and the types and numbers of anticodons and start and stop codons. The tRNAscan-SE website (http://lowelab.ucsc.edu/tRNAscan-SE/) (Lowe and Chan 2016) was used to confirm the positions of tRNA genes and predict secondary structures. PCGs were translated exactly using MEGA7.0 (Sudhir et al. 2016). Base composition and relative synonymous codon usage (RSCU) were calculated by MEGA7.0. AT:GC skew was calculated using Perna's formula (Perna and Kocher 1995), with AT Skew = (A - T)/(A + T) and GC Skew = (G - C)/(G + C). CodonW 1.4.2 software (John et al. 1999) was used to obtain the codon adaptation index (CAI), effective codon number (ENC), GC and GC3.

To study the phylogenetic position of A. barbatulus, we selected a part of the fish in the Acheilognathus and Rhodeus data provided by NCBI. Given that subfamilies Gobioninae and Acheilognathae are sister groups (Chen 2014), Pseudorasbora parva was used as the phylogenetic outgroup. A total of 33 complete mitochondrial genome sequences were selected. Table 1 shows the GenBank accession numbers and species names. Partitioned phylogenetic analysis allows the use of different nucleotide substitution models and corresponding parameters for different subsets of association data. This step helps to explore evolutionary models specific to each partition and reduce systematic errors to improve the accuracy of phylogenetic inference (Wang et al. 2012). One partitioned phylogenetic analysis method is based on Bayes principle (Ronquist and Huelsenbeck 2003) and the other is based on the ML partitioned analysis method (Stamatakis 2006). A total of 13 PCGs, two rRNA genes and 22 tRNA genes of each species were extracted using Phylosuite software (Zhang et al. 2019). MEGA7.0 (Sudhir et al. 2016) software was used to individually compare the extracted PCGs by ClustalW (John et al. 1999) (ATP6 = 684 bp, ATP8 = 165 bp, CYTB = 1141 bp, COX1 = 1551 bp, COX2 = 691 bp, COX3 = 785 bp, ND1 = 975 bp, ND2 = 1047 bp, ND3 = 351 bp, ND4 = 1381 bp, ND4L = 297 bp, ND5 = 1836 bp and ND6 = 522 bp).

The sequences were aligned, based on the mitochondrial genome by Phylosiute (Zhang et al. 2019) for the concatenation of PCGs to form a 13 PCG dataset. The best partitioning scheme and the corresponding optimal nucleotide substitution model were determined in accordance with to the Akaike Information Criterion (AIC) criteria using ModelFinder (Kalyaanamoorthy et al. 2017) and a ML tree was established by an edge-linked partitioning model based on IQ-TREE with 5000 ultrafast self-spreading values (Stamatakis 2006). Bayes were constructed, based on Bayesian Information Criterion (BIC) using PartitionFinder (Lanfear et al. 2017) to select the best partition and associated optimal base substitution model (Huelsenbeck 2001, Drummond and Rambaut 2007). Bayes inference phylogenies were inferred using MrBayes 3.2.6 (Ronquist and Huelsenbeck 2003) under an N/A model (two parallel runs, 2000000 generations), in which the initial 25% of sampled trees were discarded as burn-in.

The complete mitochondrial genome (GenBank accession number ON815031) of A. barbatulus is 16726 bp in total length and a typical covalently closed double-stranded cyclic molecule (Fig. 1). The mitochondrial genome of A. barbatulus is broadly similar to that of other vertebrates and includes 13 PCGs, namely, ATP6, ATP8, COX1, COX2, COX3, Cytb, ND1, ND2, ND3, ND4, ND4L, ND5 and ND6, 2 rRNA genes, namely, 12S rRNA and 16S rRNA genes and 22 tRNA genes, for a total of 37 genes. A non-coding region (the control region) that controls gene replication and transcription (Forst and Schulten 2001) and a light strand replication initiation region (OL) were detected. The OL is located between tRNA^Asn and tRNA^Cys and has a length of 31 bp and it can fold to form a secondary structure with a stem-loop structure. This region is highly conserved and related to the replication function of the L-strand. Except for eight tRNAs and one PCG (ND6) located on the light strand (L strand), the remaining 28 genes were located on the heavy strand (H strand) of the mitochondrial genome and the arrangement of the genes was consistent with the typical genetic composition of teleost fish (Fig. 1, Table 2).

Figure 1.  

Gene map of the Acheilognathus barbatulus mitochondrial genomes. The genome contained two rRNA genes (in yellow), 13 coding genes (in black), 22 tRNA genes (in red) and a control region (D-loop) (in brown). The outer ring corresponds to the H- (outermost) and L-strands and depicts the location of PCGs (except for ND6 which is encoded in the L-strand and is portrayed in red). The inner ring (black sliding window) denotes GC content along the genome.

Similar to other bony fish, spacers and overlap between genes were present. Gene overlapping promotes miniaturisation of mitotic genes, shortens genome replication time and offers a natural selection advantage (Boyce et al. 1989). The overlaps of coding genes were located between tRNA^Ile-tRNA^Gln, ND2-tRNA^Trp, OL-tRNA^Cys, ATP8-ATP6, ATP6-COX3, COX3-tRNA^Gly, ND3-tRNA^Arg, ND4L-ND4, ND5-ND6 and tRNA^Thr-tRNA ^Pro and ten gene intervals had a total of 29 bp, which accounted for 0.17% of the total gene length. A total of 15 spacers encoding genes with a total length of 542 bp, accounting for 3.24% of the total gene length, were obtained.

The highest base content in A. barbatulus was that of A (29.33%), followed by those of T (27.6%) > C (26.12%) > G (16.95%). The complete mitochondrial genome of A. barbatulus showed an anti-G bias and AT preference (Table 3), a phenomenon similar to the base composition of most other teleost fishes (Perna and Kocher 1995). This base preference was also present in other parts of the mitochondrial genome, with the D-loop region, tRNAs and PCGs all having considerably greater AT content than CG content and showing a strong AT preference. Of the 13 PCGs, 10 showed an anti-A bias, except for ATP8, COX2 and ND2 which showed an A bias. Except for ND6, which showed an evident G bias, 12 genes showed anti-G bias. Meanwhile, the RNAs exhibited strong A-bias and anti-G-bias. The 12S rRNA and 16S rRNA genes with larger absolute values of AT-skew and smaller absolute values of GC-skew had a strong A-bias. ND6 and tRNAs genes have the same degree of base G-bias and T-bias. From the base composition of mitochondrial genes on both sides of A. barbatulus, PCGs showed a significant anti-G bias, followed by the D-loop region, 12S rRNA and 16S rRNA. Only tRNAs revealed a G bias.

The mitochondrial genome of A. barbatulus contains 13 PCGs with a total length of 11,420 bp, which accounts for 68.28% of the entire mitochondrial genome.

The PCGs of A. barbatulus encode 3798 amino acids and, amongst the encoded amino acids, the Leu, Ser, Pro and Thr exhibit high contents. Low contents of most amino acids, such as Arg, Cys, Glu and Asp, can be observed (Table 4), with the codon encoding Leu being used the most frequently and the one encoding Cys being used the least frequently. The A. barbatulus codons with RSCU values > 1.0 have positive codon usage bias (CUB) and are defined as abundant codons, whereas those with RSCU values < 1.0 have negative CUB and are defined as less-abundant codons (Gun et al. 2018). The statistics of codon usage frequency and RSCU of PCGs (Fig. 2) showed that codons containing bases A and U, such as AAA, AAU and UAA, are used more frequently, whereas those containing bases C and G, such as GGG, GGC and CGC, are used less frequently. According to the comparison of RSCU, for the genes encoding A. barbatulus protein, the most frequently used codons were UCU (1.46%), GCC (1.38%), AAU (1.37%), and AAA (1.36%), which encode amino acids Ser, Ala, Ile and Lys, respectively, whereas codons encoding amino acids Ala (GCG, 0.36%), Ser (UCG 0.46%) and Thr (ACG, 0.52%) had the lowest codon frequency.

Figure 2.  

Results from analysis of Relative Synonymous Codon Usage (RSCU) of the mitochondrial genome of Acheilognathus barbatulus. Codon families are plotted on the x-axis. The label for the 2, 4 or 6 codons that compose each family is shown in the boxes below the x-axis and the colours correspond to those in the stacked columns. RSCU values are shown on the y-axis.

Most of the start codons of A. barbatulus are ATG. Only the COX1 gene has GTG as the start codon. ND1, COX1, ATP6, ND4L, ND5 and ND6 use TAA as the complete stop codon. ND2, ATP8 and ND3 use TAG as the complete stop codon, whereas COX2, COX3, ND4 and CYTB use TA- and T-- as incomplete stop codons.

Codon usage bias (CUB) refers to the unequal use of synonymous codons in organisms. CUB is affected by mutation and selection pressure and is an important feature of biological evolution. It not only affects gene function and expression potential, but also the accuracy and efficiency of translation. The higher the gene expression level, the stronger the CUB. To investigate codon usage preference in the mitochondrial genomes of A. barbatulus, we calculated CAI, ENC, GC and GC3 using CodonW1.4.2.

Effective number of condon (ENC) can describe the extent to which codon usage deviates from random selection, with values generally varying from 20 to 61. The larger the ENC value, the lower bias of expression genes towards the use of rare codons. The smaller the ENC value, the greater the preference for codons of genes with high expressions. The ENC values of PCGs in the mitochondrial genome of A. barbatulus ranged from 39 to 55.71 (Table 5) and were mostly concentrated in the range of 43–46, indicating a degree of codon usage preference. Codon adaptation index (CAI) values ranged from 0 to 1, with high values indicating high gene expression levels and pronounced CUB (Sharp and Li 1987). ATP8 and the three cytochrome oxidase subunits (COX1, COX2 and COX3) showed high CUB and gene expression levels. ATP6 exhibited the lowest CUB and expression levels amongst the PCGs (Table 5). GC3 refers to the GC content of the third position of all codons in a gene. In addition to methionine, tryptophan and termination codons, G and C may appear in the third codon position. The results showed that the 13 PCGs had a low GC content and a low probability of G and C in the third codon position. This finding also demonstrated that the A. barbatulus PCGs showed a strong AT preference (Table 3).

The animal mitochondrial genome has two types of ribosomal units: the large 16S and small 12S subunits. The 16S subunit is more conserved than the 12S and the secondary structure of both rRNA genes is more conserved than the sequence (Noack et al. 1996). The 12S rRNA is located between tRNA^Phe and tRNA^Val, with a total length of 956 bp, accounting for 5.72% of the complete mitochondrial genome and 50.57% of the AT content. The 16S rRNA is located between tRNA^Val and tRNA^Leu, with a total length of 1660 bp, accounting for 9.92% of the total mitochondrial genome and 56.15% of the AT content.

A. barbatulus has 22 tRNA genes (eight in the L-strand and 14 in the H-strand). Their individual gene lengths ranged from 68 bp to 76 bp. Except for Leu and Ser, which both contain two types, the other 18 tRNAs have only one type. Except for tRNA^{Ser (GCT)} (L strand), which lacks the D-arm resulting in an incomplete secondary structure, all 21 tRNAs can fold into the canonical cloverleaf secondary structure (Fig. 3), which is a condition that has been reported for the mitochondria of other fish species (Broughton et al. 2001, Hwang et al. 2013). The cloverleaf structure consists of four domains (AA stem and D, AC and T arms) and a variable loop (Fig. 3). For their normal functioning, these aberrant tRNAs may require co-evolutionary interaction factors or post-transcriptional RNA editing (Masta et al. 2004).

Figure 3.  

Secondary structure of the 22 tRNA genes of the mitochondrial genome of Acheilognathus barbatulus predicted by tRNAScan-SE 2.0.

In this study, the ML tree and MrBayes were constructed by tandemly linking 13 PCGs and using Pseudorasbora parva as an outgroup partition to accurately reveal the phylogenetic relationships of A. barbatulus. The best partitioning ML and Bayes models were based on AIC and BIC for 34 PCG tandem sequences and the most suitable nucleotide-substitution models, respectively. The best partitioning models for ML were GTR+F+I+G4: ND1, ND3, ND4, ATP6; GTR+F+I+G4: ND2; GTR+F+I+G4: ND5; GTR+F+I+G4: ND6; HKY+F+I+G4: ATP8; GTR+F+I+G4: COI; GTR+F+I+G4: COII; GTR+F+I+G4: COIII, Cytb, ND4L. The optimal partitioning models of ML and Bayesian were slightly different, with a slight difference in node support. The topological structure was consistent overall, but parts of the support rate or posterior probability of the nodes is inconsistent. This finding may be due to different ML and Bayes algorithms being different, resulting in differences in the subsequent phylogenetic trees. The results support the conclusion that Acheilognathus and Rhodeus constitute monophyletic taxa (Yang 2010) A. barbatulus is most closely related to Acheilognathus tonkinensis and Acheilognathus cf. macropterus (Chen 2011). The node support rate of ML and MrBayes is 100/100 (Figs 4, 5, respectively).

Figure 4.  

Phylogenetic trees derived from the Bayes approaches, based on concatenation of PCGs. The numbers on the nodes are the bootstrap values of Bayes. The number after the species name is the GenBank accession number. Outgroup taxa are shown, the text is bolded by the Acheilognathus barbatulus.

Figure 5.  

Phylogenetic trees derived from the Maximum-Likelihood (ML) approaches, based on concatenation of PCGs. The numbers on the nodes are the bootstrap values of ML. The number after the species name is the GenBank accession number. Outgroup taxa are shown, the text is bolded by the Acheilognathus barbatulus.