We used individual specimens identified as the following five species collected around Japan: (1) Apseudes nipponicus Shiino, 1937 (Apseudoidea: Apseudidae), collected from an open-air experimental aquarium in the Shimoda Marine Research Center (University of Tsukuba), Shizuoka, Japan (for details, see Kakui et al. 2017). (2) Apseudes ranma Matsushima and Kakui, 2024, collected from an aquarium tank in the Port of Nagoya Public Aquarium, Aichi, Japan (for details, see Kakui and Hiruta 2013); this specimen was a descendant of a single hermaphrodite isolated in March 2019 (for details, see Matsushima and Kakui 2024). (3) Carpoapseudes spinigena Bamber, 2007 (Apseudoidea: Apseudidae), collected at 1028–1075 m depth off the southeastern coast of Hokkaido, Japan (for details, see Kakui et al. 2020). (4) Neotanais sp. (Neotanaoidea: Neotanaidae), collected at 1018–1042 m depth, East China Sea (29.322500°N, 127.622900°E to 29.341983°N, 127.632633°E), Japan, on 21 November 2009 during Cruise N295 of TRV Nagasaki-maru (Nagasaki University). (5) Sinelobus kisui Hirano and Kakui, 2022 (Tanaidoidea: Tanaididae), collected from a brackish stream at Hagi, Yamaguchi, Japan (for details, see Hirano and Kakui 2022). All specimens were fixed in 80–99% ethanol and preserved in 99% ethanol.

Total DNA was extracted from the whole body of each individual by using a NucleoSpin Tissue XS Kit (Macherey–Nagel, Germany). As seed sequences, previously determined partial COI sequences were used for A. ranma (LC791470; Matsushima and Kakui 2024) and C. spinigena (LC545558; Kakui et al. 2020), and newly determined sequences for A. nipponicus, Neotanais sp., and S. kisui. Primers LCO1490 and HCO2198 (Folmer et al. 1994) were used for PCR and cycle sequencing. PCR amplification conditions for COI with TaKaRa Ex Taq DNA polymerase (TaKaRa Bio, Japan) were 94°C for 1 min; 35 cycles of 98°C for 10 s, 42–50°C for 30 s, and 72°C for 1 min; and 72°C for 2 min. Nucleotide sequences were determined with a BigDye Terminator Kit v3.1 and a 3730 DNA Analyzer (Life Technologies, USA).

For S. kisui, we designed primer pair Sk_COI_LF (5’-TCCTATAGGTGGGGGTGATCCAA TTTTATTTCAGC-3’) and Sk_COI_LR (5’-CTATATGCAGCCAAGGATAGCGGTGGGTAA ATTG-3’), by using Primer3Plus (Untergasser et al. 2007) and based on the determined COI sequence, to amplify the nearly complete mitogenome by long PCR. The long-PCR amplification conditions with KOD ONE PCR Master Mix (Toyobo, Japan) were 45 cycles of 98°C for 10 s, 68°C for 5 s, and 68°C for 3 min. Amplicons were purified with the Extractor PCR & Gel Clean Up Kit (Toyobo).

We used the purified amplicon of S. kisui and DNA extracted from each of the other four species in the whole-genome shotgun sequencing described below.

Whole-genome shotgun sequencing was performed by using 2 × 200 bp (A. nipponicus) and 2 × 150 bp (A. ranma, C. spinigena, and Neotanais sp.) paired-end reads on the DNBSEQ-G400 platform (MGI Tech, China) at Bioengineering Lab Co., Japan. Sequencing for S. kisui was performed by using 2 × 150 bp paired-end reads on the DNBSEQ-G400RS platform at Genome Read Inc., Japan. A total of 32,444,610 reads (2 × 4.8 Gbp) were obtained for A. nipponicus; 17,433,684 (2 × 3.4 Gbp) for A. ranma; 32,127,279 (2 × 4.7 Gbp) for C. spinigena; 34,957,034 (2 × 5.2 Gbp) for Neotanais sp.; and 6,228,917 (2 × 1.87 Gbp) for S. kisui.

