Convergent Adaptation of True Crabs (Decapoda: Brachyura) to a Gradient of Terrestrial Environments

Abstract For much of terrestrial biodiversity, the evolutionary pathways of adaptation from marine ancestors are poorly understood and have usually been viewed as a binary trait. True crabs, the decapod crustacean infraorder Brachyura, comprise over 7600 species representing a striking diversity of morphology and ecology, including repeated adaptation to non-marine habitats. Here, we reconstruct the evolutionary history of Brachyura using new and published sequences of 10 genes for 344 tips spanning 88 of 109 brachyuran families. Using 36 newly vetted fossil calibrations, we infer that brachyurans most likely diverged in the Triassic, with family-level splits in the late Cretaceous and early Paleogene. By contrast, the root age is underestimated with automated sampling of 328 fossil occurrences explicitly incorporated into the tree prior, suggesting such models are a poor fit under heterogeneous fossil preservation. We apply recently defined trait-by-environment associations to classify a gradient of transitions from marine to terrestrial lifestyles. We estimate that crabs left the marine environment at least 7 and up to 17 times convergently, and returned to the sea from non-marine environments at least twice. Although the most highly terrestrial- and many freshwater-adapted crabs are concentrated in Thoracotremata, Bayesian threshold models of ancestral state reconstruction fail to identify shifts to higher terrestrial grades due to the degree of underlying change required. Lineages throughout our tree inhabit intertidal and marginal marine environments, corroborating the inference that the early stages of terrestrial adaptation have a lower threshold to evolve. Our framework and extensive new fossil and natural history datasets will enable future comparisons of non-marine adaptation at the morphological and molecular level. Crabs provide an important window into the early processes of adaptation to novel environments, and different degrees of evolutionary constraint that might help predict these pathways. [Brachyura; convergent evolution; crustaceans; divergence times; fossil calibration; molecular phylogeny; terrestrialization; threshold model.]

Over 80% of estimated species comprising extant multicellular life inhabit terrestrial and freshwater ("non-marine") settings (Román-Palacios et al. 2022).Microbial life began to populate terrestrial habitats in the Precambrian, with eukaryotes potentially originating in non-marine settings around 1.6 Ga (Jamy et al. 2022), although major multicellular groups such as animals and plants were ancestrally marine.Their terrestrialization followed in the early Paleozoic (approximately 538-444 Ma), led by arthropods entering coastal and marginal marine settings (e.g., estuaries, lagoons), and plants that transformed the land and its sediments (Buatois et al. 2022).Although molecular divergence time estimates infer early Paleozoic ages for terrestrial arthropod crown groups (e.g., Benavides et al. 2023;Bernot et al. 2023), recognizable body fossils of millipedes, arachnids, and hexapods have recorded their presence on land by the onset of the .Subsequently, these groups radiated to become prominent components of terrestrial biodiversity.Fossil evidence suggests potential transitions through marginal marine settings (Edgecombe et al. 2020;Lamsdell et al. 2020) but transitions for many modern groups lack such clues (e.g., the remipede sister group of hexapods is now predominantly restricted to marine layers within anchialine caves), hinting at complex ecological pathways.Here, we examine the evolutionary history of a clade, the true crabs (Decapoda: Brachyura), which might provide insights into the early phases of adaptation from marine to non-marine environments, now obscured by extinction.
As with life, in general, crabs have an unequivocally marine ancestor (Watson-Zink 2021).The largest group of Brachyura, called Eubrachyura, which contains all non-marine members, could be as old as the mid-Jurassic  based on phylogenomic divergence time estimates (Wolfe et al. 2019).During the "Cretaceous Crab Revolution" (145-66 Ma), many now-extinct lineages appeared briefly, accompanied by the divergence of many extant superfamilies (Luque et al. 2019b;Wolfe et al. 2019;2021).Although the direct record of fossil crabs from non-marine sediments is depauperate, one well-preserved example of a completely extinct non-marine eubrachyuran lineage is known from around 100 Ma (Luque et al. 2021), and chelipeds of uncertain affinity from non-marine sediments around 74 Ma (Robin et al. 2019).Together, these fossils suggest that crabs have been entering non-marine habitats for the majority of their evolutionary history.
Complementary to direct fossil evidence, dated phylogenies and character mapping have recently been applied to investigate the evolution of crab terrestriality (Davis et al. 2022;Tsang et al. 2022).Eubrachyura has been previously divided into 2 presumed clades based on the position of the male gonopores: "Heterotremata" and Thoracotremata.In a 10-gene molecular study focused on the relationships of the clade Thoracotremata, the common ancestor of this clade was found to be "semi-terrestrial" (in Tsang et al. (2022), this referred to intertidal habitats) and Cretaceous in origin, with at least 4 transitions to terrestrial and 2 or 3 transitions to freshwater lifestyles, all within the Cenozoic (Tsang et al. 2022).In one instance, the authors estimated at least 6 returns to subtidal marine habitats and hypothesized that sesarmid crabs (specifically Geosesarma, vampire crabs) transitioned from terrestrial to freshwater habitats (Tsang et al. 2022).A separate supertree-based study across Decapoda inferred 3 transitions to terrestriality and 3 to freshwater, and 1 reversal from terrestrial to marine habitats, within all of Brachyura (Davis et al. 2022).The oldest event, encompassing the freshwater heterotreme groups Potamoidea, Gecarcinucoidea, and Pseudothelphusoidea, occurred in the upper Cretaceous, with others in the Cenozoic.Additionally, Davis et al. (2022) inferred higher rates of speciation in non-marine crabs, but habitat shifts were not found to be a significant causal factor driving crab diversity.
The aforementioned phylogenetic studies, however, treated marine, terrestrial, and freshwater lifestyles as largely discrete ecologies for crabs.Indeed, previous studies have described a gradient of terrestrial change based on independence from standing water (e.g., Bliss 1968;Powers and Bliss 1983;Hartnoll 1988), with the caveat that no known crab is completely independent from water throughout its entire life cycle.Others (e.g., Yeo et al. 2008;Cumberlidge and Ng 2009;Cumberlidge et al. 2009) focused on the 7 exclusively ("primary") freshwater crab families and their vicariant biogeography leading to high endemicity and risks of extinction but rarely drew comparisons with terrestrial crabs.Recently, Watson-Zink (2021) unified the conceptualization of the terrestrial and freshwater crab lifestyles as a series of ecological, morphological, and physiological traits describing grades of terrestriality (described in Supplementary Table S1).Crabs can leave fully marine lifestyles (Fig. 1a-e) along either of two transition pathways: through marine-associated environments (e.g., the "direct" pathway of Tsang et al. (2022), via intertidal, mangroves, beaches: Fig. 1f-j) or through freshwater environments (e.g., the "indirect" pathway via estuaries, rivers: Fig. 1k-o, akin to the transition in amphibians).Each grade of terrestriality is loosely associated with habitats: lower intertidal and estuaries (grade 1), upper intertidal and freshwater (grade 2), beaches and riverbanks (grade 3), and coastal forests and jungles, including tree climbing (grades 4-5; Watson-Zink 2021).Less terrestrial crabs (grades 1 and 2) in either pathway can tolerate fluctuating environments, with osmoregulatory ability likely playing a major role in these lifestyles (Watson-Zink 2021).Crabs of higher terrestriality (grades 3-5) possess further morphological and developmental adaptations, such as branchiostegal lungs and water-wicking setae to prevent desiccation, and increasingly abbreviated larval development and parental care (Watson-Zink 2021).Note that the grades do not represent an ultimate "goal" of terrestriality, as many groups successfully remain and diversify within lower grades.Indeed, it is evident that multiple brachyuran families have repeatedly evolved members of both transition pathways and various grades, but their distribution across crab phylogeny (within and beyond Thoracotremata) over time remains unclear.
To resolve the convergent evolution and timing of terrestriality, we present the most robust molecular taxon sampling to date for Brachyura, representing 333 species and 88 of 109 families.As our data represent only Sanger sequences of 10 loci, revision of brachyuran systematics is beyond the scope of the current study (Timm and Bracken-Grissom 2015), and efforts are currently underway to clarify deep relationships using phylogenomics.Furthermore, to partially ameliorate false confidence in the topology, we provide additional metrics describing the degree of nodal uncertainty.Using these data, we contrast divergence times inferred using 36 newly vetted detailed calibrations (Luque et al. 2023) with traditional node dating, and 328 calibrations sampled from the Paleobiology Database under the fossilized birth-death (FBD) and skyline models of tree evolution, one of only a few empirical comparisons of its type.Finally, we summarize the natural history traits of each sampled crab to assign each along a gradient of transitions from marine to terrestrial lifestyles (Watson-Zink 2021) and use Bayesian threshold models of ancestral state reconstruction to estimate convergent events.

