Systematic assessment of the Panopeidae and broader Eubrachyura (Decapoda: Brachyura) using mitochondrial genomics

This study provides a broad phylogenetic analysis for the Eubrachyura, with the inclusion of three new Panopeidae mitochondrial genomes: Eurypanopeus depressus (flatback mud crab) (15,854bp), Panopeus herbstii (Atlantic mud crab) (15,812bp) and Rhithropanopeus harrisii (Harris, or ‘white-fingered’ mud crab) (15,892bp). These new mitogenomes were analyzed alongside all available brachyuran mitochondrial genomes (n = 113), comprising 80 genera from 29 families, to provide an updated phylogenetic analysis of the infra-order Brachyura (“true crabs”). Our analyses support the subsection Potamoida within the Eubrachyura as the sister group to Thoracotremata. The family Panopeidae aligns with the family Xanthidae to form the Xanthoidea branch, which is supported by current morphological and genetic taxonomy. A unique gene arrangement termed ‘XanGO’ was identified for the panopeids and varies relative to other members of the subsection Heterotremata (within the Eubrachyura) via a transposition of the trnV gene. This gene arrangement is novel and is shared between several Xanthoidea species, including Etisus anaglyptus (hairy spooner crab), Atergatis floridus (brown egg crab), and Atergatis integerrimus (red egg crab), suggesting that it is a conserved gene arrangement within the Xanthoidea superfamily. Our study further reveals a need for taxonomic revision of some brachyuran groups, particularly the Sesarmidae. The inclusion of panopeid mitogenomes into the greater brachyuran phylogeny increases our understanding of crab evolution and higher level Eubrachyuran systematics.


Introduction
Brachyura ("true" crabs) is the largest subgroup of the Decapoda (Crustacea). It is a ubiquitous group, whose members thrive in terrestrial and aquatic habitats but are particularly prevalent in marine environments (Tsang 2009;Jia et al. 2018;Ma et al. 2019;Tan et al. 2019). Marine Brachyura boast a broad range of morphological and ecological diversity, leading to a complex taxonomy (Yong-kun et al. 2014). Historically, brachyurans were divided into two sections: the Eubrachyura and the Podotremata, with the Eubrachyura being further divided into the subsections: Heterotremata and Thoracotremata (Guinot 2013). The sub-sectioning of these Eubrachyura is based almost entirely on typological morphology (particularly the genital openings) and has been subject to debate with regards to monophyly. Due to the morphological complexity of this group, genetic tools and analytical methods are typically used to resolve systematic discrepancies (Basso et al. 2017;Jia et al. 2018;Wang et al. 2018).
High throughput sequencing (HTS) has proven effective in advancing and resolving taxonomies . Early studies on brachyuran mitogenomics relied on long PCR and primer walking techniques to read and assemble the mitogenome (Yamauchi 2003;Miller 2005). The rapid sequencing and assembly of the mitochondrial genome (mitogenome) using HTS has proven to be a powerful tool for conducting phylogenetic studies of eukaryotes ). The small size of mitogenomes (~ 14-16 kb), potential for high mutation rate, and a simple closed structure, make them an ideal marker for inferring an organism's mitogenetic phylogeny (Boore 1999). Brachyura has been classified into 93 families with over 7000 species (Ng 2008); for 111 species (representing 28 families) the complete mitogenomes are known (NCBI and Supplementary material 1). The Brachyuran genomes that have been sequenced thus far all are between 10-25 kb in length. Much of these data represent three brachyuran families [Portunidae (n = 15); Varunidae (n = 14); Sesarmidae (n = 10)]. The remaining 25 families have < 5 mitogenomes sequenced per family, many of them only having 1 (see Supplementary material). Recent publications using mitogenomics have challenged the validity of the morphologically founded Eubrachyuran subsections of Heterotremata and Thoracotremata. For example, freshwater crabs in the family Potamidae fall into the Heterotremata based on morphology but based on mitogenomics align with members of the Thoracotremata (Basso 2017;. Mitogenomes also offer insights into gene arrangement, which can have diagnostic properties at different systematic levels Boore 2000;Moret 2001;Perseke et al. 2008;Zhuang 2010;Babbucci 2014;Mindell 2016;Nakjima et al. 2016;Zhang 2020). Within the Malacostraca, mitogenome gene arrangements are conserved within certain groups (Shen 2011;Tan et al. 2019), which allows for simple comparisons at different taxonomic levels. Grouping the Crustacea with the Insecta to form the Pancrustacea has strong support based on the near-identical arrangement of the shared genes across taxa . However, gene arrangement can vary greatly within crustacean orders. Specifically, in Brachyura many species share a gene arrangement that is thought to be ancestral to Brachyura (termed BraGO). However, recent research has shown that most groups deviate from this pattern forming new arrangements at the family and subfamily levels (Basso 2017;Wang et al. 2018;Wang 2020a;Wang 2020b). BraGO differs from the Pancrustacea genomic order (PanGO) by a transposition of the gene trnH from a location between the nad5 and nad4 genes to a location between the trnD and trnF genes (Basso 2017). To date, 20 different gene arrangements have been identified within the Brachyura (Basso 2017;Zhang 2020a), but many brachyuran groups remain un-sequenced and this syntenic diversity could be much higher.
An example of an understudied brachyuran group is the superfamily Xanthoidea (Brachyura), which boasts high diversity across the world's oceans (Karasawa 2006). Species within this superfamily share a high level of morphological similarity and are often poorly described both morphologically and genetically (Ng 2008;Thoma 2014). The number of families and subfamilies within the Xanthoidea has changed drastically in recent years (Lai 2011). Two common families, the Xanthidae and Panopeidae, share several morphological features that can lead to systematic confusion and difficulty in identifying them beyond the family level (Shih 2011). Both families are found in temperate and tropical shallow intertidal and subtidal zones, but xanthid crabs have a circumtropical distribution while panopeids are only found in the Americas, excluding global invasions (Thoma 2014). To date, there are only four mitogenomes available for the Xanthidae and none for the Panopeidae, whose systematics have primarily relied upon a select number of genes or morphological keys (Williams 1984;Schubart 2000). Studies using conventional PCR to amplify and sequence mitochondrial and nuclear markers revealed that the genera Eurypanopeus and Panopeus are polyphyletic (Schubart 2000). Similarly, studies on the panopeid genus Hexapanopeus using 12S and 16S genes as markers have also suggested that this genus is polyphyletic (Thoma 2009). Later studies using three mitochondrial markers (COI, 12S and 16S) and three nuclear markers [18S, enolase (ENO) and Histone H3 (H3)] revealed that Xanthoidea is monophyletic, but its two families are not and are in need of taxonomic revision (Thoma 2014).
In this study, we enhance understanding of brachyuran systematics by adding three complete mitogenomes for the Panopeidae: Eurypanopeus depressus, Panopeus herbstii and Rhithropanopeus harrisii from their native range along the Atlantic coast of North America. The genetic composition, genetic similarity and gene ar-rangement of these three panopeid species are described relative to other brachyuran mitogenomes, allowing us to update the brachyuran mitogenomic phylogeny and explore brachyuran-wide classification. A new gene arrangement for the superfamily Xanthoidea is described as well as a renaming of previously reported gene arrangements suggested for other Brachyura.

