Morphological data were obtained from published taxonomic descriptions in which phyllosoma I was illustrated (see Supplementary Table S1 for full list of references). To avoid any taxonomic confusion, only descriptions based on larvae obtained from laboratory-reared ovigerous lobsters or collected in the field and identified by DNA barcoding, were used in this study. Moreover, those descriptions with low quality drawings, without a scale bar or without an upright dorsal/ventral view, were discarded. We restricted our study to the phyllosoma I because it is the larval stage for which more reliable taxonomic descriptions exist. As a result, a total of 37 species were retained for this study ( Figure S1 ), which cover the most ecological divergent species within Achelata.

Phyllosoma I drawings were digitized and the outlines of the cephalic shield including antennae were manually traced. We used the “polygon selection” option to trace drawings and the “spline fitting” option for smoothing outlines in ImageJ v. 1.52p (Schneider et al., 2012). Outlines were manually fitted to the nearest possible to the original drawings and saved as binary masks. Then, tpsDig v.2.31 (Rohlf, 2017) was used to measure the height and width of the cephalic shield, as well as the antenna lengths based on the original drawings.

To relate parental traits with the phyllosoma I shape, we collected data of the maximum body size and maximum distribution depth of the adult phase for each species. Data were extracted from Holthuis (1991; 2002; 2006), Chan and Yu (1993) and from published taxonomic descriptions of larval morphology, as mentioned before.

A total of 350 coordinates from the binary mask of each cephalic shield were selected for analyses. The middle point between the two antennae was considered as starting point of the outline ( Figure 1 ). The coordinates were sampled along the outline adapting the “stage2outline” and “stage2landmark” R codes of Wong et al. (2018) to our data.

We ran an elliptic Fourier analysis (EFA) (Rohlf and Archie, 1984) using the efourier function of the Momocs package (Bonhomme et al., 2014). After running Fourier power tests for each species’ shape, we decided to preserve the first twenty harmonics since they accounted for 99% of the cumulative variability of the shapes and closely fitted the original outlines ( Figure 1 ). Each harmonic has four coefficients; thus 80 coefficients were obtained from each phyllosoma I shape. However, we excluded the first three coefficients of the first harmonic to standardize the orientation, size, and rotation of outlines.

We used principal component analysis (PCA) to ordinate the Fourier coefficients. We plotted PCA scores and the major shape changes from the first two principal components (PC1 and PC2) to observe shape variation of phyllosoma I. To model the outlines on PCs, coefficients of harmonics were calculated from the product of the PC scores and their corresponding eigenvector using mean ± standard deviation (SD) of PC scores from the first two principal axes. Then, inverse Fourier transformations were used to rebuild the hypothetical shapes from the estimated coefficients. Shape differences regarding the average shape outline and the ± SD and ± 2SD of each axis were drawn within thin-plate splines

The correlations (R) and angles (θ = arccos(R)) between regression vectors of shape (PC1 and PC2) on linear measurements (height, width, aspect ratio = width/height, and relative antenna length = antenna length/height) were estimated using the rv.test function of the barnlarv package. R ranges from 0 to 1, and both represent orthogonal and parallel positions of the compared vectors, respectively.

