TREES SIFTER 1.0: an approximate method to estimate the time to the most recent common ancestor of a sample of DNA sequences

trees sifter 1.0 implements an approximate method to estimate the time to the most recent common ancestor (TMRCA) of a set of DNA sequences, using population evolution modelling. In essence, the program simulates genealogies with a user‐defined model of coalescence of lineages, and then compares each simulated genealogy to the genealogy inferred from the real data, through two summary statistics: (i) the number of mutations on the genealogy (M_n), and (ii) the number of different sequence types (alleles) observed (K_n). The simulated genealogies are then submitted to a rejection algorithm that keeps only those that are the most likely to have generated the observed sequence data. At the end of the process, the accepted genealogies can be used to estimate the posterior probability distribution of the TMRCA.     

LIKELIHOOD-BASED INFERENCE IN ISOLATION-BY-DISTANCE MODELS USING THE SPATIAL DISTRIBUTION OF LOW-FREQUENCY ALLELES

Estimating dispersal distances from population genetic data provides an important alternative to logistically taxing methods for directly observing dispersal. Although methods for estimating dispersal rates between a modest number of discrete demes are well developed, methods of inference applicable to "isolation-by-distance" models are much less established. Here, we present a method for estimating ρσ², the product of population density (ρ) and the variance of the dispersal displacement distribution (σ²). The method is based on the assumption that low-frequency alleles are identical by descent. Hence, the extent of geographic clustering of such alleles, relative to their frequency in the population, provides information about ρσ². We show that a novel likelihood-based method can infer this composite parameter with a modest bias in a lattice model of isolation-by-distance. For calculating the likelihood, we use an importance sampling approach to average over the unobserved intraallelic genealogies, where the intraallelic genealogies are modeled as a pure birth process. The approach also leads to a likelihood-ratio test of isotropy of dispersal, that is, whether dispersal distances on two axes are different. We test the performance of our methods using simulations of new mutations in a lattice model and illustrate its use with a dataset from Arabidopsis thaliana .     

Estimating a geographically explicit model of population divergence

Patterns of genetic variation can provide valuable insights for deciphering the relative roles of different evolutionary processes in species differentiation. However, population-genetic models for studying divergence in geographically structured species are generally lacking. Since these are the biogeographic settings where genetic drift is expected to predominate, not only are population-genetic tests of hypotheses in geographically structured species constrained, but generalizations about the evolutionary processes that promote species divergence may also be potentially biased. Here we estimate a population-divergence model in montane grasshoppers from the sky islands of the Rocky Mountains. Because this region was directly impacted by Pleistocene glaciation, both the displacement into glacial refugia and recolonization of montane habitats may contribute to differentiation. Building on the tradition of using information from the genealogical relationships of alleles to infer the geography of divergence, here the additional consideration of the process of gene-lineage sorting is used to obtain a quantitative estimate of population relationships and historical associations (i.e., a population tree) from the gene trees of five anonymous nuclear loci and one mitochondrial locus in the broadly distributed species Melanoplus oregonensis. Three different approaches are used to estimate a model of population divergence; this comparison allows us to evaluate specific methodological assumptions that influence the estimated history of divergence. A model of population divergence was identified that significantly fits the data better compared to the other approaches, based on per-site likelihood scores of the multiple loci, and that provides clues about how divergence proceeded in M. oregonensis during the dynamic Pleistocene. Unlike the approaches that either considered only the most recent coalescence (i.e., information from a single individual per population) or did not consider the pattern of coalescence in the gene genealogies, the population-divergence model that best fits the data was estimated by considering the pattern of gene lineage coalescence across multiple individuals, as well as loci. These results indicate that sampling of multiple individuals per population is critical to obtaining an accurate estimate of the history of divergence so that the signal of common ancestry can be separated from the confounding influence of gene flow-even though estimates suggest that gene flow is not a predominant factor structuring patterns of genetic variation across these sky island populations. They also suggest that the gene genealogies contain information about population relationships, despite the lack of complete sorting of gene lineages. What emerges from the analyses is a model of population divergence that incorporates both contemporary distributions and historical associations, and shows a latitudinal and regional structuring of populations reminiscent of population displacements into multiple glacial refugia. Because the population-divergence model itself is built upon the specific events shaping the history of M. oregonensis, it provides a framework for estimating additional population-genetic parameters relevant to understanding the processes governing differentiation in geographically structured species and avoids the problems of relying on overly simplified and inaccurate divergence models. The utility of these approaches, as well as the caveats and future improvements, for estimating population relationships and historical associations relevant to genetic analyses of geographically structured species are discussed.     