The quality of the raw FASTQ files was initially assessed by using FastQC v0.12.0 (Andrews 2010) with default parameters. Adapter sequences were subsequently trimmed from the raw reads by using fastp v1.0.1 (Chen et al. 2018). All processed reads were assembled into target genomes by using GetOrganelle v1.7.7.1 (Jin et al. 2020) on the Galaxy web platform (The Galaxy Community 2024), with a COI sequence acting as a seed for the assembly process. For S. kisui, the contig obtained and the partial COI sequence were concatenated in MEGA11 (Tamura et al. 2021). To verify the assembled genomes, the paired-end reads were mapped back to their respective assembled genomes using BWA-MEM v0.7.19 (Li 2013). The reads aligning to target genes were extracted into FASTQ format with SAMtools v1.9 (Danecek et al. 2021), and the final mapped reads were visualized and interpreted by using Integrative Genomics Viewer (IGV) v2.19.5 (Thorvaldsdóttir et al. 2013).

Based on the assembled genomes of A. nipponicus, A. ranma, and Neotanais sp., we designed the following primers by using Primer3Plus and OlvTools (https://olvtools.com/tmvalue) to determine the sequences that remained unassembled after the procedure described above: ApnF (5’-GTATCCCTTTCCCGCCCTGCATCAC-3’) and ApnR (5’-CACACATTTCGTGCTCCCCTTTGATC-3’) for A. nipponicus; AprF (5’-CAGGATACACTTTTCGTGCTCCTCTTTG-3’) and AprR (5’-CCCCAGCCTCAGG TTTTAAATGAAG-3’) for A. ranma; and NeoF (5’-GTCGTGTTTTAGTTGCTTGGTCTGCACTTC-3’) and NeoR (5’-CCACGCAACACTATAACTGCCCAAATCATC-3’) for Neotanais sp. PCR amplification conditions with KOD ONE PCR Master Mix were as follows: for A. nipponicus, 45 cycles of 98°C for 10 s, 60°C for 5 s, and 68°C for 30 s; for A. ranma, 45 cycles of 98°C for 10 s, 55°C for 5 s, and 68°C for 30 s; and for Neotanais sp., 45 cycles of 98°C for 10 s and 68°C for 5 s. The amplified sequences were determined with a BigDye Terminator Kit v3.1 and a 3730 DNA Analyzer.

In addition to annotating the sequences of the five species determined in this study, we re-annotated the previously published mitogenome of Arctotanais alascensis (LC597489; Kakui and Kano 2021), which had been annotated by using the MITOS webserver (Bernt et al. 2013). Initial annotation was performed with MITOS2 (Donath et al. 2019) on the Galaxy web platform (The Galaxy Community 2024). Protein-coding genes (PCGs) not identified by MITOS2 were searched with the Open Reading Frame Finder (ORF Finder) provided by the National Center for Biotechnology Information (https://www.ncbi.nlm.nih.gov/orffinder). All the PCG candidate regions were translated by using the invertebrate mitochondrial genetic code, and their gene identities were confirmed with reciprocal BLASTP (Altschul et al. 1990) and by searching against the Conserved Domain Database (Lu et al. 2020). Based on these alignments, start and stop codons were identified manually.

Genes that remained unidentified were identified manually through comparative alignments. The atp8 gene in A. ranma was determined by aligning the assembled genome with that of A. nipponicus, using MUSCLE (Edgar 2004) in MEGA11. After identifying putative rrnL sequences (1092 bp for A. nipponicus and 598 bp for A. ranma) by using Barrnap v0.9 (Seemann 2018), the definitive sequences were determined by aligning them with the rrnL of C. spinigena, using MUSCLE. Additional tRNA searches were conducted by using ARWEN v1.2.3 (Laslett and Canbäck 2008) and tRNAscan-SE v2.0 (Chan et al. 2021). ARWEN successfully detected trnD and trnI in A. alascensis and trnR in Neotanais sp., none of which had been identified by MITOS2. Manual curation and editing of the mitogenomes were performed by using UGENE v52.1 (Okonechnikov et al. 2012).

The annotated mitochondrial genome sequences determined in this study have been deposited in the DDBJ/EMBL/GenBank databases under accession numbers LC923624–LC923628.

The concatenated dataset of 13 PCGs and two rRNAs included sequences from seven tanaidacean species, 71 non-tanaidacean peracarid species (41 amphipods, one cumacean, 25 isopods, one lophogastrid, two mysids, and one stygiomysid), and two non-peracarid species (one decapod and one euphausiacean) (Dataset 1; Table S1). Nucleotide sequences for the PCGs from all species except Pseudotanais sp. were translated into amino acid sequences according to the invertebrate mitochondrial genetic code. For Pseudotanais sp., translation with that code was not possible (the codon TAA translates as tyrosine rather than a stop codon; Jakiel et al. 2025). We thus used the amino acid sequences provided in the GenBank CDS features. The amino acid sequences were aligned by using the L-INS-i algorithm (Katoh et al. 2005) implemented in MAFFT v7.526 (Katoh and Standley 2013) and subsequently back-translated to nucleotide sequences by using PAL2NAL v14 (Suyama et al. 2006). For Pseudotanais sp., the back-translated aligned nucleotide sequences were replaced with the original sequences, retaining the gap positions in the back-translated sequences. The rRNAs were aligned by using the Q-INS-i algorithm (Katoh and Toh 2008) in MAFFT v7.526. Poorly aligned regions were removed by using trimAl (Capella-Gutiérrez et al. 2009). After the final alignment, Dataset 1 was 11,650 positions long. The optimal substitution model for 41 partitions (three codon positions for the 13 PCGs + two rRNAs) for Dataset 1 was determined by using ModelFinder (Kalyaanamoorthy et al. 2017) with the MFP+MERGE option, based on the Bayesian information criterion (BIC), resulting in the consolidation of these 41 partitions into 22 partitions (File S1). A maximum likelihood (ML) analysis was performed by using IQ-TREE v3.0.1 (Wong et al. 2026) under the bnni option (Hoang et al. 2018). Branch support values were obtained from an ultrafast bootstrap analysis of 10,000 pseudoreplicates (UFBoot; Hoang et al. 2018). Branches with UFBoot ≥ 95% were regarded as strongly supported. The resulting phylogenetic tree was visualized by using FigTree v1.4.4 (Rambaut 2018) (File S1).

The initial analysis with IQ-TREE detected extensive substitutional saturation in Dataset 1, particularly at third-codon positions in the PCGs. We thus do not present the results from Dataset 1, but instead analyzed two new datasets derived from Dataset 1. Dataset 2 included amino-acid sequences for the 13 PCGs and nucleotide sequences for the two rRNAs; that is, the aligned amino acid sequences were not back-translated to nucleotide sequences. Dataset 3 included nucleotide sequences for the first and second codon positions of the 13 PCGs and for the two rRNAs; that is, the nucleotides in the third codon position were removed from Dataset 1.

Dataset 2 comprised 4518 positions. The initial 15 partitions were merged into 13 optimal partitions by using ModelFinder with the MFP+MERGE option based on the BIC. Table S2 presents the final partitions and optimized substitution models. Dataset 3 comprised 8182 positions. The initial 28 partitions were merged into 17 optimal partitions by using ModelFinder with the MFP+MERGE option based on the BIC. Table S3 presents the final partitions and optimized substitution models. Methods for the ML analyses, branch support estimation, and drawing of the trees were as described above for Dataset 1.

To assess whether long-branch attraction (LBA) artifacts (Bergsten 2005) affected the relationships within Tanaidomorpha, we performed parametric simulations using AliSim (Ly-Trong et al. 2022) in IQ-TREE. We simulated 100 artificial alignments under the assumption that the true evolutionary process was represented by the partitioning scheme, optimized substitution models, and parameter estimates obtained from Dataset 2 (comprising amino acid sequences of the 13 PCGs and nucleotide sequences of the two rRNAs). The simulations were performed using an ML tree constrained to recover Tanaidoidea and Neotanaoidea as sister groups. For each simulated dataset, we reconstructed an ML tree using IQ-TREE, applying these exact empirical settings of Dataset 2. Finally, the 100 resulting simulated ML trees were mapped onto the unconstrained ML tree to calculate the parametric bootstrap support values (PB) for each branch.

All sequence alignments used in this study have been deposited in the Figshare repository (Kakui 2026).

The assembled genomes for Apseudes nipponicus and Apseudes ranma were single contigs, 14,717 bp and 14,599 bp long, respectively (Fig. 1A, B), each containing an unassembled region. Ten contigs (197–15,857 bp) were obtained for Carpoapseudes spinigena (Fig. 2A), with the longest contig containing PCGs, rRNAs, and tRNAs (Fig. 2B). The assembled genome for Neotanais sp. was recovered as a single contig 14,746 bp long (Fig. 3A) containing an unassembled region. The complete circular mitochondrial genome of Sinelobus kisui was determined, which is 14,126 bp long (Fig. 3B).

For A. nipponicus and A. ranma, multiple bands were observed in electrophoresis of the PCR products amplified with specific primers, preventing identification of the target sequences. Cycle sequencing was thus not carried out for these bands. In Neotanais sp., a single band of approximately 900 bp was obtained. Sequencing of this product successfully determined only a ~600-bp portion, which completely overlapped with and was identical to the terminal region of the 14,746-bp contig mentioned above, including the primer region.

Nucleotide composition varied markedly among the species (Figs 1, 2, 3), with the following GC contents: 42.60% in A. nipponicus, 40.77% in A. ranma, 34.33% in C. spinigena, 43.30% in Neotanais sp., and 20.64% in S. kisui.

In the five species, 13 PCGs and two rRNAs were identified. The complete set of 22 tRNAs was not identified, with the following being absent: trnI and trnM in both A. nipponicus (Fig. 1A) and A. ranma (Fig. 1B); trnI, trnM, trnN, trnS1, and trnV in C. spinigena (Fig. 2B); trnI in Neotanais sp. (Fig. 3A); and trnD in S. kisui (Fig. 3B).

Our re-annotation of the mitogenome of Arctotanais alascensis identified trnI, which was not identified in Kakui and Kano (2021). The rrnL gene was 455 bp longer than annotated in Kakui and Kano (2021), and the extended portion of the gene included the whole length of the putative control region suggested by Kakui and Kano (2021).

Figure 4 shows the gene order for seven tanaidacean mitogenomes and the hypothetical ground pattern for Pancrustacea (Boore et al. 1998). Gene order was identical in the two species of Apseudes. The PCG order was the same between Carpoapseudes and Apseudes, but there were several tRNA translocations. The three apseudoid mitogenomes shared three specific groupings, trnL1 + rrnL, nad6–trnT, and trnL2–nad1. No shared gene clusters were detected among the four tanaidomorph species examined, nor between the two tanaidoid species.

While all tanaidaceans studied exhibited gene arrangements distinct from the pancrustacean ground pattern, several ancestral clusters were observed (Fig. 4). These are [cox2 + trnK], [atp6 + cox3], [trnG + nad3 + trnA + trnR], [trnF + nad5], and [trnY + cox1] in the three apseudoids; [atp8 + atp6 + cox3 + trnG + nad3 + trnA], [nad4 + nad4l], [trnM + nad2], and [trnC + trnY + cox1 + trnL2 + cox2 + trnK] in Neotanais sp.; [atp8 + atp6] in Pseudotanais sp.; and [nad3 + trnA] in the two tanaidoids.

The ML trees based on Datasets 2 (Fig. 5) and 3 (Fig. 6) strongly support the monophyly of each of the peracarid orders that include two or more species in our analyses (UFBoot ≥ 98%). The interordinal relationships derived from the two datasets were identical except for the position of the “mysidacean” clade (= Mysida + Lophogastrida + Stygiomysida; cf. Meland et al. 2015). Branch support values also differed slightly between the two datasets. The mysidacean clade was sister to a strongly supported Mancoida clade (Cumacea + Isopoda + Tanaidacea) in Dataset 2 (Fig. 5) with weak support (UFBoot = 75%), but sister to the Amphipoda clade in Dataset 3 (Fig. 6) with even weaker support (UFBoot = 68%). The monophyly of Mancoida was supported with full or high support (UFBoot = 100% for Dataset 2; 98% for Dataset 3); the Isopoda + Tanaidacea clade was weakly supported in Dataset 2 (UFBoot = 64%) but moderately in Dataset 3 (UFBoot = 88%). In Dataset 1 (File S1), relationships among the orders were largely unresolved except for those among the three mysidacean orders; the monophyly of Tanaidacea received only weak support; and the mysidacean clade was recovered outside a weakly supported clade comprising the other peracarids.

The topology within Tanaidacea is identical between the two trees, with slight differences in branch support. The fully supported Apseudoidea clade was sister to the Tanaidomorpha clade (Tanaidoidea + Paratanaoidea + Neotanaoidea) with strong to moderate support (UFBoot = 99% for Dataset 2; 92% for Dataset 3). Within the latter clade, the neotanaoid Neotanais sp. lay outside a Paratanaoidea + Tanaidoidea clade, with weak (UFBoot = 76% for Dataset 2) or moderate (UFBoot = 92% for Dataset 3) support. The two tanaidoids formed a clade with high or low support (UFBoot = 97% for Dataset 2; 75% for Dataset 3). Our parametric simulations (Fig. S1) showed that a Paratanaoidea + Tanaidoidea clade was never reproduced in the simulated datasets (PB = 0%).

The branches leading to Pseudotanais sp. and S. kisui were markedly longer than those leading to the other five tanaidacean species (Figs 5, 6).

We determined the complete, circular mitogenome for S. kisui and partial mitogenomes containing the complete set of PCGs and rRNAs and most tRNAs for the other four species. Our attempt to Sanger-sequence the unassembled regions in three species failed. Based on our observation in Neotanais sp. that the PCR amplicon sequence was identical to the terminal region of the long contig, the unassembled regions are suggested to be composed of long repeat sequences. To determine the complete mitogenomes of the species for which we were unable to obtain circular mitogenomes, long-read sequencing technologies may be necessary.

Our (re-)annotation showed that no typical control region (CR) was present in the circular mitogenomes of A. alascensis and S. kisui, and that the region previously annotated as the CR in A. alascensis corresponds to part of rrnL (Fig. 3B). The non-coding regions in these two species were less than 215 bp long. As suggested for some other crustaceans (Pons et al. 2014), their CR may be extremely short. We also failed to identify the CR in the partial mitogenomes of A. nipponicus, A. ranma, C. spinigena, and Neotanais sp. Their CRs may lie in the unassembled region. Although Jakiel et al. (2025) did not elaborate, the mitogenome of their Pseudotanais sp. has a long non-coding region that potentially contained a CR.

We found that mitogenome gene order in Tanaidacea varies markedly among superfamilies and is highly divergent from the putative pancrustacean ground pattern. In addition, gene order stability varied considerably. Gene order was highly conserved in the two studied genera in Apseudoidea, whereas extensive rearrangements were evident in the two studied genera in Tanaidoidea.

Pseudotanais sp., A. alascensis, and S. kisui shared only one short cluster with the pancrustacean ground pattern, indicating a highly divergent mitogenome gene order. These three species, especially Pseudotanais sp. and S. kisui, exhibited much longer branch lengths than other tanaidaceans in our tree, which suggests that a marked acceleration in molecular evolutionary rates may have occurred in their most recent common ancestor. Previous studies have reported that higher substitution rates are correlated with elevated rates of gene rearrangements (Shao et al. 2003; Su et al. 2025), and this could account for the highly divergent gene order observed in these three species.

Höpel et al. (2022) reported that mitogenomes in representatives of three mysidacean orders (Lophogastrida, Stygiomysida, and Mysida) lack shared derived gene orders and exhibit high levels of rearrangements. Similarly, in Tanaidacea, the fluidity of gene arrangement and changes in tRNA secondary structures are likely elevated. This high plasticity may have contributed to the difficulty of annotation and previous failures to detect trnI. While such extensive rearrangements are generally explained by the tandem duplication and random loss (TDRL) model, more complex rearrangements, such as those observed in Pseudotanais sp. and S. kisui, may stem from unknown transposition mechanisms involving recombination and double-stranded break repair (Luo et al. 2015).

In the trees from Datasets 2 and 3, Tanaidacea was embedded in Mancoida as a strongly or fully supported clade, along with a cumacean species and a fully supported Isopoda clade. The sister taxon to Tanaidacea was not well resolved; an Isopoda + Tanaidacea clade appeared in both trees, but with only weak to moderate support. This sister relationship was recovered in a recent large-scale phylogenomic analysis (Iwasa-Arai et al. 2026), but not in other studies, which placed Tanaidacea as sister to Isopoda + Cumacea (Bernot et al. 2023; Barta et al. 2025).

The relationships we detected among tanaidacean superfamilies differ from those in three previous studies (Kakui et al. 2011, 2021; Araújo-Silva 2016). In our trees, Tanaidoidea is sister to Paratanaoidea, with only weak to moderate support, whereas the three previous studies recovered a Neotanaoidea + Tanaidoidea clade with high to full support. Our support values are insufficient to confidently overturn the robust topologies of previous studies, so this discrepancy requires careful interpretation. Our parametric simulations provide little evidence that the recovered Paratanaoidea + Tanaidoidea clade is an artifact of LBA. However, they cannot exclude the possibility that LBA arising from heterotachy – a phenomenon in which functional constraints on individual sites change over time, resulting in shifts in site-specific evolutionary rates (Lopez et al. 2002; Kolaczkowski and Thornton 2004) – contributed to the inferred relationship. Alternatively, the recovered relationship may reflect mito–nuclear discordance in the evolutionary history of Tanaidacea, a phenomenon reported in many other taxa, including algae (e.g., Kao et al. 2022), ants (e.g., Rahman et al. 2021), bears (e.g., Lammers et al. 2017), and fish (e.g., Baraf et al. 2025).

The long branches in Pseudotanais sp. and S. kisui might have been caused by a unique genetic code and a low-salinity habitat, respectively. In the mitogenome of Pseudotanais sp., TAA – normally a stop codon in the invertebrate mitochondrial genetic code – translates as tyrosine (Jakiel et al. 2025), which has not been observed in other peracarids. Sinelobus kisui is a brackish water species. Conrad et al. (2021) observed unique gene arrangements in fiddler crabs adapted to brackish environments, implying that habitat disparities exert specific selective pressures on the mitochondrial genome. Similar selective pressures may also have contributed to the long branch of S. kisui in the phylogenetic analysis.

This is the first phylomitogenomic analysis incorporating all of the four superfamilies in Tanaidacea. However, the taxonomic coverage is still highly limited: our dataset contains only seven of around 1600 known species; 14 apseudoid families and 23 paratanaoid families are not represented. It is unclear whether the high evolutionary rate observed in Pseudotanais sp., evident as a remarkably long branch, is common in Paratanaoidea. Aguinaldo et al. (1997) clarified the placement of Nematoda – previously difficult to resolve due to rapid evolution – by restricting their analysis to slowly evolving lineages. Likewise, adding more slowly evolving paratanaoid taxa may alter the inferred topology. Expanded mitogenome sequencing and phylogenetic analyses across a broader taxonomic range are needed to test whether the observed contradiction with nuclear genomic data represents genuine mito–nuclear discordance or LBA artifacts.