Application of a time-dependent coalescence process for inferring the history of population size changes from DNA sequence data

Using Genomic Data to Infer Historic Population Dynamics of Nonmodel Organisms

Genome sequence data are now being routinely obtained from many nonmodel organisms. These data contain a wealth of information about the demographic history of the populations from which they originate. Many sophisticated statistical inference procedures have been developed to infer the demographic history of populations from this type of genomic data. In this review, we discuss the different statistical methods available for inference of demography, providing an overview of the underlying theory and logic behind each approach. We also discuss the types of data required and the pros and cons of each method. We then discuss how these methods have been applied to a variety of nonmodel organisms. We conclude by presenting some recommendations for researchers looking to use genomic data to infer demographic history.

Time to the MRCA of a sample in a Wright-Fisher model with variable population size

Determining the expected distribution of the time to the most recent common ancestor of a sample of individuals may deliver important information about the genetic markers and evolution of the population. In this paper, we introduce a new recursive algorithm to calculate the distribution of the time to the most recent common ancestor of the sample from a population evolved by any conditional multinomial sampling model. The most important advantage of our method is that it can be applied to a sample of any size drawn from a population regardless of its size growth pattern. We also present a very efficient method to implement and store the genealogy tree of the population evolved by the Galton-Watson process. In the final section we present results applied to a simulated population with a single bottleneck event and to real populations of known size histories.

Recovering population parameters from a single gene genealogy: an unbiased estimator of the growth rate

We show that the number of lineages ancestral to a sample, as a function of time back into the past, which we call the number of lineages as a function of time (NLFT), is a nearly deterministic property of large-sample gene genealogies. We obtain analytic expressions for the NLFT for both constant-sized and exponentially growing populations. The low level of stochastic variation associated with the NLFT of a large sample suggests using the NLFT to make estimates of population parameters. Based on this, we develop a new computational method of inferring the size and growth rate of a population from a large sample of DNA sequences at a single locus. We apply our method first to a sample of 1,212 mitochondrial DNA (mtDNA) sequences from China, confirming a pattern of recent population growth previously identified using other techniques, but with much smaller confidence intervals for past population sizes due to the low variation of the NLFT. We further analyze a set of 63 mtDNA sequences from blue whales (BWs), concluding that the population grew in the past. This calls for reevaluation of previous studies that were based on the assumption that the BW population was fixed.

Estimators of the human effective sex ratio detect sex biases on different timescales

Determining historical sex ratios throughout human evolution can provide insight into patterns of genomic variation, the structure and composition of ancient populations, and the cultural factors that influence the sex ratio (e.g., sex-specific migration rates). Although numerous studies have suggested that unequal sex ratios have existed in human evolutionary history, a coherent picture of sex-biased processes has yet to emerge. For example, two recent studies compared human X chromosome to autosomal variation to make inferences about historical sex ratios but reached seemingly contradictory conclusions, with one study finding evidence for a male bias and the other study identifying a female bias. Here, we show that a large part of this discrepancy can be explained by methodological differences. Specifically, through reanalysis of empirical data, derivation of explicit analytical formulae, and extensive simulations we demonstrate that two estimators of the effective sex ratio based on population structure and nucleotide diversity preferentially detect biases that have occurred on different timescales. Our results clarify apparently contradictory evidence on the role of sex-biased processes in human evolutionary history and show that extant patterns of human genomic variation are consistent with both a recent male bias and an earlier, persistent female bias.

Alternatives to the Wright-Fisher model: the robustness of mitochondrial Eve dating

Methods of calculating the distributions of the time to coalescence depend on the underlying model of population demography. In particular, the models assuming deterministic evolution of population size may not be applicable to populations evolving stochastically. Therefore the study of coalescence models involving stochastic demography is important for applications. One interesting approach which includes stochasticity is the O'Connell limit theory of genealogy in branching processes. Our paper explores how many generations are needed for the limiting distributions of O'Connell to become adequate approximations of exact distributions. We perform extensive simulations of slightly supercritical branching processes and compare the results to the O'Connell limits. Coalescent computations under the Wright-Fisher model are compared with limiting O'Connell results and with full genealogy-based predictions. These results are used to estimate the age of the so-called mitochondrial Eve, i.e., the root of the mitochondrial polymorphisms of the modern humans based on the DNA from humans and Neanderthal fossils.

Complex population genetic and demographic history of the Salangid, Neosalanx taihuensis, based on cytochrome b sequences

BACKGROUND

