A total of twelve species from four thoracotreme superfamilies ( Table 1 ) were sampled for sequencing. All specimens were preserved in 95% ethanol following collection and deposited at Marine Biological Museum, Chinese Academy of Sciences, Qingdao, China. Identification was performed morphologically according to the morphological information on crabs of the China Seas (Dai and Yang, 1991). Total genomic DNA was extracted from the samples using the E.Z.N.A^® Tissue DNA kit (OMEGA, Wuhan, China) following the manufacturer’s protocols. The paired-end libraries were constructed using TruSeq™ Nano DNA Sample Prep Kit (Illumina, San Diego, CA, USA) with an insert size of 450 bp. The above libraries were sequenced by an Illumina (San Diego, CA, USA) HiSeq 4000 platform (2 × 150 bp paired-end reads).

The paired-end raw sequences for the twelve samples were trimmed using Trimmomatic 0.39 (Bolger et al., 2014) with the following parameters: ILLUMINACLIP : TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:75. Resulting clean reads ( Table 1 ) were assembled de novo by NOVOPlasty software (Dierckxsens et al., 2016) with default parameters. In the seed extension algorithm achieved via NOVOPlasty, the cox1 or 16S rDNA gene fragments of the species or their closed related species were used as seed sequences ( Supplementary Table 1 ). NOVOPlasty does not try to assemble every read, but will extend the given seed until the circular mitogenome is formed. The mitogenomes will be circularized when the length is in the expected range and both ends overlap by at least 200 bp. If circularized mitogenome was not obtained, MITObim v.1. (Hahn et al., 2013) was used to achieve assemblies using short mtDNA sequences from related species as “baits”. The clean reads can be accessed from the SRA database with the accession numbers SRR18217430-SRR18217435 and SRR18217456-SRR18217461.

MITOS web server (Bernt et al., 2013) was preliminarily used to annotate the protein-coding genes (PCGs), 2 ribosomal RNA (rRNA) and transfer RNA (tRNA) genes under default settings, and the invertebrate genetic code for mitochondria. The PCGs and rRNA genes boundaries were edited manually based on similarity with the homologous genes of other reported thoracotreme mitochondrial gene sequences. The tRNA genes and their secondary structures were further identified with ARWEN 1.2.3.c (Laslett and Canbäck, 2007). Correlations between mitochondrial genome size and length of non-coding regions were measured by the Spearman correlation coefficient carried out with IBM SPSS Statistics, release 19.0.0.1. The complete mitochondrial DNA (mtDNA) sequences can be accessed from the GenBank database with accession numbers OL661260-OL661268, OL579739-OL579740, and OL657229.

