Bayesian inference of species trees from multilocus data

Abstract

Until recently, it has been common practice for a phylogenetic analysis to use a single gene sequence from a single individual organism as a proxy for an entire species. With technological advances, it is now becoming more common to collect data sets containing multiple gene loci and multiple individuals per species. These data sets often reveal the need to directly model intraspecies polymorphism and incomplete lineage sorting in phylogenetic estimation procedures. For a single species, coalescent theory is widely used in contemporary population genetics to model intraspecific gene trees. Here, we present a Bayesian Markov chain Monte Carlo method for the multispecies coalescent. Our method coestimates multiple gene trees embedded in a shared species tree along with the effective population size of both extant and ancestral species. The inference is made possible by multilocus data from multiple individuals per species. Using a multiindividual data set and a series of simulations of rapid species radiations, we demonstrate the efficacy of our new method. These simulations give some insight into the behavior of the method as a function of sampled individuals, sampled loci, and sequence length. Finally, we compare our new method to both an existing method (BEST 2.2) with similar goals and the supermatrix (concatenation) method. We demonstrate that both BEST and our method have much better estimation accuracy for species tree topology than concatenation, and our method outperforms BEST in divergence time and population size estimation.

FIG. 1.

Species tree visualization. One locus for three individuals from each of the three species giving a total of nine samples. Current population size (t = 0) of A is 2 and at time 1.5 (where it split from B) the population size is 1.

FIG. 2.

(a) A simulated species tree from a birth–death process with continuous population sizes. (b) A single gene tree embedded inside the same species tree.

FIG. 3.

(a) Species tree estimation error and (b) 95% credible interval size as a function of the number of loci. The number of individuals sampled per species is four for all experiments. Each graph point is obtained by averaging the error measure (described in the main text) over 100 analyses of simulated data sets. The “branch score” is a measure of the distance in tree space of the estimated species tree to the true tree, incorporating both topology and divergence times. The “tree score” is a measure of the distance between the estimated species tree and the true species tree incorporating information about the population size as well. For details of the tree metrics used, see main text.

FIG. 4.

(a) Relative error and (b) credible interval size for both population size and speciation time point estimates. The number of individuals sampled per species is four for all experiments. Each graph point is obtained by averaging over 100 analyses of simulated data sets (see main text for details).

FIG. 5.

Speciation time estimate error as a function of speciation time. Data are taken from the main 100 runs with two loci, four individuals per species, and the sequence length of 1,600 bp.

FIG. 6.

(a) Relative error and (b) credible interval sizes as a function of sequence length. The number of individuals sampled per species is four for all experiments. The number of independent loci are four for all experiments. Each graph point is obtained by averaging over 100 analyses of simulated data sets. The horizontal line represents the theoretical maximum when sequence length approaches infinity and is calculated by using the gene trees directly without error.

FIG. 7.

(a) Relative error and (b) credible interval sizes, as a function of number of individuals sample from each species. Each graph point is obtained by averaging over 100 analyses of simulated data sets. The analysis used the true gene trees to reduce the computational cost.

FIG. 8.

Phylogeny for seven groups of western pocket gophers (Geomyidae, Thomomys). The analysis is based on seven noncoding nuclear genes from 28 individuals. Clade posterior probability is indicated on the branch. (a) Analysis with no monophyly constraints and (b) analysis with ingroup monophyly enforced.

FIG. 9.

Western pocket gophers (Geomyidae, Thomomys) species tree with embedded gene trees, each in a different color. The species tree was generated using median estimates for the divergence times and population sizes. Note that this representation is for the purpose of visual inspection only, and any inferences should be made directly from the posterior data.