Testing the molecular clock using mechanistic models of fossil preservation and molecular evolution

Molecular sequence data provide information about relative times only, and fossil-based age constraints are the ultimate source of information about absolute times in molecular clock dating analyses. Thus, fossil calibrations are critical to molecular clock dating, but competing methods are difficult to evaluate empirically because the true evolutionary time scale is never known. Here, we combine mechanistic models of fossil preservation and sequence evolution in simulations to evaluate different approaches to constructing fossil calibrations and their impact on Bayesian molecular clock dating, and the relative impact of fossil versus molecular sampling. We show that divergence time estimation is impacted by the model of fossil preservation, sampling intensity and tree shape. The addition of sequence data may improve molecular clock estimates, but accuracy and precision is dominated by the quality of the fossil calibrations. Posterior means and medians are poor representatives of true divergence times; posterior intervals provide a much more accurate estimate of divergence times, though they may be wide and often do not have high coverage probability. Our results highlight the importance of increased fossil sampling and improved statistical approaches to generating calibrations, which should incorporate the non-uniform nature of ecological and temporal fossil species distributions.

The age of the root is fixed at 100 Ma, which is defined as one time unit. The time period between the present and the age of the root was subdivided into 50 stratigraphic intervals of equal length (2 Myr). During each interval, each extant lineage was collected with a probability that is the product of sampling probability p and the fraction of the interval during which the lineage is extant. The sampling probability p here reflects the joint effects of fossil preservation potential and sampling intensity, effects that are indistinguishable in such a model. Under the uniform model of preservation p is simply equal to the specified sampling intensity, s. We used four values: s = 0.001, 0.01, 0.1, and 1, to reflect the perceived completeness of the fossil record [1,2].
To simulate non-uniform occurrence data we used a three-parameter Gaussian model of preservation [3,4], which uses water depth as a proxy for fossil preservation or sampling potential in the marine stratigraphic record. The sampling probability is given by where d is the current water depth, PD is preferred depth, DT is depth tolerance and PA is peak abundance of the species. PD and DT are the mean and standard deviation of the distribution, respectively. Water depth (d) was simulated using a simple sine wave function d(t) = .
(2) substitution/site/unit time or 10 -8 substitutions/site/year. Then independent rates for branches on the tree were sampled from a lognormal distribution with mean rate μ i and standard deviation of the log rate σ = 0.1 (e.g. [6]). This is the independent-rates model and allows variable rates both among multiple loci and among branches at each locus. Branch lengths, in the expected number of substitutions per site, were calculated as the product of the time duration of the branch and the rate. Given the tree topology and branch lengths, the HKY85+ Γ5 substitution model was used to 'evolve' the sequences along the tree, with equal base frequencies, a transition to transversion rate ratio κ = 5, and the gamma shape parameter α = 0.25 for rate heterogeneity across sites. The number of replicate datasets is 100, and each replicate dataset consists of L loci.

(b) Minimum and maximum constraints on divergence times
The simulated fossil occurrence data was used to establish minimum (t L ) and maximum (t U ) constraints on ages of nodes on the tree, which were used as calibrations in Bayesian estimation of divergence times (figure S1). In all cases we assumed no topological uncertainty in the placement of fossils or the relationships among modern terminals. Minimum constraints were based on first appearances, using the youngest age of the time bin from which the fossil was recorded. Three approaches were used to establish maximum constraints.
First, we used a stratigraphic bracketing based approach, which uses the density of fossil occurrences over a given interval to estimate 95% confidence intervals on stratigraphic ranges ( [7], equation 14 -modified from [8]), where b is the number of lineages with a fossil record, and H is the number of fossil localities.
For a given node in each tree, H was calculated using the combined occurrence data of the descendent branches that may be ancestral to younger nodes in the phylogeny. In equation 3, t L is given by the age of the oldest fossil. Second, phylogenetic bracketing was used to emulate best-practice approaches of establishing calibrations (e.g. [9][10][11]). Maxima were established based on the oldest age interpretation of the minimum constraints of ancestral clades. For a given node j, if the immediate ancestor did not have a minimum constraint older than the minimum of j, the next nearest ancestor with a minimum older than j was used to inform the maximum of j. Maxima were not established for nodes that did not have a minimum constraint. This approach will not produce a maximum constraint for the root of the tree. Third, we generated the maximum bounds to be 110, 125, 150 and 175% the age of the minimum constraints. In all cases, the age of the fossil used to inform the maximum constraints was equal to the maximum age of the interval from which the fossil was sampled. . The minimum and maximum constraints were generated using stratigraphic bracketing, phylogenetic bracketing and arbitrary maxima (at 110% and 150%). Minimum constraints were always based on the youngest secure age interpretation of first appearances.
The fossils used to inform the maximum constraints are highlighted in the trees, including stratigraphic bracketing or phylogenetic bracketing, which uses the oldest secure age of the first appearance of the nearest ancestral clade. The age interpretation of the fossil specimens used to define minimum or maximum constraints are indicated by the dashed grey line. The minimum and maximum constraints were used to construct uniform or skew-t prior densities in Bayesian molecular clock analysis, shown right, in (b) and (c). In these plots, the solid grey line indicates the true node age.   The methods of analysis include arbitrary maxima (at 110% and 150%), stratigraphic bracketing, and phylogenetic bracketing, using skew-t calibration densities.  Figure S9-S12. Infinite-sites plots for data simulated on the unbalanced tree under the nonuniform (S8) and uniform (S9) preservation models and the balanced tree under the nonuniform (S10) and uniform (S11) models and analyzed using different calibration approaches (arbitrary maxima at 110% and 150%, stratigraphic bracketing and phylogenetic bracketing).
Non uniform pres.
Uniform pres.           Arbitrary maxima (110%) -skew-t priors  Arbitrary maxima (125%) -skew-t priors  s is the probability of sampling during each interval under the uniform preservation model. % is the proportion of calibrated nodes across all simulated replicates. CP is the coverage probability across all replicates. RMSE is the relative root mean squared error andŵ is the relative confidence interval width, both averaged across all replicates. The tree is shown in Fig. 1.   Phylogenetic bracketing -skew-t priors  Arbitrary maxima (125%) -skew-t priors  Arbitrary maxima (175%) -skew-t priors  PA is the peak abundance parameter that determines the probability of sampling during each interval under the non-uniform preservation model. % is the proportion of calibrated nodes across all simulated replicates. CP is the coverage probability across all replicates. RMSE is the relative root mean squared error andŵ is the relative confidence interval width, both averaged across all replicates. The tree is shown in Fig. 1.   Phylogenetic bracketing -skew-t priors  Arbitrary maxima (110%) -skew-t priors  Arbitrary maxima (125%) -skew-t priors  Arbitrary maxima (150%) -skew-t priors  s is the probability of sampling during each interval under the uniform preservation model. % is the proportion of calibrated nodes across all simulated replicates. CP is the coverage probability across all replicates. RMSE is the relative root mean squared error andŵ is the relative confidence interval width, both averaged across all replicates. The tree is shown in Fig. 1.     Phylogenetic bracketing -skew-t priors  Arbitrary maxima (110%) -skew-t priors  Arbitrary maxima (125%) -skew-t priors  Arbitrary maxima (150%) -skew-t priors  Arbitrary maxima (175%) -skew-t priors  PA is the peak abundance parameter that determines the probability of sampling during each interval under the non-uniform preservation model. % is the proportion of calibrated nodes across all simulated replicates. CP is the coverage probability across all replicates. RMSE is the relative root mean squared error andŵ is the relative confidence interval width, both averaged across all replicates. The tree is shown in Fig. 1.