Time Dependency of Molecular Rate Estimates and Systematic Overestimation of Recent Divergence Times

Abstract

Studies of molecular evolutionary rates have yielded a wide range of rate estimates for various genes and taxa. Recent studies based on population-level and pedigree data have produced remarkably high estimates of mutation rate, which strongly contrast with substitution rates inferred in phylogenetic (species-level) studies. Using Bayesian analysis with a relaxed-clock model, we estimated rates for three groups of mitochondrial data: avian protein-coding genes, primate protein-coding genes, and primate d-loop sequences. In all three cases, we found a measurable transition between the high, short-term (<1–2 Myr) mutation rate and the low, long-term substitution rate. The relationship between the age of the calibration and the rate of change can be described by a vertically translated exponential decay curve, which may be used for correcting molecular date estimates. The phylogenetic substitution rates in mitochondria are approximately 0.5% per million years for avian protein-coding sequences and 1.5% per million years for primate protein-coding and d-loop sequences. Further analyses showed that purifying selection offers the most convincing explanation for the observed relationship between the estimated rate and the depth of the calibration. We rule out the possibility that it is a spurious result arising from sequence errors, and find it unlikely that the apparent decline in rates over time is caused by mutational saturation. Using a rate curve estimated from the d-loop data, several dates for last common ancestors were calculated: modern humans and Neandertals (354 ka; 222–705 ka), Neandertals (108 ka; 70–156 ka), and modern humans (76 ka; 47–110 ka). If the rate curve for a particular taxonomic group can be accurately estimated, it can be a useful tool for correcting divergence date estimates by taking the rate decay into account. Our results show that it is invalid to extrapolate molecular rates of change across different evolutionary timescales, which has important consequences for studies of populations, domestication, conservation genetics, and human evolution.

mutation rate, substitution rate, rate calibration, purifying selection

Introduction

Rates of mutation and substitution show considerable variation among genes, among taxa, and among different sites in DNA sequences. A number of recent studies of mitochondrial DNA extracted from subfossils and from detailed human pedigrees have reported exceptionally high estimates of mutation rate. Analyses of the mitochondrial control region, for example, have yielded mutation rates as high as 32%–260% per million years in humans (Parsons et al. 1997; Sigurdardóttir et al. 2000; Howell et al. 2003) and 95% per million years in Adélie penguins (Lambert et al. 2002). These estimates (of the mutation rate in the control region) vastly exceed the traditionally recognized substitution rate of 1% per million years for protein-coding mitochondrial DNA, which was originally derived from studies of various metazoan groups (Brown, George, and Wilson 1979) and has been supported by a number of comparable estimates from subsequent studies (e.g., Randi 1996; Fleischer, McIntosh, and Tarr 1998).

The most marked differences have been seen between rate estimates directly measured from pedigree or population studies and those inferred in phylogenetic (species-level) studies. The large rate disparities present a direct challenge to the neutral hypothesis of molecular evolution (Kimura 1983), which postulates that the overwhelming majority of mutations are selectively neutral. The substantial difference between mutation and substitution rates cannot be satisfactorily explained by the slow fixation rate of neutral mutations. Similarly, the slow purification rate of slightly deleterious mutations suggests that this rate disparity is not entirely explained by the nearly neutral hypothesis (Ohta and Kimura 1971; Ohta 1973). Although mutation rates are, by definition, higher than (or equal to) substitution rates, the two are not always distinguishable in practice. Henceforth, in the absence of a clear dichotomy, we refer to them together as the “rate of change.”

To investigate the transition between the short-term mutation rate and long-term substitution rate, we estimate rates of change from mitochondrial sequences of avian and primate taxa and compare these rates in the context of the timescales on which they were calibrated. These quantitative analyses are followed by a discussion of the underlying processes and an investigation of the factors that could be responsible for the observed variation in rate estimates. A novel method for correcting divergence date estimates is also proposed and is applied to human and Neandertal data.

Methods

Data sets and calibration points (based on paleontological or biogeographical data) were obtained from a number of previous studies that performed rate estimates based on the mitochondrial DNA of avian or primate taxa (see caption of fig. 1 for a list of the studies). Six additional data sets comprising primate taxa, one of which consisted only of Neandertal and modern human sequences, were constructed for the present study from mitochondrial sequences available in GenBank and manually aligned (these sequence alignments are available as Supplementary Material online). Three separate groups of mitochondrial data were analyzed: seven alignments of avian protein-coding genes, four alignments of primate protein-coding genes, and seven alignments of primate d-loop sequences.

