Mitochondrial Genomes of Kinorhyncha: trnM Duplication and New Gene Orders within Animals

Many features of mitochondrial genomes of animals, such as patterns of gene arrangement, nucleotide content and substitution rate variation are extensively used in evolutionary and phylogenetic studies. Nearly 6,000 mitochondrial genomes of animals have already been sequenced, covering the majority of animal phyla. One of the groups that escaped mitogenome sequencing is phylum Kinorhyncha—an isolated taxon of microscopic worm-like ecdysozoans. The kinorhynchs are thought to be one of the early-branching lineages of Ecdysozoa, and their mitochondrial genomes may be important for resolving evolutionary relations between major animal taxa. Here we present the results of sequencing and analysis of mitochondrial genomes from two members of Kinorhyncha, Echinoderes svetlanae (Cyclorhagida) and Pycnophyes kielensis (Allomalorhagida). Their mitochondrial genomes are circular molecules approximately 15 Kbp in size. The kinorhynch mitochondrial gene sequences are highly divergent, which precludes accurate phylogenetic inference. The mitogenomes of both species encode a typical metazoan complement of 37 genes, which are all positioned on the major strand, but the gene order is distinct and unique among Ecdysozoa or animals as a whole. We predict four types of start codons for protein-coding genes in E. svetlanae and five in P. kielensis with a consensus DTD in single letter code. The mitochondrial genomes of E. svetlanae and P. kielensis encode duplicated methionine tRNA genes that display compensatory nucleotide substitutions. Two distant species of Kinorhyncha demonstrate similar patterns of gene arrangements in their mitogenomes. Both genomes have duplicated methionine tRNA genes; the duplication predates the divergence of two species. The kinorhynchs share a few features pertaining to gene order that align them with Priapulida. Gene order analysis reveals that gene arrangement specific of Priapulida may be ancestral for Scalidophora, Ecdysozoa, and even Protostomia.


Introduction
Mitochondrial genomes provide a set of important tools for evolutionary studies of animals owing to their accessibility and higher evolutionary rate in comparison with the nuclear genomes. The animal mitochondrial genomes are compact and typically encode 13 proteins of the respiratory chain (nad1-6, nad4L, cox1-3, cytb, atp6, and atp8), two subunits of ribosomal RNA and a variable complement of transfer RNAs on a circular DNA molecule with an average size of about 15 Kbp [1,2]. In addition to the gene sequence data, the mitochondrial genome displays other features that are also applicable for evolutionary studies. Although the gene content of animal mitogenomes is nearly constant, the arrangement of these genes on the DNA molecule varies among major taxonomic groups. It was suggested that changes in the mitochondrial gene order can be used as cladistic characters for studying relationships among higher-level taxa [3]. In some taxonomic groups such as insects and vertebrates the gene order tends to be conservative, while in some other groups such as Mollusca [4], Bryozoa [5][6][7], Acari [8], and Tunicata [9,10] it is highly variable. Generally, the taxon-specific gene order remains identical over long periods of time [11][12][13]. Furthermore, the analyses of gene order in bilaterian mitogenomes revealed conservation of specific gene blocks encompassing both protein-coding and ribosomal RNA genes [14]. The gene order within blocks remains uniform in most species, while the blocks themselves experience transpositions. The arrangements of these blocks allow to deduce putative ground patterns for some bilaterian taxa, including Ecdysozoa [15], Lophotrochozoa [16], and Deuterostomia [17]. However, the ground patterns for taxa with highly variable gene order (such as Nematoda or Chaetognatha) or for higher-level taxa (such as Bilateria or Metazoa) remains intractable. Within Ecdysozoa the mitochondrial gene order is generally more variable than in other taxa of the same rank, and it is potentially a valuable resource for phylogenetic information in the group [18][19][20].
The phylogenetic relationships within Ecdysozoa remain contradictory [21]. One of the important for evolutionary studies but relatively poorly studied ecdysozoan groups is Kinorhyncha. It is a small phylum of free-living, meiobenthic segmented pseudocoelomate wormlike invertebrates that accommodates 222 described species [22,23] distributed worldwide. Usually they are much smaller than 1 mm, but a few Arctic species reach the size of 1.1 mm. They inhabit the upper layers of marine sediment with densities of 1-10 animals per 10 cm 2 in the deep sea regions and 45 animals per 10 cm 2 in the shallow regions [24,25]. Taxonomically, Kinorhyncha is classified with Priapulida and Loricifera into a clade Scalidophora on the basis of morphological similarities [26]. The resent fossil findings point to early divergence of these three phyla-about 535 million years ago in the Cambrian Period [27]. Currently, the monophyly of a group uniting Kinorhyncha and Priapulida is confirmed by ribosomal phylogenies [28][29][30] and multigene phylogenetic analyses [31,21]. However, only two mitochondrial genomes from the Scalidophora lineage have been sequenced to date, both from representatives of phylum Priapulida [15,32].

