Three species of panopeid mud crabs were collected for this study. First, an individual Eurypanopeus depressus was sampled on December 1, 2018 from Hoop Pole Creek, a polyhaline site located in Atlantic beach, North Carolina (NC), USA. The individual was hand-collected at low tide from an intertidal oyster reef and then brought back to the lab for dissection. Second, an individual Panopeus herbstii was sampled on August 12, 2019 from Middle Marsh (Beaufort, NC), another polyhaline site, using a passive sampler attached to a wooden stake that had been driven into the sediment. The sampler design is a small plastic milk crate (19×22×16 cm) filled with autoclaved oyster shell (Roche 2007). Third, an individual Rhithropanopeus harrisii was sampled on February 5, 2020 from Mallard Creek (Washington, NC), a mesohaline site, using the same passive sampling design as above, but this time attached to a small fishing dock. Crabs were brought back to the lab and anesthetized prior to dissection in a –20°C freezer. Dissections for all three species were carried out using a sterilized razor blade, and part of the hepatopancreas and gills were removed and placed into separate tubes for later DNA extraction.

The DNA extractions were conducted on the hepatopancreas and gill tissue of E. depressus, P. herbstii, and R. harrisii using a Zymo DNA extraction kit, according to manufacturer’s protocols. The DNA samples were shipped on dry ice to Novogene, California, who conducted library preparation using the NEBNext Ultra DNA Library Prep Kit. The library was loaded on to a NovaSeq 6000 (Illumina) system using the 150 bp NovaSeq 600 SP reagent kit (300 cycles) for paired end metagenomic sequencing for each individual sample. The resulting data were delivered to the University of Florida for bioinformatic analysis. The data were quality checked and trimmed using Trimmomatic v.0.36 (Bolger 2014) using default parameters. The paired and unpaired reads were assembled using SPAdes v.3.13.0 (Bankevich et al. 2012) with default parameters and k-mer lengths: 21, 33, 55, 77 and 99. The resulting datasets provided a series of contigs that were compared to the NCBI nr database using BLASTx. The mitochondrial genomes of E. depressus (574.088X coverage), P. herbstii (122.084X coverage) and R. harrisii (400.520X coverage) were each identified and circularized. Confirmation of their sequence coverage was conducted using CLC genomics workbench v.12.