FIG. 1.—

Estimated rate of change (in average nucleotide changes per site per million years) plotted against the date used to calibrate the estimate. Data points represent the mean value, and error bars denote 95% credibility intervals. (a) Protein-coding genes of avian taxa, with data (from left to right) taken from Krajewski and King (1996), Fleischer, McIntosh, and Tarr (1998), Randi (1996), Nunn and Stanley (1998), and Warren et al. (2003). (b) Protein-coding genes of primate taxa, with data (from left to right) taken from Telfer et al. (2003), new cytochrome b data set, Horai et al. (1995), and new ND4 data set. (c) D-loop sequences of primate taxa, with data (from left to right) taken from Howell et al. (2003), Kolman et al. (1995), and Ward et al. (1991), new Neandertal-human data set, new chimpanzee data set, new hominid data set, Vigilant et al. (1991), and new catarrhine data set. Note the different scale on the vertical axis of figure 1c. The curves were fitted to the data using a least-squares criterion and are described by equations (4, 5, and 6) in the text.

We used the following fossil calibration points for the six new primate data sets: human-gorilla split, 6.5–8 Myr; human-orangutan split, 12–16 Myr; anthropoid-cercopithecoid split, 24–27 Myr; and platyrrhine-catarrhine split, 35 Myr (C. Soligo, personal communication). Due to lack of adequate fossil evidence, we used estimates of the chimpanzee-bonobo split (1.4–1.8 Myr) and human-chimpanzee split (4.5– 6 Myr) derived from a conspectus of molecular and paleontological estimates (e.g., Hartwig 2002; Glazko and Nei 2003; Schrago and Russo 2003). For the data set containing hypervariable region 1 sequences from Neandertals and modern humans, the rate estimate was calibrated using the radiocarbon dates of the four Neandertal sequences in the alignment (GenBank accession numbers AY149291, AF282971, AF254446, and AF011222).

Rates were estimated using Bayesian analysis, as implemented by the program BEAST (Drummond et al. 2002; Drummond and Rambaut 2003). The molecular clock assumption was relaxed by allowing the rate to vary throughout the tree in an autocorrelated manner, with the rate in each branch being drawn from an exponential distribution whose mean was equal to the rate in its parent branch. This model was also implemented by Aris-Brosou and Yang (2002) in their program PhyBayes and was recently shown to be superior to other models of rate change. Unlike previous relaxed-clock methods for rate estimation (Thorne, Kishino, and Painter 1998; Aris-Brosou and Yang 2002; Sanderson 2002), however, BEAST does not require a user-specified tree topology. For alignments containing multiple intraspecific sequences (i.e., population-level data sets), a coalescent prior was used for the tree. For each of these coalescent priors, two population models (constant size and exponential growth) were tested, and the final rate estimates were obtained using the model that yielded the highest posterior probability. Following a burn-in of 1,000,000 cycles, rates were sampled once every 500 cycles from 10,000,000 Markov Chain Monte Carlo (MCMC) steps. Rates were estimated under an HKY + Γ + I model of sequence evolution with six gamma rate categories (Hasegawa, Kishino, and Yano 1985; Yang 1994). Convergence of the chains to the stationary distribution was checked by visual inspection of plotted posterior estimates using the program Tracer (Rambaut and Drummond 2004), and the effective sample size for each parameter sampled from the MCMC analysis was almost always found to exceed 100, usually by an order of magnitude. The only exception to this occurred during the analysis of the alignment of Vigilant et al. (1991), which contained 134 d-loop sequences; for this data set, the MCMC was run for 20,000,000 steps following a burn-in period of 2,000,000 cycles.

There have been a number of rate estimates derived from transmission data in human pedigree analyses (e.g., Howell, Kubacka, and Mackey 1996; Parsons et al. 1997; Jazin et al. 1998; Sigurdardóttir et al. 2000). These data were pooled together by Howell et al. (2003); we use the published value estimated from the pooled data to represent the output of these studies (see fig. 1).

Results and Discussion

