Specimens of P. rufescens, S. maculata, Euapta sp., P. bidentata, P. ooplax and S. reticulata were collected along the China coast ( Table 1 ). All specimens were preserved in 95% ethanol prior to DNA extraction. Total genomic DNA of each single specimen was extracted using E.Z.N.A^® Tissue DNA kit (OMEGA) according to the manufacturer’s protocol. The paired-end libraries were constructed by TruSeq™ Nano DNA Sample Prep Kit (Illumina) and sequenced to produce 2×150 bp paired-end reads on an Illumina HiSeq 4000 platform in Shanghai Lingen Biotechnology Co., Ltd. The Trimmomatic V0.33 software (Bolger et al., 2014) was used to remove low quality sequences from the raw reads. De novo assemblies were conducted with SPAdes v3.10.1 (k-mer = 21-77) (http://bioinf.spbau.ru/spades). To identify the putative mitochondrial contigs, BLASTN searches (Altschul et al., 1997) were conducted using the apodan mitogenomes available from the NCBI database. All mitogenomes sequenced in this study have been deposited to GenBank ( Table 1 ).

The protein coding genes (PCGs), rRNA and tRNA genes were identified by the MITOS web server (Bernt et al., 2013b) using the genetic code for invertebrate mitochondria. Gene boundaries of the PCGs and rRNAs were further manually examined and adjusted by comparing them with other related apodan mitogenomes. Base composition was calculated using MEGA 5 (Tamura et al., 2011). The AT and GC-skews were calculated following the formulae: AT-skew = (A-T)/(A + T); GC-skew = (G-C)/(G + C) (Perna and Kocher, 1995), where positive values in skewness indicate that the coding strand favour more A and G, while negative values infer more Ts and Cs. The codon usage was counted by MEGA 5 (Tamura et al., 2011) based on the invertebrate mitochondrial genetic code. The mitochondrial codons were partitioned into GT-rich codons (GT, TG, GG, TT codons at the first and second positions), CA-rich codons (CA, AC, CC, AA codons at the first and second positions) and other codons. All statistical analysis was performed using IBM SPSS Statistics 19.

CREx (Bernt et al., 2007) was used to conduct pairwise comparisons of mitochondrial gene order patterns among the holothurian mitogenomes. The gene rearrangement events include transpositions (T), reversals (R), reverse transpositions (RT), and tandem-duplication-random-loss (TDRL). The gene order distances were measured by common intervals parameter. MLGO web server (http://www.geneorder.org/server.php, Hu et al., 2014) was used to reconstruct the ancestral state and gene rearrangement history in Apodida. The tRNA genes were not included in the MLGO analysis, because of their variable arrangements among the Apodida mitogenomes. The reference tree used for mapping gene arrangements was modified from the phylogenetic analyses in our study.

The phylogenetic relationships were constructed by partitioning sequences of 13 PCGs using Bayesian inference (BI) and maximum likelihood (ML). The dataset comprised 53 taxa, including 49 Holothuroidea species representing six orders (Apodida, Molpadida, Synallactida, Holothuriida, Elasipodida and Dendrochirotida) and four Ophiuroidea species used as outgroups ( Supplementary Table 1 ). Each of the 13 mitochondrial genes was aligned individually by MEGA 5 (Tamura et al., 2011) using the codon based alignment algorithm. The aligned genes were concatenated into a single supermatrix by Sequencematrix (Vaidya et al., 2011). The best-fit partition model was selected by PartitionFinder 2 (Lanfear et al., 2017), respectively. The phylogenetic relationships were analysed based on two data matrixes, i.e. nucleotide sequence from all codon positions (PCG₁₂₃), and a Neutral Transitions Excluded (NTE) recoding scheme, which can ameliorate the misleading influence of strand-bias in the data set by recoding bases at neutral and quasi-neutral positions (Hassanin et al., 2005). The NTE data were obtained using Tidyverse package in R v.4.3 (Wickham, 2017). The ML analysis was conducted by IQ-TREE web server (Trifinopoulos et al., 2016) with the best-fit partition schemes, using 1000 ultrafast bootstrap (UFBoot) replicates (Minh et al., 2013). The BI analysis was performed by MrBayes 3.2.6 (Ronquist and Huelsenbeck, 2003) applying the best-fit models, with the Markov chain Monte Carlo (MCMC) runs of 10,000,000 generations. The effective sampling size (ESS) (>200) was checked by Tracer v1.7 (Rambaut et al., 2013, 2018). The first 2500 trees (25%) were burning and the remaining 2500 sampled trees were used to compute the 50% majority-rule consensus tree and the posterior probabilities (PP). Pairwise amino acid distances (p-distance) were calculated in MEGA 5 (Tamura et al., 2011). And the results are presented as boxplots, which were constructed by IBM SPSS Statistics 19.

The divergence time was estimated by BEAST 1.8.1 (uncorrelated lognormal relaxed clock, Yule prior) (Drummond et al., 2012), with two criteria (PCG₁₂₃ and NTE) set as before. For each criterion, the BEAST analyses were performed for 10,000,000 steps, with sampling every 1000 steps, following a burn-in of 25% of the tree. The convergence was also checked by Tracer v1.7 (Rambaut et al., 2013, 2018) as the above BI analyses. The TreeAnnotator (Drummond et al., 2012) was used to summarize the maximum-clade-credibility tree topology with median ages and the 95% highest posterior densities (HPD). Divergence times were constrained with three time points from the fossil record according to a previous study (Miller et al., 2017). The crown group of Holothuroidea was constrained with a median age of 435 Ma with a lognormal distribution (431-449 Ma) (Reich, 2010; Reich et al., 2015). The median age of the stem Elasipodida was constrained at 387 Ma using a lognormal distribution (384-394 Ma) (Boczarowski, 2001). The stem Holothuriidae was calibrated at 246 Ma with a lognormal distribution (244-257 Ma) (Reich, 2004).

The dN, dS and dN/dS (ω) were estimated for each 13 individual PCG and the concatenated PCGs data set by PAML v.4.9e (Yang, 2007) under a phylogenetic context. Combining the divergence time and selection tests, we analysed and compared the evolutionary rates of mitogenomes between Apodida and other holothurians. The evolutionary rate was calculated as V_subs = K/2T, where K represent the mean number of substitutions for each nucleotide site, and T is the divergence time (Pesole et al., 1999). The synonymous substitution rate (V_syn) and non-synonymous substitution rate (V_nonsyn) were also calculated for every lineage. The divergence time and substitution numbers were extracted from the result file generated by TreeAnnotator and PAML v.4.9e, respectively. The p values with significant differences were constructed by IBM SPSS Statistics 19.

The six new mitogenomes of Apodida were all assembled into a closed circular contig and varied in length from 15,896 (P. bidentata) to 16,732 bp (P. ooplax). In the six species, the mitogenome contained 13 PCGs, two rRNA genes and 22 tRNA genes ( Supplementary Table 2 ). However, the gene arrangement differed from those of other Apodida mitogenomes, as well as other Holothuroidea ones ( Figure 1 ). The A+T content was from 61.44% in P. rufescens to 66.86% in P. bidentata on the coding strand of the whole mitogenome. The third codon position (PCG₃) exhibited the highest A+T content in E. godeffroyi (67.24%), P. bidentata (76.76%), P. ooplax (69.52%) and P. rufescens (65.74%). While in S. maculata and S. reticulata, the highest AT content occurred in major non-coding region (MNR) (75.56%) and rRNA (67.15%), respectively ( Supplementary Table 3 , Figure 2A ). The AT and GC skews were calculated for 49 holothurian mitogenomes ( Supplementary Table 3 , Figure 2B ) to measure the strand asymmetry. All the holothurian species fell into two groups: the first one contained thirty-seven species, i.e., ten species in Synallactida and all species in Molpadida, Holothuriida, Elasipodida and Dendrochirotida, presenting negative GC skews and weak negative or positive AT skews; and the second one included twelve species, i.e., one species in Synallactida and all species in Apodida, which are characterized by positive GC skews and negative AT skews. These results indicated that the apodan species have strand asymmetry reversal on the entire plus strand compared to other holothuroids.

Considering the contrasting patterns of strand compositional asymmetry in apodans, we chose eleven mitochondrial protein-coding genes on the plus-strand (atp6, atp8, cox1, cox2, cox3, nad1, nad2, nad3, nad4, nad4l, nad5) to explore if DNA asymmetric strand bias affects the amino acid composition in apodans. The abundance of amino acids encoded by GT-rich codons (C, V, F, G and W) was significantly higher in Apodida (Mann-Whitney U test, p < 0.001) ( Figure 2C ), while that of the CA-rich codons was significantly higher in other holothuroids (p < 0.001) ( Figure 2D ). Moreover, this significant differences (P < 0.001) occurred at almost all the individual amino acids, except I, N and E ( Supplementary Table 4 ), indicating their clear responses to shifting nucleotide use. Especially, Apodida contains twice as many L_UUR and V as do their orthologs in other holothuroids. On the other hand, the proportions of L_CUN, T and K in Apodida are approximately half the value found in the orthologs of other holothuroids. So the amino acid composition was influenced by the reversals of strand-compositional bias.

Analysis of pairwise amino acid distances of 13 PCGs showed a significantly higher genetic divergence between Apodida and the other holothuroid mitogenomes than any other pairs ( Supplementary Table 5 , Figure 3 ). Between Apodida spp. and other holothuroid species, the nad2 gene showed the highest genetic divergence from 0.632 to 0.657 and the cox1 gene presented the lowest genetic divergence from 0.222 to 0.245.

The gene order among the eleven mitogenomes of Apodida differed from each other, with the arrangement of both PCGs and tRNAs varying significantly ( Figure 1 ). A total of nine gene order patterns were known for the eleven available apodan mtgenomes. No conserved gene block was found for the Apodida mitogenomes. In the CREx analysis, the gene order distance ranged from 4 to 24 between Apodida spp. and other holothuroids, and from 12 to 22 between Apodida spp. and the Holothuroidea putative ground pattern. The lowest gene order distance among the apodan species was 16, and the highest one was 956. The gene order of Chiridota spp. was highly different from those of other apodan species ( Table 2 ).

Only the PCGs and rRNAs were included in the ancestral gene order reconstruction. The gene orders were mapped on the phylogenetic tree to recover the evolutionary pattern of PCG and rRNA genes ( Figure 4A ). Inference of gene arrangement evolution revealed the most likely ancestral gene order for Apodida (Node A1) as cox1-cox3-nad2-rrnL-nad5-cytb-atp8-atp6-cox2-rrnS-nad4-nad1-nad4l-nad3-nad6 ( Figure 4B ), which is retained by the most recent common ancestor (MRCA) of Synaptidae (A2) and that of the Chiridotidae sp. KJ-Belize-E1_1-P. ooplax clade (A8). One reversal even and three transposition events might have occurred in Chiridotidae sp. KJ-Belize-E1_1 (Pattern 7) and P. ooplax (Pattern 8), respectively. Two transposition and one reversal events have occurred in the MRCA of the S. reticulata-S. maculata-E. godeffroyi-Euapta sp.-Protankyra sp. IDXB3-P. bidentata clade (A3) and the Protankyra sp. IDXB3-P. bidentata clade (A7). The gene order of Protankyra sp. IDXB3 featured a transposition (Pattern 5), and P. bidentata displayed a transposition and one reversal event (Pattern 6). The gene order of the S. reticulata-S. maculata-E. godeffroyi-Euapta sp. clade can be explained by a reversal of the rrnS gene (A4), which was also retained by the four extant species (Pattern 1-4). The differences between the gene orders in the MRCA of Chiridotidae (A9) and Apodida (A1) can be explained by four transpositions, one reverse transposition and one reversal event. The cox3-nad3 gene block translocated their positions in Chiridota sp. and C. heheva mitogenomes. The difference in gene order between P. rufescens and the MRCA of Chiridotidae (A9) can be explained by two transpositions, four reversals and a tandem-duplication-random-loss (TDRL) event.

The ML and BI based on nucleotide sequences from all codon positions (PCG₁₂₃) recovered that the holothuroid species grouped into five well supported clades ( Figure 5A ). The orders Holothuriida, Dendrochirotida, Elasipodida and Apodida formed monophyletic clades, respectively, while Molpadida was not recovered as a sister group to Synallactida, instead appeared as a poorly supported subgroup within Synallactida. The order Apodida was placed in a well-supported basal branch, as a sister to the well-supported lineage containing all the remaining holothuroids. Notably, Apodida exhibits a long branch length in the tree, as well as the subclades. The eleven apodan species used in this analysis were grouped into Synaptidae + Chiridotidae clade. Within Synaptidae, the species were split into four sub-groups with high support values. The Protankyra and Euapta species formed two different branches. S. reticulata and S. maculata were grouped as a well-supported clade. The analyses showed a Chiridotidae sp. KJ-Belize-E1_1 + Patinapta ooplax clade as sister to other Synaptidae groups. Based on this topology, the Chiridotidae sp. KJ-Belize-E1_1 should be put in Synaptidae, but not as marked as its name. Three species belonging to the family Chiridotidae were clustered together.

The recoded nucleotide data according to the Neutral Transitions Excluded (NTE) model yielded a phylogenetic tree that is slightly different from that of the full nucleotide data ( Figure 5B ). Unlike in the analysis with PCG₁₂₃, all the five orders with multiple taxa were recovered as monophyletic in the NTE phylogeny. The order Molpadida with single species was recovered as sister group to the Synallactida clade, albeit with low support (65/0.98). In the NTE phylogenetic tree, Dendrochirotida, instead of Holothuriidaand, showed sister relationship with Molpadida + Synallactida. Relatively, the NTE analysis revealed a low variation in branch length. However, it did not ameliorate the long branch length of Apodida.

The molecular dating analysis under PCG₁₂₃ and NTE dataset generated similar patterns of node ages, but NTE data produced younger ages for most nodes. ( Figure 6 ). The split date between Apodida and the other holothuroids occurred in the boundary of Ordovician and Silurian (PCG₁₂₃ vs NTE, 434.59 vs 434.09 Mya). The Apodida split into Chiridotidae and Synaptidae branches in the Middle Carboniferous (351.78 vs 321.46 Mya). In the Late Permian (283.23 vs 261.87 Mya), the family Synaptidae diverged into two lineages. Chiridotidae sp. KJ-Belize-E1_1 and P. ooplax shared a common ancestor, which diversified in the Lower Cretaceous (127.23 vs 113.52 Mya). Other Synaptidae species shared a common ancestor, starting their divergence in the boundary of the Triassic and Jurassic (225.62 vs 198.01 Mya). The modern Protankyra species diverged in the Early Cretaceous (125.14 vs 117.48 Mya). The split date between the Euapta species and S. maculata + S. reticulata group occurred in the boundary of the Jurassic and Cretaceous (149.01 vs 121.83 Mya). The lineage of Chiridotidae diversified in the Middle Jurassic (171.58 vs 166.11 Mya).

The nucleotide evolutionary rate (V_subs) of the 13 PCGs in the 49 Holothuroidea mitogenomes were estimated and fixed in the topology of the NTE tree ( Figure 6 ). The mean substitution rates of Holothuroidea was 0.0048 subs/s/my, and different order had different evolutionary rates. The mean substitution rate of Apodida was 0.0069 subs/s/my, which was significantly higher than that of other holothuroids (0.0043 subs/s/my) (Mann-Whitney U-test, p = 0.039) ( Supplementary Table 6 ). The order of the average V_subs among the 13 PCGs in Holothuroidea was as follows: cox1 < cox3 < cytb < cox2 < nad4l < atp6 < nad1 < nad5 < nad2 < nad3 < nad4 < atp8 < nad6.

The strand asymmetry in nucleotide composition is one of the most remarkable features in animal mitogenomes, which may be mainly generated by the reversed asymmetric mutational constraints during the mtDNA replication process (Hassanin et al., 2005). The AT and GC-skew values are used to measure the nucleotide compositional asymmetry (Frank and Lobry, 1999). A negative GC-skew and negative or positive AT-skew are found in other holothuroids and other echinoderms (Sun et al., 2022a; Sun et al., 2023). However, apodans present an opposite pattern with positive GC-skews and negative AT-skews, suggesting their nucleotide usage bias towards an excess of T relative to A and of G relative to C nucleotides. The reverse strand compositional bias of Apodida presents in all the individual gene data sets, including PCGs, thus, may affect the amino acid composition (Min and Hickey, 2007; Sun et al., 2018), and further promote the homoplasious substitutions in the apodan holothuroids. Consequently, these biases may hamper the phylogeny reconstruction from the PCGs.

In the Metazoa, the taxa with reverse strand biases are usually coupled with rearranged gene order in the mitogenome (Hassanin et al., 2005). The Apodida seems a good representative. The gene order of the eleven apodan mitogenomes diverged substantially with other available holothuroid mitogenomes, and also present high interspecific variability. This association was also found in some molluscs, echinoderms and annelids (Scouras and Smith, 2006; Sun et al., 2018; Sun et al., 2021a). Sun et al. (2018) compared the strand asymmetry across 140 mollusc species and suggested that both the presence of the strand asymmetry and rearranged gene order maybe result from the same mechanism. The inversion of the control region is hypothesised to explain the occurrence of both the nucleotide bias and gene rearrangement (Hassanin et al., 2005; Scouras and Smith, 2006).

No putative ground pattern of Holothuroidea was found in Apodida, and the gene order in the two Apodida clades, Synaptidae and Chiridotidae, differed dramatically from each other, as well as from other holothuroids. Although the tRNA genes were not considered in the reconstruction of ancestral gene order, the gene arrangement of apodans still diversified. The conserved gene block cox1-cox3-rrnL-nad2-nad5-cytb-nad4 might represent the ground pattern for Synaptidae. The estimation of ancestral gene order reveals a total of four ancestral states (as marked in Figure 4 ). Only S. reticulata, S. maculata, Euapta sp. and E. godeffroyi shared the gene order of their most recent common ancestor, and none of other apodan gene arrangement pattern can represent their ancestral states. We suggest that higher variation in Apodida mitochondrial gene order may be not restricted to the analysed species, but should be found order-wide.

In the metazoans, high-level mitochondrial gene rearrangement usually corresponds to accelerated nucleotide substitution rates. This correlation has been observed in the nematode, insect, molluscs, brachiopod and crustacean mitogenomes (Hu et al., 2003; Shao et al., 2003; Rawlings et al., 2010; Luo et al., 2015; Plazzi et al., 2016; Sun et al., 2022b). Also, the high diversity of mitochondrial gene order in Apodida coupled with the high evolutionary rates observed in our present study. Some underlying factors, like relaxed repair mechanisms, and high mutational stress in combination with lower importance of mitochondrial efficiency, have been hypothesised to affect both gene rearrangement and substitution rate (Bernt et al., 2013a; Sun et al., 2021a).

In this study, we presented the first mitochondrial phylogenetic study for apodan holothuroids. Phylogenetic analyses supported the monophyly and basal position of Apodida (Synaptidae + Chiridotidae), which agree with the previous results base on morphology and molecular phylogenetic studies of multiple gene fragments (Becher, 1909; Haude, 1992; Smith, 1997; Kerr and Kim, 2001; Lacey et al., 2005; Reich, 2010; Smirnov, 2012; Miller et al., 2017). Within the Apodida lineage, an unnamed species Chiridotidae sp. KJ-Belize-E1_1 (MT877116) was placed in a branch as a sister group to synaptid P. ooplax, instead of Chiridotidae species. So the two families, Synaptidae and Chiridotidae, were not recovered as monophyletic, which seems not consistent with that reconstructed based on morphology data, as suggested by Smirnov (1998). In the study of Smirnov (1998), Synaptidae and Chiridotidae formed well supported sister group. Using cladistic analyses of 34 morphological characters, Kerr (2001) revealed the Chiridotidae was paraphyletic at the subfamily level, in which Chiridotinae was clustered with the Synaptidae clade instead of being a sister group to species of Taeniogyrinae. Owing to the absence of taxon sampling from Taeniogyrinae, the relationship between Chiridotinae and Taeniogyrinae cannot be validated in our study. The family Synaptidae was previously recognised as a paraphyletic clade in a phylogeny based on six gene fragments (COI, 16S, 12S, H3, 18S, 28S) (Miller et al., 2017). A similar placement of Synaptidae was also revealed based on mitogenomes in our analysis. Of course, if the species Chiridotidae sp. KJ-Belize-E1_1 was from the family Synaptidae, instead of Chiridotidae, the relationship between Synaptidae and Chiridotidae needs to be revised.

It has been suggested that a reverse strand compositional bias in the mitogenome can drive convergent nucleotide substitutions, which may further mislead phylogenetic analyses (Hassanin et al., 2005; Hassanin, 2006). This scheme has been used to overcome the confounded phylogenetic inference caused by the strand bias (Hassanin, 2006; Jones et al., 2007; Sun et al., 2023). In the NTE model, all neutral and nearly-neutral sites were removed, because the substitution types of these sites are most likely to be influenced by the reverse strand bias. To ameliorate the misleading effects of reversed patterns in nucleotide bias, a Neutral Transitions Excluded (NTE) strategy was applied in our analysis. In the present study, almost all the same lineages were recovered using both codon positions and the NTE model. The most notable exception was that the order Synallactida was revealed as monophyletic in the NTE phylogeny, with Molpadida as a sister group. The monophyly of Synallactida and Molpadida has been supported in recent molecular analyses (Miller et al., 2017). On the other hand, in the NTE phylogeny, the order Holothuriida was noticed as the next group to diverge after the order Elasipodida. These results are consistent with previous views (Miller et al., 2017; Sun et al., 2021b). Thus, we believe the NTE-phylogeny may be a more accurate summary of holothuroid relationships. However, the order Apodida still showed a long branch in the NTE-phylogeny. It’s difficult to know why this is, but we see an alternative explanation. The NTE algorithm recommends RY recoding of the first codon positions of CTN (L_CUN) and TTN (F and L_UUR), and of all three codon positions of ACN (T), ATN (I and M), GCN (A), and GTN (V) (Hassanin et al., 2005). Our results indicate that, the greatest differences in amino acid abundance between Apodida and other holothuroids are in L_CUN, L_UUR, V, T and K. Thus, NTE recoding may lead to a loss of some useful information for phylogenetic relationships. In addition, the abundance of Lysine (K) is significantly different between Apodida and other holothuroids, while NTE recoding cannot account for this bias.

Former phylogenetic analyses of bilateralian mitogenomes revealed that some taxa with long-branch frequently present a negative AT skew and positive GC skew (Bernt et al., 2013a). Similar patterns of the AT and GC skews were discovered in the apodan species examined in this study. This coincides with the hypothesis that taxa featured by reverse strand bias tend to gather together owning to the long-branch attraction artifacts (Hassanin et al., 2005). The strand asymmetry and the high amino acid genetic distance between Apodida spp. and other holothuroid species indicated distinct selective pressures and faster evolutionary rates in apodan mitogenomes. In the phylogenetic analysis using mitochondrial (COI, 16S, 12S) and nuclear (H3, 18S, 28S) gene fragments, Synaptidae and Chiridotidae did not show long branches, while the branch length for Myriotrochidae (Myriotrochus spp.) was much longer than other branches, and the authors believe that this disparity resulted from COI sequence (Miller et al., 2017). Thus, the accelerated substitution rate may occur more in the mitogenomes of Apodida as indicated by (Sun et al., 2021b).

The divergence time estimated here provided insights into the understanding of the apodan evolutionary history. Apodida represents one of the most ancient clades within Holothuroidea, and the split between the Apodida and other holothuroid orders occurred between Silurian and Ordovician. Under PCG₁₂₃ data, the ancestor of Apodida originated in the Early Carboniferous (351.78 Mya) and the major diversification happened in the Cretaceous. The NTE algorithm reduced the influence of strand bias by removing the bases at neutral and nearly-neutral positions, thus generating younger origin time (321.46 Mya). Our estimated origin time is much earlier than the estimation of Miller et al. (2017) (259-282 Mya) based on gene fragments (COI, 16S, 12S, H3, 18S, 28S). However, one matter of agreement is that both two studies support that the apodans radiation transitioned the Permian-Triassic mass extinction event, which was believed to result in approximately 80% marine biodiversity loss (Alroy et al., 2008). This result suggested that apodans can survive in the extinction events that occurred in the boundary of the Permian and Triassic, which were accompanied by abrupt climate change and widespread oceanic anoxic zones.

The evolutionary rate of mitogenomes and individual genes in Holothuroidea were estimated in our study. The average evolutionary rate of Holothuroidea mitogenomes was 0.0048 subs/s/my. The order of the mean V_subs among the six holothuroid orders was as follows: Apodida (0.0069 subs/s/my) > Synallactida (0.0050 subs/s/my) > Holothuriida (0.0039 subs/s/my) > Dendrochirotida (0.0022 subs/s/my) > Elasipodida (0.0021 subs/s/my) > Molpadida (0.0018 subs/s/my). Apodida showed significantly higher V_subs than that of other holothuroids (0.0043 subs/s/my), indicating that the mitogenome of Apodida was subjected to more relaxed selection pressure. One possible explanation for the accelerated evolutionary rate discovered in apodans is an adaptation to their special motion and respiratory patterns. The two remarkable features of Apodida are the absence of tube feet and the respiratory tree (Martins and Souto, 2020). Thus, various parts of the body, such as the body wall, tentacles, papillae and/or dermal ossicles, are used to locomote (Kerr, 2000), and gas exchange is achieved through their thin body wall (Smiley et al., 1991). The locomotory capacity and respiratory efficiency can affect the evolution of animal mitogenome (Das, 2006). Furthermore, the positive correlation between the low locomotory capacity and the stronger purifying selection on the mitogenome has been confirmed in molluscs, crustaceans and insects (Mitterboeck and Adamowicz, 2013; Sun et al., 2017; Jakovlić et al., 2021). Following the evidence that mitogenomic evolution is faster in the species with low locomotory ability, we hypothesise that locomotory capacity may be one of the factors which explain the higher mitogenomic evolutionary rates. In other words, the apodans may have evolved specific adaptive mechanisms for locomotion and respiration, which is linked to the fast evolution of the mitochondrial respiratory genes.

The individual genes showed different evolutionary rates in our analysis. The nad6 gene has the fastest evolutionary rate (1.8031 subs/s/my), while cox1 has the slowest one (0.0035 subs/s/my). This pattern indicated that the nad6 gene was under relaxed selection pressures, at the same time, the cox1 gene received strict selection constrains. It is noteworthy that the evolutionary rate of the cox1 gene estimated by our molecular clock differs drastically from that of echinoderms (0.02-0.05 per million years) based on the divergence after the closure of the Isthmus of Panama (Benzie, 1999; Lessios et al., 2001). Thus, the evolutionary rate of cox1 corresponding to the closure of the Central American Seaway should be used with caution. Just as suggested by Marko (2002) and Borrero-Perez et al. (2010), the rates of molecular evolution can be greatly overestimated when the clocks of marine organisms were calibrated with final seaway closure.