Euchelicerata is a successful clade of arthropods including Xiphosura, Chasmataspidida, Eurypterida, and Arachnida. Xiphosura is a group with extant representatives, which has an extensive fossil record dating back to the Ordovician (Rudkin et al., 2008; van Roy et al., 2015; Lamsdell et al., 2023). The earliest described xiphosuran fossil remains come from the Williams member of the Stone Mountain formation of Manitoba, Canada which dates to the latest Ordovician at c. 443Ma (Rudkin et al., 2008), but recent fossil discoveries suggest they first evolved in the early Ordovician (van Roy et al., 2015). Molecular clock estimates suggest a late Cambrian origin for the group (Lozano-Fernandez et al., 2020). Representatives of Xiphosura have long been referred to as “living fossils” (Stoermer, 1952) with major morphological traits seemingly unaltered by the ravages of time over hundreds of millions of years. Even recently, they have been cited as an example of extreme morphological conservatism (Bicknell and Pates, 2020). The term “living fossils” is somewhat problematic as the subtext of the term implies a lack of evolution taking place in the group, whereas it is well established that broad-scale evolutionary stasis results from gradual evolutionary changes around a relatively static morphological average position through time (Simpson, 1944; Eldredge et al., 2005; Tëmkin and Eldredge, 2015) (Figure 1). While a certain degree of morphological conservatism is recognized in Xiphosura (Bennett et al., 2018), especially in late Mesozoic and Cenozoic forms (Avise et al., 1994; Rudkin and Young, 2009; Kin and Błażejowski, 2014; Lamsdell and McKenzie, 2015; Bicknell et al., 2019b), most late Paleozoic and early Mesozoic forms are considered to go through a much more pronounced evolutionary exploration of morphological space (Lamsdell, 2016; Bicknell, 2019; Bicknell et al., 2019a; Bicknell et al., 2020; Lamsdell, 2021a; Lamsdell, 2021b; Lustri et al., 2021; Bicknell et al., 2022). Freshwater colonization during the late Paleozoic resulted in xiphosurans adapting to many new habitats, possibly on multiple occasions, and is associated with the first record of remarkable radiation of the group in the fossil record (Lamsdell, 2016; Lamsdell, 2021a; Bicknell et al., 2022).

Heterochrony has been proposed as the main pattern to described the freshwater invasion of early xiphosuran forms (Lamsdell, 2021a; Lamsdell, 2021b). To understand these developmental patterns the following Paleozoic taxa are key: Alanops magnifica (Racheboeuf et al., 2002) and the Euproops complex, including the species nicknamed “Piesproops” (more formally Andersoniella sp.) (Lamsdell, 2020), Euproops danae and Euproops rotundatus (Haug et al., 2012; Haug and Rötzer, 2018b; Tashman et al., 2019; Haug and Haug, 2020). Insights provided by these taxa, alongside studies of the development of extant species (Scholl, 1977; Jegla and Costlow, 1982; Sekiguchi et al., 1988; Shuster and Sekiguchi, 2003; Haug and Rötzer, 2018a), has allowed the recognition of peramorphic and paedomorphic patterns in two of the four main clades of Xiphosura, Austrolimulidae and Belinuridae (Bicknell, 2019; Lamsdell, 2021a; Bicknell et al., 2021b; Lustri et al., 2021). Those studies critically enhanced our understanding of evolution and development, but less has been done to compare the developmental patterns of xiphosurans with other euchelicerates inhabiting the same aquatic realm.

Arthropods are generally characterized by an extreme specialization of their different life stages (Minelli et al., 2016). However, this is usually not the case for euchelicerates including xiphosurans (Haug and Rötzer, 2018a; Kaiser and Schoppe, 2018). Euchelicerates are characterized by a gradual, direct and usually epimorphic development, where there is no addition of segments after hatching (Sekiguchi et al., 1988; Braddy, 2001; Lamsdell and Selden, 2013; Haug, 2019; Lamsdell et al., 2019; Fusco and Minelli, 2021).

