Efficient Bayesian inference under the multispecies coalescent with migration

Abstract

Analyses of genome sequence data have revealed pervasive interspecific gene flow and enriched our understanding of the role of gene flow in speciation and adaptation. Inference of gene flow using genomic data requires powerful statistical methods. Yet current likelihood-based methods involve heavy computation and are feasible for small datasets only. Here, we implement the multispecies-coalescent-with-migration model in the Bayesian program bpp, which can be used to test for gene flow and estimate migration rates, as well as species divergence times and population sizes. We develop Markov chain Monte Carlo algorithms for efficient sampling from the posterior, enabling the analysis of genome-scale datasets with thousands of loci. Implementation of both introgression and migration models in the same program allows us to test whether gene flow occurred continuously over time or in pulses. Analyses of genomic data from Anopheles mosquitoes demonstrate rich information in typical genomic datasets about the mode and rate of gene flow.

Keywords: BPP; gene flow; genomics; migration; multispecies coalescent.

Fig. 1.

(A) A species tree for three species ($A,B,C$) with migration between species $B$ and $C$ showing model parameters, $\mathrm{\Theta }$ = (${\tau }_R$, ${\tau }_S$, ${\theta }_A$, ${\theta }_B$, ${\theta }_C$, ${\theta }_R$, ${\theta }_S$, $M_{\mathit{BC}}$, $M_{\mathit{CB}}$). There are three sets of parameters in the model: species divergence times (${\tau }_R\equiv {\tau }_{\mathit{ABC}}$, ${\tau }_S\equiv {\tau }_{\mathit{AB}}$), population sizes (${\theta }_A$, ${\theta }_B$, ${\theta }_C$, ${\theta }_R$, ${\theta }_S$), and migration rates ($M_{\mathit{BC}},M_{\mathit{CB}}$). Both $\tau$ and $\theta$ are measured in the expected number of mutations per site. The (population) migration rate is defined as $M_{\mathit{ij}}=N_j{m}_{\mathit{ij}}$, the expected number of migrants from species $i$ to $j$ per generation, where $N_j$ is the (effective) population size of species $j$ and $m_{\mathit{ij}}$ is the proportion of immigrants in population $j$ from population $i$. (B) A possible gene tree with the complete history of coalescent and migration events at a locus with two sequences from each of the three species ($a_1,a_2$ from $A$; $b_1,b_2$ from $B$; and $c_1,c_2$ from $C$). In the backward-in-time process of coalescent and migration, the five coalescent events occur at times $t_1$–$t_5$, while sequence $b_2$ experienced two migration events from $B$ to $C$ at time $s_1$ and back from $C$ to $B$ at time $s_2$, with $t_1<t_2<s_1<s_2<t_3<t_4<t_5$.

Fig. 2.

A balanced tree for four species with migration from $S\to C$ at the rate $M_{\mathit{SC}}$. The migration rate parameter $M_{\mathit{SC}}$ exists when (A) ${\tau }_S<{\tau }_T$ but disappears when (B) ${\tau }_S\ge {\tau }_T$.

Fig. 3.

(A) The saturated migration model for three species with eight migration rates, used to simulate data for analysis by bpp. The parameters used are ${\tau }_R=0.02$, ${\tau }_S=0.01$, ${\theta }_A={\theta }_S=0.015$, and ${\theta }_B={\theta }_C={\theta }_R=0.025$. The eight migration rates are $M_{\mathit{AB}}=0.12$, $M_{\mathit{BA}}=0.21$, $M_{\mathit{AC}}=0.13$, $M_{\mathit{CA}}=0.31$, $M_{\mathit{BC}}=0.23$, $M_{\mathit{CB}}=0.32$, $M_{\mathit{CS}}=0.2$, and $M_{\mathit{SC}}=0.1$. (B) Posterior means and 95% HPD CIs of parameters in 10 replicate datasets of different sizes, with $L=250$, 1,000, 4,000 loci. Horizontal lines represent the true parameter values. A large dataset of $L=$ 16,000 loci is analyzed in SI Appendix, Table S1.