The circularized mitogenomes were annotated using MITOS (Bernt 2013). Using the MITOS output, the location of the cox1 gene was determined and the sequences were re-annotated with the cox1 gene at the start of the genome. The putative amino acid and rRNA sequences determined by MITOS were checked using BLASTn and BLASTp (Tables 1–3). The completed genomes were then annotated graphically using Circa (http://omgenomics.com/circa). The genomes are deposited in GenBank under accession numbers MN399962 (E. depressus), MT024989 (R. harrisii), and MT024990 (P. herbstii).

There were 112 brachyuran mitogenomes (see Supplementary material 1) obtained from the GenBank database (NCBI) for phylogenetic comparison using the Brachyura taxonomic ID (txid6752) and filtering results to yield DNA sequences of 10,000–25,000 bp (search date: January 2020). The amino acid and nucleotide sequences were retrieved and annotated from these genomes (13 and 15 sequences, respectively) using the Mitophast pipeline (Tan et al. 2015) which downloads each gene in the mitochondrial genome as separate files. The amino acid and nucleotide sequences files were then aligned individually using MAFFT in Geneious (v10.0.2), trimmed to the smallest sequence and concatenated using Geneious. Phylogenetic analyses were conducted in IQtree (Trifinopoulos 2016), which computed the most appropriate evolutionary model (mtMet+F+I+G4) according to BIC for both the amino acid sequences and the nucleotide sequences. A maximum likelihood tree using the amino acid sequences was created using 1000 bootstrap replicates and an SH-aLRT branch test (Guindon 2010) over 3733 positions; the tree had a log score of –157295.9511. A maximum likelihood tree using the nucleotide sequences was created using 1000 bootstrap replicates and an SH-aLRT branch test over 13790 positions; the tree had a log score of –598390.2632. The resulting trees were annotated using FigTree (http://tree.bio.ed.ac.uk/software/figtree). Both trees were rooted with the mitogenomes of Coenobita brevimanus (KY352233), C. rugosus (KY352235) and C. perlatus (KY352234).

All previously reported gene orders for the Brachyura were annotated according to Basso et al. (2017) and Tan et al. (2018). Pairwise comparisons of the gene orders were performed using CREx software (Bernt et al. 2007) at common intervals. The nomenclature for the gene orders follows Basso et al. (2017). MITOS was used to determine the putative location of the control region (CoRe) for E. depressus, P. herbstii, R. harrisii, Etisus anaglyptus, Leptodius sanguineus, Atergatis floridus, and A. integerrimus through manual examination of the start and stop codons of the open reading frames to look for intergenic spacers. The CoRe for the crabs in the family Potamidae were obtained from Genbank. The CoRes for Echinoecus nipponicus and Pilumnus vespertilio were determined using MITOS following the same method as for the Panopeidae.

The mitochondrial genomes of the panopeid crabs used in this study were closed circular molecules containing 13 protein coding genes, 22 tRNA genes, 2 rRNA genes, and a single control region (CoRe) (Fig. 1). The E. depressus mitogenome was 15,854 bp in length. The P. herbstii mitogenome was 15,812 bp in length. The R. harrisii mitogenome was 15,967 bp in length. As with most brachyurans, the rrnL (16S) and rrnS (12S) genes were located on the negative strand, as are the nad5, nad4, nad4L, nad1 and 8 tRNA genes (Table 1).

The nucleotide composition of the complete E. depressus mitochondrial genome was as follows: A=5442 (34.32%), T=5509 (34.75%), G=1652 (10.42%), C=3251 (20.51%). The A+T and G+C contents were 69.07% and 30.93%, respectively. The protein coding regions include 7 NADH dehydrogenases (nad1–nad6 and nad4L), three cytochrome c oxidases (cox1–cox3), 2 ATPases (atp6 and atp8) and 1 cytochrome b (cob) and account for 10,838 bp of the mitogenome. The 22 rRNA genes present in the mitogenome range in size from 62 (trnD)–71 (trnL1) bp in length, and the ribosomal RNA genes rrnL (16S) and rrnS (12S) have a length of 1393 bp and 817 bp, respectively. The 13 protein coding genes and majority of the ncRNA sequences showed similarity among the panopeid crabs used in this study (Table 1).

The nucleotide composition of the complete P. herbstii mitochondrial genome was as follows: A=5520 (34.91%), T=5687 (35.97%), G=1627 (10.29%), C=2980 (18.85%). The A+T and the G+C contents were 70.87% and 29.13%, respectively. The protein coding region contains 7 NADH dehydrogenases (nad1–nad6 and nad4L), three cytochrome c oxidases (cox1–cox3), 2 ATPases (atp6 and atp8) and 1 cytochrome b (cob) and accounts for 10947 bp of the mitogenome of P. herbstii. The 22 rRNAs present in the mitogenome range in size from 63 (trnD, trnA, trnR) – 69 (trnL1, trnQ, trnM) bp in length, and the ribosomal RNA genes rrnL (16S) and rrnS (12S) have a length of 1392 bp and 820 bp, respectively. All 13 protein coding genes showed high similarity to the panopeid crabs used in this study. The ncRNAs all showed similarity to decapod crustaceans with the majority having high similarity with the Panopeidae (Table 2).

The nucleotide composition of the complete R. harrisii mitochondrial genome was as follows: A=5595 (34.21%), T=5873 (37.00%), G=1556 (9.82%), C=2866 (17.99%). The A+T and the G+C contents were 72.20% and 27.080%, respectively. The protein coding region contains 7 NADH dehydrogenases (nad1–nad6 and nad4L), three cytochrome c oxidases (cox1–cox3), 2 ATPases (atp6 and atp8) and 1 cytochrome b (cob) and account for 10,848 bp of the mitogenome. The 22 rRNAs present in the mitogenome range in size from 63 (trnD, trnG, trnA) – 69 (trnL1, trnQ) bp in length, and the ribosomal RNA genes rrnL (16S) and rrnS (12S) have a length of 1400 bp and 817 bp, respectively. The 13 protein coding genes showed high similarity with the panopeid crabs used in this study (Table 3).

To establish where the panopeid crabs align within the Eubrachyrua, amino acid and nucleotide sequences from 112 mitogenomes comprising 77 genera from 28 families of brachyuran crabs were used along with the three new mitogenomes (Fig. 2). Two sequences that are publicly available for brachyurans were not included in our analysis due to inconsistences with the sequences. (1) The protein sequences for Gecarcoidea natalis contained ambiguous amino acid identifications, resulting in poor alignment with other members within the superfamily Grapsoidea. (2) The protein sequences for Pyrhila pisum aligned poorly with other members of the Brachyura; however, there were no missing protein codes. When tested in BLASTp, the proteins for P. pisum yielded low identity with other brachyurans; < 60% identity in most cases.

Four distinct clades were identified (Fig. 2). One clade belongs to crabs in the subsection Heterotremata (n=40), a second belongs to crabs in the subsection Thoracotremata (n=44) and a third belongs to crabs in the section Podotremata (n=7). The fourth clade belongs to the ‘Old World’ freshwater crabs in the superfamilies Potamoidea and Gecarcinucoidea (n=20). This fourth clade forms a subsection termed Potamoida, a sister group to Thoracotremata. The split between the Heterotremata and the Potamoida/Thoracotremata clades is well supported using both amino acid sequences (Sh-aLRT/UFBoot: 100/100) as well as nucleotide sequences (Sh-aLRT/UFBoot:100/100). The Potamoida and Thoracotremata split is also well supported using both sequence types (amino acids- Sh-aLRT/UFBoot: 98.9/99; nucleotides- Sh-aLRT/UFBoot: 92.9/98).

The panopeid crab species E. depressus, P. herbstii and R. harrisii formed a branch for the family Panopeidae (Fig. 2 and Fig. 3; “β”) aligned alongside the xanthid branch to form the superfamily Xanthoidea (amino acids- Sh-aLRT/UFBoot: 100/100; nucleotides- Sh-aLRT/UFBoot: 100/100). The xanthid branch contains members of the Xanthidae family: E. anaglyptus, A. floridus and A. integerrimus (Fig. 2 and Fig. 3; “α”). When considering its amino acid sequences, the crab species Epixanthus frontalis from the family Oziidae aligns with the Xanthoidea superfamily with moderate support (Sh-aLRT/UFBoot: 89.4/94) (Fig. 2). The nucleotide sequences for E. frontalis show a similar pattern; however, Leptodius sanguineus is part of the branch with middling support (Sh-aLRT/UFBoot: 66.7/93) (Fig. 3). Based on amino acid comparison, L. sanguineus (considered a member of the Xanthidae) aligns between E. frontalis and members of the Pilumnidae, on a branch separate from other xanthid crabs (Sh-aLRT/UFBoot: 16/66) (Fig. 2). The amino acid phylogeny suggests that the hydrothermal vent crabs in the family Xenograpsidae align with the terrestrial crabs in the family Ocypodidae (Sh-aLRT/UFBoot: 83.2/72), yet the nucleotide sequences suggest that the xenograpids form their own branch alongside of the sesarmid crabs (Sh-aLRT/UFBoot: 100/99).

The family Sesarmidae (10 mitogenomes) appears to be polyphyletic. Rather than grouping together, the genus Chiromantes is split, where C. dehaani aligns with Sesarma neglectum (amino acids- Sh-aLRT/UFBoot: 99.5/100; nucleotides- Sh-aLRT/UFBoot: 100/100), and C. haematocheir aligns with Sesarmops sinensis (amino acids- Sh-aLRT/UFBoot: 99.7/100; nucleotides- Sh-aLRT/UFBoot: 100/100) (Fig. 2 and Fig. 3).

The gene arrangements for the panopeid crabs E. depressus, P. herbstii and R. harrisii (Fig. 3) correspond in synteny to other sequenced xanthid species: E. anaglyptus, A. floridus and A. integerrimus. This gene arrangement differs from both the PanGO and BraGO, where the rrnL and rrnS are adjacent to each other and the trnV is transposed past the CoRe (Table 1). The gene order for the xanthid crab species L. sanguineus reported by Tan et al. (2018) is different from the other xanthids presented in this study, following the basic BraGO pattern rather than the shared pattern of the superfamily Xanthoidea. Based on the CREx test, the new gene arrangement XanGO shares 870 common intervals with PanGO and 988 common intervals with BraGO (Fig. 3), suggesting it to be a low-level rearrangement relative to the common gene arrangements. The new XanGO is most different to the MaVaGO, sharing only 80 common intervals.

The mitogenomes of the crabs in the family Panopeidae all shared a ~600 bp long intergenic spacer between the rrnS and trnV ncRNA genes (E. depressus, 618 bp; P. herbstii, 622 bp; R. harrisii, 644 bp) representing the control region (CoRe) (Fig. 1). The CoRe in the panopeid mitogenomes are A + T skewed (78.40–80.22%) and contain the repeated motifs TA (125–107), AT (112–104), TAA (47–39), TTA (30–40), ATA (43–35) and TAT (41–37). The mitogenomes of the xanthid crabs used in this study also have similar sized intergenic spacers in this region, suggesting that this is the putative location of the CoRe for members of the superfamily Xanthoidea. The CoRe nucleotide sequence for all species within Xanthoidea were isolated and run through BLASTn, resulting in a lack of any significant similarity, suggesting high mutability.

The mitogenomes of the panopeid crabs E. depressus, P. herbstii and R. harrisii support the position of the Panopeidae within the Heterotremata, helping to build/support the branch belonging to the superfamily Xanthoidea (Ng 2008). Along this branch, the Xanthidae and Panopeidae form sister groups, additionally supported by previous genetic data using five or less mitochondrial and nuclear genes (Thoma 2009; Lai 2011; Thoma 2014). The genera within these families have been historically identified as polyphyletic (Thoma 2009) and the limited number of mitogenomes available makes it difficult to determine their validity. We acknowledge that the families Xanthidae and Panopeidae both occur in two forms: sensu stricto and sensu lato (Ng 2008). There are 4 publicly available mitogenomes for the Xanthidae (GenBank) and we provide 3 additional mitogenomes for the Panopeidae. We have treated these families in their simple form due to the lack of genetic information to split them further. As more mitogenomes become available, the validity of the two forms should be revisited.

Several taxonomic conflicts appear when considering mitogenetics surrounding the Xanthoidea. First, based on morphology and limited mitogenome availability, the genus Leptodius is considered a member of the family Xanthidae. However, despite this genus having 12 separate species, only one mitogenome (the species L. sanguineus) is available for analysis. Previous studies showed that L. sanguineus aligns closely with other members of the Xanthidae (Sung 2016; Karagozlu 2018; Xie 2018; Ma 2019), but these studies use fewer brachyuran mitogenomes in their analysis prior to our study. When considering all mitogenomes available for the Brachyura in our investigation, L. sanguineus aligns more closely with the members of the family Oziidae, rather than the Xanthidae. This interesting observation merits further exploration.

Gene arrangement changes were once thought to be rare (Boore 2000) but with greater availability of mitogenome sequencing, it appears that changes in gene arrangements can be common across groups. For example, gene order is conserved within Osteichthyes and some subgroups of Mammalia, while it varies strongly in e.g. Ctenophora (Arafat 2018), Mollusca (Guerra 2018), Hymenoptera (Dowton 1999) and Anomura (Tan et al. 2018). For Crustacea, some species within the Stomatopoda, Amphipoda and Dendrobranchiata still carry the PanGO ground pattern of Pancrustacea (Shen 2011), while no sequenced species within the Brachyura have retained this gene order. Studies on gene order rearrangement are ongoing with some hypothesizing that the evolution to living within harsh environments, such as the deep sea or hydrothermal vents, can lead to new gene synteny (Nakajima 2016; Gan 2018; Tan et al. 2019).

The brachyurans include several families found in the deep sea. Two of them are represented herein: Bythograeidae and Xenograpsidae. Bythograeidae possess the BraGO arrangement plesiomorphic for Brachyura, while Xenograpsidae have their own gene arrangement (XenGO). In contrast, the freshwater crab family Potamidae has 9 different gene arrangements (Zhang 2020). Brachyuran crabs represent both cases: the adaptation from a marine to a freshwater environment was likely harsh and may have resulted in several new gene arrangements, while in contrast, the evolution of crabs to the deep-sea benthos resulted in some retaining the ancestral gene order in the face of a new environmental extreme. Therefore, when considering crabs, living within harsh environments does not seem to be the only answer to gene arrangement plasticity, but perhaps requires consideration at the finer scale of environmental adaption. Similar findings have been reported by Tan et al. (2019) who found little evidence for linking gene order rearrangements with adaptations to extreme environments, concluding that these cues are poorly understood and merit a more detailed approach.

A comparison of the eubrachyuran subsections shows that Heterotremata has a higher diversity of gene arrangements than Thoracotremata. Both subsections share species whose gene arrangement follows the basic BraGO pattern. Aside from the BraGO, Thoracotremata only has 3 unique gene arrangements while Heterotremata has 8 unique gene arrangements (including the herein newly established XanGO). This does not include the gene arrangements for the freshwater crabs in the superfamilies Potamoidea and Gecarcinucoidea. The freshwater crabs have more unique gene arrangements than the known Heterotremata.

The panopeid crabs E. depressus, P. herbstii and R. harrisii all have the trnV gene transposed from between the rrnL and rrnS genes to a location past the CoRe. This differs from the PanGO, BraGO, SesGO, XenGO, DamGO, MajGO and DynGO, which all have the trnV gene located between the rrnL and rrnS genes, with the CoRe following the rrnS gene. The xanthid crabs E. anaglyptus, L. sanguineus, A. floridus and A. integerrimusi all share the latter gene arrangement, suggesting that it might be a conserved arrangement within Xanthoidea and thus support our interpretation of the new Xanthoidea gene arrangement (XanGO). The intergenic spacer found between the rrnS and trnV genes in panopeids appears to be the putative location of the CoRe for these species and is shared with xanthid species, E. anaglyptus, A. floridus and A. integerrimus. All have similarly sized intergenic spacers (600–750 bp long) at this location, suggesting that this may be the location of the CoRe across Xanthoidea. Apart from L. sanguineus, the Xanthidae all follow the new gene arrangement XanGO. Leptodius sanguineus follows the plesiomorphic brachyuran gene arrangement BraGO and based on its amino acid sequences, it groups more closely with the family Pilumnidae than the members of the Xanthidae or the panopeids presented here; however, nodal support is low, meriting further study and sequencing of closer relatives. Higher nodal support is offered with the nucleotide tree, where L. sanguineus groups with Epixanthus frontalis from Oziidae rather than with the xanthids. Based on the molecular taxonomy and its gene arrangement, the placement of L. sanguineus within Xanthidae appears to be invalid and in need of revision, adding to our explanation above.

The mitogenome analysis we performed also supports the renaming of two gene arrangements and confirms the correct gene sequence for another. Two mitogenomes were available for the pilumnid crabs, Echinoecus nipponicus and Pilumnus vespertilio. They follow the gene arrangement reported by Tan et al. (2018) and differ from BraGO in having the trnL gene transposed from its location between the cox1 and cox2 genes to a location between the second trnL and rrnL genes. This gene arrangement was reported by Tan et al. (2018) as number 12, but we propose Pilumnidae gene order (PilGO) to follow the original gene nomenclature determined by Basso et al. (2017). Similarly, the gene arrangement reported as number 5 by Tan et al. (2018) we rename to the Somanniathelphusa gene order (SomGO). Basso et al. (2017) report the gene arrangement GeoGO as having the trnL gene between the cox1 and cox2 genes, but based on the gene arrangement listed in Genbank, this is nonconcurrent. The correct gene arrangement was reported by Tan et al. (2019) and is supported here with the addition of the mitogenome for Geothelphusa sp. (MG674171), where the trnL gene is located between nad1 and the second trnL gene. This corrected nomenclature should be incorporated into further taxonomic assessments.

This study provides an updated mitophylogeny for the Brachyura, utilizing all available mitogenomes, along with the first mitogenomes for the Panopeidae, a highly abundant group of ecologically important estuarine crabs with a limited phylogenetic understanding. Our data support the subsection, Potamoida, within the Eubrachyura. The addition of E. depressus¸ P. herbstii and R. harrisii mitogenomes provides a greater phylogenetic understanding of a group that has been taxonomically challenging in the past. Moreover, the addition of mitogenomes from the Panopeidae further supports the split of the Xanthoidea into multiple families. The novel gene arrangement we describe within the Heterotremata, increases the total number of unique gene arrangements within this subsection to eight. Whilst our results clarify some phylogenetic relationships, they also highlight the need for further study of the genus Leptodius which appears to be incorrectly placed within the subfamily Xanthoidea. Greater sequencing efforts will provide more comparative data for these underrepresented crab groups, and should include the incorporation of nuclear genetic data where possible.