A Compound Poisson Process for Relaxing the Molecular Clock

Corresponding author: John P. Huelsenbeck, Department of Biology, University of Rochester, Rochester, NY 14627., johnh{at}brahms.biology.rochester.edu (E-mail)

Communicating editor: S. TAVARÉ

	ABSTRACT

TOP ABSTRACT METHODS RESULTS DISCUSSION APPENDIX LITERATURE CITED

The molecular clock hypothesis remains an important conceptual and analytical tool in evolutionary biology despite the repeated observation that the clock hypothesis does not perfectly explain observed DNA sequence variation. We introduce a parametric model that relaxes the molecular clock by allowing rates to vary across lineages according to a compound Poisson process. Events of substitution rate change are placed onto a phylogenetic tree according to a Poisson process. When an event of substitution rate change occurs, the current rate of substitution is modified by a gamma-distributed random variable. Parameters of the model can be estimated using Bayesian inference. We use Markov chain Monte Carlo integration to evaluate the posterior probability distribution because the posterior probability involves high dimensional integrals and summations. Specifically, we use the Metropolis-Hastings-Green algorithm with 11 different move types to evaluate the posterior distribution. We demonstrate the method by analyzing a complete mtDNA sequence data set from 23 mammals. The model presented here has several potential advantages over other models that have been proposed to relax the clock because it is parametric and does not assume that rates change only at speciation events. This model should prove useful for estimating divergence times when substitution rates vary across lineages.

THE molecular clock hypothesis states that the evolutionary rate of a gene is roughly constant among different lineages (ZUCKERKANDL and PAULING 1962 ). ZUCKERKANDL and PAULING 1965 also suggested that the substitution process is approximately Poisson. If the rate of substitution is constant across lineages, then the distances between species should be ultrametric (i.e., all tips are an equal distance from the root of the tree). Furthermore, if the substitutions follow a Poisson process, then the variance and the mean of the number of substitutions that occur on different lineages in the same amount of time should be equal. Since the early 1970s, however, neither prediction has been shown to hold true; the variance to mean ratio of the number of substitutions is generally greater than one, suggesting that the substitution process is overdispersed (OHTA and KIMURA 1971 ; LANGLEY and FITCH 1973 , LANGLEY and FITCH 1974 ). Moreover, rates of substitution have been shown to vary across lineages (see GILLESPIE 1991 ).

Despite the observation that the molecular clock hypothesis does not perfectly explain the substitution process, it remains a powerful conceptual and analytical tool in evolutionary biology. Conceptually, the molecular clock provides a timescale, albeit an imperfect one, for evolution as well as a mechanistic description of the substitution process. As an analytical tool, the molecular clock has also proven useful. Many of the more interesting applications of phylogenies assume that the molecular clock holds. For example, statistical analysis of host-parasite cospeciation often assumes a molecular clock (HUELSENBECK et al. 1997 ). Also, the estimation of divergence times using molecular data depends upon (1) an accurate calibration date for at least one speciation event on the tree and (2) a constant substitution rate among lineages. Several recent studies have attempted to date the divergence times for eubacteria/eukaryotes (DOOLITTLE et al. 1996 ), metazoa (WRAY et al. 1996 ), birds (COOPER and PENNY 1997 ), mammalian orders, and major lineages of vertebrates (KUMAR and HEDGES 1998 ) using the molecular clock assumption. To a large extent, the validity of any such analysis depends on how well the data conform to the clock assumption.

Several different approaches have been taken to accommodate rate variation among lineages. One method, commonly used for estimating phylogenetic trees, is to assign each branch of a phylogenetic tree its own rate parameter. However, this procedure does not allow estimation of divergence times of clades because rate and time are confounded. Another approach was proposed by SANDERSON 1997 , who used a nonparametric method for smoothing the rate differences across speciation events on the tree. Sanderson's method allows rates to vary on branches of the tree and also allows divergence times to be estimated. At the other extreme, THORNE et al. 1998 proposed a parametric model for relaxing the clock. Like SANDERSON 1997 , they assume that rates are autocorrelated across speciation events; their model assigns new rates to descendent lineages from a lognormal distribution with the mean of the distribution equal to the rate of the ancestral lineage. The model we present here differs from the models proposed by SANDERSON 1997 and THORNE et al. 1998 by allowing rates to change anywhere on the tree. Yet, the model introduces only two additional parameters over the strict molecular clock.

In this article, we propose a compound Poisson process for introducing rate variation across lineages on a tree. We assume that nucleotide substitutions occur along branches of the tree according to a Poisson process, as in the Markov substitution models widely used for phylogenetic inference. However, an independent Poisson process also generates events of substitution-rate change. At each of these events, the rate of substitution is changed by multiplying the current rate by a gamma-distributed random variable. We use Markov chain Monte Carlo when making Bayesian inferences.

	METHODS

TOP ABSTRACT METHODS RESULTS DISCUSSION APPENDIX LITERATURE CITED

A compound Poisson process of rate variation:
We assume that the phylogeny of a group of species can be represented by a rooted binary tree, an arbitrary example of which is shown in Fig 1. The tips are labeled n₁ to n_s, and the internal nodes are labeled n_s+1 to n_2s-1, where s is the number of sequences. The root of the tree is always labeled n_2s-1. The times of the nodes are denoted t = (t₁, t₂, ... , t_2s-1). The branches are scaled such that the tips are at time 0 (t₁ = t₂ = ... = t_s = 0) and the root is at time 1 (t_2s-1 = 1). The sum of the scaled branch lengths is referred to as the total length of the tree, T. The total length of the tree of Fig 1, for example, is T = 4.55. If the branch lengths are multiplied by a parameter m, representing the expected number of substitutions per site on a single branch reaching from the root of the tree to the tips, the branch lengths of the tree are in terms of expected number of substitutions per site. Moreover, the branch lengths on the tree conform to a molecular clock; there has been no variation in rates across lineages with the result that all of the tips can be drawn to lie on a single time line. Fig 2 shows a rescaling of the tree of Fig 1 with m = 0.17.

View larger version (18K): In this window In a new window Download PPT slide   Figure 1. We assume that a rooted binary tree describes the genealogy of the sequences n1 to ns. This arbitrarily chosen tree illustrates the labeling of the nodes used in this article. The unscaled times of the tree have the tips at time 0 and the root at time 1. All other nodes on the tree have times between 0 and 1.

View larger version (18K): In this window In a new window Download PPT slide   Figure 2. The branch lengths that result when the tree of Fig 1 has all branches multiplied by m = 0.17. The rate of substitution on the tree, m, is in terms of expected number of substitutions per site that occur along a single lineage reaching from the root of the tree (t = 1) to the tips (t = 0). Multiplying the branch lengths in Fig 1 by m converts branch lengths into expected number of substitutions per site. Note that this tree has no rate variation among lineages. The branch lengths of the tree conform to a molecular clock.

We use a compound Poisson process to introduce rate variation across lineages of the tree. Events (positions on the tree at which substitution rate changes) occur according to a Poisson process with parameter . When an event occurs, the rate of substitution just prior to the event (m) is multiplied by a gamma-distributed random variable (taking value r) to produce a new substitution rate above the event (m'; m' = mr). The gamma distribution has density

(1)

