Introduction
Genome size, the bulk amount of DNA in the genome, displays an extremely high variability across the Tree of Life (Gregory, 2021). Comparative studies of genome size variation revealed universal patterns of correlations between genome size and cytological traits across a variety of taxa, spanning from plants to animals. A positive correlation between genome size and nucleus/cell size has been reported, as well as a negative relationship between genome size and cellular division rate (Cavalier-Smith, 1985; Lynch, 2007; Gregory, 2021). When life-history and eco-physiological traits are considered, however, their correlation patterns with genome size variation often become unpredictable and universal trends are elusive. These loose correlations suggest that genome size is affected by multiple synergistic factors that do not act likewise across taxa (Grime and Mowforth, 1982; Smith and Gregory, 2009; Herben et al., 2012; Suda et al., 2015; Alfsnes et al., 2017; Lertzman-Lepofsky et al., 2019; Glazier, 2021). For this reason, comparative studies among low taxonomic levels (e.g., among families) can help to disentangle the factors guiding genome size evolution and adaptation.
Crustaceans are characterized by some of the most variable genome sizes among animals (Jeffery and Gregory, 2014), with values ranging from 0.14 picograms (pg) for the copepod Cyclops kolensis to 64.62 pg for the amphipod Ampelisca macrocephala (Gregory, 2021). Some studies focused on specific crustacean orders demonstrated a relationship between genome size and their ecological, physiological and morphological features. For example, genome size is related to body size, latitude, temperature, and water depth for the Amphipoda (Hessen and Persson, 2009; Jeffery et al., 2017b; Ritchie et al., 2017; Hultgren et al., 2018; Hancock et al., 2021) and to body size in the Ostracoda and Cladocera (Hessen and Persson, 2009; Jeffery et al., 2017a), while it relates to metabolic and developmental rates, body size and latitude in the Copepoda (Gregory et al., 2000; White and McLaren, 2000; Wyngaard et al., 2005; Hultgren et al., 2018). A consistent correlation pattern is yet to be found in one of the most specious order of crustaceans, the Decapoda.
The order Decapoda includes more than 14,750 extant species, most of which belong to the suborders Caridea (true shrimps), Anomura (hermit crabs), and Brachyura (true crabs) (De Grave et al., 2009). Decapods inhabit a broad variety of habitats, including oceanic waters, seafloors, coral reefs, mangroves and estuaries, intertidal muddy and sandy areas, freshwater streams and lakes, and even terrestrial ecosystems (Martin and Davis, 2001; De Grave et al., 2009). Their ecological plasticity is reflected in a wide range of morphological, physiological, reproductive, and behavioral traits. Indeed, the importance of decapods is unquestionable also from an economic point of view, with many species being extensively fished and commercialized worldwide (Bondad-Reantaso et al., 2012) or marketed as ornamental animals (Calado et al., 2003). These characteristics make the Decapoda an interesting model group in the study of crustacean evolution, including the evolution of genome size and complexity.
In Decapoda, genome size did not show a clear correlation with any of the life-history and eco-physiological traits investigated so far. In this group, genome size does not correlate with body size (Hessen and Persson, 2009; Hultgren et al., 2018), latitude (with the exception of caridean shrimps, Hultgren et al., 2018), or habitat (Hessen and Persson, 2009; Hultgren et al., 2018), although the categories of habitat considered thus far were limited to marine/freshwater or artic/temperate. Recent studies found that life-history traits, such as direct/indirect development and the number of larval stages, do correlate with variation in genome size (Hultgren et al., 2018, 2021). This correlation, however, was not present when the influence of phylogeny was accounted for and excluded from the analyses. An absence of correlation between genome size and life-history traits was also reported by Hessen and Persson (2009).
The high diversity in morphology, physiology, and ecology of Decapoda calls for a broad analysis to clarify the role of genome size in the diversification of this taxon. We hypothesize that multiple comparisons among species with different life-style characteristics and belonging to different families should lead to resolutive outcomes, as underlined by Rees et al. (2007). Here we tested possible relationships between the genome size and the ecological adaptations related to respiration, excretion, and reproduction of 179 species of Decapoda inhabiting marine, freshwater, intertidal and terrestrial habitats. The physiological traits and habitat categories were chosen ad hoc to represent the huge range of Decapoda environmental plasticity and were never tested before. Developmental mode and number of planktonic larvae, already tested in Hessen and Persson (2009) and Hultgren et al. (2018), were also considered to test whether a broader sample set would show any hitherto hidden pattern of correlation between these traits and genome size. We also provided new genome size data for three species of decapods using flow cytometry. To our knowledge, our dataset comprises all the species of Decapoda for which the genome size is currently available.
Materials and methods
Genome size data
Data on genome size were collected for 176 decapods from the Animal Genome Size Database (Gregory, 2021) and from literature search (Supplementary Table 1). To our knowledge, this dataset includes all the decapod species for which C-values are currently available. The species recorded belong to the suborders Dendrobranchiata (superfamily Penaeoidea) and Pleocyemata (infraorders Achelata, Anomura, Astacidea, Brachyura, and Caridea) for a total of 44 families (∼20% of all families of decapods). Genome size was expressed as haploid DNA contents (C-values, in picograms). For the few species we found multiple different genome size values, we calculated an average C-value when the difference among values was less than 2 pg. Otherwise, we selected the measures obtained most recently, through flow cytometry (FCM) or Feulgen densitometry (FIAD).
We also added new data on three brachyuran crab species, the sesarmid Chiromantes haematocheir and the ocypodids Gelasimus borealis and Paraleptuca splendida, by measuring genome size using flow cytometry. Hemolymph was collected from four individuals of C. haematocheir and three individuals of G. borealis and P. splendida in tubes containing anticoagulant buffer (0.1 M glucose, 30 mM trisodium citrate, 26 mM citric acid, 10 mM EDTA and 0.1M NaCl) in a 1:1 ratio. Erythrocyte nuclei from Gallus gallus were added to cell suspensions as an internal size standard and stained using propidium iodide at a final concentration of 50 ppm before incubation in the dark for 20 min at 4°C. The suspensions were run on a BD FACSMelody Cell Sorter (BD Biosciences, NJ, United States) and the relative fluorescence of nuclei peaks of interest were isolated using the software BD FACSDIVA v 6.1.3. C-values of the samples were estimated by comparing the fluorescence intensity of the target sample with the fluorescence intensity of the internal size standard.
Eco-physiological and life-history data
Taxonomy of all the species was checked on the World Register of Marine Species database (WoRMS Editorial Board, 2022)1, and no more valid species names or incorrect taxonomic classification were updated. For each species we collected data on habitat, respiratory, and excretory adaptations, developmental mode and, when available, number of larval stages. This information was retrieved from online databases (Legall and Poupin, 2002; WoRMS Editorial Board, 2022), both scientific and gray literature (see Supplementary Table 1) and knowledge by the authors.
We classified the habitats into five categories: freshwater, terrestrial, intertidal (i.e., the area comprises between high and low tide marks), shallow-mid marine waters (0–200 m), and deep sea waters (more than 200 m). The distinction between shallow-mid and deep waters was fixed at 200 m, i.e., the maximum depth of the photic zone (Thomas, 2019). The few species of the dataset living in brackish water were classified as shallow-mid marine water species.
In terms of respiration, species were divided into three categories based on their respiratory organs and on how they function: (a) gills only capable to breathe underwater (henceforth termed “gills”); (b) gills capable to extract oxygen also from air, for a limited amount of time, if kept wet and gills kept oxygenated through branchial chamber water recirculation (henceforth termed “wet gills”); (c) and branchiostegal lung (Supplementary Table 1). Gills are the most common respiratory organ in aquatic crustaceans and are found in representatives from all the superfamilies included in our dataset. Here haemolymph flows inside thin lamellae, located inside the branchial chamber, and exchange oxygen with the water constantly recirculated in the chamber by the movement of the scaphognathite, the leaflike exopodite of the second maxilla. “Wet gills” are typical of some brachyuran crabs that can stand out of water for relatively long periods of time provided that the gills and branchial chamber are wet. In these species, gills are smaller than the ones of marine species, with bulges on the surface of the lamellae, to prevent their collapse out of the water (Hsia et al., 2013). Among species possessing “wet gills”, some sesarmid crabs also show a reticulate pattern of setae on their branchiostegites. This pattern is used to recirculate water contained in their branchial chambers to maintain a high oxygen level of this water when the animal is emersed (Hawkins et al., 1982; Hawkins and Jones, 1982). Branchiostegal lungs are typical of terrestrial hermit crabs and brachyuran crabs: these species, in addition to the gills, developed a thick and well-vascularized epithelium, located inside an enlarged branchial chamber, through which air-breathing is possible (Farrelly and Greenaway, 1994; Paoli et al., 2015). Most of this linings tissue consists of large haemolymph-filled sinuses, bound by connective tissue. The outer cuticle and the epithelium lining this cuticle are both very thin and the blood-air diffusion distance results very short (Farrelly and Greenaway, 1994).
Regarding the excretory mode, we classified species as obligate ammonotelic or not strictly ammonotelic (Supplementary Table 1). Species belonging to the first category excrete nitrogen through ammonia while species falling in the second category are able to shift from ammonia to urea production for nitrogen excretion, depending on water availability and salinity (Weihrauch et al., 2004). Only one species included in our dataset (the coconut crab Birgus latro) is uricotelic (i.e., it excretes nitrogen mostly in the form of uric acid). In our statistical analyses this species was also classified as not strictly ammonotelic.
Finally, the developmental mode was categorized as direct (no larval stages) and indirect (presence of larval stages). For the totality of the species belonging to the three decapod infraorders of Anomura, Brachyura, and Astacidae, which develop through a series of zoea stages into the megalopa, we also retrieved the effective number of larval stages (Supplementary Table 1). This allowed us to categorize the developmental mode of these infraorders into direct (no larval stages), abbreviated (number of zoeae between 1 and 4), and extended (more than 4 zoeae).
Statistical analyses and phylogenetic correction
Firstly, we reconstructed the phylogenetic relationships among species included in our dataset retrieving taxonomic information from the Open Tree of Life2 using the package rotl v3.0.12 (Michonneau et al., 2016). The length of the branches of the tree were calculated according to Grafen’s method (Grafen, 1989), i.e., the height of each node is calculated as the number of leaves of the subtree minus one, with the package ape v5.6-1 (Paradis et al., 2004). Then we checked possible phylogenetic correlations of our dependent and independent variables (C-value, habitat, respiration, excretion, developmental mode, number of larval stages) using Pagel’s λ (Pagel, 1994). Pagel’s λ is a scaling parameter of phylogenetic dependence that varies continuously from 0 (the null expectation) to 1. λ = 0 means that the trait varies independently of phylogeny, while λ = 1 indicates that there is a strong phylogenetic signal. λ values slightly greater than one are also possible because λ is not a correlation but a scaling factor for a correlation and the value is subject to constraints determined by characteristics of the tree. The index was estimated using the phylosig function under the phytools package v1.0-1 and the fitDiscrete function under the geiger package v2.0.7 (Pennell et al., 2014) for the continuous and discrete variables, respectively.
We examined the association between the response variable (genome size, log-transformed) and the predictor categorical variables (habitat, respiratory adaptations, excretion mode, and developmental mode) running the phylogenetic generalized least squares (PGLS) model. To minimize the effects of wrong model selection, for each model, we explored the contribution of phylogeny by fitting λ simultaneously with the PGLS test (Revell, 2010). We fitted models with λ = 0 (phylogenetic independence) and λ = 1 (trait evolution according to Brownian motion). Models were fitted using the gls function of the nlme package v3.1-158 (Pinheiro et al., 2022) using corPagel as a correlation structure. Considering that the PGLS algorithm does not allow for unresolved branching, prior to run each model we removed polytomies within the dendrogram using the multi2di function under the phytools package v1.0-1. We selected the most informative models using a multimodel inference approach by means of the lowest Akaike’s Information Criterion (AIC) (Burnham and Anderson, 2004). We assessed the overall significance of each PGLS model by conducting a full vs. null model comparison. The overall P-value was obtained by comparing the full model–which is simply the model fitted–and the null model–which is a model comprising only the intercept. A significant P-value indicates that at least a predictor variable has a significant effect on the response variable (Mundry, 2014). We then ran ANOVA tests, using the anova function of the nlme package v3.1-153, to identify the effect of each specific predictor variable within the model. Post hoc pairwise comparisons were applied to uncover differences among each pair of means using estimated marginal means applying Bonferroni correction using the function emmeans in the package emmeans v1.7.4-1. PGLS models were run both using the complete Decapoda dataset and on a subset including species belonging to the infraorder Anomura, Brachyura, and Astacidea (∼57% of the complete dataset). We additionally focused our analyses on these three infraorders as they display the widest range of physiological and ecological traits within Decapoda. As specified above, for this subset we also collected the effective number of larval stages. To assess the relationship between genome size (log-transformed) and number of larval stages, transformed into log10(X + 1), we ran the phylogenetically independent contrasts (PICs) (Felsenstein, 1985). The analysis was implemented in the ape package v5.6-1. Summary statistics of the two datasets are reported in Supplementary Table 2. Dendrograms and boxplots were built using the packages ggtree v3.3.1, ggnewscale v0.4.6, and ggplot2 v3.3.5 in R v4.1.2.
Discussion
The major finding of this study is that genome size variation in decapods is mainly explained by their developmental mode. We found that genome size is significantly larger in decapods that have a direct development than in species with an indirect development. This result was confirmed both when the analyses were run on the complete dataset and on the subset including only the infraorders Anomura, Astacidea, and Brachyura. Our analyses also proved that in the latter subset genome size variation is also explained by the habitat, the respiration and excretion mode of the analyzed species.
The link between developmental modes and genome sizes shows a rather homogeneous trend across all metazoans. A significant association between large genomes and developmental complexity has been described for several species (Gregory, 2002). Organisms with complex developments, such as the ones for instance characterized by many larval stages or going through metamorphosis, have time-limited developmental windows to rapidly complete single or successive metamorphoses. This process requires rapid cellular division and differentiation, with direct implications on the amount of DNA that can be maintained (Gregory, 2002). On the other hand, developmental windows are less stringent for those species showing an abbreviated or direct development, allowing for genome expansion (Gregory, 2002). In crustaceans, this hypothesis was confirmed for copepods, where species with narrow temporal windows between successive larval stages present small genomes (Wyngaard et al., 2005). Even in insects, holometabolous species (i.e., the species with the most complex developmental mode) present smaller genomes than hemimetabolous ones (Gregory, 2002; Alfsnes et al., 2017). Similar patterns have been also shown within the Agnatha, with hagfish, which show direct development, presenting larger genome sizes than lampreys, characterized by extended development with metamorphosis (Hardie and Hebert, 2004). In amphibians, smaller genomes are found in frogs, characterized by a more complex development in comparison to salamanders, which are obligate pedomorphs (i.e., they do not undergo metamorphosis) and possess larger genomes (Bonett et al., 2020).
Interestingly, as far as the phylogenetic signal is concerned, when the dataset was reduced to just three infraorders (subset), we recorded a significant reduction in the lambda value of the dependent variable (genome size). In this case, the phylogenetic signal was mainly associated with the independent variables. With such results, ordinary least-squares (OLS) models are found to perform better than PGLS models (Revell, 2010). Our results are coherent with this evidence, since, for the subset, the best-fitted models did not include phylogeny as a covariate.
Overall, the results of our analyses on the subset show that multiple eco-physiological and life-history traits correlate with genome size in Anomura, Astacidea, and Brachyura. The best-fitted model showed that respiration mode, together with developmental mode, principally explains genome size variation in these three decapod infraorders. Specifically, respiratory adaptations typical of terrestrial and intertidal species, such as branchiostegal lungs or gills capable to extract oxygen from air, explain for a reduction of genome size.
Also the excretion mode contributes to explain the genome size variation in this subset of species. Species that are not strictly ammonotelic possess larger genomes than strictly ammonotelic species, in accordance with evidence that organisms with a broad ecological tolerance have large genomes (Bennett, 1987; Beaton and Hebert, 1989). This finding is in line with studies on genome size evolution in fish, which showed that freshwater and other eurybiotic species have larger genomes than marine and stenobiotic species (Ebeling et al., 1971; Hardie and Hebert, 2004).
Habitat has been pointed out as one of the main factors driving the evolution of genome size in many species (Grime and Mowforth, 1982; Dufresne and Jeffery, 2011; Herben et al., 2012). We found that this predictor significantly affected the genome size of Anomura, Astacidea, and Brachyura, even if less predominantly than the other variables considered. The Post hoc tests indicated that freshwater species have significantly larger genomes than species living in terrestrial and intertidal habitat and in shallow-mid waters. Previous studies aimed to compare freshwater and marine crustaceans found a weak association between genome size and these habitats, with slightly larger genomes found in marine species (Hessen and Persson, 2009; Alfsnes et al., 2017). Our results are yet in line with data on genome size evolution in fish, which indicate a clear correlation between the colonization of freshwater habitat and the possession of larger genomes (Hardie and Hebert, 2004; Smith and Gregory, 2009). The developmental mode of many crustacean species is related to habitat occupation and, for instance, the strong variations of hydrodynamic and physico-chemical parameters in freshwater habitats have led to an evolutionary trend toward direct development. On the contrary, direct and abbreviated development is very rare in marine species (Vogt, 2013).
All the chosen predictors, including habitat, showed a significant effect on the genome size of Brachyura, Astacidea, and Anomura, while the only variable accounting for genome size variation of all considered crustaceans was the developmental mode. In our opinion, this is due to the wider variation in habitat occupation and physiological adaptations occurring in these three infraorders. All these taxa include freshwater species, with Brachyura and Anomura being the only decapod infraorders able to colonize the terrestrial habitats, and they often show specific adaptations to cope with these challenging ecological transitions. Moreover, the reduction of the phylogenetic signal in the dependent variable, when only Anomura, Astacidea, and Brachyura are considered, suggests that accurate phylogenetic information (i.e., not only the one inferred from taxonomy as in this study) might lead to more solid results. This evidences point out areas for further research. We hypothesize that a more detailed study on a particular group with a well-resolved phylogeny and multiple transitions between habitats and/or developmental modes could shed further light on the evolution of genome size in decapods.
One of our expectations was that terrestrial and intertidal species could have a smaller genome size than marine species, in line with the large genomes constrained hypothesis (Knight et al., 2005). This hypothesis postulates that species with large genomes are under-represented in less-stable environments due to the cost of maintaining large amounts of DNA. Terrestrial and intertidal habitats represent de facto more variable environments to cope with than the marine realm, in terms of temperature and salinity unpredictability, weather extremes and water availability, for instance. This hypothesis has been confirmed for the most successful land dwellers among Crustacea, the Isopoda (Jeffery, 2015), whose terrestrial species have smaller genomes respect to the marine ones. In our subset, terrestrial and intertidal species display the smallest genomes. Post hoc tests, however, reveal a significant differentiation in their genome sizes with respect to freshwater species only, but not with marine ones. The lack of significance in our test can be attributed to our unbalanced dataset, as for the terrestrial and intertidal categories the number of entries was very small. This might have reduced the power of our statistical analysis. We registered, however, a significant reduction in genome size in species that possess branchiostegal lungs or gills capable to extract oxygen from air, which are both adaptations typical of terrestrial and intertidal species. The colonization of terrestrial and intertidal habitats led various taxa to evolve physiological adaptations to cope with challenging new environmental conditions (Little, 2009). Land invasion in decapods mainly concerned brachyuran and anomuran species, which evolved remarkable adaptations that involved locomotion, respiration, excretion, reproduction and, eventually, their microbiome composition (Burggren and McMahon, 1988; Greenaway, 2003; Cannicci et al., 2011, 2020; Watson-Zink, 2021). The observed trend in genome sizes could be the result of the genome restructuring that led to these adaptations and the selective pressures that was exerted by the water-to-land transition on the amount of DNA. Further studies, including more species living in intertidal and terrestrial habitats, may confirm our working hypothesis.
In conclusion, our study shows that multiple eco-physiological and life-history traits correlate with genome size in Anomura, Astacidea, and Brachyura. As suggested by Gregory (2002), a strong relationship—with causation proceeding in both directions—occurs between the genome size and the phenotype of an organism, with several ecological and physiological factors exerting synergistic effects on this interaction. Whether genome restructuring involved the acquisition or loss of functional genes, spliceosomal introns, mobile elements or other genomic regions is yet to be investigated in decapods. Future research should couple genome size data with genome sequencing ones for target species, to underpinning the role of differential genome content in this interaction and shed further light on this evolutionary process.