New perspectives on the evolutionary history of xiphosuran development through comparison with other fossil euchelicerates

the


Introduction
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łazėjowski, 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).
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 Graphical explanation of the evolutionary stasis often attributed to Xiphosurida.(A) Example of stationary morphological evolution around a middle trait parameter; (B) example of morphological evolution with a trait parameter diverging through time; (C) example of a "random walk".Edited from Tëmkin and Eldredge (2015).Lustri et al. 10.3389/fevo.2023.1270429Frontiers in Ecology and Evolution frontiersin.org 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.

Studied specimens
The specimens used in this study belong to eight different species of Xiphosura, including two extant and six extinct taxa,

Regressions of morphometric measurements during ontogeny
When the slope of the linear regression for the prosomal shield (carapace) width and length was not directly available in the literature, it was calculated.Prior to performing the linear regression, the natural log of all datasets was taken to reduce the skewness.Linear regression analyses were then conducted for prosomal shield lengths vs. prosomal shield widths for the following species: Eurypterus lacustris, Hoplitaspis hiawathai, Prolimulus woodwardi, Euproops sp., Paleolimulus kunguricus and Limulus polyphemus.The regression slopes for Paleolimulus signatus, Euproops danae, Mesolimulus walchi, and Tachypleus tridentatus were taken from Bicknell et al. (2022).Therefore, at least one representative of three of the four taxa of Xiphosura (Paleolimulidae, Bellinurina, and Limulidae) were considered, alongside two non-xiphosuran euchelicerates.Linear regressions of the measurements of Eurypterus lacustris, Hoplitaspis hiawathai, Prolimulus woodwardi, Euproops sp., Paleolimulus kunguricus and Limulus polyphemus were performed with the function "lm" in RStudio 2021.09.0 + 351 "Ghost Orchid".The plots of the linear regressions were made using the function "plot" in RStudio 2021.09.0 + 351 "Ghost Orchid" and subsequently edited with Adobe Illustrator (R script in Supplementary Datasheet S1).All the linear regression slopes are reported in Table 2.

Phylogenetic analyses
Bayesian phylogenetic analyses were performed using the matrix from Lamsdell (2020), with the addition of Hoplitaspis hiawathai and Prolimulus woodwardi.The character coding for Hoplitaspis hiawathai was based on Lamsdell et al. (2019), and character coding for Prolimulus woodwardi was based on Lustri et al. (2021).Eurypterus lacustris was not present in this matrix and has not been coded.Instead, Eurypterus tetragonophthalmus was used as a proxy representing the relative phylogenetic position of Eurypterus lacustris as the utilized matrix is expected to be coded identically for them both (in Figure 4, Eurypterus lacustris would have appeared as a sister species to Eurypterus tetragonophthalmus highlighted in red).The methods are the same as in the original work from Lamsdell (2020), using MrBayes ver.3.2.7a(Huelsenbeck and Ronquist, 2001).The final data matrix includes 162 taxa and 259 discrete characters.The analyses consisted of four independent runs of 10,000,000 generations and four chains each, under the maximum likelihood model with gamma-distributed rate variation among sites (Mkv + G:) (Lewis, 2001).Characters were unordered and given equal weighting (Congreve and Lamsdell, 2016).Trees were sampled every 100 generations.The resulting trees per run is 1,000,000 and the first 25,000 sampled trees of each run were discarded as burn-in.Extended majority rule tree obtained was used for the subsequent analyses (Figure 4).The matrix used for the phylogenetic analyses and the mrBayes code are available in Supplementary Datasheet S2.

Ancestral state reconstruction and the estimation of evolutionary rates
The phylogenetic tree obtained with the Bayesian phylogenetic analyses was pruned with the "ape" (Paradis and Schliep, 2019) function "drop.tip" in RStudio 2021.09.0 + 351 "Ghost Orchid".Two different pruned trees were obtained.The first tree retained the branch length and node positions for xiphosurans species with known growth-pattern data (prosomal shield length and width ratio along the growth), resulting in a tree with 8 tips and 7 internal nodes.The second tree retained the branch length and the node positions for all euchelicerate species with known growthpattern data (prosomal shield length and width ratio along the growth), resulting in a tree of 10 tips and 9 internal nodes.The branch lengths of the trees are based on morphological character distance.Using the packages "mvMORPH" (Clavel et al., 2015), "ape" Extended majority rule tree of the Bayesian analysis performed on the matrix modified from Lamsdell (2020).(A) Root section of the tree; (B) Xiphosura section of the tree; (C) Eurypterida and Chasmataspidida section of the tree.(D) Arachnida section of the tree.The species used in the analyses of this paper are highlighted in red.The tree is based on a matrix composed of 162 taxa and 259 characters.(Paradis and Schliep, 2019) and "phytools" (Revell, 2012), ancestral state analyses have been performed on both trees, incorporating the discrete environmental data of Lamsdell (2021), here divided between marginal (not fully marine) and marine settings following Bicknell et al. (2022).The ML function was used for a maximum likelihood estimation of the ancestral state under Brownian motion models.Subsequently, two evolutionary models were tested on the ancestral states recovered from both pruned phylogenetic trees: a Multivariate Brownian motion process (BM), and a Multi-rate Brownian motion process (BMM).The BM is a model in which a single path of evolution is simulated under Brownian motion processes while the BMM is a model in which multiple paths of evolution are simulated under Brownian motion processes.The models were used to test the null hypothesis where an absence of correlation between environmental (in our case two variables marginal and marine) and development would result in the BM model (allowing only one evolutionary path) outperforming the BMM model (allowing two different evolutionary paths).The evolutionary rates of marginal and marine species were calculated under both models using the "mvBM" command to investigate the possible correlation between the environment and the evolution of ontogenetic characters over time.
The ontogenetic variation investigated is the change in shape of the prosomal shield (carapace) during ontogeny.The fit of these different evolutionary models was assessed by calculating the Akaike weight with the command "aicw".All the aforementioned analyses were made using RStudio 2021.09.0 + 351 "Ghost Orchid" (R script in Supplementary Datasheet S3).

Results
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 alleuchelicerates 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.

Discussion
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 Linear regressions of morphometric measurements of the prosomal shield length and width of 6 different euchelicerate species.Slope values resulting from the regressions are available in Table 2 alongside slopes taken from the literature.
Ancestral state reconstruction of the prosomal shield allometric growth and environments of Xiphosura.Numbers at the nodes are the reconstructed ancestral state for the prosomal shield slope.Pie charts express the probability of a marginal or marine environment at the node.0.888 is the lowest slope value and 1.02 is the highest slope value.Results of the evolutionary models (BM and BMM) tested on the tree are available in Table 3.

FIGURE 7
Ancestral state reconstruction of the prosomal shield allometric growth and environments of Xiphosura, Eurypterida and Chasmataspidida.Numbers at the nodes are the reconstructed ancestral state for the prosomal shield slope.Pie charts express the probability of a marginal or marine environment at the node.0.478 is the lowest slope value and 1.02 is the highest slope value.Results of the evolutionary models (BM and BMM) tested on the tree available are in Table 4.
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 intraxiphosurans 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.

Conclusion
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.

TABLE 2
Linear regression slopes of all the taxa examined in the study.

TABLE 3
AIC supports and estimated rate of evolution for the BM and BMM models that include only the xiphosurans.

TABLE 4
AIC supports and estimated rate of evolution for the BM and BMM models that include Xiphosurida, Eurypterida and Chasmataspidida.