Taxon and Gene Sampling, DNA Extraction, PCR, and Sequencing
The molecular dataset includes 88 families, 263 genera, 333 species, and 338 individuals within the infraorder Brachyura (Ng et al. 2008;WoRMS 2022;Poore and Ahyong 2023), as well as 6 outgroups for a total of 344 tips.Fourty-one percent of sequence data were new, with the remainder obtained from GenBank (Supplementary Table S2).A total of 10 genes were selected based on previous phylogenetic research on decapods (Spears and Abele 1998;Schubart et al. 2000;Tsang et al. 2008Tsang et al. , 2014;;Bracken-Grissom et al. 2013).These included 2 mitochondrial ribosomal RNA (rRNA) coding genes, 2 nuclear rRNA genes, and 6 nuclear protein-coding genes (Supplementary Text S1).A minimum of 2 genes were required for each taxon included in the analysis, with an average of 7 genes available per tip.
Genomic DNA was extracted from the gills, abdomen, pereopod, or pleopod, using the Qiagen DNeasy® Blood and Tissue Kit, QIAamp DNA Mini Kit, or QIAamp DNA Micro Kit.Gene regions were amplified with polymerase chain reaction (PCR) using one or more sets of primers (Supplementary Text S1).PCR amplifications and sequencing reactions were performed as described in Supplementary Text S1.

