Results and discussion
Molecular clock analysesnext section
Test of model fit
Results of the model selection based on Bayes factors (BF, Table 1) indicated that the best fitting model employed a combination of a relaxed exponential clock and a Birth-Death model of tree diversification. We used this combination of priors for subsequent analyses (Tables 2, 3). The difference in BF to other model combinations was quite low (BF 2.4 over the second ranking model), but enough to set this model as favoured. In general, model combinations involving a relaxed exponential were “strongly” favoured over the relaxed lognormal, the latter being in turn “very strongly” favoured over the strict and the random clock models. As for the tree priors, the difference between the Coalescent and the three Birth-Death type processes we tested was very low.
Because the difference of the marginal likelihood using different priors was extremely low, we advocated great care in interpreting the result of our model selection: Table 1 is however a good indication for future tardigrade clock studies using RNA sequences, rather than as decision maker to discriminate among competing time estimates.
Table 2. Molecular dating using two different sets of calibration priors. Analyses were conducted using the complete dataset and the most fitting relaxed clock and tree prior as defined in Table 1 (relaxed exponential clock plus Birth Death process). For each of the calibration sets, we provide the mean estimates in millions of years for four nodes of interest; the heights of the 95% Highest Posterior Density are in parenthesis.
Table 3. Molecular clock analyses using reduced-datasets. Analyses were conducted using the most fitting relaxed clock and tree prior as defined in Table 1 (relaxed exponential clock plus Birth Death process). For each of the calibration-sets, the mean estimates in millions of years for four nodes of interest are provided; the heights of the 95% Highest Posterior Density are in parenthesis.
Our results for the tardigrade phylogenetic relationships, based on the best fitting priors deduced from Table 1 and using calibration-set 1 (Fig. 2), were congruent with previous studies (Jørgensen et al., 2011; Guil and Giribet, 2012; Guil et al., 2013; Vicente et al., 2013; Guidetti et al., 2014). By employing different calibration sets (Table 2), and testing the effect of site selection (Table 3), we have attempted to account for the uncertainty imbedded in molecular clock studies.
The analyses conducted using two different calibration sets (Table 2) provided different estimates. Nonetheless, from the results, it was possible to infer the root of the tree, i.e. when the Heterotardigrada and Eutardigrada lines split apart, showing to be at 659 and 547 Mya, which was before the Cambrian Period (Table 2). The two methods also produced similar mean estimates, 146 and 136 Mya (early Cretaceous), for the divergence of the lineage that led to Mopsechiniscus, though the 95% Highest Posterior Density (95% HPD) for this node was very wide (ranging from 294 to 38 Mya). Both methods placed the split between the phylogenetic lines of the Mopsechiniscus species (M. franciscae: Antarctica; Guidetti et al., 2014 and M. granulosus: Chile; Jørgensen et al., 2011) between 47.8 and 32.1 Mya (Paleogene, during late Eocene early Oligocene), with a 95% HPD of 131–3 Mya (Table 2; Fig. 2). Analyses conducted on the reduced-dataset produced similar estimates (Table 3; Fig. 3) with average data in line with: a possible pre-Cambrian origin of the Tardigrada; an early Cretaceous-Eocene origin for the lineage that led to Mopsechiniscus; a late Paleogene split between M. franciscae and M. granulosus.
Fig 3. a Molecular clock analysis of the Heterotardigrada using the rate deduced from Rota-Stabelli et al. (2013) and a Log-normal distribution clock analysis. b Molecular clock analysis of the Heterotardigrada using secondary node constraints from Regier et al. (2004). Values close to branches indicate posterior probability values (> 0.9) of BI analysis while bars indicate the 95% credibility intervals.
Although the mean estimates were concordant for various nodes, particularly those describing the origin of M. franciscae, the range of visited posterior estimates was extremely large in all the analyses (see Tables 2 and 3 for a breakdown of the 95% HPD): overall, this indicated a large uncertainty in the precise estimation of tardigrade radiation using RNA 18S and 28S makers. This may be caused by a variety of reasons. First, the paucity of calibration priors currently available for tardigrades forced us in relying mostly on prior replacement rates with relatively high standard deviations: this likely reflected in highly uncertain posterior estimates. Second, the relatively short length of the alignment and the high amount of missing data may have inserted a high stochastic effect. This possibility was reinforced by the range of the HPD being slightly smaller when using a reduced (less missing data) dataset (compare Tables 2 and 3). Finally, the fast evolving nature of tardigrade genes (Campbell et al., 2011) may have complicated the correct estimation of their divergence; this possibility is compatible with a study of Ecdysozoan evolution showing that tardigrade divergence was the most unstable within the sample, with estimates strongly varying with parameter variations (Rota-Stabelli et al., 2013). Overall, our molecular dating of the Heterotardigrada indicated that their divergence estimate is a complex issue which should be tackled in the future by employing more markers and possibly outgroups to the tardigrades in order to allow external calibrations.