In order to provide a phylogenetic context to the study of xiphosuran development, it is essential to compare them with other euchelicerate groups. Modeling of evolutionary scenarios needs to account for the phylogenetic relationships between organisms as this is the only independent way to estimate rates of evolution (Garamszegi, 2014). “Stasis is generally defined as little or no net accrued species-wide morphological change during a species-lineage’s existence up to millions of years” (Eldredge et al., 2005, p. 133), yet, it is important to define what exactly “little or no” means (Eldredge et al., 2005). The only way to do so is to compare xiphosurans with other, related groups inhabiting the same environment. It is certainly very clear that the time scale matters here as what may appear to be static over millions of years may disguise a great deal of change around a mean when viewed at higher temporal resolution. Conversely, stasis at high temporal resolution may miss larger and gradual temporal trends only observable when a longer view is taken. To understand the evolutionary dynamics of any lineage, a diversity of temporal views must be taken and then contrasted to related lineages.

The inclusion of other euchelicerate groups such as eurypterids and chasmataspidids in the analyses helps to refine not only the estimation of developmental parameters at the root of Xiphosura but also the possible correlations of different evolutionary scenarios with the paleoenvironment independently from the phylogeny. Eurypterids and chasmataspidids shared similar environments with the horseshoe crabs during the Paleozoic (Dunlop, 2010; Howard et al., 2020). The development of xiphosurans has been recently explored by meta-analyses (Lamsdell, 2016; Lamsdell, 2021a; Bicknell et al., 2022), but the development of eurypterids and chasmataspidids has never been incorporated in such analyses. Data for eurypterids and chasmataspidids are also available, and research has focused on fine detailed analyses of the development of single species such as Hoplitaspis hiawathai (Lamsdell et al., 2019) and Eurypterus lacustris (Ruebenstahl et al., 2021).

In this work, a meta-analysis is presented of morphometric developmental data from eight species of Xiphosura, combined with data from one species of Eurypterida and one species of Chasmataspidida. These data have been utilized alongside environmental and phylogenetic information to perform an ancestral state reconstruction analysis for the allometric growth patterns and environment of Xiphosura. The influence of different environments on the evolution of development are then tested within a phylogenetic framework.

The slopes of the prosomal shield length and width of Eurypterus lacustris, Hoplitaspis hiawathai, Prolimulus woodwardi, Euproops sp., Paleolimulus kunguricus and Limulus polyphemus are available together with the slopes gathered from the literature of Paleolimulus signatus, Euproops danae, Mesolimulus walchi, Tachypleus tridentatus in Table 2. Gradient values greater than 1 represent a preferential growth of length over width, with higher numbers representing a more extreme allometry. A gradient of exactly 1 represents ontogenetic isometry (inflationary growth), while a gradient of less than 1 represents width increasing quicker than length during ontogeny with lower numbers representing more extreme allometry. The prosomal shield slopes range from 0.719 in Hoplitaspis hiawathai to 1.02 in Euproops sp and Paleolimulus signatus. Between these extremes, Eurypterus lacustris, Euproops danae, Prolimulus woodwardi, Paleolimulus kunguricus, Limulus polyphemus, Mesolimulus walchi and Tachypleus tridentatus range from 0.888 in Euproops danae to 0.99 in Tachypleus tridentatus. All performed regressions are shown in Figure 5.

The results for the reconstructed ancestral state of the allometric growth of the prosomal shield, partitioned by environment, are summarized in the phylogenetic trees shown in Figures 6 and 7. Figure 6 represents the tree that includes only Xiphosura, while Figure 7 includes both Xiphosura and additional euchelicerate species (Eurypterus lacustris and Hoplitaspis hiawathai). In both trees, the lowest values of the slope (indicating width increasing quicker than length) are found in species associated with marginal environments, such as Hoplitaspis hiawathai and Euproops danae, while values closer to 1 are more commonly associated with marine settings.

The two evolutionary models tested on the two different phylogenetic trees (xiphosurans-only and euchelicerate trees) with the reconstructed character history show differences in fitting the data and different statistical support. For the xiphosurans-only tree the BM model, which does not account for environmental differences, outperforms the BMM model, which considers environmental effects (BM AICw = 0.676; and BMM AICw = 0.324, see also Table 3). The Log-likelihood Ratio Test for this model yields a p-value of 0.89. For the euchelicerates tree, the BMM model, which accounts for environmental differences, outperforms the BM model (BM AICw = 0.027; and BMM AICw = 0.973, see also Table 4). The Log-likelihood Ratio Test for this model shows a p-value of 0.00246. When accounting for environmental affinities, the evolutionary rate recovered for marginal species is higher than in marine species for both xiphosurans-only and all-euchelicerates analyses. However, this pattern is much more evident in the all-euchelicerates tree. Full data regarding the comparison of the two models alongside evolutionary rates recovered in marine and marginal environments are presented in Tables 3, 4.