Fig. 4.

(A–C) Three MSC-M models for three species with (A) $A$–$B$ migration (between sister species), (B) $B$–$C$ migration (between nonsister species), and (C) $C$–$S$ migration (between sister species and involving one ancestor) used for simulating data, for analysis using bpp and IMA3. The parameters used are ${\tau }_R=0.02$, ${\tau }_S=0.01$, ${\theta }_A={\theta }_S=0.015$, ${\theta }_B={\theta }_C={\theta }_R=0.025$, with migration rates $M$ = 0.1 in one direction and 0.2 in the opposite direction. (D–F) Posterior means and 95% CIs for the nine parameters in the model obtained using bpp (blue) and IMA3 (red) in 10 replicate datasets. Black horizontal lines represent true values. IMA3 uses mutation rate per locus, not per site; estimates from IMA3 are thus divided by the sequence length. IMA3 outputs the population migration rate $2M$; the estimates are divided by 2.

Fig. 5.

(A) Stepping-stone and (B) island models used for simulating data, for analysis by using bpp and migrate. The population genetic models (Top) are approximated by the MSC-M models with very large divergence times (Bottom). All population sizes are ${\theta }_0=0.002$ except that ${\theta }_A=10{\theta }_0$ in the island model. The migration rate was $M=0.15$ migrants per generation. (C and D) Posterior means and 95% HPD CIs of parameters in bpp analyses of 10 replicate datasets (each of $L=250$, 1,000, or 4,000 loci) simulated under the models. The horizontal lines represent the true values. (E and F) The datasets of 250 loci are also analyzed using migrate (red), in comparison with bpp (blue). migrate uses the mutation-scaled migration rate, ${\varpi }_{\mathit{ij}}=4M_{\mathit{ij}}/{\theta }_j$ in the notation here; this is transformed to ${\widehat{M}}_{\mathit{ij}}={\widehat{\varpi }}_{\mathit{ij}}{\widehat{\theta }}_j/4$ by using the posterior mean ${\widehat{\theta }}_j$.

Fig. 6.

(A) MSC-M and (B) MSC-I models for six species of African mosquitoes in the A. gambiae species complex: A. gambiae (G), A. coluzzii (C), A. arabiensis (A), A. melas (L), A. merus (R), and A. quadriannulatus (Q).

Fig. 7.

The logarithm of the Bayes factor for testing gene flow obtained from bpp analysis of the 100-loci blocks of the (A) coding and (B) noncoding data from the Anopheles mosquitoes (Fig. 6A), calculated using thermodynamic integration with Gaussian quadrature (42). Gene flow is accounted for using either the introgression model (i for MSC-I) or the migration model (m for MSC-M). Model $H_0$ is the MSC model with no gene flow. Model $H_1$ assumes the $A\to GC$ gene flow, with rate $M_{A\to GC}$ under MSC-M or ${\phi }_{A\to GC}$ under MSC-I. Model $H_2$ accommodates both gene-flow events, with rates $M_{A\to GC}$ and $M_{R\to Q}$ under MSC-M or ${\phi }_{A\to GC}$ and ${\phi }_{R\to Q}$ under MSC-I. Bayes factor $B_{20}$ measures the support for $H_2$ over $H_0$, while $B_{21}$ measures the support for $H_2$ over $H_1$. The test is significant when $|logB|>4.6$ (i.e., if $B>100$ or $B<0.01$). For example, $log{B}_{21}$(m) $<-4.6$ means that the data strongly support the one-rate model with the $A\to GC$ migration (with rate $M_{A\to GC}$) over the two-rates model with both the $A\to GC$ and the $R\to Q$ migrations. Different scales are used for the $y$-axis over the intervals $(-50,-10),(-10,10),(10,50),(50,700)$.