We constructed a time-calibrated phylogenetic tree including 112 out of the 155 Achelata species described (last checklist, Chan, 2019) ( Table S2 ). We retrieved sequences of six genetic regions (12S, 16S, 18S, and 28 S ribosomal subunits (rRNAs), cytochrome oxidase c subunit I, and histone 3) from the 112 Achelata species and four outgroups: two Astacidae (Homarus americanus and Metanephrops thomsoni), and two Glypheidea (Laurentaeglyphea neocaledonica and Neoglyphea inopinata) species since they are sister suborders of Achelata (Schram and Koenemann, 2021). For those species with multiple available sequences for one gene, we selected the longest one. When necessary, we used the reverse complement of the original retrieved sequence or concatenated partial sequences of a specific genetic region (28S rRNA) ( Table S2 ). Each genetic region was independently aligned and then alignments were processed with Gblocks 0.9.1 (Castresana, 2000) to discard poorly aligned positions, but allowing the inclusion of gaps to increase the number of phylogenetically informative variable sites. The processed alignments were concatenated to form a final alignment of 3877 bp. Then, a time-calibrated phylogenetic tree was constructed using BEAST v.2.6 (Bouckaert et al., 2019). For this tree, we included two calibration points from fossil records: (1) for the Achetala clade, Yunnanopalinura schrami (241–247 Millions of years ago, Mya) (Feldmann et al., 2012); and (2) for the Palinuridae, Archaeopalinurus (210–221 Mya) (Pinna, 1974). To construct the tree, we used the best substitution model from BEAST model test (Bouckaert and Drummond, 2017) under all reversible model search, the Yule process to build the prior tree, and a fast relaxed Lognormal clock (Zhang and Drummond, 2020). The analysis was run with a Markov Chain Monte Carlo simulation of 2x10⁸ steps, storing 10000 trees from the posterior distribution. The maximum clade credibility tree (MCC) was calculated considering a 10% of burn-in with the program TreeAnnotator included in BEAST v.2.6. Then, the MCC tree was pruned to only contain 34 out of the 37 species from which phyllosoma I shape information is available (bold species, Supplementary Table S2 ) given the lack of DNA sequences for three species. The MCC tree of 112 species was used to study lineage diversification, while the MCC tree only containing 34 species was used in comparative phylogenetic analyses of phyllosoma I shape.

We assessed whether lineage accumulation has experienced a slowdown along the evolutionary history of Achelata following an early burst of diversification. We constructed a lineage-through-time (LTT) plot and evaluated the fit of the LTT curve to a model of diversification with constant rate through a Monte Carlo Constant Rate (MCCR) test using the phytools package (Revell, 2012). This test estimates the γ statistic for incompletely sampled phylogenies by contrasting the distribution of inter-node distances between the root and its time-calibrated midpoint to the distribution of distances between this midpoint and the tips (Pybus and Harvey, 2000). Negative values of γ indicate that inter-node distances between the root and midpoint are shorter at early history than those from midpoint to the tips, suggesting that most branching events occurred in early history of Achelata. In contrast, if lineage diversification follows a constant rate process, γ has a mean value of 0. It has been demonstrated that incomplete taxon sampling in the phylogeny increases type I error rates in diversification analyses. Hence, the MCCR test calculates distributions of γ through simulations of phylogenies with taxon sizes equalling the known number of species (~155 species in Achelata) (Chan, 2019) under the null hypothesis of diversification with constant rate. Species within the simulated trees were randomly pruned to equal the number of species with retrieved sequences (112), replicating the incomplete taxon sampling. The LTT plot was built for each simulated pruned tree to construct the confidence interval.

To complement the LTT plot and MCCR test, we fitted six lineage accumulation models to branching times extracted from the time-calibrated MCC tree. Two of the models assumed constant rates of diversification: the pure-birth (Yule) model that does not incorporate extinctions, and the constant rate birth-death (crBD) model that includes extinctions but with a constant rate for speciation and extinction along evolutionary history (Nee et al., 1994). The next two models were Yule models with one (Yule2R) and two (Yule3R) shifts in speciation rate at some point along with the evolutionary history. The last two models: density-dependent logistic (DDL) and density-dependent exponential (DDE) assume diversity-dependence and estimate diversification rates in relation to the accumulation of diversity over time while accounting for extinctions. All these models were fitted using the fitdAICrc function of the laser package (Rabosky, 2006). To determine the model with best fit, we considered the bias-corrected version of the Akaike Information Criterion (AICc) approach (Burnham and Anderson, 2002; Posada and Buckley, 2004) and chose the model with the lowest AICc score.