Evolution of the human immunodeficiency virus envelope gene is dominated by purifying selection

The evolution of the human immunodeficiency virus (HIV-1) during chronic infection involves the rapid, continuous turnover of genetic diversity. However, the role of natural selection, relative to random genetic drift, in governing this process is unclear. We tested a stochastic model of genetic drift using partial envelope sequences sampled longitudinally in 28 infected children. In each case the Bayesian posterior (empirical) distribution of coalescent genealogies was estimated using Markov chain Monte Carlo methods. Posterior predictive simulation was then used to generate a null distribution of genealogies assuming neutrality, with the null and empirical distributions compared using four genealogy-based summary statistics sensitive to nonneutral evolution. Because both null and empirical distributions were generated within a coalescent framework, we were able to explicitly account for the confounding influence of demography. From the distribution of corrected P-values across patients, we conclude that empirical genealogies are more asymmetric than expected if evolution is driven by mutation and genetic drift only, with an excess of low-frequency polymorphisms in the population. This indicates that although drift may still play an important role, natural selection has a strong influence on the evolution of HIV-1 envelope. A negative relationship between effective population size and substitution rate indicates that as the efficacy of selection increases, a smaller proportion of mutations approach fixation in the population. This suggests the presence of deleterious mutations. We therefore conclude that intrahost HIV-1 evolution in envelope is dominated by purifying selection against low-frequency deleterious mutations that do not reach fixation.     

The strong-migration limit for the genealogical process in geographically structured populations

In this article we assume that the entire population is subdivided into a finite number of panmictic colonies, each of which consists of a respective number of haploid individuals. We also assume that random genetic drift occurs in each colony and migration among colonies, which is independent of time and ergodic. We study the genealogical process of sampled genes from geographically structured populations. We prove that if the actual total population number is replaced by the effective population number, the mean coalescence time converges to that in a panmictic population in the strong migration limit. We also obtain the geographical distribution of the common ancestor.     

Detecting Small Amounts of Gene Flow from Phylogenies of Alleles

The method of coalescents is used to find the probability that none of the ancestors of alleles sampled from a population are immigrants. If that is the case for samples from two or more populations, then there would be concordance between the phylogenies of those alleles and the geographic locations from which they are drawn. This type of concordance has been found in several studies of mitochondrial DNA from natural populations. It is shown that if the number of sequences sampled from each population is reasonably large (10 or more), then this type of concordance suggests that the average number of individuals migrating between populations is likely to be relatively small (Nm less than 1) but the possibility of occasional migrants cannot be excluded. The method is applied to the data of E. Bermingham and J. C. Avise on mtDNA from the bowfin, Amia calva.     

The Sampling Distribution of Disease-Associated Alleles

A theory is developed that provides the sampling distribution of low frequency alleles at a single locus under the assumption that each allele is the result of a unique mutation. The numbers of copies of each allele is assumed to follow a linear birth-death process with sampling. If the population is of constant size, standard results from theory of birth-death processes show that the distribution of numbers of copies of each allele is logarithmic and that the joint distribution of numbers of copies of k alleles found in a sample of size n follows the Ewens sampling distribution. If the population from which the sample was obtained was increasing in size, if there are different selective classes of alleles, or if there are differences in penetrance among alleles, the Ewens distribution no longer applies. Likelihood functions for a given set of observations are obtained under different alternative hypotheses. These results are applied to published data from the BRCA1 locus (associated with early onset breast cancer) and the factor VIII locus (associated with hemophilia A) in humans. In both cases, the sampling distribution of alleles allows rejection of the null hypothesis, but relatively small deviations from the null model can account for the data. In particular, roughly the same population growth rate appears consistent with both data sets.     

On the size distribution of private microsatellite alleles

Private microsatellite alleles tend to be found in the tails rather than in the interior of the allele size distribution. To explain this phenomenon, we have investigated the size distribution of private alleles in a coalescent model of two populations, assuming the symmetric stepwise mutation model as the mode of microsatellite mutation. For the case in which four alleles are sampled, two from each population, we condition on the configuration in which three distinct allele sizes are present, one of which is common to both populations, one of which is private to one population, and the third of which is private to the other population. Conditional on this configuration, we calculate the probability that the two private alleles occupy the two tails of the size distribution. This probability, which increases as a function of mutation rate and divergence time between the two populations, is seen to be greater than the value that would be predicted if there was no relationship between privacy and location in the allele size distribution. In accordance with the prediction of the model, we find that in pairs of human populations, the frequency with which private microsatellite alleles occur in the tails of the allele size distribution increases as a function of genetic differentiation between populations.     