Several patterns are present in the reconstructed ancestral state for the growth pattern of the prosomal shield (carapace) compared across different environments and systematic levels (Figures 6 and 7). The obtained regression slopes represent how much the carapace length increases in comparison to its width. A lower slope value indicates a smaller growth of the carapace length in comparison to its width. In the reconstruction of ancestral states performed on the xiphosurans-only dataset and tree, higher slope values are found on average in limulid species from marine environments. Isolated species of Paleolimulidae and Bellinuridae from marginal environments also exhibit high slope values (specifically Paleolimulus signatus and Euproops sp.) (Figure 6). On the other hand, lower slope values are found in the marginal environment with Euproops danae. This may reflect a certain degree of morphological plasticity associated with species inhabiting the marginal environment, as supported by previous research (Lamsdell, 2016; Lamsdell 2021a; Bicknell et al., 2022). It is further supported by a higher estimated evolutionary rate for the marginal environment (Table 3). However, when comparing the BMM, which includes the environmental variable as a potential correlate to evolutionary rates, with the BM that does not account for this, there is no significant support for the BMM over the BM (Table 3). Instead, a single Brownian motion model for all the xiphosurans included in the analyses fits the data better than two different Brownian models associated with the marginal and marine environments. This analysis yields a high p-value (Table 3), weakening the assumption that the evolution of xiphosuran development can be inferred solely considering the internal relationships of the group. While these results may be due to the absence of correlation between evolutionary rates and environments for the Xiphosura, other hypotheses can also explain this outcome. The sample size used for the analyses may have been too small, reflecting a lack of available fossil data. Another factor that may have contributed to the results of this analysis is the difficulty in discriminating between coastal, estuarian or freshwater environments for fossil specimens. This becomes clearer when examining the results obtained from the same analyses performed on a tree and dataset that includes the non-xiphosuran euchelicerates: Eurypterus lacustris and Hoplitaspis hiawathai. The inclusion of other euchelicerates increases the variability in carapace allometric growth, and provides a new perspective on the intra-xiphosurans differences recovered from the previous analyses (Figure 7). At this systematic scale, the differences among xiphosurans appear more subtle (Figure 7). This conclusion is also supported by a higher estimated evolutionary rate for the marginal environment in the second analysis compared to the previous one (Table 4). However, it is important to note some limitations of this second approach as well. In the second analyses, Hoplitaspis hiawathai represents a significant portion of the total variability. This may not reflect the average status in Chasmataspidida. Another important factor limiting our study is a possible error in environmental assignment introduced by the ethological aspect of several of the examined taxa. It is well known that aquatic euchelicerates possessed gregarious behaviors, often associated with group molting events (Daley and Drage, 2016; Bicknell et al., 2019b; Lamsdell et al., 2019; Lustri et al., 2021) that took place in shallow waters. This may have affected our environmental classification for species such as Hoplitaspis hiawathai (Lamsdell et al., 2019), Eurypterus lacustris (Ruebenstahl et al., 2021) and Prolimulus woodwardi (Lustri et al., 2021) for which gregarious behaviors are reported and specimens were moults, meaning we cannot exclude this to be the case for other taxa involved in the study.