Material collection and DNA extraction
Adult individuals of Pycnophyes kielensis were collected by bubbling from the upper layers of marine sediment in the vicinity of the Kartesh Biological Station, Zoological Institute, Russian Academy of Sciences (Chupinskaya Bight of Kandalaksha Bay, White Sea) and delivered live to the laboratory. Specimens of Echinoderes svetlanae were caught by trawling (Rugozerskaya Inlet, Kandalaksha Bay, White Sea) in the vicinity of the White Sea Biological Station of Moscow State University, collected by bubbling, and fixed by ethanol. Permissions to fieldwork were issued by the authority of Pertsov White Sea Biological Station of Lomonosov Moscow State University and Kartesh Biological Station of Zoological Institute of Russian Academy of Sciences. Field studies did not involve endangered or protected species.
Tissues of kinorhynchs were homogenized in 100 μl of buffer solution containing 0.01 M TrisHCl, 0.1 M EDTA, and 0.15 M NaCl (pH 8.0) and then lysed with 0.5% sodium dodecyl sulfate (SDS) and Pronase (100 μg/ml) at 37°C. DNA of both species was isolated from a few dozens of individuals using the DIAtom DNA Prep kit (Isogen, Russia) following the protocol provided by the manufacturer.  Genetic code and nucleotide composition AT and GC skew were determined for the complete mitochondrial genomes (major strand) according to the formula AT-skew = (A − T) / (A + T) and GC-skew = (G − C) / (G + C) [40], where the letters stand for the absolute number of the corresponding nucleotides in the sequences. Genetic code and codon usage were analyzed by GenDecoder v1. 6

Gene order analysis
We used CREx [44] for pairwise comparisons of kinorhynch mitogenomes between each other and with putative ground patterns for Panarthropoda, Priapulida, Lophotrochozoa and Deuterostomia constructed on the basis of conservative gene blocks described in [14]. CREx determines the most parsimonious genome rearrangement scenario given the gene order of two genomes, accounting for transpositions, reverse transpositions, reversals, and tandem-duplication random-loss events. The most parsimonious rearrangement scenarios for user trees were found by TreeREx [45]. The mitochondrial genomes from the OrganelleResource database of NCBI [46] were used for comparison with kinorhynch gene orders. The ground patterns for groups Panarthropoda, Lophotrochozoa, and Deuterostomia were obtained from previous studies [15][16][17]. Duplicated trnM genes in mitochondrial genomes were found using mitotR-NAdb [47].