Genealogy of Neutral Genes and Spreading of Selected Mutations in a Geographically Structured Population

In a geographically structured population, the interplay among gene migration, genetic drift and natural selection raises intriguing evolutionary problems, but the rigorous mathematical treatment is often very difficult. Therefore several approximate formulas were developed concerning the coalescence process of neutral genes and the fixation process of selected mutations in an island model, and their accuracy was examined by computer simulation. When migration is limited, the coalescence (or divergence) time for sampled neutral genes can be described by the convolution of exponential functions, as in a panmictic population, but it is determined mainly by migration rate and the number of demes from which the sample is taken. This time can be much longer than that in a panmictic population with the same number of breeding individuals. For a selected mutation, the spreading over the entire population was formulated as a birth and death process, in which the fixation probability within a deme plays a key role. With limited amounts of migration, even advantageous mutations take a large number of generations to spread. Furthermore, it is likely that these mutations which are temporarily fixed in some demes may be swamped out again by non-mutant immigrants from other demes unless selection is strong enough. These results are potentially useful for testing quantitatively various hypotheses that have been proposed for the origin of modern human populations.     

Estimating the age of alleles by use of intraallelic variability.

A method is presented for estimating the age of an allele by use of its frequency and the extent of variation among different copies. The method uses the joint distribution of the number of copies in a population sample and the coalescence times of the intraallelic gene genealogy conditioned on the number of copies. The linear birth-death process is used to approximate the dynamics of a rare allele in a finite population. A maximum-likelihood estimate of the age of the allele is obtained by Monte Carlo integration over the coalescence times. The method is applied to two alleles at the cystic fibrosis (CFTR) locus, deltaF508 and G542X, for which intraallelic variability at three intronic microsatellite loci has been examined. Our results indicate that G542X is somewhat older than deltaF508. Although absolute estimates depend on the mutation rates at the microsatellite loci, our results support the hypothesis that deltaF508 arose < 500 generations (approximately 10,000 years) ago.     

An application of the central limit theorem to coalescence times in the structured coalescent model with strong migration

The structured coalescent describes the ancestral relationship among sampled genes from a geographically structured population. The aim of this article is to apply the central limit theorem to functionals of the migration process to study coalescence times and population structure. An application of the law of large numbers to the migration process leads to the strong migration limit for the distributions of coalescence times. The central limit theorem enables us to obtain approximate distributions of coalescence times for strong migration. We show that approximate distributions depend on the population structure. If migration is conservative and strong, we can define a kind of effective population size N _e ^*, with which the entire population approximately behaves like a panmictic population. On the other hand, the approximate distributions for nonconservative migration are qualitatively different from those for conservative migration. And the entire population behaves unlike a panmictic population even though migration is strong.     

Estimating effective population size from samples of sequences: a bootstrap Monte Carlo integration method.

We would like to use maximum likelihood to estimate parameters such as the effective population size N(e) or, if we do not know mutation rates, the product 4N(e) mu of mutation rate per site and effective population size. To compute the likelihood for a sample of unrecombined nucleotide sequences taken from a random-mating population it is necessary to sum over all genealogies that could have led to the sequences, computing for each one the probability that it would have yielded the sequences, and weighting each one by its prior probability. The genealogies vary in tree topology and in branch lengths. Although the likelihood and the prior are straightforward to compute, the summation over all genealogies seems at first sight hopelessly difficult. This paper reports that it is possible to carry out a Monte Carlo integration to evaluate the likelihoods approximately. The method uses bootstrap sampling of sites to create data sets for each of which a maximum likelihood tree is estimated. The resulting trees are assumed to be sampled from a distribution whose height is proportional to the likelihood surface for the full data. That it will be so is dependent on a theorem which is not proven, but seems likely to be true if the sequences are not short. One can use the resulting estimated likelihood curve to make a maximum likelihood estimate of the parameter of interest, N(e) or of 4N(e) mu. The method requires at least 100 times the computational effort required for estimation of a phylogeny by maximum likelihood, but is practical on today's work stations. The method does not at present have any way of dealing with recombination.     

