The two major sources of heterochronous sequence data are rapidly evolving RNA and DNA viruses, whose high mutation rates enable the generation of phylogenetically informative sequence diversity within historical time frames (Drummond, Pybus, and Rambaut 2003), and ancient DNA data isolated from preserved remains (Hofreiter, Serre, et al. 2001) that may be up to several hundred thousands of years old (Willerslev et al. 2007). In both cases, the period over which sequences are isolated is sufficiently long relative to the mutation rate to allow estimation of the evolutionary rate (Drummond, Pybus, Rambaut, Forsberg, and Rodrigo 2003).

The temporal information associated with RNA virus sequences typically represents the day, month, or year of sample isolation and/or storage (Taubenberger et al. 1997). Heterochronous viral data sets have been used to estimate rates of mutation for specific viruses (Jenkins et al. 2002), to investigate rates of evolution and adaptation within hosts (e.g., Lemey et al. 2006, 2007), and to infer epidemic dynamics within and between populations of susceptible hosts (Rambaut 2000; Pybus and Rambaut 2009).

For ancient DNA data sets, ages of the genetic sequences are most often approximated using radiocarbon dates that are estimated from the same specimens from which the DNA sequences are amplified (e.g., Shapiro et al. 2004; Bunce et al. 2009; Campos et al. 2010). Dates from material associated with stratigraphic context have also been used as calibrating information, however (e.g., Lambert et al. 2002; Valdiosera et al. 2008). For ancient DNA sequences, leaf-ages are normally assigned some number of thousands of years before the present (ka BP). Molecular clock analyses of these data have been used to identify periods of population turnover (e.g., Hadly et al. 1998; Barnes et al. 2002; Hofreiter et al. 2004) to estimate divergence times within species (e.g., Shapiro et al. 2004; Debruyne et al. 2008; Stiller et al. 2010) and to correlate population-level changes in genetic diversity with external events, such as climate change (e.g., Hadly et al. 2004; Chan et al. 2006; Barnett et al. 2009; Campos et al. 2010). For example, there has been considerable debate about the relative roles of climate change associated with the last glacial maximum (LGM; ca 21 ka BP) and the appearance and increase in human populations (ca 14 ka BP) in the recent disappearance of the North American megafauna (Alroy 2001; Barnosky et al. 2004; Stuart et al. 2004). Heterochronous data sets comprising sequences directly dated to before, during, and after these events allow explicit tests of these alternative hypotheses.

Despite significant technical and chemical pretreatment advances in Accelerator Mass Spectrometry radiocarbon dating (Bronk Ramsey et al. 2004), the oldest samples for which finite radiocarbon dates can be routinely generated are around 50-–55 ka BP (e.g., Barnett et al. 2009). The period 0–55 ka BP includes several specific large-scale environmental events that are likely to have affected the distribution and abundance of plant and animal species. However, there are a number of circumstances under which the leaf-ages of ancient DNA or viral sequences may be unknown, or at best, highly uncertain. First, ancient mitochondrial DNA (mtDNA) sequences are routinely amplified from permafrost-preserved specimens older than the 50–55 ka BP radiocarbon limit. For example, nearly 100 of the bison sequences reported in Shapiro et al (2004) were too old to be assigned finite radiocarbon ages, and ancient DNA sequences have been reported that are perhaps as old as 800 ka BP (Willerslev et al. 2007). In this case, only censored temporal information is often available (i.e., age > 55 ka BP). Second, ancient DNA samples are routinely recovered from situations in which the stratigraphic context provides calibrating information (e.g., Lambert et al. 2002; Coolen and Overmann 2007; Valdiosera et al. 2008) but only within a wide range of uncertainty, such that assigning a specific mean or median date to such sequences is statistically inappropriate (Ho and Phillips 2009). Third, for rapidly evolving viral sequences, the date of sampling may simply be unknown due to the loss or absence of accurate archival information. Even if the viral sampling date is known to the nearest year, it may be important to know the isolation date more accurately. Fourth, it may also be important to independently verify posited sampling dates due to their extreme age (Zhu et al. 1998; Taubenberger et al. 2005; Worobey et al. 2008) or because they are in some way contentious (Sonoda et al. 2000; Coolen and Overmann 2007). Because frozen viral isolates do not accumulate mutations while in storage, a leaf-dating method also has the potential to identify transmitting viruses that, after a period of storage, have been released into the environment (Worobey 2008).

