The Impact of Modelling Rate Heterogeneity among Sites on Phylogenetic Estimates of Intraspecific Evolutionary Rates and Timescales

Phylogenetic analyses of DNA sequence data can provide estimates of evolutionary rates and timescales. Nearly all phylogenetic methods rely on accurate models of nucleotide substitution. A key feature of molecular evolution is the heterogeneity of substitution rates among sites, which is often modelled using a discrete gamma distribution. A widely used derivative of this is the gamma-invariable mixture model, which assumes that a proportion of sites in the sequence are completely resistant to change, while substitution rates at the remaining sites are gamma-distributed. For data sampled at the intraspecific level, however, biological assumptions involved in the invariable-sites model are commonly violated. We examined the use of these models in analyses of five intraspecific data sets. We show that using 6–10 rate categories for the discrete gamma distribution of rates among sites is sufficient to provide a good approximation of the marginal likelihood. Increasing the number of gamma rate categories did not have a substantial effect on estimates of the substitution rate or coalescence time, unless rates varied strongly among sites in a non-gamma-distributed manner. The assumption of a proportion of invariable sites provided a better approximation of the asymptotic marginal likelihood when the number of gamma categories was small, but had minimal impact on estimates of rates and coalescence times. However, the estimated proportion of invariable sites was highly susceptible to changes in the number of gamma rate categories. The concurrent use of gamma and invariable-site models for intraspecific data is not biologically meaningful and has been challenged on statistical grounds; here we have found that the assumption of a proportion of invariable sites has no obvious impact on Bayesian estimates of rates and timescales from intraspecific data.     

Exploring among-site rate variation models in a maximum likelihood framework using empirical data: effects of model assumptions on estimates of topology,branch lengths,and bootstrap support

We have investigated the effects of different among-site rate variation models on the estimation of substitution model parameters, branch lengths, topology, and bootstrap proportions under minimum evolution (ME) and maximum likelihood (ML). Specifically, we examined equal rates, invariable sites, gamma-distributed rates, and site-specific rates (SSR) models, using mitochondrial DNA sequence data from three protein-coding genes and one tRNA gene from species of the New Zealand cicada genus Maoricicada. Estimates of topology were relatively insensitive to the substitution model used; however, estimates of bootstrap support, branch lengths, and R-matrices (underlying relative substitution rate matrix) were strongly influenced by the assumptions of the substitution model. We identified one situation where ME and ML tree building became inaccurate when implemented with an inappropriate among-site rate variation model. Despite the fact the SSR models often have a better fit to the data than do invariable sites and gamma rates models, SSR models have some serious weaknesses. First, SSR rate parameters are not comparable across data sets, unlike the proportion of invariable sites or the alpha shape parameter of the gamma distribution. Second, the extreme among-site rate variation within codon positions is problematic for SSR models, which explicitly assume rate homogeneity within each rate class. Third, the SSR models appear to give severe underestimates of R-matrices and branch lengths relative to invariable sites and gamma rates models in this example. We recommend performing phylogenetic analyses under a range of substitution models to test the effects of model assumptions not only on estimates of topology but also on estimates of branch length and nodal support.     

Substitution rate variation among sites in mitochondrial hypervariable region I of humans and chimpanzees

Mitochondrial D-loop hypervariable region I (HVI) sequences are widely used in human molecular evolutionary studies, and therefore accurate assessment of rate heterogeneity among sites is essential. We used the maximum-likelihood method to estimate the gamma shape parameter alpha for variable substitution rates among sites for HVI from humans and chimpanzees to provide estimates for future studies. The complete data of 839 humans and 224 chimpanzees, as well as many subsets of these data, were analyzed to examine the effect of sequence sampling. The effects of the genealogical tree and the nucleotide substitution model were also examined. The transition/transversion rate ratio (kappa) is estimated to be about 25, although much larger and biased estimates were also obtained from small data sets at low divergences. Estimates of alpha were 0.28-0.39 for human data sets of different sizes and 0.20-0.39 for data sets including different chimpanzee subspecies. The combined data set of both species gave estimates of 0.42-0.45. While all those estimates suggest highly variable substitution rates among sites, smaller samples tend to give smaller estimates of alpha. Possible causes for this pattern were examined, such as biases in the estimation procedure and shifts in the rate distribution along certain lineages. Computer simulations suggest that the estimation procedure is quite reliable for large trees but can be biased for small samples at low divergences. Thus, an alpha of 0.4 appears suitable for both humans and chimpanzees. Estimates of alpha can be affected by the nucleotide sites included in the data, the overall tree length (the amount of sequence divergence), the number of rate classes used for the estimation, and to a lesser extent, the included sequences. The genealogical tree, the substitution model, and demographic processes such as population expansion do not have much effect.     