We evaluated whether phyllosoma I shapes, centroid size and antenna length are more similar between more related species, and less similar between distant species (i.e., phylogenetic signal). Phylogenetic signal was estimated in the extent of the multidimensional shape data using K_mult (Adams, 2014), which indicates whether shape variation is as arranged by phylogenetic position of species as expected under a Brownian motion (BM) model of evolution (Felsenstein, 1985) (K_mult = 1) or has higher (K_mult > 1) or lower (K_mult < 1) phylogenetic signal. K_mult and its significance were estimated using the physignal function of the geomorph package (Adams and Otárola-Castillo, 2013). In addition, we estimated phylogenetic signals of PC1, PC2, centroid size (CS), and antenna length through the λ estimator (Pagel, 1999), which indicates the level of correlation between species regarding the expected correlation under the BM model. The λ estimators and their significances were calculated utilizing the phylosig function of the phytools package. For visualization, a “phylomorphospace” was constructed using the first two PCs of the shape analysis as reference of morphological information. The phylomorphospace projects the phylogenetic tree into a space of morphological traits (Sidlauskas, 2008). For interpretation, the morphological diversity or “disparity” among species was considered: on the one hand, when a subclade occupies broad regions of the morphospace, recent divergence considerably contributes to the overall “disparity”; on the other hand, when a subclade occupies small regions, its greater “disparity” is explained by earlier divergence, supporting an early burst of trait evolution.

We assessed whether shape divergences between species pairs show a tendency to stasis/convergence or rapid diversification. Small shape divergences between largely distant species suggest convergent evolution or stasis, while large shape distances between closely related species indicate rapid diversification. Thus, we built pairwise distance-contrast plots in the same way as Muschick et al. (2012). We plotted phylogenetic divergences against shape divergences of species pairs. Phylogenetic divergences were estimated using the cophenetic function, while shape divergences were calculated using the dist function in R. Then, we compared the observed shape divergences with the ones obtained from 1000 simulated datasets under the BM model. For these simulations, we used the ratematrix and sim.char functions of the geiger package (Harmon et al., 2008). The simulated divergences under BM were extracted from the observed ones. Negative values were obtained when the observed species pairs comparisons were more similar than expected due to their phylogenetic distance, and positive values when the observed species pair comparisons were less similar in the simulated dataset. For visualization, the plots were color-scaled using the differences between the observed and simulated shape datasets. To account for variation among simulations, we counted the number of species pairs with lower observed shape divergence than the expected under BM in 95% of the simulations and plotted the observed and simulated Procrustes distances with a grey-scale colouration based on the number of simulations in which the observed Procrustes distances were lower than the simulated ones for each species pair.

We fitted phylogenetic multiple regression models to predict phyllosoma shape in relation to body size (here CS), antenna length, and maximum distribution depth and maximum size of adults, while accounting for the phylogenetic position of species through the procD.pgls function of the geomorph package. This function performs regression models in a phylogenetic context under a BM model, accommodating high-dimensional datasets. In addition, since the main trajectories of shape variation (PC1 and PC2) were correlated with body size, antenna length, and aspect ratio (see Results section), we fitted univariate regression models to evaluate whether these variables were independently correlated to maximum size and maximum distribution depth of adults while considering phylogenetic position of species. For these univariate models, we log-transformed the variables and specified an expected covariance under the BM model.

We explored how shape disparity has changed over the Acheleta’s evolutionary history. We used the ddt function of the geiger package to calculate and plot disparity-through-time (DTT). The DTT plot exhibits the contribution of subclades to the total disparity as the linage diversity of the clade increases (Harmon et al., 2003). The observed DDT is usually contrasted by the expected value under the BM model, thus the difference between the observed and median expected values is called the morphological disparity index (MDI). Negative values of MDI are indicative of an early burst of disparity in trait evolution within the clade, while positive values suggest a higher variation than expected by phylogeny. Finally, we estimated the ancestral states of PC1 and PC2 to visualize changes at specific trajectories of shape along with the evolutionary history of Achelata. We estimated the ancestral state using the fitContinuous of the geiger package and superimposed them over the time-calibrated phylogenetic tree with the contMap function of the phytools package.

The first two PCs explained 76% of the total variation of phyllosoma I shapes ( Figure 2 ). Species were divided into two clusters located in the negative and positive extremes of PC1. For PC1, species with negative values had large antennae and pear-shaped cephalic shields, while species with positive values had small antennae and wider rounded cephalic shields ( Figure 2 ). For PC2, positive values were characterized by a wider anterior cephalic shield and antennae, while negative values represented narrower anterior cephalic shields and antennae ( Figure 2 ). Few species with intermediate shapes appeared far from the two main clusters located at extremes of PC1. For instance, the shape of Panulirus argus, which was the closest one to the mean shape configuration ( Figure 2 ). PC1 was significantly correlated with relative antenna length and presented a moderate non-significant correlation with aspect ratio. PC2 only showed moderate non-significant correlation with body width ( Table 1 ).