In each of the three categories of mitochondrial sequence data, there is a distinct relationship between the estimated rate and the age of the corresponding calibration point (fig. 1). Rate estimates based on recent calibration points (<1 Myr for avian protein-coding sequences and primate d-loop sequences and <2 Myr for avian and primate protein-coding sequences) are much higher than those calibrated by older nodes. This pattern, which has been noted in a previous study of avian taxa (García-Moreno 2004), could simply be a manifestation of the difference between population- and species-level analyses, given that the average geological lifetime of a mammalian species is approximately 1–2 Myr (Raup and Stanley 1978; Martin 1986). There are large credibility intervals on the rate estimates obtained from recent calibration points, with low sequence variation in these alignments, leading to large branch length uncertainties, but most of them do not overlap with credibility intervals on rate estimates based on older calibration points. Consequently, the decrease in the estimated rate of change with increasing calibration depth is unambiguous. The rapid decline in rates suggests that the strength of selection is considerable, which is inconsistent with both the neutral and nearly neutral theories of molecular evolution.

The rate estimates from nodes beyond about 2 Myr are broadly uniform, with rates hovering around 0.01 substitutions per site per million years corresponding to the traditional sequence divergence rate of 2% per million years. This residual rate, which appears to stay constant over long periods of time, is probably due to two processes that are not mutually exclusive: (1) a neutral or near-neutral mutation rate resulting from neutral and/or slightly deleterious changes that become fixed by genetic drift and (2) adaptive evolution. It is difficult, if not impossible, to separate or quantify the contributions of each of these two processes to the long-term substitution rate.

There are several factors that could be responsible for the observed patterns in rate estimates, all of which are consistent with a measurable decline from the instantaneous mutation rate to the long-term substitution rate. We examine and evaluate these below.

Purifying Selection

Many of the sequence polymorphisms that are seen among individuals of a population and in intraspecific comparisons are removed over evolutionary time due to the action of purifying selection or by random genetic drift. This is consistent with observations that the ratio of nonsynonymous to synonymous mutations, dN:dS, is (1) higher within species than among species (e.g., Hasegawa, Cao, and Yang 1998) and (2) higher in terminal branches than internal branches (e.g., Nielsen and Weinreich 1999; Moilanen and Majamaa 2003). In order to verify this, we examined a data set of 1,137 aligned sites from the cytochrome b sequences of 5 chimpanzees, 2 bonobos, 3 gorillas, and 22 humans (this alignment was a subset of one of those used in the rate analyses above). Using the program codeml (Yang 1997), we calculated dN:dS for terminal branches (intraspecific comparisons) and internal branches (interspecific comparisons) under the M0 codon model specified by Yang et al. (2000) based on the tree (((Pan paniscus, Pan troglodytes), Homo sapiens), Gorilla gorilla). Allowing a different dN:dS model for each species was a significant improvement on a model with only a single dN:dS ratio (P = 0.018, assuming a

\(\mathrm{{\chi}}_{1}^{2}\)

-distributed log-likelihood ratio). The intraspecific and interspecific dN:dS ratios were 0.1509 and 0.0348, respectively. We suggest that the difference between the two values is indicative of purifying selection, which acts on the majority of nonsynonymous changes (Fay, Wyckoff, and Wu 2001; Holmes 2003).

Mutational Hot Spots and Saturation

It is known that rates of change can vary substantially among sites and that certain sites appear to change at remarkably high rates. Multiple changes can occur at these mutational hot spots, which leads to an underestimation of the mutation or substitution rate if the method of analysis does not account for the changes hidden by superimposed mutations. For example, it is known that the control region has a number of hot spots at which mutations are likely to be superimposed over long timescales (Sigurdardóttir et al. 2000). Allowing for rate variation among sites can alleviate or remove the impact of mutational saturation because this permits a proportion of the sites to change at a high rate and corrects for multiple replacements accordingly. However, while the more sophisticated substitution models are able to take saturation into account to some degree, their accuracy in modeling the actual saturation process is not well understood (Ho and Jermiin 2004).