Fossil Calibration and Divergence Time Inference
We compared 2 strategies for fossil calibration: (1) 36 newly vetted node calibrations, and (2) 328 fossil occurrences from the Paleobiology Database (PBDB; http:// paleobiodb.org/).For node calibration, all calibrations followed best practices regarding specimen data, morphological diagnosis, and stratigraphy (Parham et al. 2012;Wolfe et al. 2016; extensive details in Luque et al. 2023, summary in Supplementary Table S3 herein), and were assigned to a crown group node at the family level or higher.This node dating strategy used a birth-death tree prior and uniform calibration age distributions.
We downloaded fossil occurrences from the PBDB on 23 March 2022, for Brachyura at family-level taxonomic resolution (details in Supplementary Text S1).We randomly subsampled 10% of the 3276 remaining occurrences, resulting in a computationally tractable 328 occurrences.The subsample was a reasonable sample of ages (including the oldest possible age), represented 69% of fossil families, and slightly overrepresented non-marine paleoenvironments (details in Supplementary Table S4).All fossil occurrences were assigned age ranges from the PBDB, each with a uniform distribution (Barido-Sottani et al. 2019, 2020).To incorporate these fossil samples as part of the inferred evolutionary process, we used the unresolved time-homogeneous FBD tree prior (Stadler 2010;Heath et al. 2014;Bapst et al. 2016; O'Reilly and Donoghue 2020), with parameters described in Supplementary Text S1.
To reflect the complete absence of brachyuran fossils from earlier than the Jurassic, which represents a known ghost lineage when compared to the diversity and abundance of outgroup anomuran fossils, thus earlier diversification of the sister group (Hegna et al. 2020;Wolfe et al. 2021), we also analyzed the fossil occurrence calibration set using a birth-death skyline tree prior with sequential sampling (BDSS; Stadler et al. 2013;Culshaw et al. 2019).Fossil sampling proportion was modeled as time-heterogeneous with time slices before and after the oldest fossil sample (details in Supplementary Text S1) using the TreeSlicer function from the skylinetools package (https://github.com/laduplessis/skylinetools).
All divergence time analyses were conducted in BEAST2 v.2.6.7 (Bouckaert et al. 2019) using a fixed starting tree derived from our ML concatenated results (detailed parameters in Supplementary Text S1).Fossil occurrences (in FBD and BDSS analyses only) were added to the starting tree as "rogues" (able to move within pre-assigned family-level constraints following Barido-Sottani et al. [2023], as most decapod fossils are fragmentary and cannot be confidently assigned).For each calibration strategy and tree prior, we compared two clock models: relaxed lognormal (Drummond et al. 2006) and random local (Drummond and Suchard 2010).Analyses used 4-6 runs for at least 450 million generations with 25% burnin.Convergence was assessed as above.We visualized the results from different parameter sets using chronospace scripts (Mongiardino Koch et al. 2022), which provide a multidimensional representation of inferred node ages that can be broken down by different factors and models.

Ancestral State Reconstruction
To code character states representing a gradient of terrestriality, we used a modified version of the trait-by-environment associations defined by Watson-Zink (2021), additional details in Supplementary Table S1 and Supplementary Text S1.Distinct transition pathways were defined for marine and freshwater routes.For both, we added a grade 0, indicating that the ancestral state for all crabs is fully marine (Fig. 1a-e).Using this framework, we coded discrete grades of terrestriality following 2 schemes.The first scheme coded the taxa that were sequenced in our molecular phylogeny and required justification from natural history literature on adult habitat, osmoregulatory status, larval developmental strategy, primary respiratory structure, water-wicking setae, burrow type, and diurnal activity period (Supplementary Table S5, Supplementary Text S2).As our phylogeny sampled 4% of brachyuran species as tips, we also constructed a scheme to estimate grades for unsampled species for which the phylogenetic positions are unknown.For this scheme, we downloaded all taxonomic data, including non-marine taxa, from WoRMS as of 8 June 2021.For families sampled in our tree, we used WoRMS to estimate the number of species that fall into each grade and accordingly assigned prior distributions to each tip on the molecular phylogeny (Supplementary Table S6).
First, we used stochastic character mapping (Bollback 2006;Revell 2012) to infer ancestral states at each node, using a simplified single dataset (details in Supplementary Text S1).Next, we used Bayesian threshold models (Felsenstein 2012;Revell 2014) to account for gradients of change.A character coded with discrete ordered states (i.e., our grades) was assumed to evolve according to an unobserved continuous trait called "liability" (here representing the coded natural history traits combined with additional unobserved factors).Following Sallan et al. (2018), we assume thresholds that represent the amount of change in terrestriality traits that allow a habitat shift.As there are 2 independent transition pathways, these were analyzed separately for each pathway (from marine to non-marine).Threshold models for ancestral state reconstruction were implemented as the ancThresh function in phytools (Revell 2012).Each ancThresh model was run for 150 million generations with 20% burnin (Revell 2014;Sallan et al. 2018).

Phylogenetic Relationships
The concatenated alignment length comprised 7516 bp in total from 2 mitochondrial rRNA, 2 nuclear rRNA, and 6 nuclear protein-coding genes (gene trees visualized in Supplementary Figs.S1-10).Results using ML and BI were similar, with some deeper nodes (higher than family level) maintaining low to moderate support (UF bootstrap = 50-94; Fig. 2, Supplementary Fig. S11) with ML and generally stronger support (most posterior probabilities ≥ 0.98; Supplementary Fig. S12) with BI.For each node of the ML tree, the gene concordance factor (gCF) reflects the percentage of loci containing all the descendant taxa, and site concordance factor (sCF) the percentage of sites supporting the node (Minh et al. 2020a).Both concordance factors illustrate a spectrum of support across nodes that are fully supported with UF bootstraps (Supplementary Fig. S13).The average gCF was 32.73 (Supplementary Figs.S11 and S13), indicating one third of loci support the average node.However, nearly half of the sites support the average node (average sCF = 45.64),demonstrating the benefit of concatenation for small numbers of loci.
Revision of brachyuran systematics at nodes above the family level would best be undertaken with phylogenomic scale data (Wolfe et al. 2019 and phylogenetic informativeness profiles, Supplementary Fig. S14); therefore, we only briefly summarize the topology results here and in Figure 2. Podotremes are paraphyletic with respect to Eubrachyura, forming the following successive clades: Dromioidea + Homoloidea (this pairing has low support from ML, but strong from BI), Raninoidea, and Cyclodorippoidea (latter two clades with full support).Within Eubrachyura, subsection Heterotremata is paraphyletic with respect to monophyletic Thoracotremata.The so-called primary freshwater crabs are polyphyletic.The African and Eurasian groups (Potamoidea and Gecarcinucoidea) form a clade (UF bootstrap = 95, posterior probability = 0.89), and are themselves the sister clade of the Gondwanan Hymenosomatoidea (UF bootstrap = 92, posterior  S7), by Javier Luque and Harrison Mancke.probability = 1).Together, this group comprises the sister group of Thoracotremata, with moderate support from ML (albeit with low concordance factors) and full support from BI. Within Thoracotremata, some higher-level relationships are weakly supported by ML (UF bootstraps < 75%, low gCF), but most nodes are similar to BI (where they have moderate to high support).Both Grapsoidea and Ocypodoidea are polyphyletic.
Meanwhile, the American freshwater groups branch off within clades including the deepest (Pseudothelphusoidea) and second deepest (Trichodactyloidea) divergences within the main heterotreme group, although these nodes are not strongly supported by traditional metrics in either analysis (concordance factors for both are >60, some of the highest in our data).The remaining heterotremes are subdivided into Majoidea, and two large supported clades containing 24 and 23 families, respectively.Within these latter clades, the superfamilies Eriphioidea and Goneplacoidea are strongly polyphyletic.Some deep splits within both clades are poorly supported (some nodes UF bootstrap < 50, concordance factors = 0, posterior probability < 0.8).

Divergence Times
Results of divergence time inference vary depending on the parameters used, with results of the vetted calibration strategy distinct on the major axis (Fig. 3a), and the FBD and BDSS results being similar to one another.Although none of the random local clock analyses converged after extensive runtime, we plotted samples from their individual chains with burnin of 50% to reduce the effect of poor mixing; when included, the choice of clock model also differs significantly (Fig. 3b), but mostly in the same direction at the same nodes as the calibration strategy (Fig. 3c).In case the unconverged analyses were skewing the results, we plotted results from relaxed lognormal clocks alone (Supplementary Fig. S15), finding the groupings by calibration strategy were mostly upheld, although the spread along axis 2 decreased for vetted calibrations especially.
Similar to Wolfe et al. (2019), using the vetted node calibrations, the divergence of Meiura (i.e., the root node) is inferred in the Permian (mean age at 268 Ma), while crown group Brachyura diverged in the Triassic (mean age 231 Ma), and crown-group Eubrachyura in the Jurassic (mean age 172 Ma).Superfamily-level divergences were inferred in the Jurassic for podotremes, and almost entirely within the Cretaceous within Eubrachyura (Supplementary Fig. S16).
Divergence estimates using fossil occurrence sampling were both considerably younger than with the vetted calibrations, with the root estimate in the Jurassic in both cases (mean ages at 180 Ma for FBD, and 191 Ma for BDSS; Supplementary Figs.S17 and S18).The PBDB-calibrated analyses are relatively immune to the root prior (Supplementary Fig. S19a).In the BDSS analysis, the inferred parameters were: sampling proportion 0.062 (95% HPD of 2.72e−3 to 0.23), death rate 0.091 (95% HPD of 0.024-0.23),and birth rate shown in Supplementary Figure S20.Most other nodes were similarly compressed, with crown-group Brachyura in the Jurassic and superfamily level divergences pushed to the Upper Cretaceous and Paleogene, although the posterior did not follow the marginal prior at some nodes (Supplmentary Fig. S19b).The PBDB analyses often inferred the placement of "rogue" fossils within the stem groups of their families.Consequently, a number of family-level crown group ages were underestimated relative to their known vetted calibrations (e.g., Dromiidae, Dynomenidae, Lyreididae, Percnidae, Varunidae, Euryplacidae, and Panopeidae; Supplementary Table S3).These families are among the most sensitive nodes, in addition to several nodes within Majoidea (Supplementary Fig. S19b).

Evolution of Terrestriality
In the summary of stochastic character mapping, 6 total shifts from marine to non-marine were inferred (Supplementary Fig. S21).These were at the base of Pseudothelphusoidea, Trichodactylidae, Menippidae, Eriphiidae, and Oziidae, and at the base of the clade of (Thoracotremata, Hymenosomatidae, Potamoidea, Gecarcinucidae).The latter node was split, with a slightly higher posterior probability of being freshwater than terrestrial or marine, leading to a freshwater node for the base of (Hymenosomatidae, Potamoidea, Gecarcinucidae) and a terrestrial common ancestor of Thoracotremata.As individual stochastic character maps inferred shifts on branches leading to a single coded tip (e.g., Carcinus maenas), the median number of shifts to non-marine was 13.One shift from the freshwater to terrestrial pathway was found in Hymenosomatidae (note, all non-zero grades were collapsed, so the "terrestrial" members of this family are intertidal).Two shifts from terrestrial to freshwater were found in Glyptograpsidae and Varunidae, respectively.Two reversals to marine were inferred at the clades of Xenophthalmidae and Pinnotheridae, and Xenograpsidae and Cryptochiridae.
Using threshold models for both transition pathways, the best-fitting model was OU based on the lowest output Deviance Information Criterion (Supplementary Table S8).For the direct pathway, 3 shifts to non-marine grades were inferred at nodes: one in Ocypodidae, one at the base of the thoracotreme clade (Percnidae, Grapsidae, Plagusiidae, Mictyridae, Heloeciidae, Macrophthalmidae, Varunidae, Gecarcinidae, Sesarmidae, and Dotillidae), and 1 in Menippidae (Fig. 4, Supplementary Fig. S22).If we consider grades 1 and 2 to be "semi-terrestrial" (similar to a character state from Tsang et al. (2022) referring to intertidal habitats), then the number of node origins for grades 3-5 is 3 or 4: potentially twice in Ocypodidae, once in Mictyridae, and once at the base of the clade formed by (Gecarcinidae, Sesarmidae, and Dotillidae).Based on the estimated liabilities (i.e., thresholds of change required to transition to a different grade), it is 8-36 times easier to move to grades 1 and 2 than to grades 3 and above (Supplementary Table S8, Supplementary Fig. S23a).Three reversals to marine or shifts to freshwater via this pathway (i.e., scored as grade 0) are inferred in: the clade of Xenograpsidae and Cryptochiridae, Plagusiidae, and Varunidae.In the case of Varunidae, the shift is to the indirect freshwater pathway (Fig. 4).Some tips that were coded with a majority of the prior probability failed to infer a shift at any nodes, such   4. Composite of ancestral state reconstructions for the 2 transition pathways under best-fitting OU models in ancThresh, with fully marine clades (all families that are not labeled) reduced for clarity and outgroups removed.Legend for grades and colors representing each pathway at bottom left: fully marine crabs (grade 0), lower intertidal and estuaries (grade 1), upper intertidal and freshwater (grade 2), beaches and riverbanks (grade 3), and coastal forests and jungles (grades 4-5).Pies at nodes represent the estimated ancestral state with the outer ring indicating the pathway (at some nodes, both pathways are shown; when node is inferred marine, no pie is shown).Tip codings are based on estimates by family (Supplementary Table S6), with the collapsed clades showing the color that represents the largest slice of their prior probabilities (split in the case of equal probabilities for 2 grades).For clades that have a small number of taxa in a grade from the opposite pathway, a small triangle is added.Line drawings at right (numbers corresponding to taxa in Supplementary Table S7), by Javier Luque and Harrison Mancke.
For the indirect pathway, 4 shifts to non-marine grades were inferred at nodes: one at the base of the clades Potamoidea + Gecarcinucidae, and one each for Glyptograpsidae, Pseudothelphusoidea, and Trichodactylidae (Fig. 4, Supplementary Fig. S24).If we consider grades 1 and 2 to be "semi-terrestrial," then only Pseudothelphusoidea and Potamidae + Gecarcinucidae are inferred (each with nodes at grade 3).The liabilities indicate that it is extremely easy to move to grade 1, but 11 times harder to move to grade 2, and nearly 100 times harder to move to grade 3 (Supplementary Table S8, Supplementary Fig. S23b).Hymenosomatidae has 25% tip prior probabilities at grade 1, but no shifts were inferred.When analyzed as subclades, a transition was inferred from ancestrally freshwater to non-freshwater in Hymenosomatidae (Supplementary Fig. S25), but we could not infer a transition from ancestrally terrestrial associated (i.e., intertidal/mangroves) to freshwater for Sesarmidae (Supplementary Fig. S26).

Relationships and Divergence of True Crabs
Previous molecular phylogenies of Brachyura have been constructed from eight to 10 Sanger loci (Tsang et al. 2014(Tsang et al. , 2022)), mitochondrial genomes (Tan et al. 2019;Sun et al. 2022;Zhang et al. 2022), transcriptomics (Ma et al. 2019), and genomic target capture (Wolfe et al. 2019).However, the deep relationships among families and superfamilies remain uncertain, as the most extensive study (Tsang et al. 2014) sampled only 58 of 109 families with 8 genes and low support at deep nodes, and more extensive gene sampling was coupled with even lower taxon sampling (Timm and Bracken-Grissom 2015;Wolfe et al. 2019).Although it is evident that the genes we and others have used are insufficient to resolve deep relationships even with current taxon sampling (Supplementary Fig. S14), many regions of the tree are strongly supported.For example, the broadest strokes of our topological results (Fig. 2) contribute to the chorus of molecular and morphological analyses rejecting the monophyly of podotremes (e.g., Ahyong et al. 2007;Tsang et al. 2014Tsang et al. , 2022;;Luque et al. 2019aLuque et al. , 2019b;;Tan et al. 2019) andheterotremes (e.g., von Sternberg andCumberlidge 2001;Tsang et al. 2014;Ma et al. 2019;Tan et al. 2019).Our divergence time estimates are older than most previous publications, except the deeper ages inferred by Wolfe et al. (2019) and the hypotheses of Guinot et al. (2019), but see below for evaluation.
The relationships among thoracotremes were recently examined by Tsang et al. (2022).As in their study, we find polyphyly of Ocypodoidea (fiddler, ghost crabs, and relatives), although our Grapsoidea (shore crabs, land crabs, and relatives) are separated into 5 clades, as opposed to four in Tsang et al. (2022).The main nodes where our results differ are: (i) the derived position of Dotillidae (sand bubbler crabs), (ii) separation of the symbiotic groups Cryptochiridae (coral gall crabs) and Pinnotheridae (pea crabs), and (iii) the position of Plagusiidae (different between ML and BI in our data: see Supplementary Fig. S12).For points ( 2) and (3), we do recover weak support using all metrics.Ultimately, both studies use 8 of the same loci.Our data incorporates the nuclear rRNA genes, with relatively low phylogenetic informativeness above the family level (Supplementary Fig. S14), yet both studies produce strong support at the base of thoracotreme families (note that Tsang et al. (2022) designated UF bootstraps > 90% as strong, which would add several nodes to our "strongly supported" category in Fig. 2 if we followed this metric).Finally, the backbone of thoracotreme phylogeny has been briefly addressed in two phylogenomic studies (Ma et al. 2019;Wolfe et al. 2019), and our results are similar to both, including the position of Sesarmidae (e.g., mangrove and vampire crabs) under models analyzing nucleotide data.A more robust understanding of internal thoracotreme relationships may be derived with additional phylogenomic data, but our results are sufficient to infer ancestral states, with the caveat of lower confidence at the aforementioned nodes.
Polyphyly of the "primary" freshwater crab families (Deckeniidae, Epiloboceridae, Gecarcinucidae, Potamidae, Potamonautidae, Pseudothelphusidae, Trichodactylidae) has been found previously (e.g., von Sternberg and Cumberlidge 2001;Tsang et al. 2022).The inclusion of Hymenosomatidae (pillbox crabs) with the African and Eurasian freshwater groups in our results is novel, as it is the first analysis to incorporate this family in a larger tree.The grouping of African and Eurasian freshwater crabs with thoracotremes is otherwise fairly well supported in previous studies (e.g., Ma et al. 2019;Tan et al. 2019;Zhang et al. 2022).We contradict previous analyses, where the American Pseudothelphusoidea were closely related to the freshwater group (Tsang et al. 2014) and share a number of morphological synapomorphies (Cumberlidge et al. 2021).The American Trichodactylidae were more closely related to other heterotremes.Our results weakly support convergent origins of Pseudothelphusoidea (Fig. 2), leading to subsequent inference of separate transitions to freshwater.Divergence time estimates under all models push the origin of all primary freshwater groups except Trichodactylidae older than 66 Ma, with the African and Eurasian group over 100 Ma, consistent with a deeper cryptic history (Wolfe et al. 2019), particularly in non-marine environments (Tsang et al. 2014;Luque et al. 2021).Nevertheless, the 95% CIs of our divergence time estimates post-date the complete breakup of Pangaea (approximately 175 Ma), and do not support the Gondwanan origin of freshwater crabs (Klaus et al. 2011;Tsang et al. 2014) with these data.
Many of the relationships among families and superfamilies of the remaining heterotremes are not critical for the question of terrestriality.However, most of the deep relationships that we do observe, even with poor nodal support from ML, are congruent with the broad results retrieved from target capture of over 400 loci (Wolfe et al. 2019), except for the position of Menippidae (stone crabs).Our topologies for the relationships within superfamilies are largely congruent with previous Sanger data (Tsang et al. 2014), with some families exchanging places (e.g., Evans 2018;Mendoza et al. 2022) and some other groups with previous conflict still unresolved (Hultgren and Stachowicz 2008;Lai et al. 2014;Windsor and Felder 2014).The position of Dorippoidea nested well within heterotremes is consistent with previous molecular analysis (Tsang et al. 2014), but may confound hypotheses of early fossils that assume dorippoids are the earliest eubrachyuran branch (e.g., Guinot et al. 2019).One potentially novel result is the polyphyly of Bythograeoidea (hydrothermal vent crabs), which had not previously been included in a global analysis with as many genera (absent from Tsang et al. (2014); no outgroups in concatenated analyses of Mateos et al. (2012)).Several non-monophyletic families were broken up by the insertion of closely related families, perhaps representing morphological groups that could be redefined as subfamilies (e.g., Lyreididae, Iphiculidae, potentially the cylodorippid genus Tymolus).In summary, we observe many similarities with previous analyses, and some new hypotheses that await improved molecular sampling before suggesting new systematic changes.

Which Divergence Time Estimates are Reliable?
The factors that we investigated for different methodological choices were the calibration strategy and the clock model.Although divergence time estimates have only been visualized previously for only one dataset of echinoids with chronospace (Mongiardino Koch et al. 2022), we see significant differences in our data.Our results were impacted by the tested methods even more substantially than in the echinoid data, for which the results were sensitive to clock model choice, but not to substitution model or subsets of loci (Mongiardino Koch et al. 2022).In our data, a strong separation of variance of both clock model and calibration strategy was observed (Fig. 3a,b), with a similar pattern of posterior ages from both factors (Fig. 3c).However, the random local clock model, which accounts for the evolution of evolutionary rates in a clade-specific manner (Drummond and Suchard 2010;Ho and Duchêne 2014), failed to converge after months of runtime, so we cannot be certain of the ultimate effect on divergence time estimates.
We were initially interested in using FBD and BDSS because these tree models describe the processes of speciation, extinction, and fossil sampling that led to the true tree and could be more accurate (Wright et al. 2022).
A precursor method that incorporates fossil counts, but not directly into the tree model, has been previously used for lobsters (Bracken-Grissom et al. 2014).The gold standard for improved age accuracy may be to couple process-based tree models with joint inference of topology and ages, particularly total evidence tip dating (i.e., morphological character data from all incorporated fossils; Wright et al. 2022).Brachyura preserve abundant biomineralized hard parts, but the majority of fossils are dorsal carapaces and cheliped fragments.It is challenging to identify phylogenetically diagnostic characters in fragmented fossils due to substantial convergence in carapace and cheliped morphology (Guinot 2019;Luque et al. 2019bLuque et al. , 2021;;Wolfe et al. 2021), and extensive missing data that could compromise their placement within crown groups.For these reasons, it was not possible to include a morphological matrix containing over 300 extant taxa, plus numerous fossils.Nevertheless, our evaluation of the brachyuran fossil record in the context of vetted crown and stem groups identified calibration fossils that were much older than previously appreciated (Luque et al. 2023).Most crown group families have exemplars preserved in the Eocene and Paleogene (66-34 Ma), as well as an increasing number of well-preserved Cretaceous fossils (Ossó 2016;Luque et al. 2017Luque et al. , 2021Luque et al. , 2023)).
Although records are abundant in the PBDB, brachyuran fossil sampling is not as evenly distributed as in many clades used as test cases for the performance of FBD and skyline models (e.g., Gavryushkina et al. 2014;Heath et al. 2014;Renner et al. 2016;Barido-Sottani et al. 2019;O'Reilly and Donoghue 2020).Various biases are observed in PBDB for crabs, including some families and early time slices with few or no fossils, and a geographic bias towards higher latitudes (we are continually working to improve the latter; e.g., Luque et al. 2017).Subsampling a fossil record using the unresolved FBD can be accurate with an evenly sampled clade such as cetaceans (Barido-Sottani et al. 2019), but can be highly inaccurate if the fossils have limited phylogenetic information (O'Reilly and Donoghue 2020).An analysis of mammals (Luo et al. 2021), however, found that fossil sampling density does not have linear effects on divergence time estimates.Altogether, it is difficult to generalize about the best models, unless prior knowledge about node ages is contradicted.
A major issue was estimating the root age of Meiura.Brachyura has a known ghost lineage, with stemgroup members Eocarcinus and Eoprosopon from the early Jurassic (approximately 190 Ma;Haug and Haug 2014;Hegna et al. 2020;Scholtz 2020;Wolfe et al. 2021).The oldest crown group brachyuran fossil occurrences in PBDB were also from the early Jurassic, representing less than 1% of total occurrences (including one occurrence in our subsample).Meanwhile, multiple modern anomuran families were already present in the late Jurassic (164-145 Ma;e.g., Fraaije et al. 2019e.g., Fraaije et al. , 2022;;Robins and Klompmaker 2019), when their divergence likely took place (Bracken-Grissom et al. 2013;Wolfe et al. 2019).This strongly suggests the divergence of the common ancestor of Meiura, and probably of Brachyura, occurred at least by the very earliest Jurassic.Yet, the FBD and BDSS analyses, despite incorporating one Jurassic occurrence and a relatively representative subsample (Supplementary Table S4), inferred impossibly young ages, as observed by O'Reilly and Donoghue (2020).Even the BDSS analysis including a time slice with no fossil sampling before the oldest occurrence (allowing a ghost lineage prior to the Jurassic: Culshaw et al. 2019;O'Reilly and Donoghue 2020) seemed to be a poor fit, and could not avoid estimating an unreasonably young root age (evident from the BDSS marginal prior in Supplementary Fig. S19a).As such, given the data currently available for brachyurans, we recommend caution if using FBD models (and their extensions) to estimate divergences when a morphological matrix is unavailable, because such models may not perform according to expectations in all clades.