The maximum clade credibility tree of the Achelata species ( Figure 3 : 34 species, Supplementary Figure S2 : 112 species tree) was congruent with previous multi-locus phylogenies of Achelata (Palero et al., 2009; Tsang et al., 2009; Yang et al., 2012). The monophyly of Stridentes and Silentes groups within Palinuridae, and monophylies of subfamilies within Scyllaridae were also supported in these trees except for Ibacinae, which was paraphyletic. Ibacinae genera were placed next to the subfamily Arctidinae. Also, our tree supported the polyphyletic state of genera within Scyllarinae subfamily observed in a previous phylogeny (Tsang et al., 2009).

The MCCR test did not find significant evidence for bursts of lineage diversification in late evolutionary history of Achelata (γ = -2.3992, P = 0.138) as suggested by the LTT plot ( Figure 4 ), which showed two bursts of lineage diversification. These bursts occurred about mid- (ca. 160–125 Mya) and late-history (ca. 50 Mya afterwards) of the Achelata clade. However, the pattern of LTT was slightly supported by the fitting of branching times to a Yule model with two shifts in diversification rate (Yule3R) in late history over models of diversity-dependence diversification (DDL and DDX) dAICc_DDL = 2.7908 and dAICc_DDX = 5.497) and a Yule model with one shift in diversification rate (dAICc_YULE2 = 2.09).

The multivariate measure of phylogenetic signal was K_mult = 0.929, indicating that shape of phyllosoma I had a significant phylogenetic signal (P < 0.0002) under the BM expectation of 1. The K_mult value corresponded to a shape variation pattern in which closely related species’ shapes were more similar than expected under BM. However, the phylomorphospace revealed overlaps between some groups of Palinuridae and Scyllaridae species (see Figure 2 ). For instance, some Ibacinae (Parribacus antarcticus and P. caledonicus) and Arctidinae species (Arctides regalis and Scyllarides aequinoctialis) had large antennae and pear-shaped cephalic shield that resemble some palinurids. Among Scyllaridae subfamilies, Scyllarinae showed a characteristic shape (short antennae and wide rounded cephalic shield) absent in Arctidinae and Ibacine, yet like the ones present in Theninae ( Figure 2 ). Within Palinuridae, Silentes were marginally separated from Stridentes species. The broad occupation of palinurids in the morphospace suggested a more recent burst of shape variation in relation to the other subclades, while narrow occupation of Theninae and Scyllarinae indicated an early burst of morphological diversification within these subclades ( Figure 2 ). The phylogenetic signals of PC1 (λ_PC1 = 0.927, P < 0.0001), and antenna length (λ_ANT = 0.869, P < 0.0001) pointed out that correlation of each independent trait was closely linked to the expected correlation under BM. In contrast, phylogenetic signals of PC2 (λ_PC2 = 0.683, P = 0.099) and CS (λ_CS = 0.168, P = 0.413) indicated that variation of these traits may be independently determined by some external condition rather than the intrinsic phylogenetic position of species.

Pairwise-contrast plots exhibited shape convergence between phylogenetically distant species ( Figure 5 , red-scale points) and shape divergences larger than expected under the BM model ( Figure 5 , blue-scale points). Comparing our dataset of observed species-pair shape divergences with species-pair shape divergences simulated under the BM model, we found a total of 150 species pairs comparisons with lower observed divergence than expected under BM in 95% of the simulations ( Figure 5 , black points), indicating a significant shape convergence in species conforming these pair comparisons. A total of 74 out of these comparisons contained any species belonging to Scyllarinae or Panulirus species.

We found that body size and antenna length were significant predictors of multidimensional shape (body size: Z = 3.267, P < 0.002; antenna length: Z = 3.626, P < 0.002), while maximum distribution depth and maximum size at adult phase were not significant predictors. In case of the individual regressions of PCs, CS, and antenna length on adult variables, only the maximum distribution depth of adults was negatively correlated with PC2 (λ_PC2 = 0.419, P < 0.02), while the other correlations were not significant ( Supplementary Figure S3 ). This suggest that maximum depth of adults might be an evolutionary driver of shape variation in phyllosoma I.

