Ecdysozoans (Phyla Arthropoda, Kinorhyncha, Loricifera, Nematoda, Nematomorpha, Onychophora, Priapulida, Tardigrada) are invertebrates bearing a tough, periodically moulted cuticle that predisposes them to exceptional preservation. Ecdysozoans dominate the oldest exceptionally preserved bilaterian animal biotas in the early to mid-Cambrian (c. 520–508 Ma), with possible trace fossils in the latest Ediacaran (<556 Ma). The fossil record of Ecdysozoa is among the best understood of major animal clades and is believed to document their origins and evolutionary history well. Strikingly, however, molecular clock analyses have implied a considerably deeper Precambrian origin for Ecdysozoa, much older than their earliest fossils. Here, using an improved set of fossil calibrations, we performed Bayesian analyses to estimate an evolutionary time-tree for Ecdysozoa, sampling all eight phyla for the first time. Our results recover Scalidophora as the sister group to Nematoida + Panarthropoda (= Cryptovermes nov.) and suggest that the Ediacaran divergence of Ecdysozoa occurred at least 23 myr before the first potential ecdysozoan trace fossils. This finding is impervious to the use of all plausible phylogenies, fossil prior distributions, evolutionary rate models and matrix partitioning strategies. Arthropods exhibit more precision and less incongruence between fossil- and clock-based estimates of clade ages than other ecdysozoan phyla.

Supplementary material: Full methodologies used and an appendix of fossil calibration points are available at https://doi.org/10.6084/m9.figshare.c.5811381

Thematic collection: This article is part of the Advances in the Cambrian Explosion collection available at: https://www.lyellcollection.org/cc/advances-cambrian-explosion

The Cambrian ‘Explosion’ of animal (and particularly bilaterian) diversity is among the most remarkable evolutionary phenomena in the geological record. Numerous animal phyla, bodyplans and ecological guilds first appear as fossils in the early to mid-Cambrian, with many of those lineages persisting and dominating through the Phanerozoic (Erwin and Valentine 2013). The Cambrian Explosion is best captured by early to mid-Cambrian Konservat-Lagerstätten; fossil biotas with soft-tissue preservation that are often dominated by ecdysozoan worms and arthropods (Caron and Jackson 2008; Zhao et al. 2009; Vannier and Martin 2017). Ecdysozoans are bilaterians with a tough, periodically moulted cuticle that, given their inclusion of insects and nematodes, make up the vast majority of extant metazoan biodiversity. Ecdysozoa is composed of the phyla Arthropoda, Onychophora and Tardigrada (often grouped in Panarthropoda), Nematoda and Nematomorpha (often grouped in Nematoida), plus Kinorhyncha, Loricifera and Priapulida (often grouped in Scalidophora). Experimental taphonomy shows that ecdysozoan cuticles are highly recalcitrant (Murdock et al. 2014; Butler et al. 2015; Sansom 2016), and their periodic moulting means that there are many more cuticles than individuals, each with the potential for fossilization. Because of their vermiform bodyplans and/or paired appendages, ecdysozoans have also left an extensive record of bioturbation as early conspirators in the Ediacaran–Cambrian substrate revolution (Bottjer et al. 2000; Vannier et al. 2010). Thus, the fossil record of ecdysozoans is widely considered to offer a faithful record of their evolutionary history (Daley et al. 2018; Budd and Mann 2020a, b). The ecdysozoan fossil record has therefore been used as a model for assessing the veracity of the Precambrian and Cambrian record more generally, and the ability of that record to document the origins of modern bilaterian groups (Daley et al. 2018; Budd and Mann 2020a, b).