Time of the deepest root for polymorphism in human mitochondrial DNA

Summary A molecular clock analysis was carried out on the nucleotide sequences of parts of the major noncoding region of mitochondrial DNA (mtDNA) from the major geographic populations of humans. Dates of branchings in the mtDNA tree among humans were estimated with an improved maximum likelihood method. Two species of chimpanzees were used as an outgroup, and the mtDNA clock was calibrated by assuming that the chimpanzee/human split occurred 4 million years ago, following our earlier works. A model of homogeneous evolution among sites does not fit well with the data even within hypervariable segments, and hence an additional parameter that represents a proportion of variable sites was introduced. Taking account of this heterogeneity among sites, the date for the deepest root of the mtDNA tree among humans was estimated to be 280,000±50,000 years old (±1 SE), although there remains uncertainty about the constancy of the evolutionary rate among lineages. The evolutionary rate of the most rapidly evolving sites in mtDNA was estimated to be more than 100 times greater than that of a nuclear pseudogene.     

Maximum likelihood estimation of phylogenetic trees is consistent when substitution rates vary according to the invariable sites plus gamma distribution

Maximum likelihood estimation of phylogenetic trees from nucleotide sequences is completely consistent when nucleotide substitution is governed by the general time reversible (GTR) model with rates that vary over sites according to the invariable sites plus gamma (I + gamma) distribution.     

Identifying evolutionary trees and substitution parameters for the general Markov model with invariable sites

The general Markov plus invariable sites (GM+I) model of biological sequence evolution is a two-class model in which an unknown proportion of sites are not allowed to change, while the remainder undergo substitutions according to a Markov process on a tree. For statistical use it is important to know if the model is identifiable; can both the tree topology and the numerical parameters be determined from a joint distribution describing sequences only at the leaves of the tree? We establish that for generic parameters both the tree and all numerical parameter values can be recovered, up to clearly understood issues of 'label swapping'. The method of analysis is algebraic, using phylogenetic invariants to study the variety defined by the model. Simple rational formulas, expressed in terms of determinantal ratios, are found for recovering numerical parameters describing the invariable sites.     

Hierarchical analysis of variation in the mitochondrial 16S rRNA gene among Hymenoptera

Nucleotide sequences from a 434-bp region of the 16S rRNA gene were analyzed for 65 taxa of Hymenoptera (ants, bees, wasps, parasitoid wasps, sawflies) to examine the patterns of variation within the gene fragment and the taxonomic levels for which it shows maximum utility in phylogeny estimation. A hierarchical approach was adopted in the study through comparison of levels of sequence variation among taxa at different taxonomic levels. As previously reported for many holometabolous insects, the 16S data reported here for Hymenoptera are highly AT-rich and exhibit strong site-to-site variation in substitution rate. More precise estimates of the shape parameter (alpha) of the gamma distribution and the proportion of invariant sites were obtained in this study by employing a reference phylogeny and utilizing maximum-likelihood estimation. The effectiveness of this approach to recovering expected phylogenies of selected hymenopteran taxa has been tested against the use of maximum parsimony. This study finds that the 16S gene is most informative for phylogenetic analysis at two different levels: among closely related species or populations, and among tribes, subfamilies, and families. Maximization of the phylogenetic signal extracted from the 16S gene at higher taxonomic levels may require consideration of the base composition bias and the site-to-site rate variation in a maximum-likelihood framework.      