How Many Times and When Did Crabs Terrestrialize?
The number of estimated shifts to terrestriality (i.e., any shift from state 0) changes when considering the full diversity of known trait states within each family.We inferred 7 shifts to non-marine lifestyles at nodes, fewer than Tsang et al. (2022), which estimated at least 6 transitions within Thoracotremata alone.Most shifts occur at well-supported nodes (Figs. 2 and 4), although the common ancestor of (Thoracotremata, Hymenosomatidae, Potamoidea, Gecarcinucidae) has only moderate support, perhaps contributing to the uncertainty at this node with stochastic mapping (Supplementary Fig. S21) and the differences we find between stochastic mapping and ancThresh.However, at least 9 additional families have some proportion of their prior probability assigned to grades 1-2 (some all the way up to 100%, but most with lower probabilities), nested among marine sisters (Fig. 4), even though ancThresh does not infer a change at any node.Adding these brings the number of transitions to 16, distributed from Cretaceous to the last 10 myr.Although it is not included in the molecular phylogeny, there was likely another convergent transition to non-marine lifestyle in the Cretaceous fossil Cretapsara (Luque et al. 2021), resulting in a total of at least 17 terrestrialization events across at least 30 families.
The 2 or 3 losses of terrestriality that we estimate include nodes with poor support (the Cretaceous Xenograpsidae and Cryptochiridae, and the Eocene Plagusiidae).Character sequence reversals are rarely favored by threshold models (Revell 2014), so it is intriguing that we find these nodes.Across the tree of eukaryotic life, reversals to marine from non-marine lifestyles are more common than expected (Jamy et al. 2022), and more have been found in Thoracotremata (Tsang et al. 2022), so we could indeed be underestimating the phenomenon of returning to a marine environment.
Across plants and animals, more terrestrial species come from freshwater ancestors than directly from marine ancestors (Román-Palacios et al. 2022), a pathway we only observe potentially once, in Hymenosomatidae.We estimate only one instance of freshwater crabs evolving from intertidal ancestors, in Varunidae.Owing to the ordering of grades through the transition pathways, and the implementation of ancThresh, it is very challenging to infer these changes.The estimated liabilities are too high to infer freshwater Sesarmidae from an intertidal/terrestrial ancestor, although this is further complicated by all sesarmids starting from a higher base grade (i.e., starting at grade 2-3) within the 2 transition pathways, and possibly by the number of species that have convergently evolved arboreal lifestyles in mangroves (Fratini et al. 2005; Naruse and Ng 2020) and in freshwater derived habitats (Diesel 1989).These 4 clades are, however, the only brachyuran groups where the data suggest transitions between the main pathways, so the overall number of convergent events is not affected.