The Salangid icefish Neosalanx taihuensis (Salangidae) is an economically important fish, which is endemic to China, restricted to large freshwater systems (e.g. lakes, large rivers and estuaries) and typically exhibit low vagility. The continuous distribution ranges from the temperate region of the Huai and Yellow River basins to the subtropical region of the Pearl River basin. This wide ranging distribution makes the species an ideal model for the study of palaeoclimatic effects on population genetic structure and phylogeography. Here, we aim to analyze population genetic differentiation within and between river basins and demographic history in order to understand how this species responded to severe climatic oscillations, decline of the sea levels during the Pleistocene ice ages and tectonic activity.

RESULTS

We obtained the complete mtDNA cytochrome b sequences (1141 bp) of 354 individuals from 13 populations in the Pearl River, the Yangze River and the Huai River basin. Thirty-six haplotypes were detected. Haplotype frequency distributions were strongly skewed, with most haplotypes (n = 24) represented only in single samples each and thus restricted to a single population. The most common haplotype (H36) was found in 49.15% of all individuals. Analysis of molecular variance (AMOVA) revealed a random pattern in the distribution of genetic diversity, which is inconsistent with contemporary hydrological structure. Significant levels of genetic subdivision were detected among populations within basins rather than between the three basins. Demographic analysis revealed that the population size in the Pearl River basin has remained relatively constant whereas the populations in the Yangze River and the Huai River basins expanded about 221 and 190 kyr ago, respectively, with the majority of mutations occurring after the last glacial maximum (LGM).

CONCLUSION

The observed complex genetic pattern of N. taihuensis is coherent with a scenario of multiple unrelated founding events by long-distance colonization and dispersal combined with contiguous population expansion and locally restricted gene flow. We also found that this species was likely severely impacted by past glaciations. More favourable climate and the formation of large suitable habitations together facilitated population expansion after the late Quaternary (especially the LGM). We proposed that all populations should be managed and conserved separately, especially for habitat protection.

A three-stage colonization model for the peopling of the Americas

Background

We evaluate the process by which the Americas were originally colonized and propose a three-stage model that integrates current genetic, archaeological, geological, and paleoecological data. Specifically, we analyze mitochondrial and nuclear genetic data by using complementary coalescent models of demographic history and incorporating non-genetic data to enhance the anthropological relevance of the analysis.

Methodology/Findings

Bayesian skyline plots, which provide dynamic representations of population size changes over time, indicate that Amerinds went through two stages of growth ≈40,000 and ≈15,000 years ago separated by a long period of population stability. Isolation-with-migration coalescent analyses, which utilize data from sister populations to estimate a divergence date and founder population sizes, suggest an Amerind population expansion starting ≈15,000 years ago.

Conclusions/Significance

These results support a model for the peopling of the New World in which Amerind ancestors diverged from the Asian gene pool prior to 40,000 years ago and experienced a gradual population expansion as they moved into Beringia. After a long period of little change in population size in greater Beringia, Amerinds rapidly expanded into the Americas ≈15,000 years ago either through an interior ice-free corridor or along the coast. This rapid colonization of the New World was achieved by a founder group with an effective population size of ≈1,000–5,400 individuals. Our model presents a detailed scenario for the timing and scale of the initial migration to the Americas, substantially refines the estimate of New World founders, and provides a unified theory for testing with future datasets and analytic methods.

Population size changes reshape genomic patterns of diversity

Elucidating the forces responsible for genomic variation is critical for understanding evolution. Under standard conditions, X-linked diversity is expected to be three-quarters the level of autosomal diversity. Empirical data often deviate from this prediction, but the reasons for these departures are unclear. We demonstrate that population size changes can greatly alter relative levels of X-linked and autosomal variation: population size reductions lead to particularly low X-linked diversity, whereas growth elevates X-linked relative to autosomal diversity. Genetic variation from a diverse array of taxa supports an important role for this effect in accounting for population differences in the ratio of X-linked to autosomal diversity. Consideration of this effect may improve the inference of population history and other evolutionary processes.

Historical biogeography of Mediterranean trout

Complete sequencing of the mitochondrial control region was used to describe phylogenetic relationships of brown trout populations (Salmo trutta) in the Mediterranean river basins of Iberia and to review the historical biogeography of trout from the Mediterranean regions. Phylogenetic relationships among trout lineages suggested that the Danubian one is the most ancestral, in accordance with the eastern origin of most of the European freshwater fish species. Nested-clade and mismatch analyses suggested that the present distribution of haplotypes of the Adriatic and Mediterranean lineages resulted from population expansions originated, respectively, from central and western Europe, which favoured extensive secondary contacts between lineages. Reduced diversity detected within 50% of the analysed populations and large intrabasin differentiation indicated restricted gene flow in post-glacial periods.

