Complete Mitochondrial Genomes of Polydactylus macrochir (Günther, 1867): Genome Characterization and Phylogenetic Analysis

1. Introduction

Threadfins are easily identified due to their numerous, elongated and threadlike pectoral-fin rays [1]. These filamentous fin rays are necessary for threadfins inhabiting sandy or muddy bottoms of turbid shallow waters to detect food, with both tactile and gustatory functions [2]. Threadfins belong to the family Polynemidae which includes eight genera and 42 extant species [3].

For the past years, the placement of the Polynemidae among teleosts has varied according to different studies. Polynemidae had been previously considered related to Mugilidae, Sphyraenidae, Atherinidae, and Phallostethoidei, and clustered into Mugiliformes [4]. The later study insisted that Polynemidae, Mugilidae and Sphyraenidae should remain within Perciformes [5]. Other studies based on morphological features suggested that Sciaenidae and Polynemidae were sister-groups [6,7,8]. Phylogenetic placement of Polynemidae inferred from molecular analyses was also conflicted with each other. Different studies have alternatively aligned polynemids with Menidae [9], Menidae + Lactariidae [10], or Pleuronectiformes [11,12]. A study grouped Polynemidae with Menidae as the sister group of Sphyraenidae [13]. The rencent study thought that polynemids should not be included in the Order Mugiliformes or Order Perciformes and reset the Order Polynemiformes [14]. As mentioned above, different studies have advanced alternative hypotheses of relationships for the family Polynemidae and a consensual phylogenetic position is still lacking.

The mitochondrial genome, typically spanning 16-17 k bp, generally consists of 13 protein-coding, two ribosomal RNAs (rRNAs), 22 transfer RNAs (tRNAs), and two noncoding regions-the control region (CR) and the origin of L-strand replication (OL) [15]. Notably, the mitochondrial DNA evolves at a faster rate and exhibits a remarkable degree of conservation in gene organization and low recombination level [16]. In addition, the complete mitogenome has an advantage in providing comprehensive information when used for evolutionary analysis compared with partial mitogenome gene. so, the complete mitogenome has been used widely to study phylogenetic relationships and species classification among teleosts [17,18,19].

To date, only seven species within the family Polynemidae have its complete mitogenomes been sequenced worldwide. Polydactylus macrochir (Günther, 1867) is a number of family Polynemidae. In this study, we sequenced the mitochondrial genome of P. macrochir and conducted genome characterization and phylogenetic analysis. This study will provide new insights regarding classiffcation of Threadfins and would be helpful to resolve uncertainties of phylogenetic position of Polynemidae.

2. Materials and Methods

2.1. Samples, DNA Extraction and Sequencing

The samples of Polydactylus macrochir was purchased from Fangchenggang local maket, Guangxi, China, sent to Genepioneer Biotechnologies (Nanjing, China) for sequencing. Total genomic DNA was extracted using the MagPure Tissue DNA LQ Kit (Magen Biotechnology, Guangzhou, China). The DNA purity was detected with 1.0% agarose gel. After that the qualified library was sequenced using the Illumina Novaseq platform and the paired end sequencing (PE) read length was 150 bp.

2.2. Genome Assembly and Annotation