Implications for Early Phases of Arthropod Terrestrialization
An outstanding question is the lifestyle of the common ancestor of the clade containing the majority of non-marine crabs: Thoracotremata, Hymenosomatoidea, Potamoidea, Gecarcinucoidea.Stochastic character mapping suggests the state of this node is uncertain (Supplementary Fig. S20), although the most probable common ancestor may have been estuarine and likely lived in the Jurassic.It is also quite possible that the common ancestor was marine (Fig. 4), and direct and indirect pathways were established independently in the thoracotreme and freshwater clades.If estuarine, the Early Cretaceous common ancestor of Thoracotremata may have transitioned from the indirect pathway to an intertidal grade on the direct pathway, or alternatively this ancestor may have returned to fully marine life.Crabs in low grades of both transition pathways have different osmoregulatory adaptations (Watson-Zink 2021), so perhaps the thoracotreme common ancestor ecologically resembled some modern Hymenosomatidae (Supplementary Table S5) or hypothetically resembled the only non-marine Cretaceous body fossil (which likely lived in brackish or fresh waters and may have been amphibious; Luque et al. 2021).The common ancestors at these early uncertain nodes could have had some degree of osmoregulatory ability, perhaps experiencing early development in marine or estuarine environments (Watson-Zink 2021).
The hypothetical common ancestor of many terrestrialized crabs, developed above, offers some lessons for understanding the Paleozoic terrestrialization of other arthropod groups.Despite a discrepancy of 100-150 myr between divergence time estimates (Cambrian) and body fossils (Silurian-Devonian), there are examples of Cambrian and Ordovician trace fossils that reveal limited excursions into non-marine environments, and perhaps more extensive life in marginal marine settings (e.g., Collette et al. 2010;Mángano et al. 2021;Buatois et al. 2022).Therefore, the fluidity of crab transitions into lower grades of terrestriality could hint at the types of adaptations other arthropods experienced in these early periods: osmoregulation and abbreviated and/or migratory larval development (Watson-Zink 2021).Other ecologies may resemble fiddler crabs that have adapted to coastal hypersaline environments, building their burrows well inland (Thurman 1984).Of course, other arthropods such as insects and arachnids continued on to surpass crabs in their terrestrial adaptations, as they became completely independent of water.Even the most terrestrial crab grades rely on water for, at the minimum, reproduction.
Additional instances of small numbers of taxa entering grades 1-2 through either pathway (e.g., Carcinidae, Panopeidae, freshwater Varunidae in Fig. 4) could provide some insights on the early stages of terrestrial adaptation.In particular, the above examples harbor some of the most persistent introduced crab taxa: Carcinus maenas (European green crab), Hemigrapsus sanguineus (Asian shore crab), and Rhithropanopeus harrisii.These species are notable for tolerating exceptionally wide-ranging salinities as larvae and have wide and migratory habitat preferences as adults (e.g., Young and Elliott 2020).The introduced Eriocheir sinensis (Chinese mitten crab) can tolerate estuarine to full freshwater (Zhang et al. 2019).Another example is the hymenosomatid Halicarcinus planatus, with a wide salinity tolerance that could help this species adapt and move into warming Antarctic waters (López-Farrán et al. 2021).Many other hymenosomatid genera have members in both freshwater and low-salinity estuarine/mangrove habitats, and many have plastic osmoregulatory capabilities (Chuang and Ng 1994).Perhaps the ancestors of diverse non-marine groups originated with lifestyles similar to successful introduced species.
Groups sharing convergent morphological adaptations to higher grades of terrestriality, such as branchiostegal lungs in Gecarcinucidae, Gecarcinidae, Ocypodidae, and Pseudothelphusidae, or water-wicking setae in Gecarcinidae and Sesarmidae (Watson-Zink 2021 and Supplementary Table S5), are deeply separated by over 100 Ma of evolution.Convergent terrestrial morphology in crabs, with likely pathways through a habitat gradient, perhaps with some traits of introduced taxa, could illuminate the hypothesis that arachnids convergently transitioned to terrestriality (Ballesteros et al. 2022).However, it is possible that horseshoe crabs (chelicerates, not decapods) returned to marine habitats even if they are phylogenetically nested within arachnids, as our crab analyses have inferred at least 2 reversals that involved likely intertidal (grade 2) ancestors.