The DTT analysis displayed a lower disparity than the median expected under BM from early to middle (Late Jurassic) history of Achelata. Then, disparity increased over the median expected under BM along with the whole Cretaceous, and then slowed down until mid-Paleogene, from which it experienced a positive peak (see Figure 4 ). The nonsignificant MDI value (0.007, P = 0.926) indicated that the DDT curve was closely fitted to the expectations of the BM, that in this case was translated into a compensating trend of disparity along with the whole evolutionary history of Achelata. The ancestral state reconstructions of PC1 and PC2 showed that the common ancestor of Achelata might had presented PC scores (PC1 = -0.0037, PC2 = -0.0006) whose shape closely corresponded with the shape of Panulirus argus ( Figure 6 ). This species had a phyllosoma with medium-size antennae and pear-shape shield ( Figure 2 ). In addition, the ancestral reconstruction indicated that main changes in states of PC1 occurred after the divergence of Palinuridae and Scyllaridae. The main changes in PC1 were concentrated at early history of Scyllarinae and Silentes, while later variations took place within Panulirus (Stridentes, Palinuridae) and Parribacus (Ibacinae, Scyllaridae). PC2 remained more stable throughout evolutionary history, the main changes in this component of shape occurred in Puerulus (Stridentes), Ibacinae and Arctidinae ( Figure 6 ). Scyllaridae species retained PC2 scores close to the common ancestor of Achelata.

Our findings partially support an association of phylogenetic and shape divergences between Palinuridae and Scyllaridae. As expected, phyllosoma I with more rounded cephalic shields and short antennae mainly belong to Scyllaridae, whereas Palinuridae grouped those larvae with more pear-shaped shields and large antennae. These shape differences were reflected in the high phylogenetic signals detected in the multivariate shape, main shape trajectory (PC1) and antenna length. At the family level, these shapes have been traditionally used to discriminate scyllarid from palinurid phyllosomata collected in the field (Santos and González-Gordillo, 2004). These differences in phyllosoma shape are even more evident in more developed larval stages (Baisre, 1994).

Based on our shape data, in combination with other phyllosoma traits previously used (Baisre, 1994), fossil records, distribution ranges, and genetic information (George, 2006; Palero et al., 2009; Yang et al., 2012), we can mention some insights about the phylogenetic significance of the shape of phyllosoma I. The shape data could not separate Arctidinae and Ibacinae of Scyllaridae from Palinuridae species, reflecting a discordance between taxonomic and morphological diversity. It occurred because the antennae in these Scyllaridae subfamilies are as large as those present in Palinuridae and share the biramous state of this trait with Silentes species (Jasus, Palinurellus, and Sagmariasus). Interestingly, most scyllarid featuring this shape hatch in an ontogenetic state before phyllosoma, which is named naupliosoma and is also present in palinurid lobster (Webber and Booth, 2007). The absence/presence of this ontogenetic state across Achelata families might indicate that hatching as phyllosoma is a derived state. The absence of naupliosoma has been related to shortened development (Baisre, 1994), and their association with phyllosoma larvae I with short antenna and wider cephalic shield might indicate a developmental link between the state at hatching and the two main phyllosoma shapes described here.

Within Scyllaridae, the shape variation of phyllosoma I supported the sister relationship of subfamilies Theninae and Scyllarinae observed on traits from first and last phyllosoma instars (Baisre, 1994), adult traits (Yang et al., 2012), and genetic information (Palero et al., 2009; Yang et al., 2012, this study). In addition, the pear-shaped cephalic shield and large antennae found in Parribacus (Ibacinae) and Arctidinae are concordant with previous evidence to place Parribacus closer to Arctidinae rather than to other Ibacinae genus (Baisre, 1994; Booth et al., 2005). However, this closer morphological relationship between Parribacus and Arctidinae was not observed in our phylogeny including most of species from these subclades but reported in a previous phylogenetic study (Yang et al., 2012).