Inference of demographic history from genealogical trees using reversible jump Markov chain Monte Carlo

BACKGROUND

Coalescent theory is a general framework to model genetic variation in a population. Specifically, it allows inference about population parameters from sampled DNA sequences. However, most currently employed variants of coalescent theory only consider very simple demographic scenarios of population size changes, such as exponential growth.

RESULTS

Here we develop a coalescent approach that allows Bayesian non-parametric estimation of the demographic history using genealogies reconstructed from sampled DNA sequences. In this framework inference and model selection is done using reversible jump Markov chain Monte Carlo (MCMC). This method is computationally efficient and overcomes the limitations of related non-parametric approaches such as the skyline plot. We validate the approach using simulated data. Subsequently, we reanalyze HIV-1 sequence data from Central Africa and Hepatitis C virus (HCV) data from Egypt.

CONCLUSIONS

The new method provides a Bayesian procedure for non-parametric estimation of the demographic history. By construction it additionally provides confidence limits and may be used jointly with other MCMC-based coalescent approaches.

Geographic genetic architecture of pocket gopher (Thomomys bottae) populations in Baja California, Mexico

Phylogenetic analyses of complete mitochondrial cytochrome b sequences support the monophyly of pocket gopher (Thomomys bottae) populations from the 1000 km length of the Baja California peninsula of Mexico, relative to other geographical segments of the species range in western North America. The Baja California peninsula is an area that encompasses considerable ecomorphological and infraspecific diversity within this pocket gopher species. However, detailed population analyses encompassing 35 localities distributed over the southern half of the peninsula reveal only trivial phylogeographical structure. Rather, most of the 72 unique 500-base pair haplotypes examined from 142 individuals is restricted to single populations, although a few haplotypes are shared broadly across geography. Individual populations are typically comprised of haplotype sets from different branches in a network of relationships. Analysis of molecular variance (amova) indicates that approximately half of the total pool of variation is contained among individuals within local populations, and that only about 25% can be explained by the regional subdivisions of current subspecies distributions or physiographic realms. A hypothesized historical vicariant event that has been causally linked to the phylogeographical structure of other, codistributed species has had little influence on these pocket gopher populations, explaining only 13% of the total variation. The temporal depth, estimated by coalescence parameters, of the haplotype lineage in Baja California is relatively recent, approximately 300,000 generations; both the mismatch distribution of pairwise comparisons and a significantly positive exponential growth estimate support a recent history of expanding populations; but current, or recent past, migration estimates have remained small, are largely unidirectional from north to south, and weak isolation by distance is present. All data suggest that pocket gophers have relatively recently invaded the southern half of peninsular Baja California, with the genetic signature of expansion still evident but with sufficient time having lapsed to result in a weak isolation by distance pattern. The geographical assemblage of sampled populations thus appears as a meta-population, with limited gene flow contrasting with random haplotype loss due to drift in small, localized populations.

New Explicit Expressions for Relative Frequencies of Single-Nucleotide Polymorphisms With Application to Statistical Inference on Population Growth

We present new methodology for calculating sampling distributions of single-nucleotide polymorphism (SNP) frequencies in populations with time-varying size. Our approach is based on deriving analytical expressions for frequencies of SNPs. Analytical expressions allow for computations that are faster and more accurate than Monte Carlo simulations. In contrast to other articles showing analytical formulas for frequencies of SNPs, we derive expressions that contain coefficients that do not explode when the genealogy size increases. We also provide analytical formulas to describe the way in which the ascertainment procedure modifies SNP distributions. Using our methods, we study the power to test the hypothesis of exponential population expansion vs. the hypothesis of evolution with constant population size. We also analyze some of the available SNP data and we compare our results of demographic parameters estimation to those obtained in previous studies in population genetics. The analyzed data seem consistent with the hypothesis of past population growth of modern humans. The analysis of the data also shows a very strong sensitivity of estimated demographic parameters to changes of the model of the ascertainment procedure.

The demographic history of the New Zealand short-tailed bat Mystacina tuberculata inferred from modified control region sequences