To gain some insight into the location of sequence differences over short and long timescales, we surveyed the distributions of changes seen in intra- and interspecific comparisons for (1) cytochrome b data used in the dN:dS analysis above and (2) hominid d-loop data used for one of the rate estimations in figure 1. There were no clear patterns in the cytochrome b comparisons (fig. 2a), with intraspecific differences distributed quite uniformly throughout the length of the sequence; a test for Poisson-distributed uniformity was conducted by dividing the sequence into five sections, and polymorphic sites were uniformly distributed (P = 0.05). In the d-loop comparisons, however, there were several regions and isolated sites at which intraspecific polymorphisms appeared in two or more of the species (fig. 2b); a test for uniformity was not conducted because the number of polymorphisms was too low. Variation at possible hot spots only appears to form a small proportion of the total variation among d-loop sequences; consequently, failure to adequately correct for saturation at such sites is probably insufficient to explain the relationship between the estimated rates and the calibration times. Furthermore, many of the data sets analyzed in the present study exhibited relatively low sequence variation, suggesting that mutational saturation is unlikely to be at a level sufficiently high to confound the estimation of molecular rates.

FIG. 2.—

Distribution of polymorphic sites along sequences of cytochrome b and d-loop in primates. Shorter and taller bars indicate sites at which two and three different nucleotides are present, respectively. The cytochrome b results are based on 1,137 sites from 25 humans, 5 chimpanzees, 2 bonobos, and 3 gorillas; the d-loop results are based on 214 sites from 35 humans, 6 chimpanzees, 20 bonobos, and 9 gorillas.

Sequence Errors

Rates estimated from population-level data are very sensitive to the number of polymorphisms in the sequence alignment. When sequence variation is low, the presence of even a small number of sequence errors can lead to an overestimation of the mutation rate. Error rates will depend on the laboratory of origin, the gene of interest, and the fidelity of the polymerase being used, among other factors, so it is difficult to accurately estimate the frequency of errors in published sequences. Reported estimates of GenBank sequence errors have been as high as 1 error per 500–1,000 bp (Clark and Whittam 1992), but lower error rates of 1 per 100,000 bp and 1 per 10,000 bp have been claimed for the published genomes of bacteria (Read et al. 2002) and humans (International Human Genome Sequencing Consortium 2003), respectively, due to sequencing redundancy.

If we assume that the tip of each terminal branch in a tree has ε extra “mutations” per site due to sequencing errors, then the effect of sequence errors on rate estimates can be described by:

\[\mathrm{Rate}{\,}(t){=}\frac{\mathrm{{\mu}}t{+}\mathrm{{\varepsilon}}}{t},\]

(1)

where μ is the mutation rate per site and t is time. Even very high rates of sequence error (e.g., 1 × 10⁻³ errors per base pair) introduce only a small bias in rate estimates based on short timescales (fig. 3). The resultant curves assume that the errors are randomly distributed within sequences and uniformly among sequences, which is unlikely to be the case (Clark and Whittam 1992). However, the rate estimates based on sequence errors are in some instances orders of magnitude smaller than those observed in real data, allowing us to rule out the possibility that the rate variation is solely the spurious product of sequence errors. Nevertheless, PCR amplification and sequencing errors may still present problems for population genetic analyses, where multiple PCR and sequencing reactions per sample are not routine. It is possible that sequence errors may also have influenced the dN:dS estimates presented above because they tend to inflate estimates of dN:dS (Nielsen and Weinreich 1999). Until further studies provide an accurate quantification of sequence error rates (e.g., Schmutz et al. 2004), this issue must be borne in mind when interpreting the results of rate studies, particularly those based on intraspecific data.

FIG. 3.—

Estimated rates when errors are present in the sequences being analyzed. The horizontal line indicates the observed rate when there are no errors in the sequences being analyzed, assuming a base rate of 0.01 changes per site per million years. Predicted curves are given for three different error rates.

Errors in Calibration Points

The majority of calibration points obtained from the fossil record are only able to place lower bounds on divergence dates because they represent the earliest appearance of members postdating the divergence event (Hedges and Kumar 2004). The actual date of the split can substantially predate the oldest known fossil evidence of the descendants, especially in taxonomic groups with poor preservation potential, such as soft-bodied invertebrates. The poor resolving power of the morphological features of fossils can also result in taxonomic ambiguity, which precludes usage of such fossils as reliable calibration points. On shorter timescales, lineage sorting can result in coalescence events predating actual speciation or population subdivision events, which also leads to underestimation of actual divergence times if calibration points are based on biogeography or the fossil record.

The effects of using calibration dates that are underestimated by an amount ε (measured in million years) can be seen by estimating rates using the following equation:

