Comparative Analysis of Complete Mitochondrial Genome of Ariosoma meeki (Jordan and Snider, 1900), Revealing Gene Rearrangement and the Phylogenetic Relationships of Anguilliformes

Simple Summary In this study, we report the complete mitochondrial genome of Ariosoma meeki (Anguilliformes (Congridae)), the mitochondrial genome structure and composition were analyzed. We found the mitogenome of A. meeki has undergone gene rearrangement: The ND6 and the conjoint tRNA-Glu genes were translocated to the location between the tRNA-Thr and tRNA-Pro genes, and a duplicated D-loop region was translocated to move upstream of the ND6 gene. At the same time, we speculated the possible evolutionary process of gene rearrangement, and made the hypothesis of tandem repeat and random loss for explanation. The results of phylogeny also echo this inference in Anguilliformes. Abstract The mitochondrial genome structure of a teleostean group is generally considered to be conservative. However, two types of gene arrangements have been identified in the mitogenomes of Anguilliformes. In this study, we report the complete mitochondrial genome of Ariosoma meeki (Anguilliformes (Congridae)). For this research, first, the mitochondrial genome structure and composition were analyzed. As opposed to the typical gene arrangement pattern in other Anguilliformes species, the mitogenome of A. meeki has undergone gene rearrangement. The ND6 and the conjoint tRNA-Glu genes were translocated to the location between the tRNA-Thr and tRNA-Pro genes, and a duplicated D-loop region was translocated to move upstream of the ND6 gene. Second, comparative genomic analysis was carried out between the mitogenomes of A. meeki and Ariosoma shiroanago. The gene arrangement between them was found to be highly consistent, against the published A. meeki mitogenomes. Third, we reproduced the possible evolutionary process of gene rearrangement in Ariosoma mitogenomes and attributed such an occurrence to tandem repeat and random loss events. Fourth, a phylogenetic analysis of Anguilliformes was conducted, and the clustering results supported the non-monophyly hypothesis regarding the Congridae. This study is expected to provide a new perspective on the A. meeki mitogenome and lay the foundation for the further exploration of gene rearrangement mechanisms.


Introduction
The gene rearrangement patterns among vertebrate mitogenomes were initially considered to be conservative. The convincing nature of this concept was reinforced after the arrangement patterns (traditional gene arrangement and special gene rearrangement) have been reported; the latter occurs in a small percentage of the Anguilliformes population [8,9].
The conjectures regarding gene rearrangement in mitochondrial genomes that have been proposed so far can be classified into three main hypotheses [10]. In the first scenario, Poulton et al. [11] initially analyzed human mitochondrial rearrangements and found that mitochondrial genomes were involved in DNA strand-breaking and reconnection, before homologous recombination. This model of gene rearrangement has been shown in subsequent studies involving mussels, birds, and frogs, being used as a practical tool to assist interpretation and analysis [12,13]. The second model was proposed by Arndt et al. [14] and has been widely accepted; it features the interpretation of tandem replication and the random loss of mitochondrial genomes (TDRL model). In this model, the mitochondrial rearrangement was considered to be the result of structural variation caused by the random deletion of duplicates after tandem replication in some genes [15]. This hypothesis has been adopted in many gene rearrangement studies [16,17]. The third model was proposed by Lavrov et al. [18] in the mitochondrial genome report of millipedes. The model indicates that random loss is activated after full replication by mitochondrial genes, and the loss ultimately depends on the polarity and position of the genes. Unlike the TDRL model, this one emphasizes tandem repetition and nonrandom loss (TDNL model). Given the above pioneering and characteristic hypotheses, there is no conclusive interpretation of how gene rearrangements occur and what their implications are. Slowly but surely, details about facilitation in gene rearrangement are starting to trickle in with the enrichment of the genomic databases.
Comparative genomics analysis is widely regarded as an efficient approach for phylogenetic analysis, owing to the richness of molecular information that is characterized by maternal inheritance, a high rate of evolution, and a relatively low rate of intermolecular recombination [19]. Generally, most eels are snake-like in size, featuring narrow gill pores and no pelvic fins [20]. For these reasons, much fuzziness is involved in their identification and the corresponding phylogenetic research. In this context, the advantages of molecular approaches become prominent when it comes to species with insufficient morphological data. Nevertheless, genetic information for most Anguilliformes species is insufficiently disclosed, and the location of each species within the Anguilliformes phylogeny remains elusive [8], resulting in poor understanding and controversies.
In this paper, we report a new version of the complete mitochondrial genome sequence of A. meeki; the genetic composition and arrangement of A. meeki have been described in detail and a comparative genomic analysis has been conducted for A. meeki (this study) and Ariosoma shiroanago. Interestingly, gene rearrangement was detected in both Ariosoma mitogenomes, differing from the genetic features of two published A. meeki mitogenomes (KX641476 and MN616974), which have been marked as "Unverified" in the NCBI database. Therefore, we have performed the correct sequence test and publication. At the same time, for this study, we have carried out a systematic analysis of the species evolution of Anguilliformes.