Phylogenetic analysis
The protein-coding gene sequences of E. svetlanae and P. kielensis were translated using the invertebrate mitochondrial genetic code and aligned with sequences from 80 metazoans, selected from the NCBI's organelles genome database (S1 Table). We excluded Nematoda from the phylogenetic analysis to decrease the influence of the long branch attraction (LBA) artefact [14,32].
The aminoacid sequences were aligned with MAFFT v7.130 [48] using the pairwise Needleman-Wunsch algorithm (-globalpair) with the offset parameter set to 0.123 (-ep 0.123) and the maximum number of iterative refinement set to 1000 (-maxiterate 1000). The alignments were trimmed with trimAl 1.2rev57 [49] using a gap threshold of 0.9 and a similarity threshold of 0.0005 over a window of size 3 (-w 1), and concatenated using SCaFoS The nucleotide sequence alignments were constructed on the basis of the aminoacid alignments using the TranslatorX program [57], trimmed according to the mask derived from the trimmed aminoacid alignments and concatenated with SCaFoS. The inferences with the nonstationary models BP and BP+CAT were performed with NH PhyloBayes 0.2.3 [58], and the inferences with the CAT and CAT+GTR models were performed with PhyloBayes 4.1c [54]. Each analysis was run in 2 independent chains, the analyses with the CAT and CAT+GTR models were sampled across 20,000 cycles and the analyses with the nonstationary models were sampled across 200 points; the trees were summarized with a 10% burn-in.
Transfer RNA gene duplication was analyzed by phylogenetic approach. Each tRNA was divided into acceptor, anticodon, D-and T-arm regions according to the predicted secondary structure obtained by MiTFi [37]. Each stem or loop region was aligned with the corresponding region from other tRNAs and alignments were concatenated. Concatenated tRNA alignment was minimized to a mask of 52 bp. A NEXUS file was generated from the alignment by the STEMS program [59] accounting the predicted secondary structure and creating stem and loop partitions. Phylogenetic inference for tRNA genes was performed by MrBayes v.3.2.5 [60]. The evaluation of site-wise evolutionary rates was performed under the four-by-four model with covarion for the loops partition and the doublet model without covarion for the stems partition. Two independent runs of four Markov Chain Monte Carlo (MCMC) were performed for 1,000,000 generations, sampling trees every 1,000 generations. The first 500 trees were discarded as burn-in, and the remaining set was used to generate a consensus tree with posterior probability values. The resulting Bayesian tree was visualized in MEGA 5.2.2 [56].

Mitochondrial genome organization and nucleotide composition
Mitochondrial genomes of Echinoderes svetlanae (GenBank KU975552) and Pycnophyes kielensis (GenBank KU975551) are closed circular DNA molecules with lengths of 15304 bp and 14985 bp respectively. The kinorhynch mitogenomes have a rather low GC content of about 26% in both species (Table 1), although the lowest known GC content in mitochondrial genomes is 12,6% [14]. The GC content in tRNA and rRNA genes of kinorhynchs is even lower relative to the whole mitochondrial genome ( Table 2), but in protein-coding genes (PCGs) the GC content slightly exceeds the average. The non-coding regions in both species have a lower GC content than the whole genome average. The GC and AT skews characterize the asymmetry of nucleotide content between the two strands of mitochondrial DNA [40,61]. The prevalence of thymine over adenine and guanine over cytosine in the major strands provides negative AT-skew and positive GC-skew in kinorhynchs similarly to most other animals [14]. The mitogenomes of E. svetlanae and P. kielensis contain a common set of 37 metazoan mitochondrial genes (13 PCGs, 22 tRNA genes, and two rRNA genes) and one additional methionine tRNA gene (Tables 3 and 4). All genes in both species are located in the major strand.
The non-coding regions in the mitogenome of E. svetlanae are 860 bp in total and consist of 33 intergenic segments, ranging from 1 to 270 bp and include two major non-coding regions of more than 100 bp. The non-coding part of P. kielensis mitogenome consists of 32 intergenic segments ranging from 1 to 245 bp and totaling 941 bp, where three major non-coding regions have a length of more than 100 bp.
The mitogenome of P. kielensis has two gene overlaps: the first one is between nad4l and trnD (11 bp) and the second one is between trnA and trnS1 (5 bp). In the mitogenome of E. svetlanae there are also two overlaps: one between trnL1 and atp8 (18 bp) and another between trnQ and trnM1 (4bp).
Excluding all termination codons, the cumulative length of 13 PCGs of E. svetlanae is 10446 bp, encoding 3722 amino acid residues. The cumulative length of PCGs of P. kielensis is 10650 bp, encoding 3637 amino acid residues. The GC% at the first two codon positions exceeds the average GC% for the whole mitochondrial genome (Tables 2 and 5), and the GC% at the third codon position is very low-14% and 16% for E. svetlanae and P. kielensis respectively.
Transfer RNA genes in kinorhynchs are distributed throughout the circular molecule and have a total length of 1556 bp in E. svetlanae and 1431 bp in P. kielensis. Mitochondrial tRNAs range from 62 bp to 74 bp in E. svetlanae and from 56 bp to 66 bp in P. kielensis (Tables 3 and  4). The secondary structures of the E. svetlanae and P. kielensis tRNAs predicted by MiTFi are shown in additional files (S3 and S4 Figs). All tRNA genes were folded into the typical cloverleaf secondary structure, except for the E. svetlanae trnS1 gene, where the dihydrouridine arm is simplified to a loop, and the P. kielensis trnF with a loop instead of the TCC arm (the T-loop).
Two rRNA genes (rrnL and rrnS) are located on the major strands. The rrnL gene is located after the trnK in both species and before nad1 gene in E. svetlanae or nad3 gene in P. kielensis. The rrnS gene is located before trnG gene in both species and after trnL2 gene in E. svetlanae or trnV gene in P. kielensis (Tables 3 and 4). The length of the rrnL gene is 1014 bp in E. svetlanae and 982 bp in P. kielensis. The length of the rrnS gene is 766 bp in E. svetlanae and 733 bp in P. kielensis.