Specimen collection and dissection
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.

DNA extraction, sequencing and assembly
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  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).

Phylogenetic and mitochondrial gene order assessment
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 . 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.

Phylogenetics
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/UF-Boot:100/100). The Potamoida and Thoracotremata split is also well supported using both sequence types (amino  -cox3, cob, atp6, atp8, nad1, nad5, nad4, nad4-L). Some families have been collapsed for increased clarity (triangles). Black circles on nodes represent an SH-aLRT and bootstrap support of greater than 90/90. Stars (*) indicate areas on the tree with taxonomic conflicts related to previous literature. The symbol α indicates the family Xanthidae; β indicates the family Panopeidae. See Supplementary material 1 for a list of the species used and their accession numbers involved in this analysis.

Gene arrangement among the Brachyura, incorporating the Panopeidae
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  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) Green boxes indicate rRNAs. Purple boxes indicate the control region (CoRe). The red lines along the bottom of the gene orders represents areas within the gene order that are located on the negative strand. Not shown are the 9 unique gene orders for the freshwater crabs (see Zhang et al. 2020). The CREx results are listed for the different gene orders. In the associated table, gene orders with high similarity (> 1000) have red boxes while those with low similarity (< 200) have blue boxes. Intermediate similarity remains white.
nucleotide sequence for all species within Xanthoidea were isolated and run through BLASTn, resulting in a lack of any significant similarity, suggesting high mutability.

Discussion
This study provides the first mitochondrial genomes for three members of the Panopeidae and an updated concatenated mito-phylogenetic analysis for the Eubrachyura (excluding nuclear genetic data), informing upon the systematics of multiple families and higher taxonomic rankings. In addition, the mitochondrial genomes for members of the Panopeidae are identified with a consensus gene arrangement shared with other Xanthoidea (XanGO). These results advance our systematic understanding of the brachyurans through the exploration of mitochondrial genomics and gene synteny rearrangement events.

Xanthid systematics considering panopeid mitogenomic data
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.

Mitogenomic gene arrangements across the Brachyura
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  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  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  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  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.

Conclusions
This study provides an updated mitophylogeny for the Brachyura, utilizing all available mitogenomes, along with the first mitogenomes for the Panopeidae, a high-ly 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.

Author contributions
AMHB collected the crabs used in the study. JB performed the extraction and bioinformatic processing/assembly of the mitogenomes. LAJ and JB performed the phylogenetics and gene similarity assessments. Gene order analysis and annotation was performed by LAJ and JB. LAJ, AMHB, KAM, DCB and JB contributed to the writing of the manuscript.

Competing interests
The authors declare no competing interests.

Acknowledgements
Thanks to Mr. Christopher Moore, who helped to collect the crabs used in the study. Funding for the research and staff time was attained from East Carolina University (AMHB and KAM), the University of Florida (LAJ and DCB) and National Horizon Centre, Teesside University (JB).