\[\mathrm{Rate}{\,}(t){=}\frac{\mathrm{{\mu}}(t{+}\mathrm{{\varepsilon}})}{t},\]

(2)

where μ is the rate of change per site per million years and t is the estimated age of the calibration point. With a rate of change of 0.01, rates were substantially overestimated when the true age of the calibration point was 1 Myr older than the date actually used for calibrating the rate estimate, but the estimation bias is minimal for calibration points older than about 1.5 Myr (fig. 4). The magnitude of the estimation bias is large enough to explain the variation in rates estimated from primate d-loop sequences (fig. 1c) but not in protein-coding genes (fig. 1a and b); however, inaccuracies in calibration points could still be contributing to the apparent decline in rates over time.

FIG. 4.—

Estimated rates when dates used for calibration are inaccurate. The horizontal line indicates the observed rate when the dates of calibration points are accurate, assuming a base rate of 0.01 changes per site per million years. The other curves are based on calibration points that are understated by 0.1 (lowest curve), 0.5, and 1.0 Myr (highest curve).

Correction for Calibration Times

The factors responsible for the observed patterns in rate estimates are not entirely clear, and it may be difficult to quantify the contribution of each factor to the overall trend in rate estimates. However, uncertainty about the underlying mechanisms does not preclude measurement of the decay in rates with growing calibration age. The relatively predictable manner in which rates decline with increasing calibration time may allow a correction for estimates of rates and divergence dates. For this purpose, we used a least-squares criterion to fit a vertically translated exponential curve to the data:

\[\mathrm{Rate}{\,}(t){=}\mathrm{{\mu}}e^{{-}\mathrm{{\lambda}}t}{+}k,\]

(3)

where μ is the instantaneous mutation rate (minus the rate of lethal mutations) and λ is inversely proportional to the half-life of the rate decay. The constant term, k, is a finite asymptote (or plateau) that vertically translates the curve and represents the evolutionary rate over long time periods. The estimated curves for the avian mitochondrial sequences, primate protein-coding sequences, and primate d-loop sequences were, respectively (fig. 1):

\[\mathrm{Rate}_{\mathrm{Aves,protein-coding}}{\,}(t){=}0.0400e^{{-}0.445t}{+}0.0054\]

(4)

\[\mathrm{Rate}_{\mathrm{Primates,protein-coding}}{\,}(t){=}0.5204e^{{-}2.042t}{+}0.0144\]

(5)

\[\mathrm{Rate}_{\mathrm{Primates,d-loop}}{\,}(t){=}0.4535e^{{-}6.408t}{+}0.0148.\]

(6)

Uncertainty in these expressions can arise due to (1) the wide confidence intervals for the rate estimates based on recent calibration points and (2) the limited number of data points used to estimate the curve, especially in its critical portions. The terms μ and λ are strongly influenced by the points near the y axis (i.e., the rate estimates calibrated on short timescales) and should be regarded with caution. However, the plateau values (k), which are far more robust, suggest that the long-term protein-coding substitution rate is about 0.5%–1.5% per million years.

Provided that we know the function for the rate curve, we can estimate the divergence date of two sequences (t_d) given their sequence divergence (d) by solving the following equation for t_d, which follows from equation (3):

\[d{=}{{\int}_{0}^{t_{d}}}\mathrm{{\mu}}e^{{-}\mathrm{{\lambda}}t}{+}k{\,}dt{=}\frac{1}{\mathrm{{\lambda}}}({-}\mathrm{{\mu}}e^{{-}\mathrm{{\lambda}}t_{d}}{+}k\mathrm{{\lambda}}t_{d}{+}\mathrm{{\mu}}).\]

(7)

This has a unique solution if k ≥ 0, which should always be the case. While this approach requires that the rate decay curve is accurately known and assumes that there is no saturation, it is more defensible than applying a single rate of change across population- and species-level timescales.

Using an alignment of d-loop sequences from four Neandertals and four humans (African, Caucasian, Chinese, and Indian) (alignment available as Supplementary Material online), we recalculated the ages of the last common ancestors of Neandertals, humans, and of Neandertals and humans. The new date estimates were computed by solving equation (7) numerically with the parameter values estimated from the primate d-loop data (eq. 6) and with HKY85-corrected distances estimated from the alignment using the program baseml (Yang 1997). Ages of last common ancestors were as follows: Neandertals and humans 354 ka (222–705 ka), Neandertals 108 ka (70–156 ka), and humans 76 ka (47–110 ka). These date estimates are intermediate between the values obtained when either the lowest or highest rates of change from figure 1c are used for divergence date estimation (table 1).