The fastp (v0.20.0, https:// github. com/ OpenGene/ fastp) software was used to filter the original data and obtained the clean data. To reduce the complexity of subsequent sequence assembly, the software bowtie2 v2.2.4 (http:// bowtie-bio. sourceforge. net/ bowtie2/ index. shtml) was used to compare the mitochondrial genome database built by the company and the compared sequencing sequences were regarded as the mitochondrial genome sequencing sequences of project samples. Then the mt DNA sequence is assembled by SPAdes soft-ware to obtain seed sequence of mitochondrial genome [20]. Kmer iterative extend seed, if the result is a contig, the result is determined as pseudo genome sequence. The sequenced sequence is aligned to pseudo genome for genome correction. Otherwise, SSPACE v2.0 (https://www. baseclear. com/ services/bioinformatics/ basetools/ sspace-standard/) software is used to connect the previously obtained contig sequences to obtain scaffolds; Gapfiller v2.1.1 (https:// sourceforge. net/ projects/ gapfiller/) software was used to fill in the scaffolds sequence. If there was still gap after the above operation, primers were designed for PCR sequencing and assembly until a complete pseudo genome sequence was obtained, and then we aligned the sequenced sequence to the pseudo genome to determine the accuracy of the final assembly results for genome correction. Finally, according to the structure of mitochondrial, the corrected pseudo genome was rearranged to obtain a complete mitochondrial circular genome sequence. The annotations of mitochondrial genome were performed in the Mitos2 web server [21] (Parameters: E-value Exponent=5, Maximum Overlap=100, ncRNA overlap=100). The Mitos2 annotation results were compared with closely related species, and the final annotation results were obtained after manual correction. Mitochondrial gene structure maps were drawn using the OGDRAW(https:// chlor obox. mpimp-golm. mpg. de/ OGDraw. html).

2.3. Codon Usage Analysis

Relative Synonymous Codon Usage (RSCU) is thought to result from a combination of natural selection, mutation and genetic drift, the numerical value is the ratio of the actual frequency of codon usage to the theoretical frequency of codon usage.A script wrote in Perl was used to filterUniq CDS(choose one of multiple copies of CDS) and do the calculations.

2.4. Ka/Ks Value Analyses

To understand the natural selection pressure in the evolution of the family Polynemidae, homologous protein sequences between P. macrochir and other threadfins were obtained using BLASTN. Then, the shared protein-coding genes were aligned using MAFFT v7.427 [22]. The non-synonymous (Ka)and synonymous (Ks) ratios (Ka/Ks) were calculated using KaKs_Calculator v2.0 [23].

2.5. Phylogenetic Analysis TimeTree Estimation

The mitochondrial genomes of 22 species were downloaded from GenBank (Table 1). The shared CDS was aligned using MAFFT [5] procedure. Maximum likelihood (ML): Maximum likelihood (ML) phylogenetic tree was conducted by RAxML v8.2.10 (https://cme.h-its.org/exelixis/software.html)(GTRGAMMA model) estimation with 1000 bootstrap replications. Bayesian inference (BI): The optimal nucleotide substitution model was calculated by jModelTestv2.1.10 (https://github.com/ddarriba/jmodeltest2), and then MrBayes v3.2.7a (http://nbisweden.github.io/MrBayes/) was used to establish Bayesian inference (BI) phylogenetic tree, the parameters of MrBayes v3.2.7 software are based on jModelTest v2.1.10 results. After inputting the sequence data, the constructed Bayesian inference topology (.nwk format) was used as a baseline tree. Fossil records acquired from the TimeTree website (http://www.timetree.org) were used to calibrate the divergence times. Estimate divergence times using PAML mcmctree (http://abacus.gene.ucl.ac.uk/software/paml.html)

3. Results

3.1. Mitogenomic Structure and Organization

In the present study, the mitogenome of P. macrochir (GenBank Accession No. PQ675801) was determined, with a total length of 16,738 bp. The complete length of P. macrochir mitogenome is similar to that of other sequenced Polynemidae species. It comprised 13 protein-coding genes (PCGs), 22 transfer RNA genes (tRNAs), 2 ribosomal RNA genes (rRNAs), and 1 noncoding region (D-loop) (Figure 1 ). NAD6 and 8 tRNAs (trnQ, trnA, trnN, trnC, trnY, trnS2, trnE, and trnP) were located on the light (L) strand, while the remaining genes including 12 PCGs, 2 rRNAs, and 14 tRNAs were encoded by the heavy (H) strand. The P. macrochir mitogenome had six overlapping regions, with the longest overlapping region (10) detected between atp6 and cox3 genes (Table 2). Nucleotide composition analysis showed that P. macrochir mitogenome was AT biased(54.15%). PCGs, tRNAs, and rRNAs exhibited an AT content similar to that of the total mitogenomes. CR had the highest A+T content. All the AT skews in the mitogenome, tRNAs, rRNAs, and Dloop were positive, while AT skews of PCGs was negative. The GC skew of mitogenome, PCGs, rRNAs, and Dloop was negative, while GC skew of tRNAs was positive. In addition, the GC skew for the complete mitogenomes and PCGs were strongly negative, with a value of -0.286, -0.329, respectively (Table 3).

3.2. Protein-Coding Genes and Codon Usage

The total length of P. macrochir PCGs was 11,447 bp, accounting for 68.39% of the whole mitogenome. The the largest PCG was nad5, with a length of 1839 bp, the shortest PCG was atp8, with a length of 1839 bp. All the PCGs started with the regular codon ATG, except for the cox1 which started with GTG. The termination codon varied across PCGs. nad6 terminated with the TAG stop codon. nad2 terminated with the incomplete stop codon TA. The cox2, nad3, and nad4 used an incomplete T stop codon. The other PCGs shared the same complete stop codon TAA (Table 2).

The relative synonymous codon usage (RSCU) values of PCGs are revealed in Figure 2. There were 30 codons with RSCU > 1. The number of codons ending with the A base was 14, with 1 codon ending with U. Similarly, there were 14 codons ending with C and 1codon ending with G. The results also showed that serine and leucine were encoded by six codons with greater codon abundance compared to other amino acids.

(The y-axis represents the usage frequency of the corresponding amino acid codons in 13 PCGs. Different colors represent the different codons in the amino acid)

3.3. Ka/Ks Value Analyses

To evaluate the selective pressures of PCGs, we detected Nonsynonymous (Ka) and synonymous (Ks) substitution rates. It was shown that the ka/ks ratios of all PCGs were less than 1. The Ka/Ks ratio was the lowest in cox1 (Figure 3).

3.4. Phylogenetic Analysis and TimeTree Estimation

To find the position of P. macrochir among 22 other species, the phylogenetic tree was constructed using maximum likelihood (ML) and Bayesian inference (BI) based on coding sequence. It was shown that ML trees and BI trees shared a similar topological structure (. As an outgroup, D. rerio was distinct from the other fish species. The remaining fish species was grouped into two major clades. Eight polynemid species formed a well-supported monophyletic cluster. P. macrochir and the monophyletic cluster consists of Eleutheronema tetradactylum Polydactylus plebeius formed a sister group. Polynemidae formed a sister cluster with the clade of Sphyraenidae with high bootstrap support value (100%). Then, (Polynemidae + Sphyraenidae) and (Menidae + Xiphiidae + Carangidae) were sister lineage. The other representative species from different family display distinct monophyletic clustering patterns. According to TimeTree analyses, the most recent common ancestor (MRCA) of Polynemidae was traced back to about 52.81 million years ago(MYA), while 6 species within Polynemidae diverged from 11.70 MYA to 20.05 MYA.

Figure 4. ML phylogenetic trees based on the nucleotide datasets for 13 PCGs from the mitogenomes of 23 species. The posterior probabilities were overlaid with each node.

Figure 4. ML phylogenetic trees based on the nucleotide datasets for 13 PCGs from the mitogenomes of 23 species. The posterior probabilities were overlaid with each node.

Figure 5. BI phylogenetic trees based on the nucleotide datasets for 13 PCGs from the mitogenomes of 23 species. The posterior probabilities were overlaid with each node.

Figure 5. BI phylogenetic trees based on the nucleotide datasets for 13 PCGs from the mitogenomes of 23 species. The posterior probabilities were overlaid with each node.

Figure 6. Divergence time estimation of the threadfins. Number at nodes indicate estimated age.

Figure 6. Divergence time estimation of the threadfins. Number at nodes indicate estimated age.

4. Discussion

The mitochondrial genome of P. macrochir has a total length of 16,738 bp, containing 13 PCGs, 22 tRNAs, 2 rRNAs, and 1 noncoding CR, with a typical gene content of vertebrates [24]. In terms of gene arrangement, P. macrochir is consistent with other sequenced related species [14]. Both of these results indicated that P. macrochir was conserved. The nucleotide composition of the P. macrochir mitogenome was rich in A + T (54.15%), with slightly positive AT-skew (0.049), while GC-skew was moderately negative (-0.286). The overlapping and intergenic region, which has commonly been observed in other species, were both present in the mitogenome of P. macrochir. The length of overlap or intergenic region between adjacent genes varied with special lineage. All the PCGs started with the typical ATN codons, except for cox1, initiated by the GTG. In P. macrochir, TAA was the most frequent codon used as stop codon, while TAG terminating codon was used once, with the remaining incomplete T and TA codon. The atypical initiation codons and incomplete terminating codons were also reported in mitogenome of other fish species [25].

As we know, the RSCU usually reflects the preference for codon usage [26]. In P. macrochir, the PCGs were more prone to use A and C than T and G at the third codon, different from Epinephelidae [27]. The bias in the usage of codons may provide hits to understand the evolutionary history of P. macrochir.

The analysis of nonsynonymous (Ka) and synonymous (Ks) substitution rates within miogenome PCGs provides evidence for studying natural selection and local adaption [28]. In this study, the Ka/Ks ratio of all PCGs was much less than 1, indicating a strong negative selection among the selected species. Mitochondrial PCGs play a crucial role in oxygen utilization and energy metabolism, which was vital for the survival and development of organism. The low Ka/Ks ratio of all PCGs ensure their functional conservation, reflecting the evolutionary strategy for ecological niche adaptation.

Previous studies did not resolve the uncertainties of Polynemidae phylogenetic position. In this study, we reconstructed a phylogenetic tree for P. macrochir and 22 other species based on 13 PCGs from the mitochondrial genome. Same topologies were generated by both BI and ML analysis.P. macrochir another 7 species of Polynemidae were clustered into a clade with high nodal support values, validating the monophyly of the family Polynemidae. Within Polynemidae, P. macrochir and the other Polydactylus species cannot form a monophyletic group. The ability to fully understand the phylogeny of Polynemidae was restricted by limited mitogenome data. Because, there are 8 genera including 42 extant species, only 8 mitogenomes were determined. Therefore, more mitogenomes need to be sequenced to reconstruct Polynemidae phylogenetic tree. It was shown that Polynemidae and Sphyraenidae were sister groups. Furthermore, sister-group relationships between (Polynemidae + Sphyraenidae) and (Menidae + Xiphiidae + Carangidae) were robustly supported. These results conflicted with earlier studies which proposed that polynemids were closer to Mugilidae, Sphyraenidae, Atherinidae, and Phallostethoidei, Mugilidae and Sphyraenidae or Sciaenidae according to different phenotypic analyses [4-8]. Alternative hypotheses for the phylogenetic placement of Polynemidae based on molecular analyses, such as, Polynemidae + Menidae as the sister group of Sphyraenidae [13], Polynemidae were grouped with Menidae and Lactariidae [10], Polynemidae and Menidae was sisiter-lineage [9], were also in conflict with present study. Rencently, sister-group relationship between Polynemidae and Pleuronectiformes was supported by molecular and morphological data [11,12]. Our study was not invloved in the relationship between Polynemidae and Pleuronectiformes, which should be resolved by future research. According to the present results, we tend to locate Polynemidae in Perciformes and do not support independent Order Polynemiformes [14]. Polynemidae and Sphyraenidae clustered together as a monophyletic group, might form a new suborder of Perciformes needing more evidences. The estimation of divergence time displayed that Polynemidae and Sphyraenidae diverged from the most recent common ancestor in the late Cretaceous at 77.80 MYA. Threadfins may originate in the Paleogene at 52.81 MYA. Some species of Polynemida may undergo an accelerated diversification during the Neogene period. This enhances our understanding of the phylogenetic relationship and speciation of Polynemidae.

5. Conclusions

In this study, the complete mitogenome of P. macrochir was sequenced and analyzed. The mitochondrial genome had a total size of 16,738 bp and had a conserved gene contents and arrangement. The PCGs tended to use A and C than T and G at the third codon. All of the PCGs were under purifying selection. Phylogenetic analysis revealed that Polynemidae were clustered into a monophyletic group being sister group to Sphyraenidae. Threadfins may originate in the Paleogene at 52.81 MYA. Some species of Polynemida may undergo an accelerated diversification during the Neogene period. These results provide valuable information for future studies on the phylogenetic relationships and speciation of Polynemidae.