Duplicated tRNA genes
The overwhelming majority of metazoan mitogenomes include two tRNA genes for serine and leucine and only one tRNA gene for each of the other 18 amino acids [62,63]. Two serine and leucine tRNA genes are associated with two codon groups coding each of these amino acids: each tRNA anticodon is adapted for recognition of a codon group.
The mitogenomes of both E. svetlanae and P. kielensis encode two methionine tRNA genes (Fig 2). The tRNA M1 gene in both species is located after a glutamine tRNA gene, and the M2 gene-before nad6 gene (Tables 3 and 4). Priapulids Priapulus caudatus [15] and Halicryptus spinulosus, which are presumed to be the closest relatives of kinorhynchs, have only one methionine tRNA gene which is positioned after a glutamine tRNA gene like the trnM1 of E. svetlanae and P. kielensis. The predicted kinorhynch methionine tRNAs have a classic "clover leaf "secondary structure with two lateral arms, and their loops and stems have similar sizes. There are two compensatory changes in the M2 tRNA helices of E. svetlanae and P. kielensis, and in M1 tRNAs there is one compensatory change (Fig 2). The presence of these compensatory changes strongly suggests that kinorhynch methionine tRNAs are not pseudogenes. The Bayesian inference with tRNA genes of E. svetlanae and P. kielensis revealed an association between the same amino acids from both species for the majority of tRNAs including two pairs of methionine tRNA genes (Fig 3). The observed placement of trnM1 and trnM2 suggests that the additional gene originated by duplication of the methionine tRNA gene before the divergence of E. svetlanae and P. kielensis. The ancient nature of this duplication is also evidenced by the large difference in the primary structures of paralogous genes and the T-arm length difference in the predicted secondary structures. The functional significance of mitochondrial tRNA-Met gene duplication in kinorhynchs is not clear, but it might represent Note that in plastid genomes the presence of two tRNA-Met genes is relatively common [64]: one methionine tRNA functions during the elongation phase of protein synthesis, and the other is charged with formyl methionine to function as an initiator tRNA. In mitochondrial genomes the duplication of trnM was also found in some placozoans, cnidarians, insects, vertebrates. In molluscs and tunicats, each methionine tRNA, tRNA-Met(AUG) and tRNA-Met (AUA), can recognize a specific methionine codon instead of one tRNA with AUR specificity [65][66][67]. In some sponges the product of one trnM gene is post-transcriptionally edited to function as an additional isoleucine tRNA [68,69]. Besides the trnM there are other rare cases of mitochondrial tRNA duplications. For example, a duplicated trnL (CUN) was found in the Japanese freshwater crab Geothelphusa dehaani [70] and a duplicated trnI is found in three blow-fly species from genus Chrysomya [71]. In other cases one copy of the duplicated tRNA gene accumulates multiple mutations and undergoes deletions decreasing the gene length in comparison with functional tRNAs [70,71]. Thus, one of the duplicated isofunctional tRNA genes may become a pseudogene that is eventually eroded from the mitogenome. Duplicated tRNAs may also appear as a result of tRNA remolding (the change of tRNA specificity)-like in gecko Tropiocolotes tripolitanus with two trnQ, one of which has evolved from trnR [72]. However in this case the duplicated genes strongly differ in primary and secondary structures and the remolded tRNA likely remains as an inactive form.