In the all-euchelicerates analysis, stronger support for the BMM over the BM is demonstrated, suggesting the presence of two different evolutionary rates for the two different environments. This signal was not recoverable when using the dataset that relies on only xiphosurans. The all-euchelicerates analysis reinforces the results of the xiphosurans-only analysis and gives a phylogenetic perspective to the evolutionary patterns of xiphosurans. Both analyses show an increase in morphological plasticity, independent of phylogeny but associated with the colonization of new environments. Furthermore, the second analysis shows that while these changes happened and are likely related to different environments within Xiphosura, they are much less pronounced than in the absence of outgroups. Body proportions in adulthood are generally stereotypical for any given species, and ontogenetic development is the process leading to their establishment. This implies that growth and form are related, but not by a simple relationship of cause and effect, because the starting point of body proportions at hatching/birth plays an important role too. Nevertheless, in light of the highest evolutionary rates of allometric growth recovered from our analyses associated with the freshwater environment, is still important to note extreme proportions of the prosomal shield even when they occur in species known from only one or few specimens where developmental data are lacking. Extremes in the proportions of the prosomal shield are often recovered in Mesozoic freshwater taxa, several of which have been excluded by our analyses owing to the lack of data about the development. This is the case for the radiation of Austrolimulidae. Austrolimulidae such as Austrolimulus fletcher (Riek, 1955) and Dubbolimulus peetae (Pickett, 1984) for example, shows an exploration of extreme prosomal shield proportions at least at a single point in their development (Bicknell et al., 2022). Other examples are present among the grade belinurines of Belinuridae. Belinurus bellulus (König, 1825), Parabelinurus lunatus (Lamsdell, 2020) and Macrobelinurus arcuatus (Lamsdell, 2020), to name a few are all freshwater species with a prosomal shield with a relative width greater than the length resulting in a crescentic moon shape of the carapace. Even if information about the evolution of development is not available for these species, their wide prosomal shield proportion provides support for the hypothesis of freshwater environments being positively correlated with higher evolutionary rates in xiphosurans. The exploration of different prosomal shield proportions took place during an anatomical radiation as the group invaded freshwater environments (Lamsdell, 2016; Lamsdell 2021a).

The evolutionary scenario for the development of xiphosurans, depicted by the analyses accounting only for the intra-xiphosuran variability, shows a similar pattern to a random walk scenario of evolution (Figure 1C) or even a trend (Figure 1B) towards isometric growth. This is especially true in the case of Limulidae (Figure 6). However, this is not the case for Belinuridae, which appear to have had explored a wide range of allometric patterns, neither is it the case for Paleolimulidae, which, even if at a lower degree, did experience different developmental patterns (Figure 6). On the other hand, in the case of Limulidae, a broader phylogenetic perspective finds a general accordance with stasis (Figure 7), showing gradual evolutionary change around a relatively static average morphological position through time (Figure 1A). A broader phylogenetic perspective also reduces the perceived variability in the families Belinuridae and Paleolimulidae. This second analysis provides the appropriate systematic level for the study of developmental evolution in xiphosurans and their related aquatic euchelicerates. In other words, changes in the allometric growth of xiphosuran species are associated with different environments, but they are relatively minor compared to the different allometric patterns found in their closest relatives (Figure 7).

Uniting knowledge of non-Xiphosura euchelicerates with knowledge of Xiphosura development has improved the understanding of the evolution of Xiphosura allometric growth patterns. While this study supports the idea that the colonization of new environments has led to increased evolutionary rates for allometric growth in xiphosurans, the wider phylogenetic framework of our analyses suggest that the entirety of those changes were still somewhat limited when compared to changes seen more broadly in euchelicerates, as it represents a small portion of the total variability observed for euchelicerates. It appears clear that the evolution of xiphosurans cannot be pigeonholed into simplistic terminology such as “living fossils”. Less impactful but more concrete definitions such as “gradual morphological evolution around a middle trait parameter” may better explain the observed pattern, at least regarding allometric growth. Furthermore, this research compares phylogeny-based evolutionary modelling without and with outgroups, emphasizing the importance of the latter to contextualize and to properly interpret the evolution for the target group.

The results show that the evolutionary rates of development of Xiphosura undergoes significant changes throughout the evolutionary history of the group, in concert with the adaptive radiation of the group as they exploit different environments through evolutionary time, and independently from their phylogenetic position. They also highlight the importance of considering outgroups when attributing evolutionary trends to a specific group. The magnitude of allometric growth among Xiphosurais was lower than in other euchelicerates with similar environmental affinities, which flattens what might otherwise appear as an explosion in diversity based solely on the observation of Xiphosura.

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

LLu: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Writing – original draft, Writing – review & editing. JA: Methodology, Supervision, Writing – review & editing. FS: Writing – review & editing. CH: Conceptualization, Data curation, Supervision, Writing – review & editing. LLa: Conceptualization, Methodology, Writing – review & editing. RG: Conceptualization, Methodology, Writing – review & editing. JH: Supervision, Writing – review & editing, Conceptualization, Data curation. AD: Funding acquisition, Resources, Supervision, Writing – review & editing, Conceptualization.