Short-tailed bats Mystacina tuberculata were widespread throughout the forest that dominated prehuman New Zealand, but extensive deforestation has restricted them to scattered populations in forest fragments. In a previous study, the species' intraspecific phylogeny was investigated using multiple mitochondrial gene sequences. Six phylogroups were identified with estimated divergences of 0.93-0.68 Ma. In the current study, the phylogeographical structure and demographic history of the phylogroups were investigated using control region sequences modified by removing homoplasic sites. Phylogeographical structure in the North Island was generally consistent with an isolation-by-distance dispersal model. Coalescent-based analyses (i.e. mismatch distributions, skyline plots, lineage dispersal analysis and nested clade analysis) indicated that the three phylogroups found in central and southern North Island expanded before the last glacial maximum, presumably during interstadials when Nothofagus forest was most extensive. Genetic structure within a central North Island hybrid zone was consistent with range expansion from separate refugia following reforestation after catastrophic volcanic eruptions. Phylogeographical structure in the South Island was consistent with southern populations originating during rapid southward range expansion from refugia in northern South Island following postglacial reforestation of the South Island 10-9 kya.

Demographic expansion of parasitic nematodes of livestock based on mitochondrial DNA regions that conflict with the infinite-sites model

Mitochondrial ND4 sequences of populations of four species of parasitic nematodes of livestock were subjected to demographic analyses. Deviation from selective neutrality was detectable using the frequency spectrum of segregating sites and highly negative neutrality statistics. However, the mitochondrial data sets do not comply with the infinite-sites model that underlies these tests, and as a consequence, it was not established whether these features are solely a result of population expansion, or whether aspects of the molecular evolution of these mitochondrial regions are also involved. Coalescent analyses based on Fu's Fs neutrality test, which incorporated estimates of rate heterogeneity, the transition-transversion ratio and nucleotide bias, as well as analyses that are fairly robust to deviations from the infinite-sites model supported population expansion. Also analyses that do not depend on the infinite-sites model suggested historical population expansion of these nematodes. The very similar time since expansion, the absence of signatures of positive selection in ND4 and the logical association with human demography imply that selective sweeps of mitochondrial variants are less probable, and that expansion is the most likely scenario for the parasitic nematodes of livestock. The methods used to characterize the expansion have different assumptions and emphasize different aspects of expansions. The resulting restrictions on the interpretation of expansions are outlined.

Historical demography and genetic structure of sister species: deermice (Peromyscus) in the North American temperate rain forest

To investigate the evolutionary and biogeographical history of Peromyscus keeni and P. maniculatus within the coastal forest ecosystem of the Pacific Northwest of North America, we sampled 128 individuals from 43 localities from southeastern Alaska through Oregon. We analysed mitochondrial DNA variation using DNA sequence data from the mitochondrial cytochrome-b (cyt-b) gene and control region, and we found two distinct clades consistent with the morphological designation of the two species. The sequence divergence between the two clades was 0.0484 substitutions per site for cyt-b and 0.0396 for the control region, suggesting that divergence of the two clades occurred during the middle to late Pleistocene. We also examined the historical demography of the two clades using stepwise and exponential expansion models, both of which indicated recent rapid population growth. Furthermore, using the program migrate we found evidence of migration from populations north of the Fraser River (British Columbia) to the south in both clades. This study demonstrates the utility of these model-based demographic methods in illuminating the evolutionary and biogegographic history of natural systems.

Demographic history of Diadema antillarum, a keystone herbivore on Caribbean reefs

The sea urchin Diadema antillarum was the most important herbivore on Caribbean reefs until 1983, when mass mortality reduced its populations by more than 97%. Knowledge of its past demography is essential to reconstruct reef ecology as it was before human impact, which has been implicated as having caused high pre-mortality Diadema abundance. To determine the history of its population size, we sequenced the ATPase 6 and 8 region of mitochondrial DNA from populations in the Caribbean and in the eastern Atlantic (which was not affected by the mass mortality), as well as from the eastern Pacific D. mexicanum. The Caribbean population harbours an order of magnitude more molecular diversity than those of the eastern Pacific or the eastern Atlantic and, despite the recent mass mortality, its DNA sequences bear the genetic signature of a previous population expansion. By estimating mutation rates from divergence between D. antillarum and D. mexicanum, that were separated at a known time by the Isthmus of Panama, and by using estimates of effective population size derived from mismatch distributions and a maximum likelihood coalescence algorithm, we date the expansion as having occurred no more recently than 100 000 years before the present. Thus, Diadema was abundant in the Caribbean long before humans could have affected ecological processes; the genetic data contain no evidence of a recent, anthropogenically caused, population increase.

Identifying main evolutionary mechanisms shaping genetic variation of Leuciscus cephalus L. 1758 (Cyprinidae) in western Greece: discordance between methods