Genetic code and codon usage
The kinorhynch mitochondrial PCGs do not contain all of the possible codons (Table 6). There are no CTC and CTG codons for leucine and TCC codon for serine in E. svetlanae. The PCGs of P. kielensis have no arginine codon CGC.
Both GenDecoder and FACIL programs determined that the majority of kinorhynch codons code the same amino acids as the standard invertebrate mitochondrial code. However, they also found some differences from this code in both species (S2 Table). GenDecoder specified the TGC codon as proline in E. svetlanae and as alanine in P. kielensis. The ACG codon is specified by GenDecoder as serine in P. kielensis, and as threonine in E. svetlanae. In E. svetlanae FACIL specified TTC as leucine and TGT as alanine, and in P. kielensis it specified the ATC codon as leucine.
Likewise, there are some differences in specifying the ATG and ATA codons that normally code methionine in invertebrate animals. FACIL specified these codons in both species as leucine, and GenDecoder specified the ATG as leucine in E. svetlanae. Considering that kinorhynchs have two tRNA genes with a conventional methionine anticodon it is possible that one of these tRNAs may be actually charged with leucine. In addition, the long branches of kinorhynchs in the phylogenetic tree (Fig 4) may be partially due to incorrect translation of ATG codons. It is possible that kinorhynch trnMs are charged by different amino acids (methionine or leucine), and introduce them to the protein chain randomly or deterministically. However our data do not allow us to draw a clear conclusion about the exact specification of these codons.
The results of codon predictions by the both programs do not completely coincide in two kinorhynch species. Four of the seven controversial codons (TGC, ACG, TTC and ATC) analyzed by the programs are very rare in both species and consequently may be specified mistakenly. We compared three other controversial codons with tRNA anticodons found in the kinorhynch mitochondrial genomes. The TGT codon does not form a complementary pair with the alanine tRNA at all three codon positions. The ATG and ATA codons can form two of the three complementary bonds with the leucine tRNAs. At this point it is not possible to answer whether we are observing a different genetic code in kinorhynch mitogenomes or whether the deviations detected by the genetic code analysis have other explanations and we are actually dealing with a standard invertebrate mitochondrial code. Previous research on cox1 predicted the invertebrate mitochondrial code in kinorhynchs [73].
The codon usage in the mitochondrial genomes of E. svetlanae and P. kielensis shows a strong preference for synonymous codons ending with thymine or adenine (Table 6). E. svetlanae and P. kielensis favor the NNT codon from the NNN or NNY codon families and the NNA from the NNR codon families. All of the missing codons end with guanine or cytosine and consist of guanine or cytosine completely or by 2/3. We predict five types of start codons in kinorhynch species (Table 7). The majority of PCGs are predicted to start with codons ATG or ATA. However, based on the amino acid alignments it is clear that additional codons are likely employed for translation initiation in kinorhynchs. Additional start codons are specific for some genes: cox2 is likely initiated by the leucine codon TTG whereas nad1 is likely initiated by the isoleucine codon ATT in both species; cytb is likely initiated by different additional codons in both species (isoleucine codon ATT in E. svetlanae and valine codon GTG in P. kielensis). Moreover, cox1 in P. kielensis is likely initiated by a valine codon GTG. All these codons have been already mentioned to function as start codons in some mitochondrial PCGs, for example, TTG in Epinephelus coioides [74], GTG in Tylototriton taliangensis [75], Scapanulus oweni [76] and Euthynnus affinis [77], ATT in S. oweni. Additional start codons are not rare in mitochondrial genomes, but the presence of four (in E. svetlanae) or five (in P. kielensis) start codons is unusual. There are two stop codons (TAA and TAG) in kinorhynch species (Table 7). Both species prefer the stop codon TAA ending with an adenine.   The positive GC skew and negative AT skew of the coding strand affects the amino acid composition bias in PCGs. In kinorhynchs there is an excess amount of amino acids coded by GT-rich codons: Phe, Gly, Val and Trp. At the same time the number of amino acids coded by AC-rich codons (Thr, Pro, Asn, His and Gln) is noticeably lower.