In the case of Stridentes within Palinuridae, the largest genus Panulirus occupied a wide region of the morphospace, indicating broader ecological speciation. However, shape data of phyllosoma I did not separate this genus into the two main groups as previously observed in phylogenetic analyses (Ptacek et al., 2001) and morphological characters of adults and phyllosoma (Patek and Oakley, 2003; George, 2006). Interestingly, the reconstructions of the ancestral state of PC1 and PC2 suggest that the common ancestor of Achelata (PC1 = -0.0037, PC2 = -0.0006) had a similar shape to Panulirus argus (PC1 = -0.0069, PC2 = 0.0001) ( Figure 6 , Figure 1 ), the hypothetical basal species of this genus (George, 2006). This might indicate a stasis of this shape in ancestors of Panulirus and that further shape variations in this genus occurred later.

In addition, shape data provides more evidence of phyllosoma similarity between Jasus and Puerulus suggesting a phylogenetic significance of the shape. Baisre (1994) has already proposed the common ancestry of these genera. Their phyllosoma have biramous antennae, but in comparison to scyllarid species with the biramous state, in Jasus and Puerulus the outer branch is much shorter than the inner one ( Supplementary Figure S1 ). In Puerulus, the antennae are bifurcated near the basal part, while in Silentes, bifurcation occurs in the distal part. The biramous antennae are conserved within Silentes, but this state is alternated with the uniramous state within Stridentes. Thus, uniramous antennae have been observed in Panulirus, Palinustus (final stage, Palero et al., 2010), Justitia (J. longimana stage VII, Konishi et al., 2021) and Nupalirus (N. japonicus stage VIII, Konishi et al., 2021); and biramous in Puerulus, Palinurus (stage I, Palero and Abelló, 2007), and Linuparus (L. sordidus stage IV?, Inoue et al., 2001).

The alternance of both states in Palinuridae and Scyllaridae, and the divergence times of subfamilies and groups within each family indicate that both states of antennae emerged multiple times across Achelata. Given the contrasting variation of a trait such as antenna within this clade, which plays an important role in food selection and predator detection, further analysis assessing variation of antennae among species would shed light on ecological changes over the evolutionary history of Achelata.

We expected that shape of phyllosoma should be more influenced by the embryonic phase, directly related to the adult benthic habitat and mother condition. Among arthropods, larger females frequently give rise to larger offspring because of their more efficient nutrient provision than those of smaller mothers, but this relationship is often weak and is certainly not ubiquitous (Oliphant and Thatje, 2021). For instance, in the European lobster Homarus gammarus, belonging to Astacidea, mean larval size at hatching is closely linked to maternal size (Moland et al., 2010). In our study, we did not find a relationship between the maximum adult size and the size of phyllosoma I, suggesting that this relationship, if any, may occur more frequently at intraspecific level or that any existing correlation between adult and larvae size is accounted for the phylogenetic position of species. However, our analyses suggest that maximum distribution depth of adults is an ecological driver of shape variation. This adult trait has a negative correlation with PC2 ( Supplementary Figure S3 ). This correlation shows that variation of PC2 scores occurred across different groups going from shallow to deep waters. High PC2 scores were present in Parribacus in shallow waters while small PC2 scores occurred in Ibacus and Puerulus, both genera occupying waters deeper than 300 m. High PC2 scores are related to species possessing a phyllosoma with wider anterior cephalic shields and large antenna, while low PC2 scores are associated with species with narrower anterior cephalic shields and smaller antenna. Phyllosomata of deep-water species must travel longer distance upwards to reach the upper layers of the water column, being exposed to higher predation risk. Longer relative antennae should facilitate tracking a bigger volume of surrounding water, increasing the perception of mechanical stimuli as happens in other planktonic organisms (i.e., copepods, Kiørboe, 2008). Considering that predation is a major evolutionary driver in Achelata (Spanier and Weihs, 1990), environments with higher predation pressure may have favored the selection for larger antennae. This negative correlation is also concordant with the multiple transition from shallow to deep waters found in analysis of scyllarid lobsters (Tsang et al., 2009; Yang et al., 2012) since Arctidinae that is distributed in shallow waters, has the oldest fossil records within Scyllaridae (Scyllarides punctatus, 120 Mya at mid-Cretaceous, Woods, 1925; Webber and Booth, 2007). On the other hand, the correlation contradicts the deep-shallow water transition previously proposed for palinurid lobsters (George, 2006) and documented for other decapods (Rodríguez-Flores et al., 2022) and marine invertebrate groups (Jablonski et al., 1983).