Despite these obvious problems, few studies have attempted to estimate the unknown age or sampling date of heterochronous sequences using molecular clock methods, and none have investigated the statistical reliability of such methods. Perhaps the most similar analysis was that undertaken by Korber et al. (1998), who validated their molecular clock of HIV-1 by testing whether it could accurately predict the date of an “old” isolate sampled in 1959. However, in that case, the approach consists of the visual inspection of the fit to a linear regression of viral sampling date against genetic distance rather than a statistical analysis or estimation procedure.

In previous molecular clock implementations, all nodes in the tree are given dates/ages, that is, internal nodes are treated as unknown parameters, although the dates/ages of the tree external nodes (leaves) are assumed to be known. To estimate the age of an individual sequence, we extend the framework introduced above to estimate the time associated with the sequence, that is, the sequence’s leaf-age is treated as a random variable. Thus, an additional parameter for the age of the external node is introduced and is treated identically to the internal node age parameters in terms of proposals made by the MCMC kernel. See Drummond et al. (2006) for further specifications and parameterizations of the molecular clock models used here.

The results of the empirical data analyses were mixed. Our leaf-dating method performed very well when applied to the DEN-2 data set, as the known dates of sampling fell within the 95% HPD intervals of our estimates 93% of the time. This is extremely encouraging, as it is almost certainly true that the evolutionary and coalescent models used in our analyses are a gross simplification of natural processes; yet, we were able to estimate leaf-dates reliably and accurately. In addition to demonstrating the effectiveness of our leaf-dating approach, our DEN-2 virus result also serve to validate the molecular clock methods used, for if terminal nodes throughout the tree are estimated with accuracy then it logically follows the internal nodes, including the tree root, are similarly dated correctly.

For the bison data set, the leaf-ages (as determined by radiocarbon dating) were estimated incorrectly approximately 15% of the time (table 1, supplementary table S2, Supplementary Material online). There was no noticeable pattern in the age of the tips for which the leaf-dating method fails. However, in general, the same leaf-ages, which differ significantly from their priors in all cases, are incorrectly estimated under both coalescent models. This suggests two things: first, that the incorrect dates likely arise from some heterogeneity in the evolutionary process that is not accounted for in the models and second, that the leaf-date estimates depend more on the data than on the prior distributions assumed. The MCP process assumes that effective population sizes change according to an exponential Markovian process (Drummond et al. 2005). In rare cases, with limited information about the root age that occurs when the unknown dated leaf attaches near the root, this nonstationary process can introduce excess variability (fig. 3A); the GMRF, on the other hand, is stationary (Minin et al. 2008) and modestly outperforms the MCP in these situations.

Although the vast majority of bison sequence ages were estimated correctly, it is unsatisfactory that some were not and important to consider the possible causes of such error. First, and most obviously, the “true” ages of the specimens may have been incorrectly reported or estimated. For example, the radiocarbon ages we used in our analysis may be incorrect or subject to error, in particular as the age of the specimen approaches the limits of this technique.