where is the shape parameter and ß is the scale parameter. One possible parameterization of the gamma distribution that we initially considered is to set = ß = _P, so that r is (r|_P, _P). The mean rate of change, then, is E[r] = 1. A potentially unappealing property of such a model of rate change is that as the branch lengthens, substitutions stop happening along the branch. In such a model, the rate along a branch is a base rate, m, multiplied by ⁿ_ir_i, where n is the number of Poisson change points and r_i is (_P, _P). As the branch length grows, n becomes large and the product converges to 0 in probability.

A prior model that changes rates multiplicatively with independent multipliers must be carefully chosen so that the rate does not tend to either 0 or infinity in probability. To see this, let r_i, i = 1, 2, ... be independent and identically distributed positive random variables and let R_n = ⁿ_i=1 r_i. Then for every > 0,

Thus, the product of rates converges in probability to 0 if E[log r_i] < 0 and increases without bound in probability if E[log r_i] > 0.

In this article, we largely circumvent this problem by placing a prior on such that the number of change points on the tree does not become large. In addition, we define a one-dimensional family of gamma-distributed random variables with the property that E[log r_i] = 0 by letting ß = e^(_P), where is the derivative of the logarithm of the gamma function:

Here is a short derivation of this property. Suppose that r ~ (,ß). Then r is equal in distribution to X/ß where X ~ (,1):

This is 0 precisely when ß = e⁽⁾.

In this article, we multiply the rate at an event by a gamma-distributed random variable with density

(2)

The collection of events is denoted e = (, z, r), where is the number of rate change events, z is the vector of event positions, and r is the vector of rate multipliers. When = 0, z and r are empty. For > 0, z [0,T] and r . There is a mapping that takes points from [0,T] to points on the tree. Fig 3 illustrates the compound Poisson process for rate change events acting upon the tree of Fig 1. In this example, three events of rate change occur on the tree [e = (3, z, r), z = (1.354, 3.009, 4.301), r = (2.871, 1.103, 0.776)]. Associate each branch with the smaller number of its two nodes (that node that is furthest from the root of the tree) and partition the interval [0, T] into segments associated with branches ordered from the lowest index to the highest index. For example, the branch descending from node n₁ is represented by the interval [0, 0.125]. Similarly, the branch descending from node n₆ is represented by the interval (0.979, 1.416]. The first event is located on this branch at z₁ = 1.354. Events 1 and 2 increase the rate of substitution whereas event 3 decreases the rate of substitution. Fig 4 shows the relative branch lengths in terms of expected number of substitutions per site. Note that the branch lengths of the tree no longer obey the molecular clock. The molecular clock is simply a special case of our model with = 0 and/or _P = .