Conclusions
Herein, we inferred a large molecular phylogeny of true crabs, estimated divergence times that were older than previously thought, and estimated the number of transitions from marine to non-marine lifestyles.We found up to 17 convergent transitions through direct and indirect pathways, with at least 3 climbing to higher degrees of terrestrial adaptation.The most highly terrestrial clades were some of the oldest non-marine inferences in our data, with their common ancestors having diverged over 66 Ma.At least 9 more recent events throughout the Cenozoic led to crabs living in intertidal and marginal marine environments, a shift that is estimated to be much easier based on lower threshold liability and likely fewer traits required.As instances of convergent evolution provide emerging models in the form of "natural experiments," the framework we have developed to compare the gradient of adaptations will enable future research that aims to "predict" the constraints leading to repeated trait evolution and better understand the drivers of biodiversity across related groups.

Figure 2 .
Figure 2. Summary of phylogeny and divergence time estimates for Brachyura (88 brachyuran families, 263 genera, 333 species, 338 individuals plus 6 outgroups).Posterior ages were estimated in BEAST2 using a fixed topology resulting from the concatenated ML analysis in IQ-TREE, 36 vetted node calibrations, a birth-death tree prior, and relaxed lognormal clock model.Shaded circles at nodes represent ultrafast bootstraps.Pie slices are colored by superfamily, with the outermost ring colored by taxonomic section.Line drawings, one representative per superfamily (numbers corresponding to taxa in Supplementary TableS7), by Javier Luque and Harrison Mancke.