The nucleotide composition was calculated in MEGA 6.0 (Tamura et al., 2013), based on the invertebrate mitochondrial genetic code (genetic code = 5). AT and GC skews were determined for complete mitogenomes (plus strand) according to the following formula: AT-skew = (A-T)/(A+T) and GC-skew = (G-C)/(G+C) (Perna and Kocher, 1995). Mitochondrial gene arrangements were compared using the newest version of MitoPhAST v2.0 (Tan et al., 2015) and the same gene orders were clustered into groups. Genomic rearrangement analyses were performed with the CREx program (http://pacosy.informatik.uni-leipzig.de/crex/form) (Bernt et al., 2007). The rearrangement events involve transpositions, inversions, inverse transpositions, and tandem duplication random loss (TDRL) (Boore, 2000; Bernt et al., 2007; Bernt et al., 2008; Bernt and Middendorf, 2011). We include all the mitochondrial genes (PCGs, tRNAs, and rRNAs); thus, 37 genes were considered. We explored the relationship between gene rearrangements and branch length of each taxon. Gene order changes were measured using the breakpoint distance (Blanchette et al., 1999; Bernt et al., 2013). The breakpoint distance was quantified as the minimal number of breakpoints needed to get from the gene order under view to the Brachyura ground pattern. The branch length of each taxon was determined as root-to-leaf distance in the phylogenetic tree. Correlations between different gene orders and branch length were measured by the Spearman correlation coefficient carried out with IBM SPSS Statistics, release 19.0.0.1.

The phylogenetic trees were reconstructed with 75 species, including the newly determined sequences, of which 5 Heterotremata species were used as outgroups ( Supplementary Table 2 ). The relationships were reconstructed based on nucleotide sequences from 13 PCGs, which were aligned separately using MAFFT (Katoh and Standley, 2013) with default settings, then concatenated into a single supermatrix with FASconCAT (Kück and Meusemann, 2010). PartitionFinder 1.1.1 (Lanfear et al., 2012) was used to select the best-fit partitioning schemes and substitution models ( Supplementary Table 3 ). Maximum Likelihood (ML) analysis was conducted with the IQ-TREE web server (Trifinopoulos et al., 2016) with the best-fit partition schemes and models. The node reliability was assessed using 5,000 ultrafast bootstrap replicates (Minh et al., 2013). The Bayesian inference (BI) was conducted using MrBayes 3 (Ronquist and Huelsenbeck, 2003) with partition models. Runs of 10,000,000 generations (Markov chain Monte Carlo, MCMC) were conducted, with a sampling frequency of 1,000 generations to allow adequate time for convergence. In our study, the average standard deviation of split frequencies was less than 0.01 (0.000771), the estimated sample size was >200, and the potential scale reduction factor approached 1.0 already after 10 million generations, and then the run was stopped. All trees generated prior to the achievement of stationarity of the log likelihood values were discarded as burn-in (2,500 trees). The 50% majority rule consensus tree was constructed using the remaining sampled trees to estimate the Bayesian posterior probabilities. The software Tracer v1.6 (Rambaut and Drummond, 2013) was used to check all the parameters.

Our work successfully obtained complete mitogenomes of 12 thoracotreme species from 4 superfamilies, substantially increasing the mitogenome resources for superfamilies, which were poorly represented in previous mitogenomic studies. They contained a standard set of 37 genes, and varied in size from 15,422 [Hapalocarcinus marsupialis Stimpson, 1858 sensu lato (Bahr et al., 2021)] to 16,490 bp [Arcotheres sinensis (Shen, 1932)] ( Table 1 and Supplementary Table 4 ). The average coverage of mitogenomes varied from 694-fold to >3,000-fold coverage. Although the mitogenome sizes of thoracotreme were similar to that of the most metazoans (Gissi et al., 2008), the length of thoracotreme mitogenomes seemed to be family-specific, instead of superfamily-specific. The size of mtDNAs from five ocypodoid families (Camptandriidae, Xenograpsidae, Dotillidae, Mictyridae, and Ocypodidae), five grapsoid families (Grapsidae, Gecarcinidae, Sesarmidae, Plagusiidae, and Varunidae), and a cryptochiroid family (Cryptochiridae) were usually less than 16,000 bp, while the two families, Macrophthalmidae and Xenophthalmidae, from Ocypodoidea, Varunidae from Grapsoidea, and Pinnotheridae from Pinnotheroidea had mtDNAs larger than 16,000 bp ( Figure 2A ), forming their own cluster. Two species, Chasmagnathus convexus (De Haan, 1835 [in De Haan, 1833-1850]) from Varunidae (yellow dot 1 in Figure 2A ) and Macrophthalmus pacificus Dana, 1851 from Macrophthalmidae (yellow dot 2 in Figure 2A ), were 15,107 and 17,226 bp in length, respectively, which were the outliers. The length variations in thoracotreme mitogenomes is mainly due to the length heterogeneity of the non-coding regions ( Supplementary Table 5 ). A significant and strong positive correlation was observed between mitochondrial genome size and length of non-coding regions (Spearman’s rho = 0.859, p < 0.001).

Nucleotide composition of the thoracotreme mitogenomes was biased toward A and T ( Figure 2A and Supplementary Table 5 ). The A+T content ranged from 65.0% in Grapsus tenuicrustatus (Herbst, 1783 [in Herbst, 1782-1790]) to 77.7% in Nanosesarma minutum (de Man, 1887). The values of the AT-skews were mostly negative, indicating more Ts than As, while the GC-skews for all thoracotreme mitogenomes were negative, suggesting that Cs were more abundant than Gs ( Figure 2B ). We found a strong negative relationship between the AT and GC skew (Pearson = −0.409, p < 0.0001) ( Figure 2B ). Two clusters can be found, one with negative GC-skews and AT-skews, the other with negative GC-skews and positive AT-skews. Note that the skewness of the thoracotreme mitogenomes were neither superfamily-specific nor family-specific. Species from four superfamilies were mixed in an irregular skewness.

A total of nine distinct patterns of mitochondrial gene order (MGO) are known for the 70 available thoracotreme mtDNAs ( Figure 3A ). The most widespread gene order pattern is brachyuran gene order (BraGO). For some thoracotreme species, MGO patterns are shared at the family level. For example, the BraGO pattern was shared by the crabs belonging to eight families (Dotillidae, Mictyridae, Ocypodidae, Plagusiidae, Grapsidae, Gecarcinidae, Camptandriidae, and Xenophthalmidae). The mitochondrial DNA of Sesarmidae and Cryptochiridae rearranged with respect to SesGO and CryGO, respectively. The rearrangements of the two families Macrophthalmidae and Varunidae were described by the MaVaGO pattern, whereas for other species, MGO patterns are distinct among species within the same family (e.g., Varunidae and Pinnotheridae) or even genus (e.g., Xenograpsus). The Pinnotheridae species exhibit two gene orders (hereafter named Pin1GO and Pin2GO, respectively). Asthenognathus inaequipes Stimpson, 1858 was first found to have a novel gene order (AinGO) relative to other Varunidae species.

The ancestral gene order reconstruction shows that the BraGO pattern was the ancestral MGO of Thoracotremata, which occurred at the onset of thoracotreme clade ( Figure 3B ). The other Thoracotremata GOs evolved from BraGO can be separated into three groups: (i) high rearranged GO (SesGO), which shares the most common intervals (1,256) with BraGO; (ii) moderate rearranged GOs (CryGO, Pin1GO, and Pin2GO), which share 778, 882, and 882 common intervals with BraGO, respectively; and (iii) low rearranged GOs (XngGO, XteGO, MaVaGO, and AinGO), which share 216 or less common intervals with BraGO ( Figure 3C ). Through mechanism analysis, we found that MaVaGO, AinGO, CryGO, and XngGO ( Figures 4A–D ) could be produced through TDRL events from BraGO. XteGO derived from BraGO through a very complex mechanism requiring two transposition events and two TDRL events ( Figure 4D ). In Pin1GO, Pin2GO, and SesGO, only some tRNA genes have changed their location, implying the transposition events ( Figures 4E, F ).

In the present study, we newly described four MGOs, AinGO, CryGO, Pin1GO, and Pin2GO, for the subsection Thoracotremata. The evolutionary pathways generating the new arrangements are depicted in Figures 4B, C, E . The AinGO has been produced through two TDRL events that involved the genomic portion included between trnT and trnW ( Figure 4B ). In the CryGO pattern, a single TDRL event is necessary to explain the final rearrangement, and the genomic portion involved in the process contains 10 genes ( Figure 4C ). Pin1GO differs from BraGO for the transposition of trnN, trnR, trnP, and trnQ (T1–T4), which were located between trnE and trnH, trnN and trnH, trnS2 and nad1, and rrnS and trnI, instead of its placement downstream to trnA, trnT, and trnI in BraGO, respectively ( Figure 4E ). The transformational pathway producing Pin2GO implies three transpositions (T1–T3). Only trnR, trnP, and trnQ are transposed with respect to BraGO.

Some studies have mentioned that higher variation in mitochondrial gene order may be correlated with accelerated substitution rates, resulting in long branches and problems in nucleotide sequence-based analysis. Strong positive correlations between gene orders and elevated nucleotide substitution rates were found in arthropod and bivalve mitogenomes (Shao et al., 2003; Xu et al., 2006; Plazzi et al., 2016). In contrast to these findings, our Spearman correlation analysis indicated no link between variable gene arrangements (breakpoint distances) and nucleotide substitution rates (branch lengths) in thoracotreme crabs ( Supplementary Table 6 ). The value of Spearman’s rho was 0.15 (p = 0.217). However, the branch lengths of taxa with rearranged GO (SesGO, CryGO, Pin1GO, Pin2GO, XngGO, XteGO, MaVaGO, and AinGO) were significantly higher than that of BraGO (Mann–Whitney U, p = 0.031). This result is consistent with the observation of Tan et al. (2019), who found that the Spearman correlation analysis does not suggest a link between variable MGOs and nucleotide/amino acid substitution rates in Decapoda. Although our correlation analysis indicated no such association for the Thoracotremata, we found that the symbiotic groups, the cryptochiroid and pinnotheroid crabs, display variable MGOs (CryGO, Pin1GO, and Pin2GO). The crabs of Cryptochiroidea and Pinnotheroidea shared a similar lifestyle. Pea crabs from the superfamily Pinnotheroidea are obligate symbionts with some marine invertebrates, such as the mollusks, tubeworms, and echinoderms (Campos, 2016; Theil et al., 2016). Gall crabs (Cryptochiroidea) are coral-dwelling crabs, which are obligate symbionts of living scleractinian coral reefs (Wetzer et al., 2009; van der Meij, 2015; van der Meij et al., 2015). Some studies have suggested possible correlations of rearranged MGOs to the adaptations to specialized lifestyles and extreme ecological niches (also see Tan et al., 2019), such as burrowing crayfish (Gan et al., 2018), freshwater crabs (Tan et al., 2018), and the species adapted to extreme deep-sea environments (Nakajima et al., 2016). This seems not absolute. Some deep-sea vent crab and shrimp species and burrowing axiid mud shrimps exhibited only a few rearranged MGOs (Tan et al., 2019). Therefore, more mitogenomes of representation from major groups with varied lifestyles and extreme ecological niches were needed to clarify exact drivers of mitochondrial gene rearrangements.

Previous mitogenomic thoracotreme phylogenies were limited by including only Grapsoidea and Ocypodoidea taxa in studies (Chen et al., 2018; Wang et al., 2021). In this study, we presented the first expanded Thoracotremata mitochondrial phylogenetic relationship, including all the four superfamilies (Grapsoidea, Ocypodoidea, Pinnotheroidea, and Cryptochiroidea).

The ML phylogeny shared most similarities in their topologies with that obtained from Bayesian inference ( Figure 5 ), specifically in the recovery of Pinnotheroidea as basal clade. However, several nodes were still in contention. In the BI tree, the soldier crab Mictyris longicarpus Latreille, 1806 closely grouped with the ocypodoid family Macrophthalmidae with strong nodal support (BPP = 0.96). However, in the ML tree, M. longicarpus and the grapsoid family Xenograpsidae formed a sister lineage with low support value (BS = 56). The monophyly of the grapsoid family Gecarcinidae was not supported in the BI tree (BPP = 1.00). Gecarcoidea lalandii H. Milne Edwards, 1837 was more closely related to the family Sesarmidae than its Gecarcinidae relatives Cardisoma carnifex (Herbst, 1796 [in Herbst, 1791-1796]) and Cardisoma armatum (Herklots, 1851), while the three Gecarcinidae species formed a monophyletic group in the ML tree, but with low support value (BS = 65). Although there is some discrepancy in the internal relationships among Varunidae between the ML and BI trees, the position of the newly presented grapsid species Tritodynamia horvathi Nobili, 1905 was identical. It is apparent that T. horvathi, a member of Macrophthalmidae, was shown to be more closely related to the varunid crab A. inaequipes (BS = 100, BPP = 1.00), than to other macrophthalmid crabs. This seems not surprising. The genera Tritodynamia and Asthenognathus were once placed in the family Pinnotheridae based on the morphological characteristics (Schmitt et al., 1973). Števčić (2005) raised them to family level and placed Tritodynamia in the Macrophthalmidae as a separate subfamily, the Tritodynamiinae. Cuesta et al. (2005) suggested a close relationship with the Varunidae after looking at the molecular data for several species. Later, Ng et al. (2008) placed the Tritodynamia and Asthenognathus into the family Macrophthalmidae and Varunidae, respectively. Our results questioned the validity of the subfamily Tritodynamiinae and support the inclusion of Tritodynamica in the family Varunidae.

The crabs of Pinnotheroidea are usually symbiotic with mollusks, tubeworms, echinoids, and other invertebrates. Its phylogeny was poorly understood, with the monophyly of this group still under question (Theil et al., 2016). In our study, the Pinnotheroidea, including only one family Pinnotheridae, formed a monophyletic clade. However, previous studies have pointed out the non-monophyly of the pinnotheroid crabs (Palacios-Theil et al., 2009; Theil et al., 2016). This discrepancy may due to the limited pinnotheroid sample in our study and the heterogeneity of data. Pinnotheridae was the earliest branching in Thoracotremata, located at the base of the phylogenetic tree, which was consistent with the result of Schubart et al. (2006). In the work of Schubart et al. (2006) (using 16S and 12S), only one Pinnotheridae species, Pinnotheres pisum (Linnaeus, 1767), was sampled, which was shown to hold a basal position compared to the clade of “Grapsoidea & Ocypodoidea”. However, this presents an opposite view compared with some other previous studies (Wetzer et al., 2009; Tsang et al., 2014; van der Meij and Schubart, 2014; Tsang et al., 2018). Using a fragment of the 16S gene, Wetzer et al. (2009) found that the three Pinnotheridae species constituted a monophyletic clade, nested within a clade that includes the ocypodoid families Mictyridae, Camptandridae, and Ocypodidae. An analysis of more than 140 brachyuran species (only two Pinnotheridae species) using six nuclear protein-coding genes and two mitochondrial rRNA genes suggested that the two pea crabs clustered with Dotillidae from Ocypodoidea, forming a basal lineage of Thoracotremata (Tsang et al., 2014). Using molecular data from three markers (mitochondrial 12S and 16S rRNAs, and nuclear Histone H3), Tsang et al. (2018) revealed Pinnotheridae to be monophyletic, but widely distant from its pinnotheroid relatives (Aphanodactylidae) and instead clustered with Ocypodidae (Ocypodoidea), forming two basal thoracotreme clades. Based on the analysis of the 16S gene of 82 thoracotreme crabs, van der Meij and Schubart (2014) found that the monophyletic Pinnotheridae clade was closely related to the family Ocypodidae from Ocypodoidea.

Cryptochiroidea was represented in the present study by the species H. marsupialis from family Cryptochiridae, which was grouped with an ocypodoid lineage (Dotillidae/Xenophthalmidae/Camptandriidae). Previous studies have suggested that Cryptochiridae was a sister lineage to Xenograpsidae (Grapsoidea) (Tsang et al., 2014; Tsang et al., 2018). Close phylogenetic relationships between the Cryptochiridae and Grapsidae were proposed by Wetzer et al. (2009), who, based on 16SmtDNA, considered H. marsupialis as a “highly modified Grapsidae”, because H. marsupialis evolved within the family Grapsidae (Grapsoidea), while van der Meij and Schubart (2014) recommended that gall crabs should not be considered “highly modified Grapsidae”, but an independent lineage of grapsid crabs. As the Cryptochiroidea was represented by single species in the present study and the ML analysis only has 77 as a support value in this node, no overall conclusions about its phylogenetic position in the Thoracotremata should be made.

Although the results on the polyphyly of Grapsoidea and Ocypodoidea are consistent between the present and previous studies (Tsang et al., 2014; Basso et al., 2017; Tsang et al., 2018; Chen et al., 2019; Kobayashi et al., 2020; Zhang et al., 2020; Wang et al., 2021), the internal relationships of the families within these two groups were different. Our research provides some new insights into the phylogenetic relationships of some families with a wider taxonomic coverage. For example, Dotillidae (Ocypodoidea) and the ocypodoid families Xenophthalmidae and Camptandriidae formed a sister lineage, and Sesarmidae (Grapsoidea) grouped with its grapsoid relatives Xenograpsidae and Gecarcinidae. In contrast, in the studies of Basso et al. (2017); Tan et al. (2018), and Wang et al. (2021), which were also based on mitogenomes, Dotillidae (Ocypodoidea) and Sesarmidae (Grapsoidea) formed a sister group. Certainly, when the material from more species were obtained and the mitogenomic data were substantial, the relationships of families from Grapsoidea and Ocypodoidea can be well interpreted.