Adaptive radiations are characterized by a pattern of early lineage, morphological, and ecological diversification mainly triggered by ecological opportunities (Simpson, 1949) and followed by subsequent asymptotic declines in diversification rates over time (Schluter, 2000; Harmon et al., 2010). Our evidence suggests that Achelata did not experience an early burst of diversification and neither present a decline in lineage diversification at recent evolutionary history, rather it passed through two bursts of lineage and morphological diversification at mid and late evolutionary history ( Figure 4 ). A similar pattern of nonearly radiation was also found in analyses of body plans of adults including Achelata and other groups within Decapoda (Schweitzer and Feldmann, 2014). The bursts of lineage diversification were concordant with the fit of branching times to a Yule model with shifts in speciation rate across time (Yule3rate) over models of diversity-dependence diversification (DDL, DDE), which are expected in clades diversifying through adaptive radiation processes (Schluter, 2000).

Nonelection of a diversity-dependence model is supported by low competition for resources between species as suggested by: (1) the scarce or absent geographical overlap and asymmetric depth distribution of Palinuridae species (George and Main, 1967; George, 2006) and some Scyllaridae genera; and (2) morphological specialization of some scyllarid clades that allow them to exploit a variety of different environments (Webber and Booth, 2007). Allopatric distribution and reduced competition joined with the two bursts of lineage diversification precluding long-temporal bursts of morphological disparity (over millions of years) indicate that they were produced through geographic isolation after punctual episodes of colonization followed by subsequent subtle ecological speciation (Rundell and Price, 2009; Czekanski-Moir and Rundell, 2019). Hence, lineage diversification of Achelata has been primarily triggered by nonadaptive radiation with minimal ecological diversification (Gittenberger, 1991; Gillespie, 2004; Rundell and Price, 2009). This mode of radiation explains the presence of some groups of closely related species with subtle ecological divergences but clear reproductive isolation within Palinuridae (see cases of Jasus (Brasher et al., 1992), and Panulirus (Chan and Chu, 1996) and Scyllaridae (e.g. Thenus (Burton and Davie, 2007), Galearctus and Petrarctus (Yang et al., 2012). Nonadaptive radiation is concordant with the absence of a diversity-dependence slowdown in lineage diversification at the end of the evolutionary history. This last finding also indicates that Achelata has not reached equilibrium yet and continue diversifying, in agreement with the highest diversity levels of reported at Holocene for this clade (Schweitzer and Feldmann, 2014).

Going into detail, the first burst of lineage diversification was dated around the transition from Jurassic to Cretaceous (ca. 160–80 Mya) as supported by previous estimates from fossil records (Schweitzer and Feldmann, 2014). It occurred posterior to the beginning of the major phase of Gondwana breakup (around 180 Mya) that separated Africa, South America, India, Antarctica, and Australia and allowed the widening of the Atlantic Ocean (Mcloughlin, 2001). As expected in a clade whose diversification follows a nonadaptive radiation mode, the formation of isolated lineages preluded a posterior increase of morphological disparity over millions of years as a slow response to the environmental conditions presented for each separate lineage (Gittenberger, 1991). The DDT plot indicates that morphological disparity remained higher than expected under the BM model across the whole Cretaceous, in which an episode of radiation occurred as suggested by diversity estimations from body plans (Schweitzer and Feldmann, 2014) and fossil records of several extinct genera (Astacodes, Cancrinos, Palinurina, and Palaeopalinurus) succeeding Yannanopalinurus schrani, the earliest known Achelata (Schram and Koenemann, 2021). Interestingly, this episode of radiation is coincident with changes in early ontogenetic stages of Achelata visualized in fossil records of Palinurina tenera possessing additional instars that are absent in extant species (Lavalli and Spanier, 2010), and variations in ancestral states of the main trajectories of shape principally occurring in ancestors of Stridentes, Silentes, Scyllarinae and Theninae ( Figure 6 ).

The second burst of lineage diversification started in Late Cretaceous towards present ( Figure 4 ) and is concordant with a major radiation previously reported (Schweitzer and Feldmann, 2014). Following the allopatric speciation hypothesis proposed by George (2006), this burst of lineage diversification occurred as consequence of tectonic changes that separated marine basins and then gave rise islands, providing new available environments for colonization and the subsequent second burst of morphological diversification ( Figure 4 ). This burst might be related to changes in distribution depth as shown by variations in ancestral reconstruction of PC2 scores ( Figure 6 ) and their negative significant correlation with the maximum distribution depth of adults ( Supplementary Figure S3 ). Hence, colonization of different depths (i.e., new habitats) by ancestors might have represented an ecological opportunity that allowed lineage and ecological diversification (Simpson, 1949; Yoder et al., 2010). In addition, for this second burst, we cannot discard an ecological diversification of Achelata prompted by the presence of ecological opportunities due to prolonged extinction events around the Cretaceous-Paleogene limit that affected several decapod groups (Schweitzer and Feldmann, 2014). This is because one of the extinct groups that left available environments for colonization is Cancrinidae, whose members were present from Early Jurassic to Late Cretaceous (Schweitzer and Feldmann, 2014; Schram and Koenemann, 2021) and are considered the intermediate between palinurids and scyllarids (Haug et al., 2016; Schram and Koenemann, 2021).

Given the low shape divergences between largely distant species ( Figure 5 ), their overlap in the phylomorphospace (Palinuridae with Arctidinae and Ibacinae) ( Figure 2 ), and the conserved shape (phyllosoma I with biramous large antennae and pear-shaped cephalic shields) present in the old groups from each respective family (Scyllaridae: Scyllarides punctatus, 110–120 Mya (Woods, 1925); Palinuridae: Silentes (163 Mya, median estimated age)) ( Figure 3 , Supplementary Figure S2 ), it is likely that this shape might had also been present in the extinct Cancrinidae lobsters despite they possessed intermediate states of adult antenna forms between scyllarids and palinurids (Haug et al., 2016). However, they were intermingled with an allopatric speciation mode and subtle morphological differences (i.e., are attributes of a nonadaptive radiation, Gittenberger, 1991), not causing a substantial increase of the slow lineage diversification occurred over millions of years ( Figure 3 ). This highlights the possibility that components of both types of radiation processes could have contributed to diversification of Achelata as reported in other invertebrate groups (Gillespie, 2004; Cotoras et al., 2018).

Nevertheless, as explained above some strong arguments for the dominance of nonadaptive radiation in diversification of Achelata are the degree of specialization of scyllarid species (Scyllarinae and Theninae), the overlap within the phylomorphospace ( Figure 2 ), the pairwise-distance plots showing lower shape divergence than expected under the BM model ( Figure 5 , black points); the shape convergence in 150 species pairs ( Figure 5 , dark red points), the absence of a slowdown in lineage diversification and rather the highest diversifying trend at recent times (Schweitzer and Feldmann, 2014). However, under a hypothetical scenario of complete adaptive radiation, the possibility of geographic overlap between some species of Scyllaridae (Webber and Booth, 2007) and their shape similarities, indicate low levels of character displacement (Brown and Wilson, 1956; Schluter, 2000). This refers to accentuated phenotypic differences as response of competition for resources between sympatric species (Brown and Wilson, 1956), being a clear signal of slowdown in lineage diversification as consequence of diversity-dependence accumulation (Monroe and Bokma, 2017). However, the low level of character displacement detected here, and absent signature of lineage-diversification slowdown may indicate only a marginal contribution of some components of adaptive radiation in the dominant nonadaptive radiation processes occurred along with the diversification of Achelata lobsters.

It is a fact that diversification of early life stages (e.g. larval stages) is rarely considered in evolutionary studies. Like our study in slipper and spiny lobsters, others in echinoderms (Raff and Byrne, 2006), barnacles (Wong et al., 2018), and sparid fishes (Colangelo et al., 2019) indicated that morphological novelties present in adult stages are related to conspicuous morphologies originated since early life stages that allow species to exploit different environments and to reduce interspecific competition. Thus, these morphologies not only play a crucial role within a lifespan of a single specimen, but rather along with the evolutionary history of a clade. In this sense, we encourage the scientific community to further investigate the diversification of larval forms to increase our understanding of the evolution of target groups.