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.