Allelic genealogies in sporophytic self-incompatibility systems in plants.

Expectations for the time scale and structure of allelic genealogies in finite populations are formed under three models of sporophytic self-incompatibility. The models differ in the dominance interactions among the alleles that determine the self-incompatibility phenotype: In the SSIcod model, alleles act codominantly in both pollen and style, in the SSIdom model, alleles form a dominance hierarchy, and in SSIdomcod, alleles are codominant in the style and show a dominance hierarchy in the pollen. Coalescence times of alleles rarely differ more than threefold from those under gametophytic self-incompatibility, and transspecific polymorphism is therefore expected to be equally common. The previously reported directional turnover process of alleles in the SSIdomcod model results in coalescence times lower and substitution rates higher than those in the other models. The SSIdom model assumes strong asymmetries in allelic action, and the most recessive extant allele is likely to be the most recent common ancestor. Despite these asymmetries, the expected shape of the allele genealogies does not deviate markedly from the shape of a neutral gene genealogy. The application of the results to sequence surveys of alleles, including interspecific comparisons, is discussed.     

Effective population size and population subdivision in demographically structured populations

A fast-timescale approximation is applied to the coalescent process in a single population, which is demographically structured by sex and/or age. This provides a general expression for the probability that a pair of alleles sampled from the population coalesce in the previous time interval. The effective population size is defined as the reciprocal of twice the product of generation time and the coalescence probability. Biologically explicit formulas for effective population size with discrete generations and separate sexes are derived for a variety of different modes of inheritance. The method is also applied to a nuclear gene in a population of partially self-fertilizing hermaphrodites. The effects of population subdivision on a demographically structured population are analyzed, using a matrix of net rates of movement of genes between different local populations. This involves weighting the migration probabilities of individuals of a given age/sex class by the contribution of this class to the leading left eigenvector of the matrix describing the movements of genes between age/sex classes. The effects of sex-specific migration and nonrandom distributions of offspring number on levels of genetic variability and among-population differentiation are described for different modes of inheritance in an island model. Data on DNA sequence variability in human and plant populations are discussed in the light of the results.     

Gene and Allelic Genealogies at a Gametophytic Self-Incompatibility Locus

The properties of gene and allelic genealogies at a gametophytic self-incompatibility locus in plants have been investigated analytically and checked against extensive numerical simulations. It is found that, as with overdominant loci, there are two genealogical processes with markedly different time scales. First, functionally distinct allelic lines diverge on an extremely long time scale which is inversely related to the mutation rate to new alleles. These alleles show a genealogical structure which is similar, after an appropriate rescaling of time, to that described by the coalescent process for genes at a neutral locus. Second, gene copies sampled within the same functional allelic line show genealogical relationships similar to neutral gene genealogies but on a much shorter time scale, which is on the same order of magnitude as the harmonic mean of the number of gene copies within an allelic line. These results are discussed in relation to data showing trans-specific polymorphisms for alleles at the gametophytic self-incompatibility locus in the Solanaceae. It is shown that population sizes on the order of 4 X 10(5) and a mutation rate per locus per generation as high as 10(-6) could account for estimated allelic divergence times in this family.     

A Continuous Method for Gene Flow

Most modern population genetics inference methods are based on the coalescence framework. Methods that allow estimating parameters of structured populations commonly insert migration events into the genealogies. For these methods the calculation of the coalescence probability density of a genealogy requires a product over all time periods between events. Data sets that contain populations with high rates of gene flow among them require an enormous number of calculations. A new method, transition probability-structured coalescence (TPSC), replaces the discrete migration events with probability statements. Because the speed of calculation is independent of the amount of gene flow, this method allows calculating the coalescence densities efficiently. The current implementation of TPSC uses an approximation simplifying the interaction among lineages. Simulations and coverage comparisons of TPSC vs. MIGRATE show that TPSC allows estimation of high migration rates more precisely, but because of the approximation the estimation of low migration rates is biased. The implementation of TPSC into programs that calculate quantities on phylogenetic tree structures is straightforward, so the TPSC approach will facilitate more general inferences in many computer programs.     

Phylogenetic discordance at the species boundary: comparative gene genealogies among rapidly radiating Heliconius butterflies