Second, the sequences themselves may be incorrect or be erroneous in some way. This is particularly problematic for ancient DNA data, where the degraded nature of the samples results in predictable patterns of DNA damage and base misincorporation by the polymerases used in polymerase chain reaction (Hofreiter, Jaenicke, et al. 2001; Gilbert et al. 2003, 2005; Binladen et al. 2006). Previously, an analysis of the bison data set incorporating the postmortem damage (PMD) model showed that this particular aDNA data set contained very little damage (Rambaut et al. 2009). The PMD model assumes that the probability that any given nucleotide remained undamaged decays exponentially with age; the oldest sequences in the bison data set demonstrated, on average, only 0.74 damaged sites, a level that produced no qualitative change in the demographic reconstructions. In its current implementation, the PMD model assumes a common decay rate for every site of every sequence (such that the probability of a site being undamaged decays exponentially with age). A useful extension of the PMD model may be to provide a probabilistic expectation that each individual sequence is damaged, thereby making it possible to identify problematic sequences for additional assessment. Although incorporating both PMD models and leaf-dating in a single analysis is possible, the common decay rate PMD model and leaf-dating are expected to be only weakly identifiable in the sampling density of the observed sequence data, serving primarily to increase variance on the leaf-age estimate without extracting much additional information from the data. Furthermore, the individual decay rate PMD model and a random leaf-age are not identifiable. Consequently, prior assumptions will dominate inference as to whether a sequence is damaged or incorrectly dated.

Third, some aspect of the evolutionary models used may be unrealistic and a possible source of error. We consider each model component in turn:

Note that the discussion above does not directly address the biological assumptions underlying model components, such as the coalescent and the molecular clock. This is deliberate, as we wish to highlight a common misconception. In population-level analyses, it is often assumed that it is necessary to assume neutral evolution in order to accurately estimate divergence times using a molecular clock. If divergence times are the primary parameters of interest, then the clock model used can be thought of as a statistical or phenomenological description of the relationship between genetic distance and time rather than an explicit model of a biological process. Viewed this way, all that is important is that the model “fits” the data well, statistically speaking. The same argument can be applied to the coalescent model, which, when being treated as a nuisance parameter, need only represent a suitable range of tree shapes and sizes. In such cases, the assumptions of the coalescent model (e.g., random sampling or panmixis) are not necessary conditions for accurate date estimation. As a result, estimation of leaf-dates is likely to be highly robust provided that the clock model used incorporates sufficient rate heterogeneity.

Although it is often not feasible to discriminate among the sources of error outlined above, it may be possible to identify problematic sequences by evaluating the output of the leaf-dating analysis. For example, a common problem among the small number of analyses that fail to recover the true leaf-age is that the sequences appear to be equally or nearly equally likely to fall in two locations in the genealogy, with each of these resulting in very different leaf-age estimates. This pattern is seen clearly by plotting the estimated marginal posterior distribution leaf-ages, which has an unusually high variance (fig. 3A). This result may be due to errors within the sequence itself, which could potentially be resolved by resequencing. Alternately, very precise leaf-age estimates not containing the true age may indicate a problem with the true age (fig. 3B). For ancient DNA data, recovering an additional radiocarbon date or confirming information about the stratigraphic context of the sample can be useful to rule out this potential source of error. For viral sequences, a reexamination of the documentation associated with the isolate may reveal an annotation or transcription error.

Although the estimated ages recovered by the leaf-dating method are often associated with wide credible intervals, our leaf-dating method provides a means to include in molecular clock analyses data for which little or no temporal information is known. Incorporating additional sequence data can improve the resolution of the phylogenetic, demographic, and geographic history of the sampled sequences and can extend significantly the temporal range of the analysis. Additionally, estimating leaf-ages can provide an independent means of assessing both the authenticity of heterochronous sequences and the ages to which the sequences have been ascribed, which is often a significant concern in ancient DNA research.

The possibility of treating leaf-ages as random variables also allows uncertainty to be modeled explicitly. This enables the incorporation of the uncertainty associated with layer dating, for example, in the form of a uniform prior across the age range of the source stratum. In addition, the error in isotopic dating can be reflected by choosing an appropriate prior distribution for the corresponding leaf-age (Ho and Phillips 2009).