Phylogenetic analysis
The Bayesian inference with concatenated mitochondrial protein alignments recovers wellsupported conventional monophyletic groups of Deuterostomia, Lophotrochozoa, and Ecdysozoa (S5 Fig). Two kinorhynch species are grouped together forming a long branch, which is placed within arthropods contradicting the modern conceptions of ecdysozoan taxa. To decrease the possible impact of long-branch attraction (LBA) caused by divergent kinorhynch sequences, we applied the removal of fast-evolving sites from the alignment. After removal of the fastest evolving sites from the protein alignment, the conventional taxa have retained their monophyly and the long branch of kinorhynchs was positioned near the root of Ecdysozoa (Fig 4). However, we did not observe the kinorhynchs group with Priapulida (Priapulus caudatus and Halicryptus spinulosus) in what would constitute the monophyletic Scalidophora, contrary to the results obtained with rRNA or nuclear protein coding gene datasets [28-31, 21]. The LBA artifacts in mitochondrial PCGs were noted in other ecdysozoans [32], and, perhaps, the observed placement of kinorhynchs in the mitochondrial protein tree is associated with their accelerated evolutionary rate. It was observed that the basal branching in Ecdysozoa is highly sensitive to taxon sampling and the choice of model for the mitogenome data [32]. In our analysis the first branch of Ecdysozoa is Onychophora with posterior probability 0.41. Generally the basal position of onychophores is not supported by molecular phylogeny, but an onychophore-like ancestor was suggested for Ecdysozoa by paleontological data [78].
We tested whether the alternative positions of Kinorhyncha in the mitochondrial PCGs tree are significantly worse than the reconstructed topology. The majority of tested alternative topologies differ insignificantly from the Bayesian tree. The tree with the monophyletic Scalidophora clade has the AU-test p-value of 0.563, and the tree with Scalidophora as the basal branch within Ecdysozoa has the AU-test p-value of 0.378. The results do not contradict the hypothesis of Scalidophora monophyly or the basal position of Scalidophora within Ecdysozoa. The low values of posterior probabilities and results of the AU-test indicate poor resolution of the mitochondrial ecdysozoan tree, which agrees with previous observations [32].
Because the poor resolution of the mitochondrial tree may be associated with the GC bias, we analyzed the nucleotide sequences that underlie the protein alignment using nonstationary models implemented by NH PhyloBayes as well as stationary PhyloBayes models. The Bayesian inference with nucleotide data revealed some differences with the aminoacid data inference depending on the employed model (S6-S9 Figs). The trees built with stationary models CAT or CAT+GTR (S6 and S7 Figs) are broadly similar and recover monophyletic Ecdysozoa, while Deuterostomia and Lophotrochozoa are split into several polyphyletic groups. Kinorhynchs are grouped with onychophores into a clade sister to priapulids (S6 Fig)  or within arthropods (S7 Fig) with low posterior probability. Trees under the nonstationary models BP and CAT+BP are more similar to the aminoacid tree (S8 and S9 Figs): in the first tree Ecdysozoa and Lophotrochozoa monophyly are destroyed by the kinirhynchs position only, while in the second one all the major taxa are monophyletic. We conclude that the inference with our mitochondrial nucleotide data is facilitated by using the nonstationary models of sequence evolution, presumably by alleviating the negative effect of the compositional bias.