Figure 3 .
Figure 3. Sensitivity of divergence time estimates to inference strategy plotted with chronospace, with outgroup taxa removed from these analyses.a) Between-group principal component analysis (bgPCA) separating chronograms by calibration strategy and b) by clock model (note: analyses with the random local clock model did not converge, so individual chains are sampled here with 50% burnin).c) Theoretical topologies representing change in branch lengths from the mean along the major bgPCA axis, discriminating based on calibration strategy.Negative extreme (−1 SD branch length) above, positive extreme (+1 SD) below.

Figure
Figure4.Composite of ancestral state reconstructions for the 2 transition pathways under best-fitting OU models in ancThresh, with fully marine clades (all families that are not labeled) reduced for clarity and outgroups removed.Legend for grades and colors representing each pathway at bottom left: fully marine crabs (grade 0), lower intertidal and estuaries (grade 1), upper intertidal and freshwater (grade 2), beaches and riverbanks (grade 3), and coastal forests and jungles (grades 4-5).Pies at nodes represent the estimated ancestral state with the outer ring indicating the pathway (at some nodes, both pathways are shown; when node is inferred marine, no pie is shown).Tip codings are based on estimates by family (Supplementary TableS6), with the collapsed clades showing the color that represents the largest slice of their prior probabilities (split in the case of equal probabilities for 2 grades).For clades that have a small number of taxa in a grade from the opposite pathway, a small triangle is added.Line drawings at right (numbers corresponding to taxa in Supplementary TableS7), by Javier Luque and Harrison Mancke.