The earliest unequivocal ecdysozoan body fossils are of Fortunian (Terreneuvian) age (c. 535 Ma). These occur as worm-like phosphatic microfossils with radial proboscis spines interpreted as total-group Scalidophora (Liu et al. 2014, 2018, 2019; Zhang et al. 2015, 2018; Shao et al. 2016, 2020a, 2020b; Wang et al. 2019, 2020). However, trace fossil evidence of ecdysozoan affinity extends to the latest Ediacaran (Buatois 2018), with a maximum age of more equivocal examples ranging from 551 to 555 Ma (Parry et al. 2017; Chen et al. 2018, 2019). Indeed, the trace fossil Treptichnus pedum, long attributed to a priapulan trace-maker (Vannier et al. 2010; Kesidis et al. 2019), defines the Ediacaran–Cambrian boundary. Some molecular clock estimates for the emergence of crown-group Arthropoda are in close accord with the fossil record (dos Reis et al. 2015; Cunningham et al. 2017; Daley et al. 2018), but there is a much greater discrepancy between molecular estimates and the oldest fossil evidence of Ecdysozoa as a whole (Erwin et al. 2011; Rehm et al. 2011; Rota-Stabelli et al. 2013; dos Reis et al. 2015). This discrepancy may simply reflect incompleteness of the early ecdysozoan fossil record (see Cunningham et al. 2017) or result from the systematic overestimation of clade ages by molecular clock methods (Budd and Mann 2020a, b). There are good reasons, however, to expect a discrepancy between first fossil occurrences and clock estimates of divergence. This is not only because of the usual concerns about the completeness of the rock and fossil records (Holland 2016) but also because lineage divergence events are genomic phenomena that are invisible to the fossil record until descendent lineages have acquired distinguishable and fossilizable characteristics (Donoghue and Yang 2016). Even the largest and most species-rich animal clades, of which Ecdysozoa is the ultimate example, originate with a single lineage at a centre of endemism, with no guarantee that this particular environment will be preserved and yield fossils. Low species richness, potentially low population sizes implicit in some models of rapid speciation (Mayr 1963; Eldredge and Gould 1972) and restricted geographical distributions all militate against the preservation of early clade congeners. The arthropod fossil record illustrates this point, with some studies indicating poor congruence between phylogenetic branching order and first fossil occurrence dates (Wills 2001; O'Connor and Wills 2016). Crucially, however, to produce realistic evolutionary timescales, clock-based methods rely on accurate phylogenetic hypotheses and integration of fossil evidence within accurate stratigraphic, geochronological and phylogenetic contexts.

To explore the Neoproterozoic extent of ecdysozoan evolutionary history, we assembled a phylogenomic super-matrix containing new genomic and transcriptomic data, sampling all eight ecdysozoan phyla for the first time in a dating study. Our Bayesian phylogenetic analyses support a monophyletic Scalidophora, and a sister group relationship between Nematoida and Panarthropoda. A grouping of nematoids and panarthropods has regularly emerged in previous phylogenomic studies (Hejnol et al. 2009; Campbell et al. 2011; Borner et al. 2014; Laumer et al. 2015, 2019; Yoshida et al. 2017), which is, in turn, sister group to Scalidophora here. Our Bayesian analysis was dated using 55 fossil calibrations, which are revised and improved significantly from previous studies in three ways. First, we include a number of new fully justified clade-age calibrations following the best practice principles of Parham et al. (2012). Second, our approach eschews fossil calibrations where crown-group lineage assignment is unclear (e.g. Priapulida, Nematoda and Tardigrada). Third, we substantially improve coverage of non-arthropod lineages, including new data for Nematomorpha and Tardigrada. The robustness of our inferences of ecdysozoan clade ages was tested by varying (1) evolutionary rate heterogeneity and heritability, (2) calibration strategy (i.e. alternating the prior distributions to reflect different interpretations of the fossil record) and (3) the hypothesis of ecdysozoan relationships. Overall, whereas we found that each of these experiments influenced the posterior distributions (hereafter, ‘posteriors’), the estimates of the age of crown-group Ecdysozoa were always Ediacaran, in many cases post-dating the iconic Ediacaran non-biomineralized macrobiotas. Our study indicates that a palaeontologically cryptic origin of bilaterian animal evolution before a geologically conspicuous diversification in the terminal Ediacaran to Cambrian is plausible.

We performed Bayesian evolutionary analyses to estimate the timescale of ecdysozoan evolution from molecular sequences and the fossil record, attempting to constrain sources of uncertainty surrounding clade age estimates. First, we assembled a super-matrix containing 228 protein-coding genes (with 43 852 amino acid positions) for 65 ecdysozoan taxa from all eight phyla (Supplementary material Methods 1 and 2; including three newly sequenced ecdysozoan genomes or transcriptomes, those of the nematomorph Nectonema munidae, and the heterotardigrades Actinarctus doryphorus and Echiniscus bisetosus), as well as spiralian, deuterostome and non-bilaterian metazoan outgroups. Second, phylogenies were inferred using PhyloBayes MPI v.15.a (Lartillot et al. 2013) with site-heterogeneous CAT models (Lartillot and Philippe 2004) (Supplementary material Methods 3). We employed or devised 55 fossil calibration points compatible with this phylogenetic hypothesis (Supplementary material Appendix), with a 2.5% soft tail on upper bounds and a hard minimum on lower bounds. This follows best practice in that the calibration minima are established based on such a conservative interpretation of the fossil record, eschewing older–riskier records in favour of younger–safer clade records, that we can be certain that the minimum bounds underestimate the true clade age, all other factors considered (e.g. phylogeny, stratigraphy, geochronology; Parham et al. 2012). No such assurances can be made of the soft maximum, which, as its name suggests, is uncertain; we implement an arbitrary but objective 2.5% (of the minimum–maximum span) tail distribution to the maximum bound to allow for the non-zero probability that the true clade age could lie beyond the maximum (Parham et al. 2012). These bounds augmented existing calibrations (Benton et al. 2015; Wolfe et al. 2016) formulated following best practice principles (Parham et al. 2012) (see Supplementary material Appendix). We estimated ecdysozoan clade divergence times using relaxed clock methods in MCMCTree, part of the PAML4.9 package (Yang 2007), with substitutions modelled using the LG + G model (Le and Gascuel 2008) (Supplementary material Methods 4).

To explore the uncertainty associated with our clade age estimates, we considered four key variables when generating time-trees: (1) two alternative relaxed clock models, to compare the impact of model assumptions of rate heritability (dos Reis et al. 2015); (2) four alternative matrix partitioning schemes, to compare the effect of alternative substitution modelling (Telford et al. 2014; Egger et al. 2015); (3) four alternative calibration densities, modelling interpretations of the fossil record as a good or poor approximation of the history of Ecdysozoa; (4) five alternative phylogenetic hypotheses, to explore the impact of phylogenetic uncertainty around ecdysozoan relationships (see Giribet and Edgecombe 2017). Although we were able to reject some of these variables, this was not possible for all and so in an attempt to constrain methodological uncertainty we derived a summary integrative timescale that combines the posterior time estimates from across all conditions. Full details of our methods are presented in detail in Supplementary material Methods 1–4.

Phylogenetic analyses

We undertook a sequential molecular clock analysis and, therefore, it was necessary for us first to estimate ecdysozoan phylogeny. We recovered monophyly of Ecdysozoa with full support: posterior probability (PP) = 1. Within Ecdysozoa (see Fig. 1; Supplementary material Figs S4 and S5), we resolved three major subclades that are supported on morphological grounds: (1) Scalidophora Lemburg (1995), worms with a retractable introvert bearing radially arranged spines known as scalids (Loricifera + (Priapulida + Kinorhyncha)); (2) Nematoida Schmidt-Rhaesa (1996), elongate worms with no circular body wall musculature but with layers of crossed collagenous fibrils in the cuticle (Nematoda + Nematomorpha); (3) Panarthropoda Nielsen (1995), taxa with paired and segmental ventro-lateral appendages and a supraesophageal ganglionar brain (Tardigrada + (Onychophora + Arthropoda)). Nematoida and Panarthropoda were recovered with full support, whereas Scalidophora was recovered with a PP of 0.89. Nematoida was recovered as the sister to Panarthropoda (PP = 1), as in another study (Campbell et al. 2011). The putative clade Cycloneuralia Ahlrichs (1995) (Scalidophora + Nematoida) is therefore rendered paraphyletic with respect to Panarthropoda. Although their internal relationships vary, nematoids plus panarthropods are commonly recovered as a clade in ecdysozoan phylogenomic analyses (Hejnol et al. 2009; Campbell et al. 2011; Borner et al. 2014; Laumer et al. 2015, 2019; Yoshida et al. 2017). As such, we recognize the grouping formally, proposing the name Cryptovermes taxon nov., meaning ‘hidden/secret worms’ reflecting the unknown characteristics of the (presumably) worm-like ancestor of this clade. Nematoda and Nematomorpha were recovered as sister taxa as in a number of previous studies (Schmidt-Rhaesa 1998; Bleidorn et al. 2002; Campbell et al. 2011). Among Panarthropoda, Arthropoda and Onychophora were found to be sister taxa, with Tardigrada as sister to that clade. We found no support for a putative tardigrade–nematode/nematoid relationship, which has elsewhere been attributed to long-branch attraction (Campbell et al. 2011).

Multiple phylogenomic analyses have now recovered a clade comprising Nematoda, Tardigrada, Onychophora and Arthropoda (plus Nematomorpha where these are sampled; see citations above). As such, Cryptovermes is defined using phylogenetic nomenclatural principles as the clade composed of Arthropoda, Onychophora, Tardigrada, Nematoda and Nematomorpha, their common ancestor and all of its living and extinct descendants; members of the total-group Cryptovermes would also include all extinct taxa that are more closely related to the Cryptovermes crown ancestor than to Scalidophora (or other Ecdysozoa if Scalidophora is paraphyletic, as in Laumer et al. 2019). We retrieved Cryptovermes as the sister-group to Scalidophora, but enhanced taxon-sampling (particularly of loriciferans) would be required to reject the possibility that Scalidophora is paraphyletic. Regardless, Ecdysozoa can be defined as the crown-clade composed of Cryptovermes, Priapulida, Kinorhyncha and Loricifera, their last common ancestor, and all its living and extinct descendants. Extinct taxa more closely related to crown-Ecdysozoa than to its sister-clade, Spiralia (= Lophotrochozoa), belong to stem-group Ecdysozoa; stem- and crown-group Ecdysozoa comprise total-group Ecdysozoa.

Molecular clock analyses; accounting for rate heritability

Previous studies have shown that the different assumptions of rate heritability in alternative relaxed clock models can significantly influence posterior age estimates (dos Reis et al. 2015). As such, all dating analyses were conducted in two iterations, using both autocorrelated (AR) and independent rates (IR) models in MCMCTree (Yang 2007). In the AR model, rates are autocorrelated between ancestor and descendant lineages, whereas in the IR model independent branch-specific rates are drawn from a common distribution (Rannala and Yang 2007). We therefore present the results from both models throughout (see Table 1). Across both clock models and all other parameters tested (see following sections), crown-group Ecdysozoa diverged between 636 and 578 Ma; that is, between the Cryogenian–Ediacaran boundary and mid- to late Ediacaran (see Table 1). Estimates for the divergence of crown-groups Scalidophora, Cryptovermes, Panarthropoda and Arthropoda were largely confined to the Ediacaran. Crown-group Scalidophora ranges from early Ediacaran to early Cambrian (617–534 Ma). Crown-group Arthropoda ranges from mid-Ediacaran to just prior to the Ediacaran–Cambrian boundary (589–540 Ma), the Ediacaran–Cambrian boundary being constrained to 538.8–538.6 Ma by recent radiometric dating (Linneman et al. 2019). Crown-group Nematoida yielded the least precise range of estimates for a supra-phylum level clade, ranging from the mid- to late Ediacaran to the Silurian, and similar imprecise estimates beginning in the late Ediacaran and ranging well into the Phanerozoic were recovered for crown-groups Tardigrada (up to late Devonian), Priapulida (up to late Carboniferous) and Kinorhyncha (up to early Triassic). Wholly Phanerozoic distributions were recovered with poor precision for crown-groups Nematoda (early Cambrian to early Permian), Nematomorpha (mid-Cambrian to mid-Cretaceous) and Onychophora (mid-Cambrian to early Cretaceous).

Molecular clock analyses; exploring rate heterogeneity through partition strategies

To determine the influence of evolutionary rate heterogeneity among lineages on clade age estimates, we partitioned our dataset according to the gene-wise evolutionary rate (see Figs 2 and 3). We used the approach of previous researchers (Telford et al. 2014; Egger et al. 2015) to rank genes from slowest to fastest evolving (see Supplementary material for details). We then divided the dataset into two, four and five subgroups (with each of the partitions having an equal number of genes) and analysed the sequence alignment in these alternative partitions (see methods in Supplementary material). Increasing the number of partitions led to narrower posterior confidence intervals (i.e. more precise results). Figure 2 shows the width of the 95% posterior confidence intervals plotted against the posterior mean ages for the corresponding nodes (the ‘infinite sites’ plot). The decreasing regression coefficient (slope of the regression line; see Fig. 2) in the infinite sites plots is a measure of this reduction in the posterior interval width, indicating the uncertainty resulting from the fossil calibrations (Inoue et al. 2010). On average, the amount of uncertainty added per 100 myr is relatively small (c. 8 or 9 myr), but this is not evenly distributed, with posterior interval widths remaining shorter in highly calibrated time periods; that is, those intervals occurring within or close to the Cambrian, where many node calibrations occur. We found that the mean age for most nodes became older as the number of partitions increased. This was consistent for all divergences under the IR model, whereas the AR model was variable, with some decreases as well as increases in divergence time estimates as the number of partitions increased. In addition, using the binned dataset as a reference, we performed four further divergence time estimation experiments: (1) including only the slowest evolving genes (11 883 sites); (2) including only the fastest evolving genes (7442 sites); (3) with both the slowest and the fastest genes excluded (12 784 sites); (4) including only those genes for which at least 70% of species provided data (29 796 sites). We plotted the results against our time estimates from the complete dataset (Supplementary material Fig. S7). The strong (R2 ≈ 1) correlation between these analyses indicates that our posterior estimates were not greatly influenced by these different approaches.

Molecular clock analyses; exploring the impact of competing calibration strategies

Time priors in Bayesian molecular dating analyses are usually based on minimum age constraints (often the age of the oldest fossil representative of the clade) combined with a probability distribution that reflects a prior view of the relationship between the minimum age constraint and clade age (Yang and Rannala 2006). It is well established that these choices significantly affect posterior age estimates (Inoue et al. 2010; Warnock et al. 2012). As such, it is possible to model the uncertainty introduced by the incomplete nature of the fossil record by selecting alternative probability distributions that reflect whether the stratigraphic age of the calibrating fossil is likely to be a close or a poor approximation of the true clade age (dos Reis et al. 2015; Betts et al. 2018). For example, scalidophorans and panarthropods are abundant in Cambrian Konservat-Lagerstätten, but total-group nematoids are conspicuously absent, even though the phylogeny implies that they must also have diverged by this point as their sibling lineages, like Panarthropoda, are already manifest in the fossil record (see Fig. 1). This stratigraphic gap between clade divergence and fossil first occurrence is known as a ghost lineage (Norell and Novacek 1992; Novacek 1992). Nematoida therefore has a cladistically implied ghost lineage of over 100 myr, from the oldest confirmed nematoid (a probable total-group nematode from the Lower Devonian Rhynie Chert; Poinar et al. 2008) to biomineralized and soft-bodied early Cambrian faunas containing its probable sister-clade, Panarthropoda. Likewise, crown-group representatives of Nematomorpha are first recorded from the Cretaceous (Voigt 1938; Poinar 1999; Poinar and Buckley 2006) and crown-group Nematoda from the Jurassic (Poinar et al. 1994). Similarly, the earliest crown-group tardigrade is not known until the Cretaceous (Bertolani and Grimaldi 2000), but fossil evidence of at least total-group Tardigrada is recognized in the Cambrian (Maas and Waloszek 2001; Maas et al. 2007). For these reasons, there is a prior expectation that the fossil calibration minima for Nematoida, Nematoda, Nematomorpha and Tardigrada are a poor approximation of their respective crown-clade ages. We therefore estimated divergence times across the tree with four alternative calibration strategies for 15 key nodes to represent different palaeobiological interpretations of fossil minima (see Fig. 4). All calibrations employed a hard minimum age constraint and a soft maximum age constraint (Donoghue and Benton 2007), calculated using MCMCTreeR (Puttick 2019). The four alternatives were: (1) uniform probability distribution between the minimum and maximum age, with a 2.5% tail on the maximum, to represent agnosticism concerning the true time of divergence between fossil calibration minima and maxima; (2) truncated Cauchy distribution, where the mode of the distribution was shifted towards the minimum age (i.e. 5% of the uncertainty older than the minimum age), and with a long tail extending into the past, to represent an ‘optimistic’ interpretation of the fossil record, implying that the minimum constraint is a close approximation of the true age of divergence; (3) truncated Cauchy distribution where the mode of distribution was at the midpoint (i.e. 50%) between minimum and maximum age constraints, to represent an intermediate view of the efficacy of the fossil minima; (4) truncated Cauchy distribution where the mode of the distribution was shifted towards the maximum age (i.e. 95%), with a long tail from the minimum, to represent a ‘pessimistic’ view in which the calibration minima are considered a poor approximation of the true age of divergence. Figure 4 and Table 1 show how these alternative calibration densities affected clade age estimates.

Using the uniform (agnostic) calibration strategy, both the autocorrelated (AR) and independent (IR) models estimated crown-group Ecdysozoa originating in a c. 30 myr interval of the Ediacaran (AR: 609–583 Ma; IR: 621–593 Ma). For the AR model, the ‘intermediate’ and ‘pessimistic’ calibration strategies yielded similar posteriors for the age of Ecdysozoa, but the ‘optimistic’ strategy yielded much younger posteriors (599–578 Ma). For the IR model, the ‘optimistic’, ‘intermediate’ and ‘pessimistic’ calibrations all yielded similar posteriors for the age of crown-group Ecdysozoa that were much younger than the uniform strategy (612–583 Ma).

Molecular clock analyses; exploring the impact of phylogenetic uncertainty

In addition to the uncertainty introduced by rate heterogeneity and calibration strategy, different phylogenetic hypotheses change the distribution of this uncertainty by implying different branch lengths. The internal relationships of Ecdysozoa are still debated (Giribet and Edgecombe 2017). We therefore estimated divergence times according to additional alternative sister-group relationships to the relationships we inferred in this study. Altogether, the topologies tested include the following: (1) the results of our CAT–Poisson + G analysis (see Fig. 1 and Supplementary material Fig. S4); (2) the results of our CAT–GTR + G analysis (see Fig. S5), which were almost identical to those of CAT–Poisson + G but yielded a different sister-group relationship for Hexapoda; under CAT–GTR + G, Hexapoda is sister-group to Branchiopoda, whereas under CAT–Poisson + G it is sister-group to Remipedia (Labiocarida), which is supported by several phylogenomic studies (von Reumont et al. 2011; Oakley et al. 2012; Misof et al. 2014; Schwentner et al. 2017; Lozano-Fernandez et al. 2019); (3) Nematoda + Tardigrada (see Fig. S10), which is supported by some phylogenomic studies (Hejnol et al. 2009; Borner et al. 2014; Laumer et al. 2015; Yoshida et al. 2017), although expressly refuted by others (Rota-Stabelli et al. 2010; Campbell et al. 2011); (4) Arthropoda + Tardigrada (Tactopoda; see Fig. S10), which is endorsed by some studies of fossil lobopodians that include characters of the ventral nerve cord supporting this group (e.g. Mayer et al. 2013; Smith and Ortega-Hernández 2014); (5) Priapulida + Loricifera (Vinctiplicata; see Fig. S10), which is endorsed by some morphological studies (Kristensen 1983, 2017; Wills 1998; Lemburg 1999; Ax 2003; Maas et al. 2007; Peel et al. 2013; Harvey and Butterfield 2017). All alternative topologies yielded a posterior estimate for the age of crown-group Ecdysozoa within the Ediacaran (Figs 5 and 6). For both models (autocorrelated rates and independent rates), the Labiocarida and Vinctiplicata topologies yielded posterior age estimates for Ecdysozoa somewhat older than the Hexapoda + Branchiopoda, Tactopoda and Nematoda + Tardigrada topologies (Fig. 5). Crown-group Tardigrada was estimated, with more precision, to have diverged in the Ediacaran under the Labiocarida and Vinctiplicata topologies, whereas other conditions allowed Tardigrada to be considerably younger and with more widely distributed 95% higher posterior density (HPD) (see Fig. 6). The crown-clades of Priapulida and Nematoda were not influenced as greatly by alternative topologies.

Establishing an integrated timescale for ecdysozoan evolution

Clade age estimates derived from clock methods are vanishingly unlikely to be the same as the first appearance datum of those clades in the fossil record, as clocks and rocks record different events. An offset between the two should be expected (Donoghue and Yang 2016), with molecular clock analyses yielding older clade age estimates. The first appearance of a clade in the fossil record reflects the accumulation of recognizable, phylogenetically constrained anatomical characteristics in descendent lineages, coupled with favourable conditions for preservation, discovery and interpretation. The first appearance of a clade in the fossil record cannot occur before the point of genetic divergence; the degree to which the oldest fossil approximates the age of its clade is therefore contingent on the timing of the origin of fossilizable apomorphies, preservation and sampling. A number of the fossil-based node calibrations for Ecdysozoa are not very informative, post-dating clade ages by hundreds of million years, as can be demonstrated by ghost lineage analysis of fossil evidence alone (Fortey et al. 1997). This is clearly the case for Cretaceous amber-derived calibrations of clades whose successive sister-lineages have Cambrian fossil records (e.g. the crown-clades of Tardigrada, Nematoda and Nematomorpha). As such, although these calibrations preclude very young estimates of clade ages, they introduce much uncertainty to our overall estimates. Furthermore, the methodological uncertainty, such as the choice of clock model, partition strategy and calibration density, can lead to different clade age estimates. Where it was not possible to reject competing methodological choices, we compiled the age estimates (always expressed as the 95% HPD in millions of years before present) for all clades, deriving a single integrated time-tree in which the uncertainty associated with each node reflects a compounding of the minimum and maximum bounds on the 95% HPD from all such experiments (Figs 6 and 7). As such, the clade age estimates in our integrated time-tree must be read in terms of the span of the uncertainty associated with the nodes. Thus, the position of the nodes (i.e. the mean) within our integrated evolutionary time-tree is arbitrary and uninformative; it is not a statistical mean. This follows recommendations from simulation-based analyses of the performance of molecular clocks, which showed that whereas the mean of the 95% HPD does not recover the true clade age, this is normally encompassed by the span of the 95% HPD (Warnock et al. 2017).

Our integrated time-tree estimates that crown-group Ecdysozoa diverged in the interval 636–578 Ma. This is slightly older than some previous studies that integrated sources of uncertainty in divergence time estimation (619–546 Ma, Rota-Stabelli et al. 2013; 629–567 Ma, dos Reis et al. 2015), which probably results from our denser sampling of ecdysozoan phyla. The estimated age for crown-group Cryptovermes also spans a similar early interval of 620–571 Ma, whereas crown-group Scalidophora is estimated to have diverged within 617–534 Ma (early Ediacaran–early Cambrian). Crown-group Panarthropoda exhibits a joint posterior estimate of 616–562 Ma (early Ediacaran–late Ediacaran), and crown-Nematoida ranges from 595 to 419 Ma (mid-Ediacaran–Early Devonian). The joint posterior estimates for the crown-clades of Arthropoda (589–540 Ma), Chelicerata (565–523 Ma) and Mandibulata (571–529 Ma) are notable in that they are shorter ranging than other phylum or subphylum level clades (i.e. the results are more precise). This shows that the estimates for crown-groups Arthropoda, Chelicerata and Mandibulata were more resilient to competing parameter choices than comparable clades, which generally presented much larger joint posterior intervals with greater variation in estimates from particular parameter combinations (e.g. Priapulida, Kinorhyncha, Nematoda, Nematomorpha, Tardigrada, Onychophora; see Fig. 6). This reflects the relative informativeness of the calibrations and, ultimately, the fidelity of the fossil record, which is consistent with prior studies on the effect of both fossil and tip taxon sampling on molecular clock analyses (Duchêne et al. 2015). We would anticipate that increased sampling in either dimension would increase the precision of these intervals.

A factor contributing to some early divergence dates in our study is the use of conservative soft maxima. Deep nodes in Panarthropoda and within the Arthropoda followed Warnock et al. (2012) and Wolfe et al. (2016) in using the Lantian biota (636.1 Ma) as a soft maximum, whereas some other recent studies have used younger maxima. For example, Lozano-Fernandez et al. (2020) and Howard et al. (2020a) considered the trace and body fossil records to constrain the divergence between Onychophora and Arthropoda to a soft maximum of 559 Ma (White Sea Ediacaran age) and the crown-groups of Chelicerata and Mandibulata to 538.8 Ma (based on a radiometric date for the base of the Cambrian; Linneman et al. 2019). These younger maxima inevitably led to shallower divergence times estimates for some arthropod groups such as Chelicerata and Arachnida (although mostly overlapping with estimates here) but did not alter the broad picture of an Ediacaran origin of ecdysozoans, yielding similar late to terminal Ediacaran estimates for crown-group Arthropoda. However, it remains unclear which of these two strategies reflects the best interpretation of the fossil record and, indeed, some researchers have argued that the older, Lantian-based soft maximum constraint is itself too young (Battistuzzi et al. 2015, 2018). Hence, we used the Lantian-based soft maximum constraint on deep nodes within the phylogeny, in combination with conservative hard minimum constraints, to encompass the range of interpretations of the degree to which the fossil record approximates the true timescale of animal evolutionary history.

The rock–clock mismatch in Ecdysozoa

As anticipated, the composite interval shown in Figure 6 for the divergence of ecdysozoans precedes the first appearance of ecdysozoan fossils. The oldest possible Ecdysozoan fossils are traces from the terminal portion of the Ediacaran (younger than 556 Ma) equivocally attributed to nematoids (Parry et al. 2017) and panarthropods (Chen et al. 2018, 2019), and small carbonaceous cuticular fragments compared with scalids also date from a similar time (Moczydłowska et al. 2015). The possible nematoid traces are meiofaunal worm burrows showing undulating nematoid-like movement from the Corumbá Group of Brazil, constrained to a maximum age of c. 555 Ma (Parry et al. 2017). The possible panarthropod traces are from the Shibantan Lagerstätte of South China (Xiao et al. 2021) showing locomotion perhaps associated paired appendages (Chen et al. 2018) and body segmentation (Chen et al. 2019), constrained to a maximum age of c. 551 Ma. The scalid-like cuticular fragments are from the Eastern European Platform and date to the late Ediacaran, younger than c. 552 Ma (Moczydlowska et al. 2015). The Ediacaran nematoid-like traces of the Corumbá Group are within our joint posterior interval for crown-group Nematoida (595–419 Ma), and the oldest certain nematoid (Palaeonema phyticum Poinar et al. 2008) is Early Devonian in age (c. 405 Ma). As such, our joint posterior for crown-group Nematoida is reasonably concordant with the fossil record. The joint posterior interval for crown-group Panarthropoda is at least 12 myr older (minimum estimate 562 Ma) than the late Ediacaran Shibantan panarthropod-like traces. The oldest certain panarthropod fossil taxon is the trace Rusophycus, which is geochronologically constrained to a minimum of 528.82 Ma (Wolfe et al. 2016) and is estimated to occur up to a maximum of c. 537 Ma, in the Fortunian Stage of the Cambrian (Daley et al. 2018). As such, our joint posterior estimate for crown-group Panarthropoda suggests that they diverged well before the fossil record indicates, as for the crown-clades of Ecdysozoa and Cryptovermes.

The direct fossil evidence for Ediacaran ecdysozoan taxa discussed above is equivocal. Although clocks are expected to exceed fossil first occurrences in age, the absence of Ediacaran ecdysozoans may be construed against the geobiological palatability of the offset between clock-based estimates obtained here for the age of crown-Ecdysozoa (636–578 Ma) and the beginning of the ecdysozoan fossil record around the Ediacaran–Cambrian boundary (538.8–538.6 Ma, Linneman et al. 2019). However, this mismatch is reconciled by a phylogenetic perspective on the trace fossil record, which suggests that the body and trace fossil records underestimate ecdysozoan clade ages. Specifically, the cosmopolitan ichnospecies Treptichnus pedum (= Phycodes pedum, = Trichophycus pedum), the first appearance of which defines the Neoproterozoic–Phanerozoic boundary (Narbonne et al. 1987; Buatois et al. 2013; Buatois 2018), is widely interpreted as the burrows of macrofaunal priapulan-grade organisms, conservatively indicating the presence of at least total-group Scalidophora (Vannier et al. 2010; Kesidis et al. 2019). Thus, if at least total-group Scalidophora or even crown-group Priapulida had originated before the Cambrian, our phylogeny for Ecdysozoa (see Figs 1 and 6; Supplementary material Figs S4 and S5) requires that all successive sister-lineages of Priapulida and/or Scalidophora (namely, total-groups Kinorhyncha, Loricifera and Cryptovermes, and therefore Ecdysozoa) must have diverged even earlier, by definition. Therefore, whether directly through up to >556 Ma trace fossils (Parry et al. 2017; Chen et al. 2018, 2019), or indirectly through phylogenetic inference from the scalidophoran affinity of T. pedum, as well as treptichnid traces in the terminal Ediacaran (Buatois 2018), there is evidence from the fossil record for an Ediacaran history to Ecdysozoa.

Despite the deep Ediacaran divergences estimated for the crown-clades of Cryptovermes and Panarthropoda, our joint posterior intervals for the crown-clades of Arthropoda, Mandibulata and Chelicerata are considerably shorter than for other ecdysozoan phyla (Fig. 6), and are in closer accord with the fossil record, all with minimum estimates in or approaching the early Cambrian. This close accordance of clock-based age estimates and rock-based records for arthropods has been interpreted by other researchers as evidence that the arthropod fossil record closely approximates the origin and diversification of this phylum at the beginning of the Phanerozoic (Daley et al. 2018; Budd and Mann 2020a, b). Finally, by around 535–532 Ma, unequivocal ecdysozoan body fossils are well represented (selected as hard minimum constraints for crown-group Ecdysozoa in this study), in the form of total-group scalidophorans from the Kuanchuanpu and Dengying formations of China (e.g. Liu et al. 2014; Zhang et al. 2015). These are constrained by small shelly fossil biostratigraphy to the Anabarites trisulcatusProtohertzina anabarica Assemblage Zone, Fortunian Stage, Terreneuvian Series (Steiner et al. 2007). These records are within our joint posteriors for Priapulida and Scalidophora.

Our study supports a scenario whereby the ecdysozoan crown-group diverged in the early to mid-Ediacaran. This divergence precedes the first appearance of ecdysozoans in the fossil record, as expected, but the palatability of the offset is dependent on our expectations of the nature of the ancestral organisms and their likely imprint on the fossil record. Our molecular clock estimate allows for an offset that could be as much as 98 myr (636 Ma (upper boundary of combined HPD) minus 538 Ma (approximate minimum age of T. pedum)) or as little as 23 myr (578 Ma (lower boundary of combined HPD) minus 555 Ma (approximate maximum age of Corumbá Group traces; Parry et al. 2017)). The youngest interpretation for the divergence of crown-group Ecdysozoa approximates the age of the oldest well-constrained Ediacaran soft-bodied macrofauna (Matthews et al. 2021) and post-dates the Gaskiers glaciation (c. 580 Ma; Pu et al. 2016). Thus, this younger estimate may be considered to be more geobiologically palatable than the oldest extremes. However, the reason there is so much uncertainty concerning the clade age estimates is because there is no prior objective fossil, phylogenetic or geological evidence to preclude older age estimates, or else it could have been factored into molecular clock calibrations a priori. This considerable range is largely reflective of the quality of the fossil calibrations, with many clades constrained by fossils that are obviously far younger than the clades they are age-constraining. Few fossils are reliably constrained to the stems of the higher-level clades within Ecdysozoa outside Arthropoda and Priapulida, and more broadly in Protostomia and Bilateria. Identifying representatives of such inclusive clades is challenging, as this must assume the presence of a subset of the characters diagnostic of the crown-group, yet there is a lack of consensus regarding what those diagnostic characters are in such large and diverse groups (Cunningham et al. 2017). This equates to a further uncertainty in what our expectations of ecdysozoan fossils from the Ediacaran should look like.

Evaluating molecular clock estimates using reconstructions of ancestral ecdysozoans

An obstacle to assessing the palatability of the clock–rock offset in Ecdysozoa, and other similarly large and diverse clades, is the difficulty in estimating their ancestral anatomical, ecological and behavioural traits. An understanding of the palaeobiology of ancestral crown-group ecdysozoans would be useful in determining whether clock estimates are unrealistically old or not, given our understanding of what sort of organisms the biosphere contained at the supposed time of divergence. However, this requires a well-resolved phylogeny, populated with many fossil stem-representatives of crown-clades to inform the sequence of character acquisitions across lineages. Given its exceptional fossil record, and the increasingly stable phylogeny of the extant phyla, Ecdysozoa may represent the best candidate to scrutinize molecular clock–fossil record offsets among Neoproterozoic metazoans based on ancestral reconstruction.

From current morphology-based phylogenies including Cambrian ecdysozoans, the ancestral crown-group ecdysozoan has been reconstructed as an annulated, vermiform animal with an anterior terminal mouth, a through gut (complete digestive system) and circumoral–pharyngeal armature (Smith and Caron 2015; Ortega-Hernández et al. 2019; Howard et al. 2020b; Yang et al. 2020). This approximates the morphology of the tracemaker of Treptichnus pedum to corroborate an Ediacaran phase of ecdysozoan evolution further, but our more ancient joint posterior estimates (i.e. early Ediacaran) could be interpreted as unpalatable with this reconstruction owing to lack of any bilaterian-grade fossils from earlier than the Ediacaran soft-bodied macrobiotas (Matthews et al. 2021). However, the study by Howard et al. (2020b), which used ancestral character reconstruction to constrain the appearance of ancestral crown-group Ecdysozoa, did not attempt to constrain other important characteristics such as body-size, ecology or biogeographical origin of crown-group Ecdysozoa. Furthermore, the apparent morphological complexity of the ancestral crown-ecdysozoan inferred by Howard et al. (2020b) is subject to the caveats of uncertain homologies across Ecdysozoa for several characters interpreted as ancestral (e.g. circumoral–pharyngeal armature; see Smith and Caron 2015), and their analyses were limited by the relative paucity of fossils constrained to the ecdysozoan stem-group. Future palaeontological studies will reveal more about the ecdysozoan stem-group, which may significantly influence reconstructions of the ancestral crown-group ecdysozoan. Overall, if the ancestral ecdysozoan reconstruction of Howard et al. (2020b) is corroborated by future palaeontological discoveries, then it is at odds with the older end of our joint posterior estimates (636 Ma) where evidence for metazoans is more equivocal; molecular clock calibrations can then be revised in this light. By contrast, the palatability of this reconstruction is less contentious for the younger end of our joint posteriors (578 Ma), which is considerably closer to Ediacaran fossils of likely metazoan, eumetazoan and bilaterian affinity (Cunningham et al. 2017).

Phylogenetic considerations

Our analyses recovered the monophyly of Scalidophora Lemburg (1995), placed as the sister taxon to a clade comprising Nematoida Schmidt-Rhaesa (1996) and Panarthropoda Nielsen (1995). Only one other phylogenetic study has sampled all eight ecdysozoan phyla (Laumer et al. 2019), which was largely consistent with our topology. However, Laumer et al. (2019) recovered scalidophorans as paraphyletic, with Priapulida and Kinorhyncha as successive sister lineages to the rest of Ecdysozoa, and Loricifera as sister to Nematoida (Sørensen et al. 2008), or as sister to Nematoda when Nematomorpha was excluded. Support for Scalidophora was low in our analyses (PP = 0.89) and, as such, the possibility of scalidophoran non-monophyly should not be rejected. Other phylogenomic studies, using a range of data and models, have incompletely sampled Ecdysozoa, typically lacking data for loriciferans, kinorhynchs and/or nematomorphs. However, they have also usually supported a clade comprising a selection of panarthropods and nematoids, although frequently with Tardigrada as sister group to Nematoda (Hejnol et al. 2009; Borner et al. 2014; Laumer et al. 2015), which Laumer et al. (2019) also recovered under some parameters. Of these studies, that by Campbell et al. (2011) (which did not include Loricifera) recovered topologies otherwise consistent with our own from two independent datasets.

Although our topology supports clades originally defined from morphology (Panarthropoda, Nematoida and Scalidophora), phylogenetic analyses of morphological character data typically ally Nematoida to Scalidophora instead of Panarthropoda (Ahlrichs 1995; Dong et al. 2004, 2005, 2010; Harvey et al. 2010; Nielsen 2012; Wills et al. 2012; Liu et al. 2014; Ma et al. 2014; Zhang et al. 2015; Shao et al. 2016). This Nematoida + Scalidophora clade is known as Cycloneuralia Ahlrichs (1995), and is consistently paraphyletic in molecular phylogenies (Hejnol et al. 2009; Campbell et al. 2011; Borner et al. 2014; Laumer et al. 2015, 2019). Only one phylogenomic study has recovered Cycloneuralia and Panarthropoda as sister clades (Dunn et al. 2008, fig. 2), but it did not include Loricifera and in some analyses instead recovered Tardigrada as the sister group to Nematoida (Dunn et al. 2008, fig. 1). However, morphology-based phylogenetic analyses have not adequately sampled the diversity of cycloneuralian phyla, with Loricifera, Kinorhyncha, Nematoda and Nematomorpha all usually represented by only a single terminal with few or no specific characters and states pertaining to these groups (but see Zhang et al. 2015; Shao et al. 2016). This is largely because morphology-based phylogenetic analyses have mostly focused on the position of Cambrian ecdysozoan worms relative to the well-resolved priapulan crown-group. By contrast, phylogenomic studies typically have better taxon sampling and usually recover Cycloneuralia as paraphyletic, suggesting that the scalidophoran-like larva of nematomorphs may reflect the plesiomorphic bodyplan of Nematoida (Sørensen et al. 2008), rather than being indicative of a closer relationship of nematoids to scalidophorans (Nielsen 2001).

Although our study does not support the monophyly of Cycloneuralia as do morphology-based studies, Cryptovermes is not incompatible with palaeontological studies that have attempted to infer ancestral characteristics of the ecdysozoan crown-group detailed in the previous section. These characters include a vermiform body with trunk annulations (Howard et al. 2020b), an anterior terminal mouth (Ortega-Hernández et al. 2019), radially arranged pharyngeal teeth (Smith and Caron 2015; Yang et al. 2020) and an armoured proboscis adorned with circumoral structures (Smith and Caron 2015; Yang et al. 2020). Each of these characters is present in at least the total groups of a subset of the phyla comprising both Scalidophora and Cryptovermes in our tree and, therefore, can be inferred as ancestral for crown-group Ecdysozoa if their homology is accepted, but that is not always clear and is subject to further palaeontological interrogation.

Molecular clock analyses and a phylogenetic interpretation of the fossil record each predict an Ediacaran phase of ecdysozoan evolution, preceding the onset of the Cambrian Explosion. Constraining the extent of this interval is challenging, first, because there is a high degree of uncertainty associated with molecular clock estimates, largely reflecting uncertainties in the interpretation of the fossil record, but also because competing methodological variables yield concomitantly different clade age estimates. Rather than asserting the superiority of one interpretation of the fossil record and one suite of methodological options, we have consolidated the uncertainty associated with all of these experimental variables into one synthetic evolutionary timescale for Ecdysozoa (Fig. 6). This identifies a 68 myr interval for the age of crown-group Ecdysozoa (636–578 Ma). This considerable time span reflects the fact that many calibrations are unlikely to closely approximate the true ages of the clades they are employed to calibrate and the challenge of interpreting absence of fossil evidence as definitive evidence that a clade has not yet evolved (our maximum bounds). Ancestral state estimates of the nature of the ecdysozoan crown-ancestor are tentative but, if accurate, they imply a macroscopic priapulan-grade ancestor (Howard et al. 2020b) that, based on the trace fossil record, is more compatible with the younger end of our joint posterior estimate, which is much closer to Ediacaran metazoan, eumetazoan and bilaterian fossils. This suggests that this lower end of our joint posteriors represents a more plausible age for the divergence of crown-group Ecdysozoa. However, we accept that attempts to formally estimate the nature of the ancestral ecdysozoan are nascent and insufficiently robust to preclude older ages within the range of age estimates for the ecdysozoan crown-ancestor.

The divergence of crown-group Arthropoda exhibits a more precise joint posterior interval (c. 49 myr) compared with other ecdysozoan phyla (>100 myr) and the younger end of the crown-group arthropod divergence interval (c. 540 Ma) precedes the earliest arthropod fossils Rusophycus by only a few million years. This is indicative of the perceived fidelity of the arthropod fossil record. Although wider than that for Arthropoda, the joint posterior interval for crown-group Scalidophora (617–534 Ma) is also less broad than those for the phylum-level clades, and although the older end is at odds with the fossil record, the younger end is virtually concordant with the first total-group scalidophoran body fossils, and actually younger than treptichnids. However, the post-Cambrian record of scalidophorans is limited to one Carboniferous priapulan species (Schram 1973) and so the efficacy of their fossil record cannot be considered comparable with that of arthropods, although their earliest evolutionary history does appear to be well constrained. Finally, our phylogenomic analyses recovered with full support (i.e. posterior probability = 1) the paraphyly of Cycloneuralia, the monophyly of Cryptovermes and the first molecular phylogenetic evidence, albeit with low support (PP = 0.89), for the monophyly of Scalidophora.

We are grateful for the technical support provided by the staff from Station Biologique de Roscoff, France and especially for the crucial help of sea captain G. Maron (R.V. Neomysis) during samplings of marine invertebrates from Trezen ar Skoden, Roscoff. We thank J. Vinther for discussion and help with the early drafts of this paper.

RJH: conceptualization (supporting), formal analysis (supporting), investigation (supporting), methodology (supporting), visualization (equal), writing – original draft (lead), writing – review & editing (equal); MG: conceptualization (supporting), data curation (supporting), formal analysis (lead), investigation (supporting), methodology (supporting), resources (supporting), software (equal), visualization (supporting), writing – original draft (supporting), writing – review & editing (supporting); JL-F: conceptualization (supporting), data curation (lead), formal analysis (supporting), investigation (supporting), methodology (supporting), resources (equal), software (equal), supervision (equal), visualization (supporting), writing – original draft (supporting), writing – review & editing (supporting); GDE: investigation (supporting), supervision (supporting), writing – original draft (supporting), writing – review & editing (supporting); JFF: data curation (supporting), resources (supporting); RMK: resources (supporting), writing – review & editing (supporting); XM: supervision (supporting), writing – original draft (supporting), writing – review & editing (supporting); JO: resources (supporting), writing – review & editing (supporting); MVS: resources (supporting), writing – review & editing (supporting); PFT: resources (supporting), writing – review & editing (supporting); MAW: supervision (supporting), writing – original draft (supporting), writing – review & editing (supporting); PCJD: conceptualization (supporting), formal analysis (supporting), funding acquisition (supporting), investigation (equal), methodology (equal), project administration (equal), software (supporting), supervision (lead), writing – original draft (supporting), writing – review & editing (supporting); DP: conceptualization (lead), data curation (supporting), formal analysis (supporting), funding acquisition (lead), investigation (equal), methodology (equal), project administration (lead), resources (equal), software (supporting), supervision (supporting), writing – original draft (supporting), writing – review & editing (supporting)

The NERC GW4 + Doctoral Training Partnership provided stipend and research expenditure to R.J.H. D.P. and M.G. were supported by the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement 764840. D.P. and P.C.J.D. were funded by Natural Environment Research Council (NERC) grant NE/P013678/1, part of the Biosphere Evolution, Transitions and Resilience (BETR) programme, which is co-funded by the Natural Science Foundation of China (NSFC). P.C.J.D. is funded by Biotechnology and Biological Sciences Research Council (BBSRC) grants BB/N000919/1 and BB/T012773/1. R.M.K. was funded by the Danish Natural Science Research Council (grant FNU 272-08-0576) and Carlsberg Foundation (grants CF 970345/30-488, CF 2009_01_0063, CF 2012_01_0123 and CF 16-0236). J.O. was supported by a grant from the Danish Agency for Science, Technology and Innovation (0601-12345B). M.A.W acknowledges BBSRC grant BB/K006754/1 and JTF Grant 61408. NERC Independent Research Fellowship (Grant 680 NE/L011751/1) provided salary and research expenditure for X.-Y.M. J.F.F was funded by Royal Society–Japan Society for the Promotion of Science KAKENHI grant 18F18788.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

The data associated with this paper can be found in the following link: https://doi.org/10.6084/m9.figshare.c.5811381. The genomes and transcriptomes generated as part of our study are available at NCBI through the following accession numbers: SAMN24271371, SAMN24271372 and SAMN24271374.

Scientific editing by Maoyan Zhu