The rate heterogeneity of nonsynonymous substitutions in mammalian mitochondrial genes

Substitution rates at the three codon positions (r1, r2, and r3) of mammalian mitochondrial genes are in the order of r3 > r1 > r2, and the rate heterogeneity at the three positions, as measured by the shape parameter of the gamma distribution (alpha 1, alpha 2, and alpha 3), is in the order of alpha 3 > alpha 1 > alpha 2. The causes for the rate heterogeneity at the three codon positions remain unclear and, in particular, there has been no satisfactory explanation for the observation of alpha 1 > alpha 2. I attempted to dissect the causes of rate heterogeneity by studying the pattern of nonsynonymous substitutions with respect to codon positions in 10 mitochondrial genes from 19 mammalian species. Nonsynonymous substitutions involve more different amino acid replacements at the second than at the first codon position, which results in r1 > r2. The difference between r1 and r2 increases with the intensity of purifying selection, and so does the rate heterogeneity in nonsynonymous substitutions among sites at the same codon position. All mitochondrial genes appear to have functionally important and unimportant codons, with the latter having all three codon positions prone to nonsynonymous substitutions. Within the functionally important codons, the second codon position is much more conservative than the codon position. This explains why alpha 1 > alpha 2. The result suggests that overweighting of the second codon position in phylogenetic analysis may be a misguided practice.      

Maximum-Likelihood Analysis of the Tethythere Hypothesis Based on a Multigene Data Set and a Comparison of Different Models of Sequence Evolution

Hyracoids have been allied with either perissodactyls or tethytheres (i.e., Proboscidea + Sirenia) based on morphological data. The latter hypothesis, termed Paenungulata, is corroborated by numerous molecular studies. However, molecular studies have failed to support Tethytheria, a group that is supported by morphological data. We examined relationships among living paenungulate orders using a multigene data set that included sequences from four mitochondrial genes (12S rRNA, tRNA valine, 16S rRNA, cytochrome b) and four nuclear genes (aquaporin, A2AB, IRBP, vWF). Nineteen maximum-likelihood models were employed, including models with process partitions for base composition and substitution parameterizations. With the inclusion of partitions with a heterogeneous base composition, 18 of 19 models favored Hyracoidea + Sirenia. All 19 models favored Hyracoidea + Sirenia after excluding heterogeneous base composition partitions. Most of the support for Hyracoidea + Sirenia derived from the mitochondrial genes (bootstrap support ranged from 51 to 99%); Tethytheria, in turn, received 0 to 19% support in different analyses. Bootstrap support deriving from the nuclear genes was more evenly split among the competing hypotheses (3 to 45% for Tethytheria; 17.5 to 62% for Hyracoidea + Sirenia). Lineage-specific rate variation among both mitochondrial and nuclear genes may contribute to the different results that were obtained with mitochondrial versus nuclear data. Whether Tethytheria or a competing hypothesis is correct, short internodes on the molecular phylogenies suggest that paenungulate orders diverged from each other over a 5- to 8-million-year time window extending from the late Paleocene into the early Eocene. We also used likelihood-ratio tests to compare different models of sequence evolution. A gamma distribution of rates results in a greater improvement in likelihood scores than does an allowance for invariant sites. Twenty-one rate partitions corresponding to stems, loops, and codon positions of different genes result in higher likelihood scores than a gamma distribution of rates and/or an allowance for invariant sites. Process partitions of the data that incorporate base composition and substitution parameterizations result in significant improvements in likelihood scores in comparison to models that allow only for relative rate differences among partitions.     

Testing the covarion hypothesis of molecular evolution