Mitochondrial gene order and rearrangements
The gene orders of E. svetlanae and P. kielensis mitochondrial PCGs are shown in Fig 5. These two arrangements differ by one transposition between two gene groups: nad1-nad2-nad4l-cox2-atp6-cox3-nad6 and nad3-cytb-nad5. Both kinorhynch gene orders are unique: there are no similar arrangements in the Organelle Resource database not only among the Ecdysozoa (1007 species in the database at the time of comparison), but also among Bilateria (5165 species).
Unlike most bilaterians, E. svetlanae and P. kielensis have gene orders that lack any of the previously described conservative gene blocks [14] (Fig 6). However, E. svetlanae has three pairs of adjacent genes from one conservative gene block-cox2-atp6, nad4-nad5, and rrnL- nad1, while P. kielensis has only one such pair of genes-cox2-atp6. This suggests that E. svetlanae has a more plesiomorphic gene order than P. kielensis. While all of the PCG blocks have been eroded in kinorhynchs, the conservative pairs atp6-cox3 and nad1-nad2, which represent block boundaries, are preserved as common features in kinorhynchs and a majority of taxa within protostomes and deuterostomes (Fig 7). PCG location in the same chain is another widespread, although not universal character. We can assume that some of these features are inherited by kinorhynchs from a common ancestor of Bilateria.
We compared the gene orders and putative ground patterns of genes between several ecdysozoan groups and Deuterostomia as an outgroup (Fig 7). The most conservative gene clusters (cox2-atp6, nad4-nad5, and rrnL-nad1-nad2) are present in ecdysozoans and deuterostomes. As shown in Fig 7, Priapulida has the highest similarity to deuterostomes in their gene order. Some gene clusters in Panarthropoda are positioned in the minor strand, and E. svetlanae has the rrnL-nad1-nad2 cluster, which is shared with deuterostomes but is exclusive among ecdysozoans.
We reconstructed the gene order evolution in protostomian mitochondrial genomes using the TreeREx program and the previously suggested ground patterns of Panarthropoda [15], Lophotrochozoa [16], and Deuterostomia [17] (Fig 8). Deuterostomia was selected as an outgroup because their mitogenomes generally have more conservative gene orders than the mitogenomes of lophotrochozoans. The evolutionary scenario reconstructed by the TreeREx suggests that the ancestors of Scalidophora, Ecdysozoa, and Protostomia all share the same gene order, which coincides with the gene order seen in Priapulida, and proposes this pattern as a plesiomorphic trait in the group. Considering that the priapulid gene order is reconstructed as ancestral for both Scalidophora and Ecdysozoa, the kinorhynch gene order rearrangements are equally parsimonious under the monophyly of Kinorhyncha and Priapulida or under any alternative positions of Kinorhyncha around the base of Ecdysozoa, despite the fact the gene order of Panarthropoda can be converted to E. svetlanae in five steps while the gene order of Priapulida can be converted to E. svetlanae in four.
In the ancestors of some major clades, the genome rearrangements reshuffled the mitochondrial genes of the bilaterian ancestor generating patterns that have been described in terms of conservative gene blocks [14], which presumably correspond to fragments of the mitochondrial genome of the bilaterian ancestor. In other clades, such as kinorhynchs, the rearrangements reshuffle the ancestral genome to a point where the "block boundaries" seem to be more conservative than the "blocks". The presence of taxa with highly altered gene orders in comparison to the presumed ancestral pattern, such as Nematoda [79], Chaetognata [80,81], Acoela [82], and Urochordata [67,83,84] demonstrates the lack of strong selective restrictions on the patterns of the gene order. At the same time, the conservation of mitochondrial gene order patterns in various taxa indicates that the rearrangements are relatively rare events: the rate of rearrangements appears to be comparable to or lower than the rate of formation of the highest level taxa such as classes and phyla. Presumably, the gene order rearrangements are frequently deleterious rather than neutral. The gene arrangements differ greatly in the taxa with long branches in both mitochondrial and nuclear PCGs trees, such as nematodes, chaetognats, tunicates, and kinorhynchs. One of the possible reasons for the high rate fixation of unusual mutations simultaneously in coding sequences and gene pattern is the weakening of the purifying selection caused by the low effective population size [85][86][87][88] during the ancestor's history. Meiobenthic kinorhynchs have high populations and wide geographic distributions, but their patchy habitat may decrease the effective population size. Further studies are required to elucidate the reasons for the increased rate of evolution of mitochondrial sequences and gene orders in kinorhynchs and other taxa with divergent mitochondrial genomes.

Conclusions
The complete mitochondrial genomes of two distant species of Kinorhyncha, Echinoderes svetlanae (Cyclorhagida) and Pycnophyes kielensis (Allomalorhagida), demonstrate similarity in the nucleotide composition, patterns of gene arrangements, and genome architecture. Both mitogenomes have duplicated methionine tRNA genes. The closest relatives of Kinorhyncha within Ecdysozoa are not clearly established by the mitochondrial PCGs phylogeny due to their highly divergent sequences, however the reconstructed scenario of gene order evolution does not contradict to the monophyly of Scalidophora. According to gene order analysis, Priapulida gene arrangement may be ancestral for Scalidophora, Ecdysozoa, and Protostomia.