Numerous methods can be used in intraspecific phylogeographic studies to infer the evolutionary mechanisms that shaped observed genetic variation in populations. However, these methods are scarcely used jointly, and the evolutionary outcomes they could propose are not fully compared. In this study, using a chub (Leuciscus cephalus; Cyprinidae) mitochondrial DNA data set (13 populations in Western Greece, 14 haplotypes), we compare three distinct 'historical' methods that could possibly infer relative importance of basic evolutionary mechanisms (isolation vs migration) shaping genetic variation: the nested clade analysis, the psi-test and the 'mismatch distributions'. Taking together, interpretations of these analyses allow to draw a picture of the evolutionary history of chub in Western Greece based on isolation and genetic drift for higher clades. However, results issued each method can differ for low differentiated clades. We discuss such differences and suggest that methods should be used jointly in phylogeographic studies for a better evaluation of the evolutionary mechanisms that shaped genetic variation.

A power analysis of microsatellite-based statistics for inferring past population growth

We present results concerning the power to detect past population growth using three microsatellite-based statistics available in the current literature: (1) that based on between-locus variability, (2) that based on the shape of allele size distribution, and (3) that based on the imbalance between variance and heterozygosity at a locus. The analysis is based on the single-step stepwise mutation model. The power of the statistics is evaluated for constant, as well as variable, mutation rates across loci. The latter case is important, since it is a standard procedure to pool data collected at a number of loci, and mutation rates at microsatellite loci are known to be different. Our analysis indicates that the statistic based on the imbalance between allele size variance and heterozygosity at a locus has the highest power for detection of population growth, particularly when mutation rates vary across loci.

An integrated framework for the inference of viral population history from reconstructed genealogies

We describe a unified set of methods for the inference of demographic history using genealogies reconstructed from gene sequence data. We introduce the skyline plot, a graphical, nonparametric estimate of demographic history. We discuss both maximum-likelihood parameter estimation and demographic hypothesis testing. Simulations are carried out to investigate the statistical properties of maximum-likelihood estimates of demographic parameters. The simulations reveal that (i) the performance of exponential growth model estimates is determined by a simple function of the true parameter values and (ii) under some conditions, estimates from reconstructed trees perform as well as estimates from perfect trees. We apply our methods to HIV-1 sequence data and find strong evidence that subtypes A and B have different demographic histories. We also provide the first (albeit tentative) genetic evidence for a recent decrease in the growth rate of subtype B.

Estimation of past demographic parameters from the distribution of pairwise differences when the mutation rates vary among sites: application to human mitochondrial DNA

Distributions of pairwise differences often called "mismatch distributions" have been extensively used to estimate the demographic parameters of past population expansions. However, these estimations relied on the assumption that all mutations occurring in the ancestry of a pair of genes lead to observable differences (the infinite-sites model). This mutation model may not be very realistic, especially in the case of the control region of mitochondrial DNA, where this methodology has been mostly applied. In this article, we show how to infer past demographic parameters by explicitly taking into account a finite-sites model with heterogeneity of mutation rates. We also propose an alternative way to derive confidence intervals around the estimated parameters, based on a bootstrap approach. By checking the validity of these confidence intervals by simulations, we find that only those associated with the timing of the expansion are approximately correctly estimated, while those around the population sizes are overly large. We also propose a test of the validity of the estimated demographic expansion scenario, whose proper behavior is verified by simulation. We illustrate our method with human mitochondrial DNA, where estimates of expansion times are found to be 10-20% larger when taking into account heterogeneity of mutation rates than under the infinite-sites model.

Using rare mutations to estimate population divergence times: a maximum likelihood approach

In this paper we propose a method to estimate by maximum likelihood the divergence time between two populations, specifically designed for the analysis of nonrecurrent rare mutations. Given the rapidly growing amount of data, rare disease mutations affecting humans seem the most suitable candidates for this method. The estimator RD, and its conditional version RDc, were derived, assuming that the population dynamics of rare alleles can be described by using a birth-death process approximation and that each mutation arose before the split of a common ancestral population into the two diverging populations. The RD estimator seems more suitable for large sample sizes and few alleles, whose age can be approximated, whereas the RDc estimator appears preferable when this is not the case. When applied to three cystic fibrosis mutations, the estimator RD could not exclude a very recent time of divergence among three Mediterranean populations. On the other hand, the divergence time between these populations and the Danish population was estimated to be, on the average, 4,500 or 15,000 years, assuming or not a selective advantage for cystic fibrosis carriers, respectively. Confidence intervals are large, however, and can probably be reduced only by analyzing more alleles or loci.