Recent adaptive radiations provide excellent model systems for understanding speciation, but rapid diversification can cause problems for phylogenetic inference. Here we use gene genealogies to investigate the phylogeny of recent speciation in the heliconiine butterflies. We sequenced three gene regions, intron 3 ( approximately 550 bp) of sex-linked triose-phosphate isomerase (Tpi), intron 3 ( approximately 450 bp) of autosomal mannose-phosphate isomerase (Mpi), and 1,603 bp of mitochondrial cytochrome oxidase subunits I and II (COI and COII), for 37 individuals from 25 species of Heliconius and related genera. The nuclear intron sequences evolved at rates similar to those of mitochondrial coding sequences, but the phylogenetic utility of introns was restricted to closely related geographic populations and species due to high levels of indel variation. For two sister species pairs, Heliconius erato-Heliconius himera and Heliconius melpomene-Heliconius cydno, there was highly significant discordance between the three genes. At mtDNA and Tpi, the hypotheses of reciprocal monophyly and paraphyly of at least one species with respect to its sister could not be distinguished. In contrast alleles sampled from the third locus, Mpi, showed polyphyletic relationships between both species pairs. In all cases, recent coalescence of mtDNA lineages within species suggests that polyphyly of nuclear genes is not unexpected. In addition, very similar alleles were shared between melpomene and cydno, implying recent gene flow. Our finding of discordant genealogies between genes is consistent with models of adaptive speciation with ongoing gene flow and highlights the need for multiple locus comparisons to resolve phylogeny among closely related species.     

Genealogy of neutral genes in two partially isolated populations

Gene genealogy in two partially isolated populations which diverged at a given time t in the past and have since been exchanging individuals at a constant rate m is studied based upon an analytic method for large t and a simulation method for any t. Particular attention is paid to the conditions under which neutral genes sampled from populations are mono-, para-, and polyphyletic in terms of coalescence (divergence) times of genes. It is shown tha the probability of monophyly is high if M = 2Nm less than 0.5 and T = t/(2N) greater than 1, where N is the size of ancestral and descendant haploid populations, in which case most gene genealogies are likely to be concordant with the population relatedness. This probbility decreases as the sample size of genes increases. On the other hand, the case where the probability of monophyly is low will be either that of M greater than 1 and any T or that of M less than 1 and T less than 1, but the clear distinction between these conditions appears very difficult to make. These results are also examined if the gene genealogy is reconstructed from nucleotide differences. It is then shown that the results based upon coalescence times remain valid if the number of nucleotide differences between any pair of genes is not much smaller than 10. To observe such large nucleotide differences in small populations and therefore infer a reliable gene genealogy, we must examine a fairly long stretch of DNA sequences.     

On the Genealogy of a Rare Allele

The gene genealogy is derived for a rare allele that is descended from a mutant ancestor that arose at a fixed time in the past. Following Thompson (1976,Amer. J. Human Genet.28, 442–452), the fractional linear branching process is used as a model of the demography of a rare allele. The model does not require the total population size to be constant or the mutant class to be neutral; so long as individuals in the class are selectively equivalent, the class as a whole may have a selective advantage, or disadvantage, relative to other alleles in the population. An exact result is given for the joint probability distribution of the coalescence times among a sample of alleles descended from the mutant. A method is described for rapidly simulating these coalescence times. The relationship between the genealogical structure of a discrete generation branching process and a continuous generation birth–death process is elucidated. The theory may be applied to the problem of estimating the ages of rare nonrecurrent mutations.     

Impact of ancestral populations on postzygotic isolation in allopatric speciation

Postzygotic isolation evolves due to an accumulation of substitutions (potentially deleterious alleles in hybrids) in populations that have become geographically isolated. These potentially deleterious alleles might also be maintained in ancestral populations before geographic isolation. We used an individual-based model to examine the effect of the genetic state of an ancestral population on the evolution of postzygotic isolation after geographic isolation of a population. The results showed that the number of loci at which degenerative alleles are fixed in an ancestral population at equilibrium significantly affects the evolutionary rates of postzygotic isolation between descendant allopatric populations. Our results suggest that: (1) a severe decrease in population size (e.g., less than ten individuals) is not necessary for the rapid evolution of postzygotic isolation (e.g., <10,000 generation); (2) rapid speciation can occur when there is a large difference in the equilibrium number of accumulated degenerative alleles between ancestral and descendant populations; and (3) in an ancestral population maintained at a small effective population size for a long period of time, postzygotic isolation rarely evolves if back mutations that restore the function of degenerative alleles are limited.