Sample Collection, DNA Extraction, and PCR Amplification and Sequencing
Individual A. meeki samples were collected using traditional trawling methods that have been approved in Zhoushan City, Zhejiang Province, China (30 • 40 30 N, 121 • 20 28 E). After the samples were collected, we preserved them with 95% ethanol and kept them permanently in a freezer at −80 • C in the museum of the National Marine Facility's aquaculture engineering technology research center. Total DNA was extracted from the tissue of 6 individual specimens' muscles, using an Aidlab Genomic DNA Extraction Kit (Beijing, China), following the manufacturer's instructions. The complete mitochondrial genome of A. meeki was sequenced by Sangon Biotech (Shanghai, China), and the design primers were based on the published mitochondrial genome and then amplified in strict accordance with the requirements of the kit (Takara, China) [21].

Sequence Analysis and Assembly and Mitochondrial Genome Annotation
The complete mitogenome of the invertebrate genetic sequence was annotated using the MITOS web server (http://mitos2.bioinf.uni-leipzig.de/index.py, (accessed on 15 July 2022) [22,23] and was then manually corrected. After correcting the sequenced DNA fragment results, the CodonCode Aligner 5.1.5 (CodonCode Corporation, Dedham, MA, USA) was applied to splice the fragments and complete the creation of the mitotic genome. The Sequin software (version 15.10) was used to mark and annotate the complete mitochondrion genome sequence, to determine the location and boundary of protein-coding and the ribosomal RNA gene; the results were verified via NCBI-BLAST (http://blast.ncbi.nlm.nih.gov, accessed on 4 July 2022) comparison. MitoFish (version 3.69, http://mitofish.aori.u-tokyo.ac.jp/annotation/input.html, (accessed on 4 July 2022) was used to predict the sequence characteristics of the A. meeki mitochondrial annular genome and to map the cyclical gene. The gene rearrangements were defined using the CREx program (http://pacosy.informatik.uni-leipzig.de/crex, (accessed on 28 July 2022) [24]. The relative synonymous codon usage (RSCU) values were analyzed with MEGA 11.09 [25]. Composition skew values were calculated according to the following formulas: AT-skew = (A − T)/(A + T); GC-skew = (G − C)/(G + C) [26].

Phylogenetic Analyses
Thirty-one complete Anguilliformes mitochondrial genomes were downloaded from GenBank (https://www.ncbi.nlm.nih.gov/genbank/, (accessed on 28 July 2022) for phylogenetic studies (Table 1). Saccopharyngiformes have been thought to be closely related to Anguilliformes [27]. Therefore, two Saccopharyngiformes species, Eurypharynx pelecanoides and Saccopharynx lavenbergi, were selected as the outgroup. In the analysis, 12 PCGs sequences were selected in order to construct a phylogenetic tree using DAMBE, version 7.2.3 [28]. These PCGs did not include ND6 since the base composition was more irregular than other sequences and led to poor phylogenetic performance [29]. Sequences were aligned with default parameters, using Clustal X 2.0 [30], and were manually checked using BioEdit [31]. Ambiguous sequences were eliminated using Gblock [32]. Substitution vs. the Tamura-Nei (TN93) genetic distance in pairwise comparisons was used to test for substitution saturation, using DAMBE (version5.3.19) [33]. The third codon positions showed significant saturation, which, as a result, were defined only as purines and pyrimidines (3RY) [34]. The phylogenetic analyses were conducted utilizing the MrBayes 3.2.6 and PhyML80 software, based on Bayesian inference (BI) and maximum likelihood (ML), respectively [35,36]. The best-fit models of nucleotide substitution for each of the sequences were selected, using MrModelTest 2.2, from 33 models [37]. ML analysis uses bootstrap analysis (1000 repetitions) to verify the relative support levels [38]. Bayesian analysis was performed using default settings over four independent sets; the average standard deviation of split frequencies was < 0.01, the estimated sample size was >200, the potential scale reduction factor approached 1.0, and all parameters were checked with Tracer v. 1.6 [39]. The resulting phylogenetic trees were visualized using FigTree v. 1.4.4 and its tool (https://itol.embl.de, (accessed on 29 July 2022) [40]. Note: * the length excludes the control region.

Genome Structure and Composition
Compared to the published mitochondrial genome of bony fish, the complete mitogenome of A. meeki (17,695 bp) is within the normative range (16,099-18,247 bp) ( Table 1). In all, 37 typical coding regions (13 PCGs, 22 tRNAs, 2 rRNAs) and 1 origin of replication between the WANCY structures (tRNA-Trp, tRNA-Ala, tRNA-Asn, tRNA-Cys, and tRNA-Tyr) were found in the mitogenome of A. meeki. However, two control regions (D-loops) were detected next to the tRNA-Thr and tRNA-Pro genes, respectively; the features of duplicated D-loops were distinct from most vertebrate mitogenomes [1]. The gene order of the A. meeki mitogenome was also distinctive. It indicated that the ND6 and tRNA-Glu genes were translocated upstream of tRNA-Pro and downstream of tRNA-Thr, respectively. The additional D-loop was translocated to the site ahead of ND6 as well; such structural heterogeneity is regarded as an uncommon trait among vertebrate mitogenomes [6,44].
The mitogenomic compositions and corresponding gene order for most bony fish are considered to be relatively conserved. Previous research revealed that the mitochondrial gene structure of homologous species could be analogous. For instance, many Anguilliformes species, such as the Anguilla, Simenchelys, and Synaphobranchus populations, revealed typical mitogenomic compositions (37 coding regions and two non-coding regions), and the gene arrangement among them showed a high degree of consistency, while it would be another story when it comes to other Anguilliformes species, such as Facciolella, Ariosoma, and Muraenesox, of which the species in the Ariosoma genus revealed a distinct gene rearrangement pattern [8]. In the mitogenome of A. meeki, an additional control region was detected, while one coding protein, ND6, underwent translocation ( Figure 1). Compared with other non-rearranged species, such a genetic phenomenon between Cytb and ND6 can be found across the published Ariosoma mitogenome (GenBank accession number: AP010861). Nevertheless, there is only one control region in the reported mitogenome of A. shiroanago. In this study, we performed a comparison and analysis of other genetic locations between A. meeki and A. shiroanago (Tables 2 and 3) (Figures 2-4). The overall base composition was 28.93% (A), 25.39% (C), 26.14% (T), and 19.54% (G) for A. meeki, and the overall base composition was 32.54% (A), 23.90% (C), 26.85% (T), and 16.71% (G) for A. shiroanago (Tables 2 and 3). Typically, most of the molecular components in mitochondrial genomes were closely linked with each other; 11 intergenic spacers consisting of 107 base pairs were observed in the A. meeki mitogenome. The largest region contained 49 pairs of nucleotides and was located between ND5 and Cytb; the smallest area was formed by a single 1 bp between Cytb and tRNA-Thr. As for the A. shiroanago mitogenome, the corresponding length thresholds for non-coding regions were 154 bp and 41 bp, which were situated between ND6 and tRNA-Thr, Cytb, and ND5, respectively.
The initiation codons for most PCGs were ATG, except for the COX1 gene, which was initiated with GTG ( Table 2). The utilization of such an extraordinary initiation codon for the COX1 gene is observable in most teleostean mitogenomes [47]. In the A. meeki mitogenome, five PCGs (ATP8, ATP6, COX3, ND4L, and Cybt) ended by TAA, while 5 PCGs (ND1, ND2, ND3, ND5, and ND6) were terminated by TAG, COX2 and ND4 terminated with a single T, and COX1 was stopped by AGA. The phenomenon of the coding proteins ended by incomplete codons is widespread in invertebrate as well as vertebrate mitogenomes [48,49]. The probable mechanism adopted to explain this occurrence is that the codon TAA has generated the subsequent polyadenylation process through transcription [50]. Three PCGs in the A. shiroanago mitogenome owned different types of stop codons, which showed ND5 and ND6 terminated with TAA and COX1 terminated with AGG.
The modes for the codon usage of 13 PCGs between the A. meeki and A. shiroanago mitogenomes are shown in Figure 2. The amino acids primarily used in PCGs for the previous species could be Leu1, Ser2, Pro, and Thr, with Leu1, Leu2, Ser1, and Ser2 for the latter. The relative synonymous codon usage (RSCU) analysis in two Ariosoma species, in terms of the third position, is depicted in Figure 2. The utilization among two-and four-fold degenerate codons presented an overall bias toward those codons that are abundant in A.

Transfer RNAs, Ribosomal RNAs, and D-Loops
Twenty-two tRNAs in A. meeki and A. shiroanago were detected, based on their unique codons. Fourteen tRNAs were located on the heavy chain, and the rest were tRNAs (tRNA-Gln, tRNA-Ala, tRNA-Asn, tRNA-Cys, tRNA-Tyr, tRNA-Ser, tRNA-Glu, and tRNA-Pro), which were found on the light chain ( Figure 1). A typical clover structure was detected among most tRNAs in both species, with the exception of tRNA-Ser1 (Figures 3 and 4), which was unable to form a stable clover structure due to the lack of a complete dihydrouridine arm [51]. Twenty-two tRNAs contained 1560 base pairs in the A. meeki mitogenome, with 1562 in the A. shiroanago mitogenome. The length of each tRNA gene ranged from 64 to 76 bp in the mitogenome of A. meeki and from 66 to 76 bp in the mitogenome of A. shiroanago. The base composition was 29.03% (A), 27.12% (T), 20.45% (C), and 23.40% (G) for the former, and 29.71% (A), 28.17% (T), 20.23% (C), and 21.90% (G) for the latter. The anticodons of the mitogenomes in two Ariosoma species reflected the same usage pattern ( Table 2).
The genes 12S and 16S were found to be surrounded by tRNA-Phe and tRNA-Leu1 and were separated by tRNA-Val in the A. meeki and A. shiroanago mitogenomes ( Table 3). The total length of rRNAs was 2652 bp and 2687 bp in both Ariosoma mitogenomes, respectively. The base compositions of the rRNAs were 33.78% (A), 20.06% (T), 22.25% (C), and 23.91% (G) for the mitogenome of A. meeki, and 36.55% (A), 20.54% (T), 22.29% (C), and 20.62% (G) for the mitogenome of A. meek. The AT-skew and GC-skew values were 0.04 and 0.26 for the former and 0.28 and −0.04 for the latter. Such a base bias revealed the higher percentages of the adenine and cytosine nucleotides of two rRNA genes in both mitogenomes.
In the A. meeki mitogenome, two D-loops were detected downstream of tRNA-Thr and tRNA-Pro, respectively. The total length of the D-loops was 1929 bp (969 bp for D-loop1 and 960 bp for D-loop2), with an identical AT bias of 62.00% (Tables 2 and 3). Compared to other components, both D-loops presented a higher percentage in terms of AT bias. For this reason, the D-loop region, as the main non-coding area, has been called the "AT-rich region". The AT-and GC-skew values were 0.17 and −0.17 in both D-loops, reflecting the abundance of adenine and cytosine and the thymine deficiency and guanine. The termination of heavy chain replication could be located in two D-loops, based on the palindrome element motifs, such as "TACAT" and "ATGTA" [52] (Figure 5). based on the palindrome element motifs, such as "TACAT" and "ATGTA" [52] ( Figure  5).
The exploration of a substitution saturation index was implemented, utilizing the aligned PCGs of 33 Anguilliformes mitogenomes. Compared to the third codon, substitution in the first and second codons was relatively low (Figure 6). The impact on the whole amino acid released by the mutation of the third codons is rather vulnerable, in comparison with that of the first two codons, which possess a constrained evolving freedom [53,54]. The substitution in the third codons, even in insect mitogenomes, remains observable, as the amino acid products were relatively insensitive to the variation in the third codons [55].  The exploration of a substitution saturation index was implemented, utilizing the aligned PCGs of 33 Anguilliformes mitogenomes. Compared to the third codon, substitution in the first and second codons was relatively low (Figure 6). The impact on the whole amino acid released by the mutation of the third codons is rather vulnerable, in comparison with that of the first two codons, which possess a constrained evolving freedom [53,54]. The substitution in the third codons, even in insect mitogenomes, remains observable, as the amino acid products were relatively insensitive to the variation in the third codons [55].

Gene Rearrangement
Since one of the published A. meeki mitogenomes (MN616974) was not annotated, another A. meeki mitogenome (KX641476) was selected to conduct a comparative genomic analysis. The results indicated that there were not any other genes behind the tRNA-Thr

Gene Rearrangement
Since one of the published A. meeki mitogenomes (MN616974) was not annotated, another A. meeki mitogenome (KX641476) was selected to conduct a comparative genomic analysis. The results indicated that there were not any other genes behind the tRNA-Thr in KX641476; three typical genes, ND6, tRNA-Glu, and tRNA-Pro, cannot be identified, as in most mitochondrial genomes of vertebrates. On the contrary, both A. meeki (this study) and A. shiroanago have experienced gene rearrangement, with the rearranged protein-coding ND6, transfer RNA Glu, and a duplicated non-coding region, D-loop ( Figure 7A).

Phylogenetic Analysis
Fourteen highly related families (Eurypharynx pelecanoides and Saccopharynx lavenbergi, functioning as outgroups) were selected to explore the phylogenetic location of A. meeki. The phylogenetic analysis was implemented with the information revealed by phylogenetic trees (BI and ML), based on 12PCGs (the highly varied sections were Generally, a teleostean group possesses unique or highly similar gene arrangements. Nevertheless, two types of gene arrangements have been detected in the mitogenomes of Anguilliformes. As observed in the mitogenomes of Anguilla marmorata (Anguillidae), Ilyophis brunneus (Synaphobranchidae), Serrivomer sector (Serrivomeridae), and Gymnothorax minor (Muraenidae), which possess the typical mitochondrial structure, the genes ND6 and the conjoint tRNA-Glu were found to be typically situated between ND5 and Cytb, among which the Cytb was sandwiched by tRNA-Glu and tRNA-Pro, without the insertion of a D-loop between Cytb and ND6 [56] ( Figure 7C). In the mitogenomes of A. meeki and A. shiroanago, ND6 and tRNA-Glu were translocated upstream of tRNA-Pro and downstream of tRNA-Thr, respectively, accompanied by the insertion of a larger non-coding region. Nucleotide composition and the base arrangement of two control regions display similarities to a high and low degree in the mitogenomes of A. meeki and A. shiroanago, respectively.
Three typical models are expected to explain the gene rearrangement features in Ariosoma mitogenomes. The model, based on the recombination hypothesis, was initially proposed for gene rearrangement in the nuclear genome and was generally adopted to explain small-fragment exchanges and inversion events in mitogenomes. However, our comparative analysis shows that the Ariosoma mitogenomes have undergone gene translocation and genome-scale expansion. Therefore, the discovery of mitochondrial gene rearrangement in Ariosoma mitogenomes is too far-fetched to be explained by this model. The TDNL model placed an emphasis on non-random loss. In this case, the duplication followed by the loss of genes is a predesigned resorting to the corresponding transcriptional polarity and position in the genome; the genes are clustered in the same polarity (light-or heavy-strand coding) and the gene order remains unchanged (for example, the GCT cannot be produced by gene loss alone during the duplication from TCG to TCGTCG) [18,57]. This hypothesis does not apply to the findings, as the structural variation in Ariosoma mitogenomes is not the result of alterations in the genes' transcriptional polarity. The TDRL model underlined the rearrangement, based on the incomplete deletion of repeated genes [58]. In the A. meeki mitogenome, an extra D-loop was found downstream of tRNA-Thr; both D-loops revealed a high degree of resemblance in terms of nucleotide composition and bas-e arrangement. In addition, ND6, combined with the tRNA-Glu genes, was translocated from the stream downstream of ND5 to the upstream of tRNA-Pro. In the TDRL model, the intergenic spacers or pseudogenes were commonly detected in the rearrangement region [59,60]. In the mitochondrial genome of A. meeki, the existence of a considerable interval among ND5 and Cytb further indicated the bias to this model ( Table 2). In light of that, the TDRL model was suggested for the interpretation of rearrangement events in the mitochondrial genomes of A. meeki and A. shiroanago; among them, more detections of intergenic spacers or pseudogenes may occur in the A. shiroanago mitogenome during the random loss process.

Phylogenetic Analysis
Fourteen highly related families (Eurypharynx pelecanoides and Saccopharynx lavenbergi, functioning as outgroups) were selected to explore the phylogenetic location of A. meeki. The phylogenetic analysis was implemented with the information revealed by phylogenetic trees (BI and ML), based on 12PCGs (the highly varied sections were manually corrected). The topology of both trees presented a high degree of consistency ( Figure 8) and indicated that families (e.g., the Nettastomatidae, Muraenesocidae, and Colocongridae) with gene rearrangement grouped together and formed an independent clade, while the rest of the families (e.g., the Anguillidae, Synaphobranchidae, and Muraenidae) with the traditional gene order formed another clade. The result implied the evolution and origin of the diversified Anguilliformes species, as this pattern of the presence/absence of gene arrangement may be a potential phylogenetic marker for the identification of this morphologically conservative group of eels at the suborder level.
Both the BI and ML phylogenetic trees revealed that A. meeki and A. shiroanago were strongly related to each other; they formed a separate clade of Congridae (BI posterior probabilities [PP] = 1; ML bootstrap [BP] = 100). Consistent with previous studies, Congridae was closely related to the Muraenidae, Myrocongridae, and Synaphobranchidae [8]. Except for the family Congridae, which was deemed to be the non-monophyletic taxa [10], the rest of the families unanimously revealed the monophyletic clade, while two phylogenetic trees highly supported the non-monophyletic Congridae, and two clades were observed; such clustering results were compatible with those of previously published studies [61]. Regarding the clustering results released by an increasing number of phylogenetic studies on Anguilliformes, an independent subfamily (including Ariosoma) may be forming or may even have formed. This may offer a new perspective on the Anguilliformes taxonomy, but more available mitogenomes are needed to test this hypothesis.

Conclusions
The arrangement of the mitochondrial genome is consistent in most vertebrate mitogenomes; two kinds of gene arrangements have been revealed in the order Anguilliformes, which leads us to some confusion about the mechanism involved. In particular, some Anguilliformes populations were known to converge into non-monophyletic results, such as species in the Congridae and Nettastomatidae [9]. As of now, there is no comprehensive elaboration for this, since the relevant research is in its infancy [63].
In this study, we reported the complete mitochondrial genome of A. meeki, analyzed the corresponding genomic information, compared it with the reported mitogenomes of A. meeki and congeneric species, and depicted the phylogenetic relationship among Anguilliformes species. The complete mitogenome of A. meeki exhibited a rearrangement of gene order, in contrast to the typical vertebrate mitogenomes; ND6 and tRNA-Glu genes were translocated upstream of tRNA-Pro and downstream of tRNA-Thr, respectively, accompanied by the insertion of a D-loop. Both D-loops exhibited a high degree of resemblance in terms of nucleotide composition and base arrangement. The phylogenetic analysis indicated that A. meeki was closely related to its congeneric, A. shiroanago. The two species clustered together as a separate clade of the Congridae family. According to previous phylogenetic studies on Anguilliformes, the non-monophyly of Congridae was reflected by utilizing a comparative genomic analysis, based on 33 The phylogeny within Anguilliformes can be further routed into two major clades: one consists of eels that presented non-gene rearrangement patterns (Synaphobranchidae, Nemichthyidae, and Serrivomeridae), while another one comprised eels demonstrating gene rearrangement features (Nettastomatidae, Congridae, Ophichthidae, Muraenesocidae, and Derichthyidae). Similar consequences can be identified in other phylogenetic research on Anguilliformes [10]. The phylogenetic tree implied the evolutionary relationship of diversified Anguilliformes species and revealed the non-monophyly of the Congridae. According to previous studies, gene rearrangements may have originated from a common ancestral species for the Congroidei population. These features, which are found on the mitochondrial genomes of the partial population, were viewed as cryptic evolutionary population traces to distinguish them from other related species that possess a relatively conserved morphology [20,62]. However, with few available Anguilliformes species with gene rearrangements, a thorough confirmation of this hypothesis or clarification of the phylogenetic relationships among Anguilliformes species could be unattainable; depictions of the phylogenetic relationships resort to abundant genomic resources among this taxon would be more authoritative.
Given the remarkable resemblance in terms of the composition (Tables 2 and 3) and the codon usage (Table 2) (Figure 2) of PCGs, the structure of tRNA (Figures 3 and 4) and the results of phylogenetic analysis (Figure 8) between A. meeki and A. shiroanago, it may be a proof of a congener and two species that belong to the Congridae genetically. The variations in mitogenome length may be attributed to incomplete sequencing or annotation. One of the D-loops that should have been located between tRNA-Thr and ND6 was found to be missing.

Conclusions
The arrangement of the mitochondrial genome is consistent in most vertebrate mitogenomes; two kinds of gene arrangements have been revealed in the order Anguilliformes, which leads us to some confusion about the mechanism involved. In particular, some Anguilliformes populations were known to converge into non-monophyletic results, such as species in the Congridae and Nettastomatidae [9]. As of now, there is no comprehensive elaboration for this, since the relevant research is in its infancy [63].
In this study, we reported the complete mitochondrial genome of A. meeki, analyzed the corresponding genomic information, compared it with the reported mitogenomes of A. meeki and congeneric species, and depicted the phylogenetic relationship among Anguilliformes species. The complete mitogenome of A. meeki exhibited a rearrangement of gene order, in contrast to the typical vertebrate mitogenomes; ND6 and tRNA-Glu genes were translocated upstream of tRNA-Pro and downstream of tRNA-Thr, respectively, accompanied by the insertion of a D-loop. Both D-loops exhibited a high degree of resemblance in terms of nucleotide composition and base arrangement. The phylogenetic analysis indicated that A. meeki was closely related to its congeneric, A. shiroanago. The two species clustered together as a separate clade of the Congridae family. According to previous phylogenetic studies on Anguilliformes, the non-monophyly of Congridae was reflected by utilizing a comparative genomic analysis, based on 33 mitochondrial genomes of the Anguilliformes species. Both phylogenetic trees strongly support the nonmonophyly of Congridae, offering a theoretical basis for more advanced phylogenetic studies of Anguilliformes. In addition, our results provide new insight into the evolution of A. meeki and an in-depth understanding of the mechanism of gene rearrangement in the A. meeki mitogenome.
Author Contributions: In this research report, Y.H.; conceived and designed the experiments and performed sample collection, performed the experiments, analyzed the data, prepared figures and tables, wrote the original draft preparation, and reviewed and edited drafts of the paper. Z.L., J.Y. and Y.Y. analyzed the data, prepared the figures and tables, and reviewed and edited the paper. L.F., C.J. and J.C. conceived and designed the experiments, analyzed the data, contributed reagents and materials, and reviewed and edited drafts of the paper. K.Z. and H.J. conducted sample collection, performed the experiments, analyzed the data, prepared the figures and tables, wrote the original draft preparation, reviewed and edited drafts of the paper, and approved the final draft. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement:
The study was conducted in accordance with the Declaration of Helsinki, and approved by the Institutional Review Board of Zhejiang University (protocol code GBT 35823-2018 and 2 July 2021).

Informed Consent Statement: Not applicable.
Data Availability Statement: Due to ethical restrictions, the data presented in this study are available to researchers eligible under the Research Ethics Board rules on request from the corresponding author.