The covarion hypothesis of molecular evolution states that the fixation of mutations may alter the probability that any given position will fix the next change. Tests of this hypothesis using the divergence of real sequences are compromised because models of rate variation among sites (e.g., the gamma version of the one-parameter equation) predict sequence divergence values similar to those for the covarion process. This study therefore focuses on the extent to which the varied and unvaried codons of two well-diverged taxa are the same, because fewer are expected by the covarion hypothesis than by the gamma model. The data for these tests are the protein sequences of Cu, Zn superoxide dismutase (SOD) for mammals and plants. Simulation analyses show that the covarion hypothesis makes better predictions about the frequencies of varied and unhit positions in common between these two taxa than does the gamma version of the one-parameter model. Furthermore, the analysis of SOD tertiary structure demonstrates that mammal and plant variabilities are distributed differently on the protein. These results support the conclusions that the variable and invariable codons of mammal and plant SODs are different and that the covarion model explains the evolution of this protein better than the gamma version of the one-parameter process. Unlike other models, the covarion hypothesis accounts for rate fluctuations among positions over time, which is an important parameter of molecular evolution.      

HLA-D region antigen-associated invariant polypeptides as revealed by two-dimensional gel analysis. Glycosylation and structural inter-relationships

Two-dimensional polyacrylamide gel analyses of immunoprecipitates of HLA-D region antigens prepared from [35S]methionine-labeled B lymphoblastoid cells revealed a number of invariant polypeptides (Ii and theta) that co-precipitate with the alpha and beta polypeptides of the class II (Ia) antigens. The invariant polypeptides comprised at least three Ii spots of Mr = 31,000 (Ii1-Ii3) and a series of six theta spots of Mr = 34,000 (theta 1-theta 6). The structural inter-relationships of these polypeptides have been investigated. Tryptic peptide fingerprints showed that Ii and theta have closely related amino acid sequences. In contrast, the fingerprints of the HLA-DR alpha and beta polypeptides clearly differed from those of theta and Ii as well as from each other. Analyses of immunoprecipitates prepared from cells cultured in the presence of tunicamycin revealed the presence of two N-linked oligosaccharides on each invariant polypeptide and suggested that the more acidic theta polypeptides (theta 1 and theta 2) differed from the other invariant polypeptides by the presence of sialic acid on one or both N-linked oligosaccharides. Removal of sialic acid by neuraminidase simplified the pattern of theta spots into three distinct Ii-related polypeptides. Endo-beta-N-acetylglycosaminidase H digestion indicated that the individual theta polypeptides represent stages in carbohydrate processing whereby Ii with two N-linked immature oligosaccharides are converted initially to theta 6-theta 3 with one immature and one complex, but nonsialylated, oligosaccharide and finally to theta 2-theta 1 with two complex oligosaccharides. Digestion of the theta polypeptides with N-acetylgalactosamine oligosaccharidase indicated that the theta spots are also derived by O-glycosylation from the Ii polypeptides. This assignment is supported by results obtained using monensin to block glycosylation within the Golgi. At least three spots persisted after complete removal of the N- and O-linked oligosaccharides, suggesting the presence of a family of invariant polypeptides differing in amino acid sequence.     

Invariable sites models and their use in phylogeny reconstruction

Phylogenetic inference is well known to be problematic if both long and short branches occur together in the underlying tree. With biological data, correcting for this problem may require simultaneous consideration for both substitution biases and rate heterogeneity between lineages and across sequence positions. A particular form of the latter is the presence of invariable sites, which are well known to mislead estimation of genetic divergences. Here we describe a capture-recapture method to estimate the proportion of invariable sites in an alignment of amino acids or nucleotides. We use it to investigate phylogenetic signals in 18S ribosomal DNA sequences from Holometabolus insects. Our results suggest that, as taxa diverged, their 18S rDNA sequences have altered in both their distribution of sites that can vary as well as in their base compositions.     

Patterns of Nucleotide Substitution in Mitochondrial Protein Coding Genes of Vertebrates

Maximum likelihood methods were used to study the differences in substitution rates among the four nucleotides and among different nucleotide sites in mitochondrial protein-coding genes of vertebrates. In the 1st+2nd codon position data, the frequency of nucleotide G is negatively correlated with evolutionary rates of genes, substitution rates vary substantially among sites, and the transition/transversion rate bias (R) is two to five times larger than that expected at random. Generally, largest transition biases and greatest differences in substitution rates among sites are found in the highly conserved genes. The 3rd positions in placental mammal genes exhibit strong nucleotide composition biases and the transitional rates exceed transversional rates by one to two orders of magnitude. Tamura-Nei and Hasegawa-Kishino-Yano models with gamma distributed variable rates among sites (gamma parameter, α) adequately describe the nucleotide substitution process in 1st+2nd position data. In these data, ignoring differences in substitution rates among sites leads to largest biases while estimating substitution rates. Kimura's two-parameter model with variable-rates among sites performs satisfactorily in likelihood estimation of R, α, and overall amount of evolution for 1st+2nd position data. It can also be used to estimate pairwise distances with appropriate values of α for a majority of genes.     