Table 1

Divergence Date Estimates for Humans and Neanderthals Using Three Different Methods

Evolutionary Event	Date (ka) Estimate Using Lowest Rate (0.2% per million years)a	Date (ka) Estimate Using Highest Rate (48% per million years)a	Date (ka) Estimate Using Rate Curvea
Human-Neandertal split	2,606 (2,160–3,052)	145 (120–169)	354 (222–705)
Common ancestor of Neandertals	1,399 (1,012–1,786)	78 (56–99)	108 (70–156)
Common ancestor of humans	1,074 (725–1,423)	60 (40–79)	76 (47–110)

a

Upper and lower bounds on the error range of each date estimate correspond to the date estimates made from sequence divergence plus and minus one standard deviation, respectively.

The three new date estimates were considerably younger than those estimated in previous studies, which gave ranges of 365–853 ka (Ovchinnikov et al. 2000), 550–690 ka (Krings et al. 1997), and 317–741 ka (Krings et al. 1999) for the Neandertal-human divergence; 151–352 ka (Ovchinnikov et al. 2000) for the last common ancestor of Neandertals; and 106–246 ka (Ovchinnikov et al. 2000), 120–150 ka (Krings et al. 1997), and 111–260 ka (Krings et al. 1999) for the last common ancestor of humans (fig. 5). This discrepancy arises because high, short-term rates of change were taken into account by our approach. The standard error ranges on the new estimates do not take into account uncertainty associated with equation (7). Furthermore, the revised date estimates are sensitive to μ and λ, which are difficult to estimate accurately unless abundant data are available. However, our results clearly show the effect of correcting date estimates of geologically recent divergence events using a rate decay curve estimated from sequence data and demonstrate the unsuitability of using a single short-term or long-term rate to date evolutionary events on an intermediate timescale.

FIG. 5.—

Trees showing ages of various last common ancestors for Neandertals and humans. (a) Dates obtained from previous studies (from left to right: Ovchinnikov et al. [2000], Krings et al. [1997], and Krings et al. [1999]) estimated by applying a single, unchanging substitution rate. (b) Dates corrected using equation (7), using values taken from the rate curve described by equation (6), which accounts for the transition from a high, short-term rate to a low, long-term rate.

Conclusions

Our study has shown that the relationship between mutation and substitution rates can be described by an exponential curve. The curve can be estimated from rates obtained from sequence data and offers a new method for correcting divergence date estimates. The immediate implications of our results are that molecular rates should be interpreted in the context of calibration point age and that short-term mutation rates can only be extrapolated to older times (or vice versa) after accounting for the relationship between short-term and long-term rates of change. Taking rate variation into account is particularly important for analyses of sequences on timescales of less than about 1–2 Myr before the present, such as studies on populations, domestication, and conservation genetics, which often incorrectly apply phylogenetic substitution rates to population-level analyses. Consequently, the results of previous studies need to be reevaluated.

1

Present address: School of Environmental Sciences, University of Adelaide, Adelaide, Australia.

Mark Springer, Associate Editor

The authors thank C. Soligo for advice regarding primate fossil calibration points; D. Rubenstein for reading a draft of the manuscript; E. C. Holmes, D. Penny, E. Willerslev, A. Rambaut, R. A. Fortey, and M. Bunce for advice; M. Spencer for pointing out an error in an earlier version of the paper; and M. S. Springer and two anonymous reviewers for helpful comments that led to improvements in the paper. S.Y.W.H. was supported by a Commonwealth (Oxford) Scholarship and a Domus Research Studentship and Edward Penley Abraham Cephalosporin Scholarship from Linacre College, Oxford. A.C. was supported by the Biotechnology and Biological Sciences Research Council and by the Wellcome and Leverhulme Trusts. A.J.D. was supported by the Wellcome Trust.

Author notes

*Henry Wellcome Ancient Biomolecules Centre, Department of Zoology, University of Oxford, Oxford, United Kingdom; and †Evolutionary Biology Group, Department of Zoology, University of Oxford, Oxford, United Kingdom