Introduction
One crucial aspect of molecular evolution is the exploration of variation in the rate of substitution between different lineages or under different evolutionary pressures. Investigating this variation can lead to interesting insights into the adaptation or neutral evolution of a gene.
1–3 Many papers have presented statistical methods that distinguish different types of rate variation, such as variation in substitution rate of a single site among different lineages (heterotachy),
4–6 or variation between sites.
7,8
Molecular rate variation with a consistent pattern can be indicative of underlying biological processes.
3 Using different phylogenetic models of molecular evolution, some life history traits have been shown to correlate with substitution rate including body mass, time to sexual maturity, and lifespan.
9–12 Other traits that have shown to be correlated with substitution rate include sperm competition
13–16 and sexual dimorphism, either in mass
17 or coloration.
18 A set of abiotic variables that may be correlated with rate are environmental conditions such as salinity, temperature, latitude, and elevation.
19–21
The mitochondrial genome (mtDNA) is known to exhibit substantial rate variation,
22 and two prominent hypotheses to explain this are the longevity hypothesis and generation time hypothesis. The longevity hypothesis proposes that there is selection for greater fidelity in the replication mechanism of the mtDNA between an organism's conception and production of its gametes, and thus predicts that shorter-lived species will have a faster rate of neutral mtDNA substitution.
23 The longevity hypothesis is also consistent with a prediction for selective advantage for decreased reactive oxygen species (ROS) production and/or more resistant macromolecular components.
24
The generation time hypothesis predicts a faster substitution rate for both mtDNA and nuclear DNA in species with shorter time to sexual maturity and thus shorter generation length because it assumes a constant error rate and more replications per geological time.
25 Generation time was proposed as an exception to the assumption of the ‘molecular clock’ hypothesis
26 of a single rate of molecular evolution between species.
27 It originated from the results of DNA hybridization experiments between rodents and artiodactyls, and rodents and primates which found a higher rate of nucleotide substitution measured in time, but not generations.
28,29 This was later supported in phylogenetic studies where molecular evolution was faster in rodents than primates, consistent with their generation time differences.
30,31
The pattern and causes of mtDNA rate variation in primates are of particular interest, and have been the subject of longstanding debate.
32 Much early work was concerned with correcting for primate mtDNA rate variation to estimate divergence times.
33,34 The generation time hypothesis is of particular interest for primates as it has been invoked to explain the difference in neutral rate in both nuclear and mtDNA genomes between Old World monkeys and apes,
31,33,35 but conflicts with genome-wide estimates that include strepsirhines and New World monkeys.
36 In addition, lemurs are known to have a low rate of transitions but a similar rate of transversions to apes and monkeys.
37
One set of studies compared the three major hypotheses for the third codon position of cytochrome b (
CYTB) from the mtDNA and a few nuclear loci across mammals and birds
10–12 and was later extended to include an amino acid analysis of the protein-coding regions of the mtDNA for the longevity hypothesis.
38 Welch et al
9 examined the nine longest protein-coding genes of the mtDNA (at the codon level) in mammals and found similar negative correlations with the synonymous rate and body mass, lifespan, fecundity, and time to sexual maturity. In both approaches, primates have been identified as an outlier group for the longevity and body mass phenotypes, where they either reduce the association with lifespan
9,38 or cause the association to be positive with body mass.
9 Recently, a Bayesian approach which simultaneously estimates phenotype state, divergence dates, and evolutionary rate also found a negative correlation between the evolution of
CYTB, and longevity and body mass, but no association with generation time.
39 Considering the results described above, the generation time hypothesis was rejected because of lack of a significant association with nuclear loci, whereas the longevity hypothesis was favored over the body mass association because of a stronger genotype-phenotype correlation. Non-protein-coding portions of mtDNA are predicted to show similar signals to the protein-coding regions, but have not been examined.
For these analyses, the main methods used are phylogenetically corrected correlations between the substitution rate and the relevant variable value (eg, life history, external condition, or phenotype) of the extant taxa. One difficulty with this type of approach is that the two values being compared are not properties of the same part of the tree: the substitution rate is reconstructed on the branches, while information about the phenotype is generally located only at the tips of the tree or reconstructed on internal nodes. This type of pairing may be statistically problematic.
In both cases, these methods require a broad comparison between taxa (ie, mammals) and/or concatenation of genes (or gene groups) in order to to achieve sufficient power.
10,40 Sufficient evolutionary time is a requirement in forming comparisons in order to correctly estimate the correlations as the full phylogeny is not used. Thus, these papers did not examine the primate clade on its own because the methods are not well-suited to scaling down to a relatively small number of taxa. Some studies have analyzed body mass and generation time in primates, mostly anthropoids,
31,33,34,41 but so far primates have not been specifically tested for the longevity hypothesis.
Recently, a model was implemented that can overcome many of these problems by directly estimating the association between phenotype and substitution rate in a relatively small number of taxa.
16 It achieves this by comparing extant discrete phenotypes to extant sequences across a phylogeny and, through the maximum likelihood framework, both genotype and phenotype substitutions are allowed on the branches. Thus, the phenotype can evolve across the tree. The rate of genotypic substitution while under a given phenotype can then be estimated throughout the whole tree as opposed to just the terminal branches giving increased power with limited numbers in a taxonomic group. The method uses a mixture model approach that allows for different sites (potentially codon positions) to associate with either a foreground-associated model or a background-unassociated model.
16,42
In this paper, we extend the method and examine its robustness through extensive simulations. We then use the method to explore the longevity and generation time hypotheses as well as body mass in primates. We found that the longevity hypothesis is supported for the third codon position of protein-coding genes, consistent with mammal wide estimates,
9–12 and find a similar result for the non-coding control region. We find no consistent association between the phenotypes tested and nuclear genes, including those that encode mitochondrial proteins. In light of our results, we discuss the effect of the three important phenotypes on rate variation applied to primate evolution.
Methods
Theory
The hybrid substitution models and likelihood ratio test (LRT) first introduced in O'Connor and Mundy
16 are used in the present study with some modification. This pair of models uses the coevolutionary model of Pagel
43 as inspiration to combine other substitution models for both genotype (for example General Time Reversible (GTR))
7,44,45 and discrete phenotype models (such as Mk).
46 This allows us to detect differences in genotypic evolutionary rates associated with a phenotype.
The background evolutionary rate is estimated by the independent matrix which assumes a single rate of genotypic substitution irrespective of the phenotypic state. The instantaneous independent
Q matrix is given as:
For every site in the alignment the Q matrix defines the instantaneous rate of change. In this case there are two inputs and outputs, pi is the input phenotype and pj is the output phenotype, with similar notation for the genotype state of gi changing to gj. Here, we test with nucleotide data so gi and gj can take any of four character states (A, C, G, T). The phenotypic characters can take on the characters 0, 1, … Np − 1 depending on the number of phenotypes tested (Np).
For single changes the rate is calculated based on the underlying rate matrix for the type of change; genotypes use the GTR model and phenotypes use the Mk. For example, if the phenotypic character is the one to change, then the rate of change is defined by phenotypic rate matrix (
Qp). Similarly, for a genotypic change
Qg, the genotypic rate matrix, is used. Double substitutions are fixed to zero following the arguments of Pagel.
43
To construct the null model for our likelihood ratio statistic (LRS), we expand from the basic structure of the independent matrix to include rate heterogeneity among sites. Previous studies have shown that real data can be highly susceptible to this type of variation.
8 We use a mixture model approach
42 with two different rate matrices, otherwise the time and rate are mathematically confounded. Using the independent matrix to detect the background signal creates a contrast for a different matrix to find the associated or heterogeneous signal.
16,42,47
In our null model, we combine an independent rate matrix and an independent Q matrix with a scaling weight parameter (W) that multiplies each of the genotypic rate parameters. All of the branch lengths and rate parameters are equal in the two Q matrices, except for the addition of the scaling parameter in the scaled Q matrix. Also, we add a parameter to estimate the proportion of sites (Pr) coming from the independent Q matrix, with 1 - Pr coming from the scaled Q matrix. Thus, only two additional parameters (the proportion of sites and the scaling weight parameter) are added to allow for a scaled rate heterogeneity of sites.
The dependent matrix is designed similarly to the scaled matrix but takes on scaling parameters corresponding to different phenotypes (
Wpi). The
Q matrix for the dependent matrix is defined as:
When Wpi is set to one for all phenotypic states, it reduces to the independent model. As described before, the null model is a mixture model between the independent matrix and the scaled matrix. The alternative model is constructed as a mixture model with a proportion of sites from the independent matrix (the background unassociated signal) and the dependent matrix (ie, replaces the scaled Q matrix in the previous formulation). The dependent model is a likelihood approximation as the true likelihood would be calculated conditional on all phenotypic intermediate substitutions. This mathematical simplification does not seem to bias the test for any scenario for which we are aware (see additional file 1).
As a result of the approximation, our LRT is constructed by calculating the approximate LRS (ALRS) (two times the difference in log likelihood from the alternative and null models), which is distributed as a Ç
2 distribution with the degrees of freedom being the difference in the number of parameters in the hierarchically related models.
48–50 In the general case, the degrees of freedom are equal to the number of phenotypes minus one; in the binary phenotype case, the degrees of freedom equals one. The simplification of the dependent matrix does not seem to adversely affect the asymptotic properties of our test statistic (the ALRS), even with heterogeneity among sites (see Fig. S1 of additional file 1).
We constrain all weight parameters (W in the scale matrix and all Wpi in the dependent matrix) to be greater than or equal to one.
Simulations
To test the models we developed a simulation program as used in O'Connor and Mundy
16 with one minor modification. Previously, the program normalized the input parameters, which caused a number of difficulties. First, the parameters explicitly set were not the ones used for simulation and thus not the values tested. Second, when calculating the rate matrix values for simulation, a normalized matrix with
W0 of 1 was not comparable to the probabilities for the independent matrix. Thus, the background rate simulated in part of the data (
W = 1) was not the same in a different portion under the alternative model even though both values of input were the same (
W0 = 1). Therefore, we modified the java library PAL
51 to prevent it from normalizing the rate matrix before calculating the probabilities used in simulation, since the parameters are chosen for particular purposes.
We ran two main sets of simulations to reevaluate the performance of the method.
16 First, we used a binary phenotype with 16 taxa, a strictly bifurcating tree, and concatenated 500 sites for both the independent and dependent matrices for a total of 1000 sites. We set the scale value for phenotype 0 (
W0) to one and varied the scale value for phenotype 1 (
W1) from 1 (false positive (FP)) to 10 (extreme true-positive (TP)) at every integer. We increased the sum of all branches (tree length) from 1 to 10 at every integer, creating 100 different scenarios with 50 tests each. The alternative model used in the ALRS was constrained with
W0 <
W1. Second, we tested this specification of the alternative model by simulating data sets with greater numbers of phenotypes. This served the dual purpose of evaluating the use of the ALRS beyond the binary case and testing its performance with more constrained alternative models.
We simulated data with 16 or 32 taxa, a tree length of three, and number of phenotypes ranging from two to four. Further, we tested the difference between a stricter hypothesis (W0 < ··· < W3) and a relaxed hypothesis (W0 ≠ … ≠ W3) under both the null model (weights all equal to 1) and an alternative model (Wx = 2x + 1, for x from 0 to Np − 1 with Np set as two, three, or four). This created 24 scenarios (two numbers of taxa by two hypotheses by three different numbers of phenotypes by a null and alternative model) each tested with 200 simulated data sets.
Primate genotype and phenotype data
We applied the method just described to the question of association between the rate of mtDNA evolution in primates and key life history variables. We collected two forms of data, phenotypic and genotypic. We used the primate phylogenies of Goodman et al
52 (Old and New World monkeys and apes) and Horvath et al
53 (lemurs and lorisoids).
We collected data on maximum lifespan, time to female sexual maturity, and female body mass to test hypotheses of longevity and generation time and general life history associations, respectively.
10,23 These continuous trait values were discretized because our method is based on discrete characters. We split the values into three equal groups based on the distribution of all readily available values for primates and then mapped them back to the taxa we used (
Table 1). We used three categories so that they would roughly correspond to lower than average, average, and higher than average in the distribution to account for the two tails of the distributions. Also, as will be shown from our simulation results, this is the greatest number of categories with appreciable power. The data were collected from Smith and Junders
54 and Isaac et al
55 for female body mass and An Age for time to female sexual maturity and maximum lifespan for most of the species.
56,57 Less data were available for time to female sexual maturity and thus we used more genus level assignments. The phenotype data and assignments used are reported in additional file 2.
We retrieved 26 primate mitochondrial genomes from Genbank (see additional file 3) and aligned these sequences using MultiZ genome aligner
58,59 with a human sequence as reference and extracted individual coding genes based on their annotation in the human sequence (accession number: NC_001807). Again, based on the human reference, we extracted all the tRNAs and concatenated them into a single alignment. We also retrieved the control region, which is known for its high level of divergence.
60,61 We then realigned the control region using the alignment program PRANK with the parameter -F, following the recommendation of the authors.
62,63 PRANK was designed to align rapidly divergent and gap-inclined alignments and has been shown to perform well with simulations of insertion-deletion models.
64 All alignments were manually evaluated for consistency.
In addition to the 17 mtDNA alignments, we collected 24 alignments of nuclear protein-coding genes where data were available for at least 8 taxa. This data included 6 nuclear-encoded mitochondrial complex IV genes which interact with three mtDNA genes (
COX1, COX2, COX3) and should be functionally coevolving.
65,66 Accession numbers can be found in additional file 3.
For all three hypotheses (including general life history) we constrained the alternative model to consider only W0 > W1 > W2 as described previously in the simulation section as this formulation is consistent with the directionality predictions from the hypotheses.
Nuclear genes
In contrast to the mtDNA, the nuclear genome does not have a consistent pattern of association with the phenotypes (
Table 4). 61.5% of the mtDNA protein-coding genes were significantly associated with body mass whereas only 12.5% of the nuclear genes tested were significant. Similarly, maximum lifespan is significantly associated with 69.2% of the mtDNA protein-coding genes and 4.2% of the nuclear genes. Time to sexual maturity shows a low but similar quantity of signal in both the mtDNA and nuclear genes, with 15.3% and 8.3% respectively.
Two specific molecular functions predominate the nuclear genes we tested: nuclear encoded genes that function in the electron transport chain of the mitochondria and male associated genes.
In the first set, we tested six of the mitochondrial complex IV genes and found no significant associations. Thus, the associations found with mitochondrially-encoded proteins do not extend to mitochondrial proteins encoded in the nuclear genome. Even after all six loci are concatenated, they are not significant for any phenotype tested (uncorrected P-values; body mass = 0.980; maximum lifespan = 0.850; and time to female maturity = 0.743).
The second main functional group of nuclear genes, male associated genes, show the most significant associations with phenotypes for the nuclear genome. These are genes either known to be expressed in the testis, found in the ejaculate, or involved in sperm function. In this case, however, in contrast to the mtDNA genes, the strength of signal from the third codon position is not generally greater than that from the first two codon positions (
Table 5). This indicates that some of the associations in male genes may be due to positive selection, which is consistent with the frequent occurrence of adaptive evolution at these loci.
14,15,77,78 However, the cause of the specific associations seen with different phenotypes is unclear. Several of these loci are related to sperm competition, and there is a possibility of an association arising indirectly from a correlation between the traits tested and the degree of male competition, but this requires further analysis.
Assessing the three phenotypes and their association with mtDNA rate variation in primates
Two of the three phenotypes we tested, body mass and maximum lifespan, are associated with mtDNA rate variation in primates. It is difficult to differentiate between the two significant phenotypes (corresponding to general life history and the longevity hypotheses) because the two traits are highly correlated.
25 In our results, the proportion of mtDNA loci supporting the longevity hypothesis (10/17) was very similar to that for the association with body mass (8/17), providing no discrimination between the two. When examining the distribution of ALRSs for each hypothesis, we are able to differentiate between the two. Using a two-tailed paired Wilcoxon Rank Test for the 17 regions of the mtDNA, we find that maximum lifespan is significantly different from body mass (
P-value = 6.63 × 10
−3; median maximum lifespan 38.518; median body mass 6.657), thus providing more support for a longevity association than of general life history for primates. We now discuss how our results relate to each phenotype in more detail.
Generation time hypothesis
The generation time hypothesis predicts that if a species has shorter intervals between generations then the mutation rate per year will increases with decreased generation time and there will be be greater number of substitutions fixed because of the increased number of genomic replications.
35,79 Some studies have found evidence in support of this hypothesis across a variety of taxonomic groups.
25,35,79–82
The generation time hypothesis predicts that both nuclear and mitochondrial substitution rates should be affected. In contrast, when examining
CYTB in mammals, Nabholz et al
10 found a correlation with time to female maturity (though not as the main predictor), but not in the the nuclear gene interphotoreceptor retinoid-binding protein exon 1.
Our results provide little support for the generation time hypothesis for primates in either the mtDNA or nuclear genomes. Initially, this may seem surprising as previous studies examining primate rate variation have invoked the generation time hypothesis as an explanation for the rate variation between Old World monkeys and apes (also known as the hominoid rate-slowdown hypothesis).31,35
However, these findings are consistent with what we know about substitution rates in primates. Based on a genome-wide analysis, strepsirhines and New World monkeys are known to have a significantly slower rate, about half that of apes and Old World monkeys.
36 Similarly, the mtDNA specifically is known to have a slower rate in lemurs than anthropoids, especially the rate of transitions.
37 In contrast, strepsirhines have a faster generation time than apes and monkeys (
t-test of log transformed maturity time:
P-value = 7.18 × 10
−14 for females only;
56P-value = 8.42 × 10
−16 for a maturity time function of both male and female
83) and New World monkeys also have a faster generation time than apes and Old World monkeys (
P-value = 1.98 × 10
−6;
56P-value = 1.85 × 10
−5 83). Thus, the generation time hypothesis predicts the following order for rates, from lowest to highest: apes, Old World monkeys, New World monkeys, strepsirhines, but the actual order based on
36 is: New World monkeys and strepsirhines, apes, Old World monkeys. This implies some additional trait may be involved in the differences between apes and Old World monkeys, and that the hominoid slowdown hypothesis may have alternative explanations other than generation time. Similarly, if the hominoid slowdown hypothesis is explained by generation time, then some other feature is obscuring the general signal in primates.
Another potential explanation for our lack of significant association is that the generation time signal is associated with longer evolutionary divergence that is too noisy when examining closely related species such as primates. One study examining generation time in invertebrates found that it was the best predictor of rate variation, but over a considerably longer evolutionary time than that considered here.
79 Even among mammals, the time may not be sufficient to detect a significant association other than that found in common with other life history variables.
9,10,12,39,79
Some caution is necessary in evaluating the generation time hypothesis for primates, because the phenotypic data are less complete with age at female sexual maturity and we had to make a greater number of genus level assignments to generate the taxon mapping of the discretization (see Methods). Further, only two of the 26 taxa were in the smallest phenotypic category. This may have effectively reduced our test from a three category to a two category test but significance still evaluated by a Ç
2 distribution with degrees of freedom equal to two. We reexamined this phenotype by reducing the number of categories to two (based on either the median or the
P0.66 threshold) and compared it to a Ç
2 distribution with one degree of freedom. In both cases, we still failed to achieve any level of significance greater than that reported in
Table 3.
Longevity hypothesis
The longevity hypothesis predicts that species where long life is adaptive will show selection on the molecular pathways that mitigate or reduce the production of radical oxidative species and as well as the molecular machinery involved in mtDNA replication. As a result, the selective effect should reduce the rate of mutation in the mtDNA.
23
Previous studies have supported this hypothesis in one of two ways: biochemistry and molecular phylogenetic rate variation. Biochemically, it has been shown that long-lived species have reduced rates of reactive oxygen species (ROS) generation and increased resistance to damage in structural components (for reviews see Pamplona and Barja
24 and Galtier et al
23). In terms of rate variation, the third codon position of
CYTB in mammals and birds was shown to be correlated with substitution rate.
10,12 Further, maximum lifespan was the strongest predictor of synonymous rate for the nine longest protein-coding genes in the mtDNA and some nuclear encoded genes, even with multiple regressions.
9 This association was found across all Mammalia using 160 taxa and concatenated genes (all nine coding genes as a single alignment) in order to obtain sufficient power.
9 Here, we have included a larger number of loci (including non-protein-coding mtDNA genes), but a smaller number of taxa and we were still able to detect an association. We only required 26 taxa within a single order and used separate genes to detect the association. Moreover, we were able to examine the control region because of the close evolutionary relationship of primates and because our model uses nucleotide models for the genotype substitution.
Further, these results in primates are important as some studies have found either a lack of association
9 (when only examining primates and rodents) or an insignificant mammalian association with primates and a negative association when primates are excluded.
38 One possibility is the difference in power between tests, as mentioned previously.
We found a significant association (
Table 3) with at least one representative alignment for all gene groups of Galtier et al.
38 As previously stated, we analyzed gene by gene not by gene groups. Our results also imply that the primary signal comes from neutral changes rather than positive selection or relaxation of negative selection. This third codon signal would have been missed by Galtier et al
38 as that study only examined amino acids, but our results for primates are consistent with their overall conclusion of a negative association of neutral substitution and longevity across mammals for the mtDNA.
General life history and metabolic rate
Body mass is associated with many aspects of life history.
84–88 We have so far interpreted the body mass association as one of general life history, but many studies have used body mass as a surrogate for metabolic rate.
The metabolic rate hypothesis predicts that organisms with a higher metabolic rate should have increased production of ROS and thus higher mutation rates in the mtDNA.
23 Some have hypothesized that the correlation between metabolic rate and substitution rate is general enough to recalibrate the molecular clock.
19,89 The metabolic rate hypothesis has had some support within birds and mammals,
41,90 but these studies have been criticized for their small data sets.
91 A few investigations using larger data sets have failed to find a correlation between mass-specific metabolic and substitution rates, again within birds and mammals,
25 and also in invertebrates.
91,92 Other papers with large data sets have found a negative correlation between body mass and substitution rate,
10,12 but for mammals this is removed when regressed with other life history variables.
10
Our results support a connection between body mass and substitution rate that may be separate from the metabolic rate hypothesis. Body mass is a known correlate of many life history and ecological traits
84–88 and it is plausible that an unknown phenotype, or longevity alone as indicated by the Wilcoxon Test between their ALRSs, is the root correlate.
Conclusions
We have successfully extended the method of
16 to account for a greater number of phenotypic categories and tested it on coding and non-coding sequence data from primate mitochondrial and nuclear genomes. Simulation studies showed that the method performs well, with a false-positive rate within acceptable values for a 5% cutoff as well as reasonable power (>80%) when there is greater than a three-fold change in substitution rate and greater than a tree length of three. Further, we have shown that a stricter alternative hypothesis (
W1 >
W0) is preferred over a more relaxed hypothesis for two or three phenotype categories. In simulation, the number of phenotypic categories was tested up to four and shown to have a low FP rate, but lower power with the four phenotype category case.
When applied to primate data with three phenotypic categories, the method identified many of the genes from the mtDNA as having an association with maximum lifespan and body mass consistent with the longevity hypothesis and general life history respectively. The significant signal predominantly originates from the third codon position for both hypotheses in the mtDNA protein-coding genes, implying a neutral signal on the mtDNA. The control region, which is non-coding and highly divergent, was also significantly associated with maximum lifespan and again supports the prediction of the longevity hypothesis. Those few nuclear encoded genes that were significant have a different pattern of contribution from the codon positions that indicate a separate and unknown explanation for the associations.
Primates had previously been identified as an outlier group among mammals for the maximum lifespan and body mass associations as it was inconsistent with the general mammalian pattern.
9,38 Also, some signal of adaptive evolution has been identified in the primate mtDNA.
22,93 Examining this pattern with our method, we show that primate mtDNA is evolving in a manner consistent with previous predictors for the phenotypes and that the signal coming from the third codon position is indicative of the neutral rate, which may have been overlooked in some of the previous studies. We find little support for the generation time hypothesis, which is inconsistent with the typical rationale for the hominoid slowdown,
31,33,35 but consistent with what we know of substitution rates between apes and Old World monkeys, and New World monkeys and strepsirhine primates.
36
Our method is a useful tool in identifying variation in the neutral rate of evolution associated with phenotype. Here, it has detected an association with longevity and body size, and mtDNA rate variation in primates missed by other methods, showing that primates are consistent with other mammals in this respect.