Toward a more accurate time scale for the human mitochondrial DNA tree

Several estimates of the time of occurrence of the most recent common mitochondrial DNA (mtDNA) ancestor of modern humans have been made. Estimates derived from noncoding regions based on a model that classifies sites into two categories (variable and invariable) have been consistently older than those derived from the third positions of codons. This discrepancy can be attributed to a violation of the assumption of rate homogeneity among variable sites when analyzing the noncoding regions. Additional data from the partial control region sequences allow us to take into account some of this further heterogeneity. By assigning the sites to three classes (highly variable, moderately variable, and invariable) and by assuming that the last common mtDNA ancestor of humans and chimpanzees lived 4 million years ago, the most recent common mtDNA ancestor of humans is estimated to have occurred 211,000 ±111,000 years ago (±1 SE), consistent with the estimate, 101,000 ± 52,000 years, made from third positions of codons and also with those proposed previously. We used the same technique to estimate when a putative expansion of modern humans out of Africa took place and estimated a time of 89,000 ± 69,000 years ago. Even though the standard errors of these estimates are large, they allow us to reject the multiregional hypothesis of modern human origin.Deceased July 21, 1991 Correspondence to: M. Hasegawa     

SSU rRNA-based phylogenetic position of the genera Amoeba and Chaos (Lobosea, Gymnamoebia): the origin of gymnamoebae revisited.

Naked lobose amoebae (gymnamoebae) are among the most abundant group of protists present in all aquatic and terrestrial biotopes. Yet, because of lack of informative morphological characters, the origin and evolutionary history of gymnamoebae are poorly known. The first molecular studies revealed multiple origins for the amoeboid lineages and an extraordinary diversity of amoebae species. Molecular data, however, exist only for a few species of the numerous taxa belonging to this group. Here, we present the small-subunit (SSU) rDNA sequences of four species of typical large gymnamoebae: Amoeba proteus, Amoeba leningradensis, Chaos nobile, and Chaos carolinense. Sequence analysis suggests that the four species are closely related to the species of genera Saccamoeba, Leptomyxa, Rhizamoeba, Paraflabellula, Hartmannella, and Echinamoeba. All of them form a relatively well-supported clade, which corresponds to the subclass Gymnamoebia, in agreement with morphology-based taxonomy. The other gymnamoebae cluster in small groups or branch separately. Their relationships change depending on the type of analysis and the model of nucleotide substitution. All gymnamoebae branch together in Neighbor-Joining analysis with corrections for among-site rate heterogeneity and proportion of invariable sites. This clade, however, is not statistically supported by SSU rRNA gene sequences and further analysis of protein sequence data will be necessary to test the monophyly of gymnamoebae.     

Analysis of elicited EEG synchronization and desynchronization in emotional activation: temporal and topographic characteristics

Event-related synchronization (ERS) and desynchronization (ERD) in delta, theta1, theta2, alpha1, alpha2, beta1, beta2, beta3, and gamma were measured in 20 healthy right-handed subjects in response to IAPS stimuli with low, moderate, and high arousal reactions. The 62-channel EEG was simultaneously recorded while subjects viewed sequentially presented pictures and subjectively rated them after each presentation. The results show that emotionally loaded stimuli induced higher ERS in the delta, theta1, theta2, beta1, beta3, and gamma bands along with combined ERD and ERS effects in alpha2 band. As to hemispheric asymmetries, the effects of emotional arousal were restricted not only to right parietal (theta1 and theta2 ERS, alpha2 ERD) but also to left frontal (theta2 ERS) regions. In terms of affective chronometry, lower theta was the first to catch the affective salience of incoming stimuli (time window 0-600 ms after the stimulus input). For theta2, alpha2, and gamma bands this process was delayed to 600-1000 ms.     