View larger version (18K): In this window In a new window Download PPT slide   Figure 3. The compound Poisson process discussed in this article places events of substitution-rate change on the tree according to a Poisson process. Associated with each event is a rate multiplier that multiplies the rate of substitution up to the event (m) by a gamma-distributed random variable (taking value r) to produce a new rate of substitution above the event (m'). In this example, which illustrates the process for the tree from Fig 1, three events of substitution-rate change occur. The rate multipliers are 2.871, 1.103, and 0.776 for events 1, 2, and 3, respectively. The zi specify the event locations on the tree.

View larger version (17K): In this window In a new window Download PPT slide   Figure 4. The branch lengths in terms of expected number of substitutions per site when the substitution rates are modified using the compound Poisson process. The rates were modified using the three events from Fig 3 with a starting substitution rate of m = 0.17 (at the root of the tree). Note that the branch lengths of the tree no longer conform to the molecular clock.

The prior distribution of e = (, z, r) given the tree and parameters _P and is described by the probability measure

(3)

where ₀ is the point mass measure at (0, , ). This corresponds to picking a Poisson (T) number of points, choosing their locations from mutually independent uniform [0,T] distributions, and attaching mutually independent gamma-distributed rate multipliers.

Model of DNA substitution:
We assume that DNA sequences from homologous regions are available for species n₁ to n_s. Let X = {x_ij} be the aligned nucleotide sequences, where i = 1, 2, ... , s; j = 1, 2, ... ,c; and c is the number of nucleotide sites per sequence. Each column of the data matrix x_j = {x_1j, ... , x_sj}' specifies the nucleotides for the s sequences at the jth site.

As is usual for DNA sequences, we assume that substitutions occur according to a Poisson process with rate matrix Q. A general reversible rate matrix allows a different stationary frequency for the four nucleotides = (_A, _C, _G, _T) constrained to sum to one and six rates for the 12 substitution types

(4)

(YANG 1994B ). The diagonal of the rate matrix is specified such that the row sums are equal to zero. We add an additional constraint by rescaling so that -q_iii = 1; this means that branch lengths of the tree are in terms of expected number of substitutions per site, v. This model is reversible because it fulfills the reversibility criterion that _iq_ij = _jq_ji for all i and j. Most commonly used models of DNA substitution are constrained to be reversible and are simply special cases of the model described here. In the analyses that follow, we use the HKY85 model of DNA substitution (HASEGAWA et al. 1984 , HASEGAWA et al. 1985 ). This model allows different base frequencies and for a transition/transversion rate bias (r_AG = r_CT; r_AC = r_AT = r_CG = r_GT; = r_AG/r_AC). The transition probabilities are calculated as P(v) = {p_ij(v)} = e^Qv.

Among-site rate variation can be accommodated in several different ways. One potential method partitions the DNA sequence into different regions (e.g., first, second, and third codon positions) and then estimates the rate separately for each partition, assuming that the rate within a partition is homogeneous [i.e., instead of a single substitution rate (m), applying to all sites, there are multiple substitution rates (m₁, m₂, ... )]. Another method assumes that the rate assigned to a site is a random variable. Typically, the gamma distribution, parameterized with the shape parameter equal to the scale parameter, is used to model rate variation across sites (YANG 1993 , YANG 1994A ). The shape parameter of the gamma distribution for among-site rate variation is denoted _R. In this article, we assume equal rates across sites or gamma-distributed rates across sites. Models that assume gamma-distribution rate variation are denoted "+ ." Our implementation of gamma-distributed rate variation uses the discrete gamma approximation with four rate categories (YANG 1994A ).

The probability of observing the states present at the ith site (x_i) is a sum over all possible assignments of nucleotides to the internal nodes of the tree. Suppose y = {y_k} for k = s + 1, ... , 2s - 1 is a generic data vector at the internal nodes. Branch k of the tree has length v_k expected substitutions per site and ancestral node (k). The transition probability from state i to state j along a branch with v expected substitutions is p_ij(v). The initial substitution rate at the root is m. Then, the conditional probability of observing the data at the ith site given the tree and rate events is

(5)

The summation is over all possible combinations of nucleotide states that can be assigned to the internal nodes of the tree. The expected number of substitutions on branch k (v_k) is found by integrating the rate over the length of the branch, v_k = ^tk_t(k) r_k(u)du, where r_k(u) is the rate along branch k at time u. The rate is a step function for the compound Poisson process considered in this article. Fig 5 shows a single branch starting at t_B = 0.40 and ending at t_E = 0.20. There are two events of substitution-rate change along this branch, occurring at times 0.30 and 0.25. The expected number of substitutions per site along this branch, then, is

(6)

View larger version (13K): In this window In a new window Download PPT slide   Figure 5. An example of how the average number of substitutions per site is calculated for a branch when substitution rate changes through time. Here, the substitution rate changes according to the compound Poisson process used in this article. There are two events of substitution rate change along this branch. The first doubles the current substitution rate and the second halves the current substitution rate.

The substitution rate at the base of the branch is m = 0.10.

Assuming independence of the substitutions across sites, the conditional probability of observing the full sequence data set given the tree and rate events is the product of the probabilities of observing the sites:

(7)

Estimating m, , and _P using Markov chain Monte Carlo:
We wish to estimate the rate of molecular evolution at each branch through Bayesian estimation of the parameters , _P, and m (where m is the rate of substitution at the base of the tree). The likelihood function for , _P, and m is

(8)

where the single integral denotes integration over all rate events and branch lengths consistent with the tree topology . Integration with respect to the probability measure is used to denote a summation for discrete random variables (such as the number of events on a tree or the topology of a tree) as well as integration for continuous random variables (such as the position of the events, the speciation times, and the gamma-distributed rates associated with events). The tree topology is considered to be fixed in this study, but the speciation times and the position and rates of events are treated as random variables. In a Bayesian analysis, we specify a joint prior probability distribution for the parameters , _P, and m along with the rate events and branch lengths.

All Bayesian inference arises from the joint posterior distribution of the parameters of interest (in this case, the posterior distribution of , _P, and m). The posterior probability density of , _P, and m is

(9)

where

and is the space for , _P, and m. We assume that , _P, and m have independent priors. The parameters _P and m have uniform priors on [0, B_P] and [0, B_m], respectively. We place an exponential prior on . The relative speciation times, t, which are scaled to be between 0 and 1 (see Fig 1), are distributed as the order statistics drawn from a uniform (0, 1) distribution, conditional on agreement with the tree topology . There is no biological meaning to the uniformly distributed priors. However, in Bayesian analysis, such priors are often used in cases where there is little or no prior knowledge about the parameters. Using uninformative priors is a way to avoid biasing the results of the analysis; the posterior probability distribution will mainly be determined by the likelihood function. An exponential prior on decreases the probability of substitution rate histories with a large number of small rate-change events.

Calculation of the posterior probability density (Equation 8) involves evaluating high-dimensional integrals and summations. We use Markov chain Monte Carlo (MCMC) integration to estimate the posterior distributions of interest. Specifically, we used the Metropolis-Hastings-Green (MHG) algorithm (GREEN 1995 ; see also GEYER 1999 ), an extension of the Metropolis-Hastings algorithm for problems in which the dimension of the sample space changes. The MHG algorithm constructs a Markov chain by first proposing a new state and then moving to that state with probability R. The steps of the algorithm are as follows: (1) the current state of the chain is ; (2)with probability density f(*|), a new state (*) is nominated; (3) the acceptance probability of the nominated state is calculated,

(10)

where f() is the target distribution (i.e., Equation 7 or Equation 8); and (4) a uniformly distributed random variable on the interval [0, 1] is generated. If this random variable is <R, then the nominated state is accepted and becomes the current state of the chain ( = *). Otherwise, the chain remains in state . Steps 1–4 are repeated many times.

Several different mechanisms were used to update the state of the chain. These mechanisms included (1) adding an event to the tree, (2) deleting an event from the tree, (3) changing the position of an event on the tree, (4) changing the gamma-distributed random variable associated with an event on the tree, (5) changing the time of an internal node on the tree, (6) changing the substitution rate (m) at the base of the tree, (7) changing the Poisson parameter (), (8) changing the gamma parameter (_P), (9) changing the transition/transversion rate ratio (), (10) changing the gamma shape parameter for among-site rate variation (_R), and (11) changing the base frequencies ().

The acceptance probabilities for the different moves take the form

(11)

(GREEN 1995 ). For move types 3–8, the standard Markov chain theory employed in the Metropolis-Hastings algorithm (METROPOLIS et al. 1953 ; HASTINGS 1970 ) applies. However, move types 1 and 2 add and delete an event, respectively. These move types involve a change in the dimensionality of the sample space. Hence, we construct reversible Markov chains for move types 1 and 2 that jump between parameter subspaces of different dimensionality using the methodology developed by GREEN 1995 .

Move type 1: Adding an event to the tree: The prior ratio for the addition of a single point (z*, r*) to the current state e = (, z, r) is

(12)

and the proposal ratio is

(13)

where d₊₁ and b are the probabilities of making a move that deletes one of + 1 events or adds an event when there are currently events, respectively. The Jacobian is 1. The acceptance probability, then, is

(14)

Move type 2: Deleting an event from the tree: The acceptance probability for the reverse step, deleting an event from the tree, has the same form as Equation 13, but with the ratio terms inverted. Detailed balance between the move types that add and delete an event on the tree is demonstrated in the Appendix

Move type 3: Changing the position of an event: Two different move types changed the position of an event on the tree. The first move type (3a) was chosen with probability _3a and picked at random one of the events on the tree. A new position was chosen uniformly on the interval [0, T]. The second move type (3b) was chosen with probability _3b and picked one of the events on the tree at random and moved its current position by a small amount randomly drawn from the interval [-_3b, +_3b]. The acceptance probability for both of these moves is

(15)

Move type 4: Changing the rate multiplier associated with an event: Two different move types changed the rate associated with an event on the tree. Both move types pick at random one of the events. This randomly chosen event will have its rate multiplier changed from r to r*. For the first move type (4a) a change to r* is proposed with probability _4a such that log_e(r*/r) is uniformly distributed on the interval [-¹/₂, +¹/₂]. The acceptance probability for this move is then

(16)

(GREEN 1995 ). The other move type (4b) is chosen with probability _4b and picks a new rate by drawing a random variable from the distribution g(r*|_P). The acceptance probability for this move type is simply the likelihood ratio (Equation 15).

Move type 5: Changing the time of an internal node on the tree: With probability ₅, the time of an internal node was changed. The times of the internal nodes were updated as follows. First, an internal node of the tree was chosen at random (excluding the root node). The time of this node was increased or decreased by adding a uniformly distributed random variable drawn from the interval [-₅, +₅]. The acceptance probability for this move is

(17)

where T* is the total tree length after the node time has been adjusted. A change to the speciation time on an internal node affects the length of a total of three branches (the ancestral branch and the two descendant branches). Often, events of rate change occur along these branches. The times of the rate-change events are maintained proportionally along the same branch they started on.

Move types 6, 7, 8, 9, and 10: Changing m, , _P, , and _R: The substitution rate (m), Poisson parameter (), gamma shape parameter (_P), transition/transversion rate ratio (), and gamma shape parameter for among-site rate variation (_R) were updated with probabilities ₆, ₇, ₈, ₉, and ₁₀, respectively, by adding to the current value a uniformly distributed random variable on the interval [-₆, +₆], [-₇, +₇], [-₈, +₈], [-₉, +₉], and [-₁₀, +₁₀], respectively. The acceptance probabilities are

(18)

for the update of m, , and _R,

(19)

for the update of , where * is the proposed state and is the current state, and

(20)

for the update of _p, where ^*_P is the proposed state and _P is the current state. We also considered gamma priors for m, , and _P. The gamma prior has parameters a and b. The exponential prior for used in the analyses of this article has a = 1. A change to m* is proposed such that log_e(m*/m) is uniformly distributed on the interval [-¹/₂, +¹/₂]. The acceptance probability for this move is then

(21)

The acceptance probabilities for changing and _p are the same as Equation 21, with m and m* replaced by and * and _p and ^*_p, respectively.

Move type 11: Changing the base frequencies: With probability ₁₁ a move was attempted that changed the equilibrium base frequencies, . The sum of the base frequencies is constrained to equal 1 and new values are proposed from a Dirichlet distribution with expected values at the current values. The Dirichlet distribution is the natural conjugate prior for a multinomial distribution and has probability density

(22)

where S is the state space (A, C, G, or T), _i is the Dirichlet parameter of the ith nucleotide, ₀ = _iSi, and _i is the frequency of the ith nucleotide. New base frequencies are drawn from the Dirichlet distribution with _i = _i0. We set ₀ = 100.0 in all of the analyses of this study.

Changing parameters near the boundary of the parameter space: Note that for move types 3b, 5, 6, 7, 8, 9, and 10 that the parameter is changed by adding a uniformly distributed random variable from an interval [-, +]. When the parameter is restricted to an interval (l, h) and the proposed value is outside this interval, we reflect the excess back into the interval. Namely, if ' = + U, where is the current parameter value (either z, t, m, , or _P) and U is the uniformly distributed change, the proposed parameter value * is

(23)

The proposal ratio [in (10)] is one for this type of proposal.

Estimating parameters using Bayesian inference: The chain was sampled every g_S generations. The posterior distributions of the parameters are obtained directly by noting the position of the chain and recording the values of m, , and _P for each sampled state. The proportion of the time the chain stays in different intervals is an approximation of the posterior distribution.

The chain was burned in by discarding the first g_B generations of the Markov chain. The burn-in was performed to allow the chain to approach stationarity before states are sampled.

Validation of computer program:
A computer program implementing the compound Poisson process for the HKY85 + (HASEGAWA et al. 1984 , HASEGAWA et al. 1985 ; YANG 1993 , YANG 1994A ) model of DNA substitution was written in C by one of us (J.P.H.; available via anonymous ftp to brahms.biology.rochester.edu or via the WWW at http://brahms.biology.rochester.edu). The advantage of using MCMC for Bayesian inference is that the sampled points are a valid (albeit dependent) sample from the posterior distribution: the Markov chain law of large numbers (theorem 3; TIERNEY 1994 ) states that posterior probabilities can be validly estimated from long-run sample frequencies. However, it is impossible to guarantee that an implementation of MCMC will converge for any particular problem (GEYER 1999 ). As GEYER 1999 (p. 80) states, "MCMC is a complex mixture of computer programming, statistical theory, and practical experience. When it works, it does things that cannot be done any other way, but it is good to remember that it is not foolproof."

The computer program was validated in several ways: (1) likelihoods were checked against existing computer programs (e.g., PAUP*; SWOFFORD 1998 ); (2) the acceptance probabilities for the various move types were worked out independently by two of us; and (3) the program was run without any data. When the program is run without any data, then the likelihood ratio equals one and the chain should target the prior distribution. When the chain is run without data, the expectation and variance of the number of events on the tree should both be equal to T. Moreover, the average rate associated with the events should equal _P/e^(P) and the variance should equal _P/e^(P). When the chain was run without data, the results satisfied these expectations. Also, multiple independent chains were started from different starting values of the parameters. The chains were examined to see if they converged to the same posterior probability distribution.

Analysis of DNA sequence data:
We applied the method to a single DNA sequence data set. The data set included complete mitochondrial sequences from 23 mammals (ARNASON et al. 1997 ). The tree topology was considered fixed in the analysis; the tree topology estimated using UPGMA (SOKAL and MICHENER 1958 ) with the Jukes-Cantor distance (JUKES and CANTOR 1969 ) was used in the analysis. All analyses were performed assuming the HKY85 model of DNA substitution (HASEGAWA et al. 1984 , HASEGAWA et al. 1985 ). This model allows different base frequencies as well as a transition/transversion rate bias. Parameters of the substitution model were estimated on the UPGMA tree using maximum likelihood (as implemented in PAUP*; SWOFFORD 1998 ).

The Markov chain was run for at least 500,000 generations. Every 50th state of the chain was sampled (i.e., g_S = 50). The sampled states were used to construct the posterior distribution for the variables m, , and _p.

	RESULTS

TOP ABSTRACT METHODS RESULTS DISCUSSION APPENDIX LITERATURE CITED

The molecular clock hypothesis of equal rates across lineages could be rejected for the mammalian mtDNA sequence data set. We used a likelihood-ratio test to examine the molecular clock hypothesis (FELSENSTEIN 1981 ). The HKY85 + model was assumed in the analysis. The null hypothesis is that the molecular clock holds; the likelihood is maximized under the constraint of equal rates across lineages (L₀). The alternative hypothesis relaxes the clock constraint by assigning a different rate to each of the lineages; the likelihood under the alternative hypothesis is L₁. The likelihood-ratio test statistic (twice the difference in the log_e likelihood, -2 log_e; = L₀/L₁) between the null and alternative models is approximately ² distributed with s - 2 degrees of freedom. The molecular clock hypothesis was rejected at the 5% level (mammalian mitochondrial sequences: log_eL₀ = -114,514.23, log_eL₁ = -114,431.21, -2 log_e = 166.03, P < 0.00001).

We assumed two different priors for for the mtDNA data. Fig 6 and Fig 7 show the results of the analysis of the mammalian mitochondrial sequences assuming an exponential prior for with mean 1.0. Fig 6 shows the change in the log_e likelihood through time. The probability of observing the data increased for the first 10,000 generations of the chain. The log_e likelihood then stabilized, fluctuating around a value of ~ -130,240. We discarded the first g_B = 50,000 generations as the burn-in time of the chain. Fig 7 shows the posterior distributions for the parameters m and _P. The estimates were m = 0.33 (0.32, 0.34) and _P = 5.70 (2.35, 10.82). The credibility interval for each variable was determined by taking the 2.5% tails of the posterior distributions.

View larger version (16K): In this window In a new window Download PPT slide   Figure 6. The change in the loge probability of observing the data (Equation 6) for the analysis of the mammalian mitochondrial sequences. The prior for was exponential with mean 1.

View larger version (11K): In this window In a new window Download PPT slide   Figure 7. The posterior probability distributions for the parameters m and P for the mammalian mitochondrial sequences. The prior for was exponential with mean 1.

Fig 8 and Fig 9 show the results of analyses in which an exponential prior with mean 10.0 was assumed for . Estimates of _P (the shape parameter of the compound Poisson process) and m changed when different parameter values of the prior were used [m = 0.33 (0.32, 0.33); _P = 63.89 (30.00, 104.24)]. The posterior distribution of m became broader and the posterior distribution of _P increased (there were more events, each with a smaller effect). However, the branch lengths were very similar regardless of the exact value of the exponential prior on . Fig 10 shows the relationship between the maximum-likelihood estimates of the branch lengths obtained assuming the HKY85 model of DNA substitution and the mean of the posterior distribution of each branch obtained under the compound Poisson process. In both cases, the compound Poisson process was able to accommodate rate variation across lineages of the tree (the slope was 1.00, with r² = 0.999 for both priors). We also examined the relationship between the compound Poisson process with exponential priors of means 1 and 10 (Fig 11). Again, there is a close relationship between the branch lengths regardless of the prior used for . The posterior distributions of branch lengths are quite robust to the two different priors used for .

View larger version (16K): In this window In a new window Download PPT slide   Figure 8. The change in the loge probability of observing the data (Equation 6) for the analysis of the mammalian mitochondrial sequences. The prior for was exponential with mean 10.

View larger version (12K): In this window In a new window Download PPT slide   Figure 9. The posterior probability distributions for the parameters m and P for the mammalian mitochondrial sequences. The prior for was exponential with mean 10.

View larger version (20K): In this window In a new window Download PPT slide   Figure 10. The Bayesian and maximum-likelihood estimates of the branch lengths for the mammalian mitochondrial sequences are highly correlated (slope = 1.0; r2 = 0.999). Bayesian estimates of branch lengths were taken as the mean of the posterior distribution. Both sets of branch lengths were obtained under the HKY85 model of DNA substitution.

View larger version (21K): In this window In a new window Download PPT slide   Figure 11. The Bayesian estimates of branch lengths for the two exponential priors examined (1 and 10). The branch lengths are highly correlated (slope = 1.00; r2 = 0.999). Bayesian estimates of branch lengths were taken as the mean of the posterior distribution. Both sets of branch lengths were obtained under the HKY85 model of DNA substitution.

The compound Poisson process for relaxing the molecular clock should prove practically useful for several applications. For example, estimation of divergence times for clades depends upon a calibration time for at least one speciation event on the tree and equal rates across lineages (though see SANDERSON 1997 ). The compound Poisson process can be used to relax the molecular clock while at the same time allow estimation of divergence times of clades. An application of the method applied to the mammalian mitochondrial data is shown in Fig 12. The general approach is the same as that taken by THORNE et al. 1998 , but we used the compound Poisson process discussed in this article instead of a lognormal distribution to relax the clock assumption. The estimation of divergence times required only one modification to the notation described in this article; instead of rescaling branch times between 0 (at the tips) and 1 (at the root), the tips of the tree were considered to be time 0 and the times of other branches on the tree were in terms of millions of years before present. The HKY85 + model of DNA substitution was assumed in the analysis; the gamma model fit the observed sequences better than a model that pooled sites according to codon position and allowed a different rate for each partition (log_eL = -114,431.21 for HKY85 + , log_eL = -117,091.83 for HKY85 + SS). An exponential prior with mean 0.006 was placed on . The estimates of the substitution model parameters were m = 0.010 (0.008, 0.012), _P = 7.58 (2.60, 17.13), = 13.67 (13.00, 14.37), and _R = 0.230 (0.228, 0.240). The posterior distribution of the divergence time of ferungulates (carnivores, perissodactyls, artiodactyls, and cetaceans) and Xenartha (Edentata) is shown in Fig 12. ARNASON et al. 1997 obtained an estimate of 86 mya for this speciation event using a calibration date of 60 million years for the cow-whale divergence. We used a calibration date of 56.5 million years for the cow-whale divergence that is based on the earliest occurrence of either artiodactyls or cetaceans. The estimate of the Ferungulate-Xenartha divergence is 92.5 mya (83.7, 99.6). The credibility interval overlaps the estimate obtained by ARNASON et al. 1997 . In Fig 12, we also indicate the divergence time estimates of the other clades on the tree, with the 95% credibility intervals.

View larger version (28K): In this window In a new window Download PPT slide   Figure 12. The estimates of divergence times for 22 mammalian species. Estimates were made under the HKY85 + model of DNA substitution. Numbers at internal nodes represent the point estimate of the divergence time (top) and the upper (middle) and lower (bottom) limit of the 95% credibility interval. An exponential prior with mean 0.006 was placed on .

Two different priors were used in the analysis of the speciation times for the mammals. We used exponential priors with means 0.06 and 0.006 per million years. These numbers were chosen so that the means of the exponential priors were similar to the priors used earlier (0.06 x 150 mya = 10; 0.006 x 150 = 1). Fig 13 shows the relationship between estimates of speciation times for two different priors on for the mammal mtDNA data. We used exponential priors with means 0.06 (x) and 0.006 (y). The x-axis shows that speciation times differed by at most 0.2 for the two priors. The estimates of speciation times appear robust to the choice of prior used for .

View larger version (14K): In this window In a new window Download PPT slide   Figure 13. The estimates of speciation times were similar regardless of the exponential prior placed on . Two exponential priors with means 0.06 and 0.006 were examined. The x-axis shows the proportion difference between the two estimates of speciation time. x and y represent the speciation times for the exponential priors with means 0.06 and 0.006, respectively.

	DISCUSSION

TOP ABSTRACT METHODS RESULTS DISCUSSION APPENDIX LITERATURE CITED

The molecular clock hypothesis, although useful in evolutionary biology, is inaccurate in two important details; rates across lineages appear to vary and the substitution process deviates significantly from a Poisson process (OHTA and KIMURA 1971 ; LANGLEY and FITCH 1973 , LANGLEY and FITCH 1974 ; also see GILLESPIE 1991 ). This article considers only lineage effects on rate variation; that is, we only consider changes in the overall substitution rate that occur along different lineages through processes such as a change in the generation time. We do not relax the assumption that the substitutions on a tree follow a Poisson process. It is important to relax the assumption that rates are constant across lineages because almost any study of the rate of substitution in DNA or amino acid sequences rejects the null hypothesis that rates across lineages are equal through time (LANGLEY and FITCH 1974 ; GOODMAN et al. 1975 ; WU and LI 1985 ). This is especially true now that biologists are using recently available software that allows easy testing of the molecular clock hypothesis (e.g., PAML, YANG 1997; PAUP*, SWOFFORD 1998).

Several approaches have been taken to relax the molecular clock. The most frequently used method for relaxing the rate-constancy assumption is to assign a fixed rate parameter to each branch of a phylogeny. This solution works well for the phylogeny problem where rates of substitution are often treated as nuisance parameters and the parameter of interest is the phylogeny of a group of organisms. However, allowing a rate of substitution for each lineage on a phylogenetic tree does not allow estimation of divergence times or allow comparison of speciation times for different taxonomic groups. SANDERSON 1997 made an important advance in the estimation of divergence times using molecular data by (1) showing how all of the sequences of a study can be used simultaneously to estimate the divergence time of a clade and (2) relaxing the molecular clock by smoothing the rate differences across speciation times on a tree. Similarly, THORNE et al. 1998 relaxed the molecular clock by assuming that rates change across speciation events. They assigned rates to descendant lineages by sampling rates from a lognormal distribution. They evaluated the posterior distribution of speciation times using MCMC.

The compound Poisson process considered here models rate variation across lineages on a tree as a step function. There were several reasons we chose to model rate variation as a step function, with a major reason being analytical tractability. However, there are biologically plausible mechanisms that could cause rates to change in a nearly step-like manner. For example, any process that rapidly changes the substitution rate (such as a change in the proofreading mechanisms, a change in selective constraints, or a rapid change in the generation time) can be approximately described using a step function. One of the main advantages of treating rate variation as a compound Poisson process is that the model can introduce rate variation across lineages at any point on the phylogenetic tree; other methods that have been considered assume that substitution rate changes at speciation events (e.g., SANDERSON 1997 ; THORNE et al. 1998 ). Such models may be sensitive to taxon sampling. On the other hand, these solutions may prove to be more practical for data analysis because there are fewer parameters to estimate. The gamma distribution used to modify ancestral rates, however, is not so easily justified using as an argument biological plausibility. We chose to use a gamma distribution to modify rates because of its analytical tractability and flexibility. Other methods for modifying rates, however, may work just as well.

We uncovered significant rate variation across lineages for the mammalian mitochondrial DNA sequence data sets. The clock hypothesis could be rejected using a likelihood-ratio test (FELSENSTEIN 1981 ). Importantly, the posterior distribution of the gamma parameter (_P) was close to 0, which also suggests that the data were not consistent with the clock hypothesis.

Other models may also prove useful for relaxing the molecular clock assumption. For example, a doubly stochastic (or Cox) process can be used to allow rates to vary through time (GILLESPIE 1991 ). In this process, the rate of substitution is drawn from some stationary distribution r (u) (the rate at time u) with the transition probabilities calculated as they were in this article. Importantly, this process has a variance to mean ratio (index of dispersion) greater than one. One of the main difficulties of this process is the choice of the rate distribution. It is not clear what distribution would be most appropriate for modeling rate variation through time. Also, if the rate distribution has too many parameters, it may be difficult to estimate the parameters of the rate distribution while simultaneously allowing estimation of other interesting parameters (such as the divergence time of a clade).

We were able to estimate the divergence times of ferungulates (and other mammalian speciation events) using Bayesian inference under a relaxed molecular clock. In this study, we assumed one calibration date; the date of the speciation event leading to cows and whales was assumed to be exactly 56.5 mya. However, the Bayesian approach makes it easy to accommodate multiple calibration dates each with error. The error in a calibration date might be assumed to be uniformly distributed on some interval. However, this would only provide a starting point for characterizing the error in the calibration dates. Ideally, information on the density of fossil taxa around the calibration date would be used to construct a more realistic probability distribution describing the divergence time of a clade.

One important result of this article is that we demonstrate that events can be placed onto a phylogenetic tree under a Poisson process and estimation performed using MCMC. The Poisson process is one of the most useful models in biology, and there are several possible extensions of the approach we describe in phylogenetics. For example, one could place events onto a tree according to a Poisson process with each event changing the synonymous/nonsynonymous rate ratio, the base frequencies, or whether the site has a positive rate or a rate of zero (i.e., the covarion model). The advantage of this approach is that it allows the rate of nucleotide substitution to change on a tree under a biologically reasonable model with the addition of a small number of parameters.

	ACKNOWLEDGMENTS

This article benefited from discussions with Rasmus Nielsen and Michael Newton. Mary Kuhner provided suggestions for checking the results of the computer program. The suggestions of two anonymous reviewers also improved this article. J.P.H. and B.L. thank the Isaac Newton Institute for Mathematical Sciences and the National Science Foundation (NSF grant DMS-9804739) for their support in Cambridge during the summer of 1998. B.L.'s research was supported by NSF grant DBI-9723799.

Manuscript received March 22, 1999; Accepted for publication December 17, 1999.

	APPENDIX

TOP ABSTRACT METHODS RESULTS DISCUSSION APPENDIX LITERATURE CITED

Here we demonstrate detailed balance between the move types that add and delete an event from the tree. Following the notation of GREEN 1995 , consider the move m "insert a point into position + 1 in the sequence" and its dual return move "delete the point from position + 1 in the sequence." We show that (from GREEN 1995 ; p. 714)

(A1)

holds for general Borel sets A and B. It is sufficient to show it holds when A and B are the products of open intervals in the spaces and + 1 events on the tree, respectively. When = 0, an open interval is just the point (0, , ).

Let a generic point in A be x = (, z, r), where z_i [0,T] and r_i (0,) for i = 1, 2, ... , . The point x' = ( +1, (z, z*), (r, r*)) is a generic point in B. Note that q_m(x, dx') = 0 if any of the arguments for the first u's or r's in x' differ from the corresponding arguments in x (because m adds a point at the end of the sequence). This means that the double integrals will measure 0 over any points in A that are not in B restricted to the smaller dimension or points in B that are not in A expanded to the larger dimension. Therefore, it suffices when specifying A and B to assume that the open intervals for the first dimensions of the z's and r's agree.

The purported acceptance probabilities for traversing the prior are

(A2)

and

(A3)

for adding and deleting an event, respectively. Note that the expressions depend on only and may be factored out of the integral expressions.

Let A = {} x _i=1(_i,c_i) x _i=1(s_i,t_i) and B = { + 1} x ⁺¹_i=1(_i,c_i) x ⁺¹_i=1(s_i,t_i). This means that A is the set where a_i < z_i < c_i and s_i < r_i < t_i for generic x A and there is a similar description for B. We have

(A4)

and

(A5)

We also have

(A6)

and

(A7)

The left-hand side of Equation 22 simplifies to

where G(·|_P) is the cumulative distribution function of g(·|_P). The right-hand side of Equation A1 simplifies to

(A8)

The left-hand side over the right-hand side is

(A9)

so detailed balance is satisfied for the jump moves.

Note that

(A10)

whether Td₊₁ > ( + 1)b or Td₊₁ < ( + 1)b. A similar calculation for inserting/deleting a point in position i in the list of events would be similarly derived, differing only in labeling of subscripts.

	LITERATURE CITED

TOP ABSTRACT METHODS RESULTS DISCUSSION APPENDIX LITERATURE CITED

ARNASON, U., A. GULLBERG, and A. JANKE, 1997  Phylogenetic analyses of mitochondrial DNA suggest a sister group relationship between Xenartha (Edentata) and Ferungulates. Mol. Biol. Evol. 14:762-768[Abstract].

COOPER, A. and D. PENNY, 1997  Mass survival of birds across the Cretaceous-Tertiary boundary: molecular evidence. Science 275:1109-1113[Abstract/Free Full Text].

DOOLITTLE, F. F., D.-F. FENG, S. TSAR, G. CHO, and E. LITTLE, 1996  Determining divergence times of the major kingdoms of living organisms with a protein clock. Science 271:470-477[Abstract].

FELSENSTEIN, J., 1981  Evolutionary trees from DNA sequences: a maximum likelihood approach. J. Mol. Evol. 17:368-376[Medline].

GEYER, C. J., 1999 Likelihood inference for spatial point processes, pp. 79–140 in Stochastic Geometry: Likelihood and Computation, edited by O. E. BARNDORFF-NIELSEN, W. S. KENDALL and M. N. M. VAN LIESHOUT. Chapman and Hall/CRC, London.

GILLESPIE, J. H., 1991 The Causes of Molecular Evolution. Oxford University Press, Oxford.

GOODMAN, M., G. W. MOORE, and G. MATSUDA, 1975  Darwinian evolution in the genealogy of hemoglobin. Nature 253:603-608[Medline].

GREEN, P. J., 1995  Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika 82:711-732[Abstract/Free Full Text].

HASEGAWA, M., T. YANO, and H. KISHINO, 1984  A new molecular clock of mitochondrial DNA and the evolution of Hominoids. Proc. Japan Acad. Ser. B 60:95-98.

HASEGAWA, M., H. KISHINO, and T. YANO, 1985  Dating the human-ape split by a molecular clock of mitochondrial DNA. J. Mol. Evol. 22:160-174[Medline].

HASTINGS, W. K., 1970  Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57:97-109[Abstract/Free Full Text].

HUELSENBECK, J. P., B. RANNALA, and Z. YANG, 1997  Statistical test of host-parasite cospeciation. Evolution 51:410-419.

JUKES, T. H., and C. R. CANTOR, 1969 Evolution of protein molecules, pp. 21–132 in Mammalian Protein Metabolism, edited by H. N. MUNRO. Academic Press, New York.

KUMAR, S. and S. B. HEDGES, 1998  A molecular timescale for vertebrate evolution. Nature 392:917-920.

LANGLEY, C. H., and W. M. FITCH, 1973 The constancy of evolution: a statistical analysis of a and b haemoglobins, cytochrome c, and fibrinopeptide A, pp. 246–262 in Genetic Structure of Populations, edited by N. E. MORTON. University of Hawaii Press, Honolulu.

LANGLEY, C. H. and W. M. FITCH, 1974  An estimation of the constancy of the rate of molecular evolution. J. Mol. Evol. 3:161-177[Medline].

METROPOLIS, N., A. W. ROSENBLUTH, M. N. ROSENBLUTH, A. H. TELLER, and E. TELLER, 1953  Equations of state calculations by fast computing machines. J. Chem. Phys. 21:1087-1091.

OHTA, T. and M. KIMURA, 1971  On the constancy of the evolutionary rate of cistrons. J. Mol. Evol. 1:18-25[Medline].

SANDERSON, M., 1997  A nonparametric approach to estimating divergence times in the absence of rate constancy. Mol. Biol. Evol. 14:1218-1231.

SOKAL, R. R. and C. D. MICHENER, 1958  A statistical method for evaluating systematic relationships. Univ. Kans. Sci. Bull. 28:1409-1438.

SWOFFORD, D. L., 1998 Phylogenetic Analysis Using Parsimony (* and other Methods). Sinauer Associates, Sunderland, MA.

THORNE, J., H. KISHINO, and I. S. PAINTER, 1998  Estimating the rate of evolution of the rate of molecular evolution. Mol. Biol. Evol. 15:1647-1657[Abstract].

TIERNEY, L., 1994  Markov chains for exploring posterior distributions. Ann. Stat. 22:1701-1762.

WRAY, G., J. S. LEVINTON, and L. SHAPIRO, 1996  Molecular evidence for deep Precambrian divergences among metazoan phyla. Science 274:568-573[Abstract/Free Full Text].

WU, C.-I. and W.-H. LI, 1985  Evidence for higher rates of nucleotide substitution in rodents than in man. Proc. Natl. Acad. Sci. USA 82:1741-1745[Abstract/Free Full Text].

YANG, Z., 1993  Maximum likelihood estimation of phylogeny from DNA sequences when substitution rates differ over sites. Mol. Biol. Evol. 10:1396-1401[Abstract].

YANG, Z., 1994a  Maximum likelihood phylogenetic estimation from DNA sequences with variable rates over sites: approximate methods. J. Mol. Evol. 39:306-314[Medline].

YANG, Z., 1994b  Estimating the pattern of nucleotide substitution. J. Mol. Evol. 39:105-111[Medline].

YANG, Z., 1997  PAML: a program package for phylogenetic analysis by maximum likelihood. CABIOS 15:555-556.

ZUCKERKANDL, E., and L. PAULING, 1962 Molecular disease, evolution, and genetic heterogeneity, pp. 189–225 in Horizons in Biochemistry, edited by M. KASHA and B. PULLMAN. Academic Press, New York.

ZUCKERKANDL, E., and L. PAULING, 1965 Evolutionary divergence and convergence in proteins, pp. 97–166 in Evolving Genes and Proteins, edited by V. BRYSON and H. J. VOGEL. Academic Press, New York.

This article has been cited by other articles:

				L. Himmelmann and D. Metzler TreeTime: an extensible C++ software package for Bayesian phylogeny reconstruction with time-calibration Bioinformatics, September 15, 2009; 25(18): 2440 - 2441. [Abstract] [Full Text] [PDF]

				R. Schultheiss, B. Van Bocxlaer, T. Wilke, and C. Albrecht Old fossils-young species: evolutionary history of an endemic gastropod assemblage in Lake Malawi Proc R Soc B, August 7, 2009; 276(1668): 2837 - 2846. [Abstract] [Full Text] [PDF]

				S. Y.W. Ho An examination of phylogenetic models of substitution rate variation among lineages Biol Lett, June 23, 2009; 5(3): 421 - 424. [Abstract] [Full Text] [PDF]

				D. Ford, F. A. Matsen, and T. Stadler A Method for Investigating Relative Timing Information on Phylogenetic Trees Syst Biol, June 12, 2009; (2009) syp018v1. [Abstract] [Full Text] [PDF]

				O. Akerborg, B. Sennblad, L. Arvestad, and J. Lagergren Simultaneous Bayesian gene tree reconstruction and reconciliation analysis PNAS, April 7, 2009; 106(14): 5714 - 5719. [Abstract] [Full Text] [PDF]

				H. Wang, M. J. Moore, P. S. Soltis, C. D. Bell, S. F. Brockington, R. Alexandre, C. C. Davis, M. Latvis, S. R. Manchester, and D. E. Soltis Rosid radiation and the rapid rise of angiosperm-dominated forests PNAS, March 10, 2009; 106(10): 3853 - 3858. [Abstract] [Full Text] [PDF]

				X. Didelot, A. Darling, and D. Falush Inferring genomic flux in bacteria Genome Res., February 1, 2009; 19(2): 306 - 317. [Abstract] [Full Text] [PDF]

				A. D. Cutter Divergence Times in Caenorhabditis and Drosophila Inferred from Direct Estimates of the Neutral Mutation Rate Mol. Biol. Evol., April 1, 2008; 25(4): 778 - 786. [Abstract] [Full Text] [PDF]

				T. A. Heath, D. J. Zwickl, J. Kim, and D. M. Hillis Taxon Sampling Affects Inferences of Macroevolutionary Processes from Phylogenetic Trees Syst Biol, February 1, 2008; 57(1): 160 - 166. [Full Text] [PDF]

				T. Lepage, D. Bryant, H. Philippe, and N. Lartillot A General Comparison of Relaxed Molecular Clock Models Mol. Biol. Evol., December 1, 2007; 24(12): 2669 - 2680. [Abstract] [Full Text] [PDF]

				T. Britton, C. L. Anderson, D. Jacquet, S. Lundqvist, and K. Bremer Estimating Divergence Times in Large Phylogenetic Trees Syst Biol, October 1, 2007; 56(5): 741 - 752. [Abstract] [Full Text] [PDF]

				N. Galtier A Model of Horizontal Gene Transfer and the Bacterial Phylogeny Problem Syst Biol, August 1, 2007; 56(4): 633 - 642. [Abstract] [Full Text] [PDF]

				B. Rannala and Z. Yang Inferring Speciation Times under an Episodic Molecular Clock Syst Biol, June 1, 2007; 56(3): 453 - 466. [Abstract] [Full Text] [PDF]

				C. L. Anderson, K. Bremer, and E. M. Friis Dating phylogenetically basal eudicots using rbcL sequences and multiple fossil reference points Am. J. Botany, October 1, 2005; 92(10): 1737 - 1748. [Abstract] [Full Text] [PDF]

				F. A. Reed, J. M. Akey, and C. F. Aquadro Fitting background-selection predictions to levels of nucleotide variation and divergence along the human autosomes Genome Res., September 1, 2005; 15(9): 1211 - 1221. [Abstract] [Full Text] [PDF]

				J. J. Welch, E. Fontanillas, and L. Bromham Molecular Dates for the "Cambrian Explosion": The Influence of Prior Assumptions Syst Biol, August 1, 2005; 54(4): 672 - 678. [Full Text] [PDF]

				S. Y. W. Ho, M. J. Phillips, A. J. Drummond, and A. Cooper Accuracy of Rate Estimation Using Relaxed-Clock Models with a Critical Focus on the Early Metazoan Radiation Mol. Biol. Evol., May 1, 2005; 22(5): 1355 - 1363. [Abstract] [Full Text] [PDF]

				M. M. Alba and J. Castresana Inverse Relationship Between Evolutionary Rate and Age of Mammalian Genes Mol. Biol. Evol., March 1, 2005; 22(3): 598 - 606. [Abstract] [Full Text] [PDF]

				C. D. Bell and M. J. Donoghue Dating the Dipsacales: comparing models, genes, and evolutionary implications Am. J. Botany, February 1, 2005; 92(2): 284 - 296. [Abstract] [Full Text] [PDF]

				J. P. Huelsenbeck and B. Rannala Frequentist Properties of Bayesian Posterior Probabilities of Phylogenetic Trees Under Simple and Complex Substitution Models Syst Biol, December 1, 2004; 53(6): 904 - 913. [Abstract] [Full Text] [PDF]

				M. J. Sanderson, J. L. Thorne, N. Wikstrom, and K. Bremer Molecular evidence on plant divergence times Am. J. Botany, October 1, 2004; 91(10): 1656 - 1665. [Abstract] [Full Text] [PDF]

				T.-K. Seo, H. Kishino, and J. L. Thorne Estimating Absolute Rates of Synonymous and Nonsynonymous Nucleotide Substitution in Order to Characterize Natural Selection and Date Species Divergences Mol. Biol. Evol., July 1, 2004; 21(7): 1201 - 1213. [Abstract] [Full Text] [PDF]

				P. G. Foster Modeling Compositional Heterogeneity Syst Biol, June 1, 2004; 53(3): 485 - 495. [Abstract] [Full Text] [PDF]

				K. Bremer, e. Friis, and b. Bremer Molecular Phylogenetic Dating of Asterid Flowering Plants Shows Early Cretaceous Diversification Syst Biol, June 1, 2004; 53(3): 496 - 505. [Abstract] [Full Text] [PDF]

				A. R. Lemmon and E. C. Moriarty The Importance of Proper Model Assumption in Bayesian Phylogenetics Syst Biol, April 1, 2004; 53(2): 265 - 277. [Abstract] [Full Text] [PDF]

				PATTERNS OF CALIBRATION AGE SENSITIVITY WITH QUARTET DATING METHODS Journal of Paleontology, January 1, 2004; 78(1): 7 - 30.

				REDUCTIO AD ABSURDUM: TESTING THE EVOLUTIONARY RELATIONSHIPS OF EDIACARAN AND PALEOZOIC PROBLEMATIC FOSSILS USING MOLECULAR DIVERGENCE DATES Journal of Paleontology, January 1, 2004; 78(1): 51 - 61.

				M. J. Donoghue and B. R. Moore Toward an Integrative Historical Biogeography Integr. Comp. Biol., April 1, 2003; 43(2): 261 - 270. [Abstract] [Full Text] [PDF]

				M. S. Springer, W. J. Murphy, E. Eizirik, and S. J. O'Brien Placental mammal diversification and the Cretaceous-Tertiary boundary PNAS, February 4, 2003; 100(3): 1056 - 1061. [Abstract] [Full Text] [PDF]

				T. G. Barraclough and A. P. Vogler Recent Diversification Rates in North American Tiger Beetles Estimated from a Dated mtDNA Phylogenetic Tree Mol. Biol. Evol., October 1, 2002; 19(10): 1706 - 1716. [Abstract] [Full Text] [PDF]

				A. J. Drummond, G. K. Nicholls, A. G. Rodrigo, and W. Solomon Estimating Mutation Parameters, Population History and Genealogy Simultaneously From Temporally Spaced Sequence Data Genetics, July 1, 2002; 161(3): 1307 - 1320. [Abstract] [Full Text] [PDF]

				J. B. Whitfield Estimating the age of the polydnavirus/braconid wasp symbiosis PNAS, May 28, 2002; 99(11): 7508 - 7513. [Abstract] [Full Text] [PDF]

				M. W. Gaunt and M. A. Miles An Insect Molecular Clock Dates the Origin of the Insects and Accords with Palaeontological and Biogeographic Landmarks Mol. Biol. Evol., May 1, 2002; 19(5): 748 - 761. [Abstract] [Full Text] [PDF]

				P. S. Soltis, D. E. Soltis, V. Savolainen, P. R. Crane, and T. G. Barraclough Rate heterogeneity among lineages of tracheophytes: Integration of molecular and fossil data and evidence for molecular living fossils PNAS, April 2, 2002; 99(7): 4430 - 4435. [Abstract] [Full Text] [PDF]

				M. J. Sanderson Estimating Absolute Rates of Molecular Evolution and Divergence Times: A Penalized Likelihood Approach Mol. Biol. Evol., January 1, 2002; 19(1): 101 - 109. [Abstract] [Full Text] [PDF]

				K. Strimmer and O. G. Pybus Exploring the Demographic History of DNA Sequences Using the Generalized Skyline Plot Mol. Biol. Evol., December 1, 2001; 18(12): 2298 - 2305. [Abstract] [Full Text] [PDF]

				M. J. Sanderson and J. A. Doyle Sources of error and confidence intervals in estimating the age of angiosperms from rbcL and 18S rDNA data Am. J. Botany, August 1, 2001; 88(8): 1499 - 1516. [Full Text]

				A. Drummond, R. Forsberg, and A. G. Rodrigo The Inference of Stepwise Changes in Substitution Rates Using Serial Sequence Samples Mol. Biol. Evol., July 1, 2001; 18(7): 1365 - 1371. [Abstract] [Full Text] [PDF]

				M. A. Suchard, R. E. Weiss, and J. S. Sinsheimer Bayesian Selection of Continuous-Time Markov Chain Evolutionary Models Mol. Biol. Evol., June 1, 2001; 18(6): 1001 - 1013. [Abstract] [Full Text]

				A. D. Yoder, R. M. Rasoloarison, S. M. Goodman, J. A. Irwin, S. Atsalis, M. J. Ravosa, and J. U. Ganzhorn Remarkable species diversity in Malagasy mouse lemurs (primates, Microcebus) PNAS, September 22, 2000; (2000) 200121897. [Abstract] [Full Text]

				P. S. Soltis, D. E. Soltis, V. Savolainen, P. R. Crane, and T. G. Barraclough Rate heterogeneity among lineages of tracheophytes: Integration of molecular and fossil data and evidence for molecular living fossils PNAS, April 2, 2002; 99(7): 4430 - 4435. [Abstract] [Full Text] [PDF]

				A. D. Yoder, R. M. Rasoloarison, S. M. Goodman, J. A. Irwin, S. Atsalis, M. J. Ravosa, and J. U. Ganzhorn Remarkable species diversity in Malagasy mouse lemurs (primates, Microcebus) PNAS, October 10, 2000; 97(21): 11325 - 11330. [Abstract] [Full Text] [PDF]