Estimating the contribution of mutation, recombination and gene conversion in the generation of haplotypic diversity

Recombination occurs through both homologous crossing over and homologous gene conversion during meiosis. The contribution of recombination relative to mutation is expected to be dramatically reduced in inbreeding organisms. We report coalescent-based estimates of the recombination parameter (rho) relative to estimates of the mutation parameter (theta) for 18 genes from the highly self-fertilizing grass, wild barley, Hordeum vulgare ssp. spontaneum. Estimates of rho/theta are much greater than expected, with a mean rho/theta approximately 1.5, similar to estimates from outcrossing species. We also estimate rho with and without the contribution of gene conversion. Genotyping errors can mimic the effect of gene conversion, upwardly biasing estimates of the role of conversion. Thus we report a novel method for identifying genotyping errors in nucleotide sequence data sets. We show that there is evidence for gene conversion in many large nucleotide sequence data sets including our data that have been purged of all detectable sequencing errors and in data sets from Drosophila melanogaster, D. simulans, and Zea mays. In total, 13 of 27 loci show evidence of gene conversion. For these loci, gene conversion is estimated to contribute an average of twice as much as crossing over to total recombination.     

Effects of Mg2+ and the beta gamma-subunit complex on the interactions of guanine nucleotides with G proteins

Mg2+ interacts with the alpha subunits of guanine nucleotide-binding regulatory proteins (G proteins) in the presence of guanosine-5'-[gamma-thio]triphosphate (GTP-gamma S) to form a highly fluorescent complex from which nucleotide dissociates very slowly. The apparent Kd for interaction of G alpha X GTP gamma S with Mg2+ is approximately 5 nM, similar to the Km for G protein GTPase activity X G beta gamma increases the rate of dissociation of GTP gamma S from G alpha X GTP gamma S or G alpha X GTP gamma S X Mg2+ at low concentrations of Mg2+. When the concentration of Mg2+ exceeds 1 mM, G beta gamma dissociates from G beta gamma X G alpha X GTP gamma S X Mg2+. Compared with the dramatic effect of Mg2+ on binding of GTP gamma S to G alpha, the metal has relatively little effect on the binding of GDP. However, G beta gamma increases the affinity of G alpha for GDP by more than 100-fold. High concentrations of Mg2+ promote the dissociation of GDP from G beta gamma X G alpha X GDP, apparently without causing subunit dissociation. The steady-state rate of GTP hydrolysis is strictly correlated with the rate of dissociation of GDP from G alpha under all conditions examined. Thus, there are at least two sites for interaction of Mg2+ with G protein-nucleotide complexes. Furthermore, binding of G beta gamma and GTP gamma S to G alpha is negatively cooperative, while the binding interaction between G beta gamma and GDP is strongly positive.     

gamma delta T cells in the human intestine express surface markers of activation and are preferentially located in the epithelium

We have investigated the expression of the alpha beta and the gamma delta T cell receptor (TCR) in the human intestine. By immunohistology we found that 39% of CD3+ intraepithelial lymphocytes (IEL) expressed the gamma delta TCR compared to 3% of CD3+ lamina propria lymphocytes (LPL). Cytofluorometric analysis of isolated cells revealed a significantly higher proportion of gamma delta T cells among CD3+ IEL compared to LPL and peripheral blood lymphocytes. This was due to an increase in both CD8+ (low density) and CD4-CD8- gamma delta T cells in IEL. Most alpha beta IEL expressed high-density CD8. About 30% of both IEL and LPL expressed CD25 (alpha-chain of the IL-2 receptor). HML-1 expression was detected on nearly all IEL and on 27% of LPL. CD25 and HML-1 were preferentially expressed on intestinal alpha beta and gamma delta T cells, respectively. Thus, human gamma delta T cells are located preferentially in the gut epithelium and are phenotypically different from alpha beta T cells, which constitute the majority of both LPL and IEL.