Scaling up Bayesian population phylogenomics through virtual dimension reduction

Domain: neuroscience
Published: 2026-04-08
Editorial outcome: minor_revision

Abstract
Population phylogenomics uses sampled genomes to jointly infer population genetic processes (ancestral and contemporary population sizes, historical gene flow) and a phylogenetic tree relating species or populations including species split times. This challenging problem has been tackled most successfully in the Bayesian framework under the multispecies coalescent (MSC) model via Markov chain Monte Carlo (MCMC) computational algorithms. However, MCMC methods suffer from two serious problems: (i) mixing difficulties due to the high-dimensional state space with complex constraints, and (ii) the intrinsically serial nature of MCMC algorithms that defies parallelisation. To deal with both issues, we develop a new method, called Virtual Dimension Reduction allowing Parallelisation (VDRoP), that achieves the same MCMC mixing efficiency as dimension reduction through analytical integration of parameters, but without sacrificing parallel computation and without the restriction to conjugate priors. We implement the new method in the Bayesian program bpp and apply it to genomic datasets from Adansonia baobab trees, Anopheles mosquitoes, and Heliconius butterflies. The new algorithms reduce the run-time of MCMC analyses by 3 to 8 fold and improve the mixing efficiency by up to 50 fold for representative empirical datasets.

Introduction
In population phylogenomics, genomic samples are analysed to infer phylogenetic relationships among species or populations accounting for local population genetic processes (such as genetic drift and selection) and potential gene flow. This is a challenging problem of fundamental importance in understanding the evolution of all organisms. Methods for inferring key parameters such as population split times, rates of migration, and effective population sizes under the multispecies coalescent [1, 2] rely heavily on Bayesian Markov chain Monte Carlo (MCMC) algorithms [3, 4]. Those algorithms generate a Markov chain of samples from the posterior distribution of parameters, combining information from the genomic sequence data and prior distributions.

Two significant challenges have slow progress in scaling population phylogenomic methods to handle large genomes and complex patterns of gene flow between populations [5–7]. First, MCMC algorithms are difficult to parallelise; because each MCMC step can start only after the previous step finishes, the algorithm is intrinsically serial. Second, the massive size of genomic datasets and large number of highly correlated parameters makes it difficult to achieve good mixing in MCMC algorithms; poor mixing occurs when sampled parameter values in the Markov chain are highly correlated reducing the effective number of samples from the posterior distribution. Under the MSC model, gene trees (at potentially thousands of loci) are unobserved latent variables and are part of the state space for the Markov chain. They constitute a large highly complex space with both discrete structures (the gene tree topology) and continuous variables (the branch lengths or coalescent times), and places stringent constraints on the species tree and MSC model parameters.

The first challenge, parallelisation, can be addressed by exploiting conditional independence. In a Bayesian hierarchical model, conditional independence allows a subset of variables in the state space of the Markov chain to be updated in parallel when other variables are fixed. Most high-dimensional MCMC algorithms use multiple proposal steps, each updating a subset of parameters, while collectively the steps ensure that every dimension of the state space is updated in each MCMC iteration. While the MCMC iterations are sequential, proposals updating subsets of parameters within each MCMC iteration can be applied in parallel by relying on conditional independence.

In algorithms under the MSC model, gene trees at different loci are independent given the species tree and parameters and can thus be updated in parallel. Parallelising MCMC updates of gene trees is particularly important because these are computationally the most expensive proposals in the algorithm. Indeed under the MSC model (with and without gene flow), sequenced data at different loci are independently and identically distributed (i.i.d.), so that the number of loci is the sample size in the inference problem [1, 2]. Among many factors that affect the information content in the data (such as the number of loci, the number of sequences sampled per species, the sequence length, and the mutation rate), the number of loci is often the most important [e.g., 8–10]. For example, to estimate the rates of gene flow reliably, thousands of loci are often needed and it is important to develop algorithms that can accommodate many loci even if for a relatively small phylogeny.

The second challenge, poor mixing in high dimensions, can be addressed by analytically integrating some parameters out of the model, reducing the dimension of the parameter space. Analytical integration of population size parameters and migration rates have been used to reduce dimensionality and improve mixing in population phylogenomic models [5, 11, 12]. However, this approach has two serious limitations. First, it is tractable only under conjugate priors, but conjugate priors may not be the most appropriate and may cause inference problems as well as MCMC mixing problems. Second, integrating out parameters using conjugate priors often destroys conditional independence. For example, gene trees for different loci are independent given the parameters on the species tree, but they are not independent when those shared parameters are integrated out.

In this paper we propose a solution to the integration versus parallelisation trade-off that allows their simultaneous use in an MCMC analysis of a hierarchical Bayesian model reaping the benefits of both. The method, which we refer to as Virtual Dimension Reduction allowing Parallelisation (VDRoP), proposes certain parameters from their conditional distribution (achieving the same mixing efficiency as analytical integration) but allows the updates of gene trees to be executed in parallel. Here, efficiency is measured by the variance of an estimate from an MCMC sample relative to the variance based on an independent sample.

We solve the second problem of the integration approach, the need to use a conjugate prior, by using approximations to the conditional distribution for non-conjugate priors. We show that the approximation improves with increasing data size. Thus, the VDRoP algorithm removes the limitations of both the integration and parallelisation strategies for improving MCMC efficiency.

We describe MCMC algorithms under the MSC models emphasizing the conditional-independence structure of the hierarchical model. We implement new algorithms that leverage the conditional independence and exploit the use of (approximate) conditional distributions of parameters to design MCMC moves that achieve the same efficiency as analytical integration. We develop new methods to generate posterior summaries for parameters which are integrated out in the MCMC algorithm or whose conditionals are tractable or can be approximated. A simple example involving a bivariate Gaussian target is presented to illustrate the new algorithms. We then analyse three genomic datasets to assess the performance of the new algorithms in real-world applications.

Results
The multispecies coalescent model has been extended to account for gene flow between species [2]. Gene flow is modelled as either a discrete introgression/hybridization event which occurred at a time point in the past in the MSC-introgression (MSC-I) model [13–15] or as a continuous migration process in the MSC-migration (MSC-M) model [5, 7, 16]. Both models are implemented in bpp [7, 15], and have been applied to data of more than 10,000 loci (albeit for a small number of species), providing a powerful framework for testing gene flow and estimating its rate from genomic data [e.g., 18–20]. Previous studies have demonstrated the correctness and efficiency of our implementations [7, 15]. Our focus here is on improvement of computational and mixing efficiency of MCMC algorithms. Theories and algorithms developed in this study are implemented in bpp under the MSC, MSC-I (fig. 1a) and MSC-M (fig. 1b) models [21], detailed in Supplementary Note 1. Here we briefly describe our implementation of the MSC-M model, while noting the differences from the simpler MSC and MSC-I models.

The gene-tree density under the MSC-M model

We consider an MSC-M model for two species (fig. 1b) to introduce the parameters in the model and the calculation of the density for gene trees. The probability density of the gene tree at each locus under MSC-M is given by the backwards-in-time process of coalescent and migration [22–25] and has been described previously [e.g., 5, 12]. Here we follow eqs. 2 & 3 in ref. [7].

There are three types of parameters in the MSC-M model: species/population split times (𝝉), population sizes (𝜽), and the mutation-scaled migration rates (𝝕) (fig. 1b). For example, the MSC-M model for two species of figure 1b involves six parameters, Θ = (𝝉,𝜽,𝝕) = (𝜏,𝜃 𝐴 ,𝜃 𝐵 ,𝜃 𝐴𝐵 ,𝜛 𝐴𝐵 ,𝜛 𝐵𝐴). In analysis of sequence data, we measure time by mutations so that one time unit is the expected time to accumulate one mutation per site. Then both 𝜏 and 𝜃 are measured in the expected number of mutations per site. The mutation-scaled migration rate is then 𝜛 𝐴𝐵 = 𝑚 𝐴𝐵/𝜇, where 𝜇 is the mutation rate per site per generation, and 𝑚 𝐴𝐵 is the expected proportion of immigrants in the recipient population 𝐵 from the donor population 𝐴 in each generation. We use the real-world view with time running forward to define parameters (𝑚,𝜛).

Let the multi-locus sequence data be 𝑋 = (𝑋(𝑖)), where 𝑋(𝑖) denotes the sequence alignment at locus 𝑖. Let 𝑮 = (𝐺(𝑖)) be the gene trees, where 𝐺(𝑖) includes the rooted tree, the coalescent times, and the migration history at the locus (including the number, directions and timings of migration events). Given the species tree and parameters, the gene trees are assumed to be independent, while different sites in the sequence at the same locus are assumed to share the same genealogy. These assumptions hold if there is free recombination between loci and no recombination within a locus. Simulation suggests that inference under the MSC such as estimation of the rate of gene flow is robust to moderate levels of within-locus recombination, with rates at about ten times the human rate [10, 26, 27].

The MCMC algorithm samples from the joint conditional distribution of the parameters and the gene trees.

𝑝(Θ,𝑮|𝿿) ∝ 𝑝(Θ)𝑝(𝑮|Θ)𝑝(𝿿|𝑮), (1)

where 𝑝(Θ) is the prior, 𝑝(𝑮|Θ) is the probability density of the gene trees given the parameters in the MSC-M model, and 𝑝(𝿿|𝑮) is the probability of the sequence data given the gene trees or the phylogenetic likelihood [28]. The hierarchical model is illustrated in figure 2.

When we trace the history of the sample backwards in time in any population 𝑗 at locus 𝑖, two kinds of events may occur as competing Poisson processes. Coalescence between any pair of sequences occurs at the rate of 𝜃𝑗/2ℎ𝑖, where ℎ 𝑖 is the heredity (ploidy) scalar for locus 𝑖 (e.g., 1 for a diploid autosomal locus, 3 for an X locus, and 1/4 for a Y or mitochondrial locus). Migration for each sequence in population 𝑗 from population 𝑠 occurs at the rate of 𝜛 𝑠𝑗 (here an 𝑠 → 𝑗 migration means that a sequence in population 𝑗 is traced back to population 𝑠).

The gene-tree density, 𝑝(𝑮|Θ), is a product over populations and over loci. For each population 𝑗 and locus 𝑖 we break the time period into 𝐾 𝑖𝑗 time segments (defined by species divergence, coalescent and migration events) during which the coalescent and migration rates are constant (fig. 1b). Segment 𝑘 (𝑘 = 1,··· ,𝐾 𝑖𝑗) has time duration 𝑡 𝑖𝑗𝑘 and 𝑛 𝑖𝑗𝑘 lineages ancestral to the sample. Let I 𝑠𝑗𝑘 be an indicator which takes the value 1 if migration from 𝑠 to 𝑗 is possible during time segment 𝑘 (that is, if both species 𝑠 and 𝑗 exist during time segment 𝑘 and are permitted to exchange migrants in the model) and 0 otherwise. Let 𝑐 𝑗𝑖 be the number of coalescent events in population 𝑗 at locus 𝑖 and 𝑤 𝑠𝑗𝑖 be the number of migration events from population 𝑠 to population 𝑗 at locus 𝑖. The contribution to the gene tree density from population 𝑗 and locus 𝑖 is equal to the product of the Poisson rates for events that have occurred and the probability of no events during the time duration of the population, given by the variable-rate Poisson process [see, e.g. 29, p. 322].

𝑝(𝑮|Θ) = ∏︁ 𝑗 ∏︁ 𝑖 (cid:20) ∏︁ 𝑠 (cid:16) 2 𝜛 𝑠𝑗 ℎ 𝑖 𝜃 𝑗 (cid:17) 𝑐 𝑗𝑖 𝑤 𝑠𝑗𝑖 exp (cid:26) − ∑︁ 𝑘 𝐾 𝑖𝑗 (cid:18) 𝑛 𝑖𝑗𝑘 2 𝜃 𝑗 ℎ 𝑖 + ∑︁ 𝑠 I 𝑠𝑗𝑘 𝑛 𝑖𝑗𝑘 𝜛 𝑠𝑗 (cid:19) 𝑡 𝑖𝑗𝑘 (cid:27) (cid:21)

= ∏︁ 𝑗 ∏︁ 𝑖 (cid:20) (cid:16) 2 𝜃 𝑗 ℎ 𝑖 (cid:17) 𝑐 𝑗𝑖 exp (cid:26) − 𝐶 𝑗𝑖 𝜃 𝑗 ℎ 𝑖 (cid:27) × ∏︁ 𝑠 (cid:16) 𝜛 𝑠𝑗 (cid:17) 𝑤 𝑠𝑗𝑖 exp (cid:8) −𝑊 𝑠𝑗𝑖 𝜛 𝑠𝑗 (cid:9) (cid:21)

= ∏︁ 𝑗 (cid:20) (cid:16) 2 𝜃 𝑗 (cid:17) 𝑐 𝑗 exp (cid:26) − 𝐶∗ 𝑗 𝜃 𝑗 (cid:27) × ∏︁ 𝑠 (cid:16) 𝜛 𝑠𝑗 (cid:17) 𝑤 𝑠𝑗 exp (cid:8) −𝑊 𝑠𝑗 𝜛 𝑠𝑗 (cid:9) (cid:21) , (2)

where 𝐶 𝑗𝑖 = ∑︁ 𝑘 𝐾 𝑖𝑗 (cid:0)𝑛 𝑖𝑗𝑘 2 (cid:1)𝑡 𝑖𝑗𝑘 is the total coalescent waiting time in population 𝑗 at locus 𝑖, summed over time segments and lineage pairs, 𝑊 𝑠𝑗 = ∑︁ 𝑖 ∑︁ 𝑘 𝐾 𝑖𝑗 I 𝑠𝑗𝑘 𝑛 𝑖𝑗𝑘 𝑡 𝑖𝑗𝑘 is the total migration waiting time for 𝑠 → 𝑗 migration in population 𝑗, summed over time segments, lineages and loci, and 𝑤 𝑠𝑗 = ∑︁ 𝑖 𝑤 𝑠𝑗𝑖 is the total number of 𝑠 → 𝑗 migration events over all loci.

In ref. [7] we used the population migration rate, 𝑀 𝐴𝐵 = 𝑚 𝐴𝐵 𝑁 𝐵 = 𝜛 𝐴𝐵 𝑁 𝐵 𝜇 = 𝜛 𝐴𝐵 𝜃 𝐵/4, which is the expected number of 𝐴 → 𝐵 migrants per generation. Here we reparametrize the model using the mutation-scaled migration rate (𝜛) as this leads to simpler algorithms. The change of parametrization from 𝑀 to 𝜛 means a change to the prior on migration rates while the likelihood model stays the same. Also there is an error in eq. 2 of ref. [7]: the mutation-scaled migration rate from species 𝑠 to 𝑗 should be 𝜛 𝑠𝑗 instead of 4𝑀 𝑠𝑗. The error affects analyses in which loci with different heredity scalars are included in the dataset (e.g., a mixture of autosomal and X-linked loci), and has no effect when all loci have the same heredity scalar (for example when all loci are autosomal).

Let 𝐺 𝑗 be the part of the gene trees in species 𝑗 over all loci, with 𝑮 = (𝐺 𝑗). We rewrite the gene-tree density of eq. 2 as

𝑝(𝑮|Θ) = ∏︁ 𝑗 𝑝(𝐺 𝑗|Θ) = ∏︁ 𝑗 𝑝( 𝑗 𝑐)𝑝( 𝑗 𝑚), (3)

This is the same factorisation as in equation 2 of ref. [12] [see also 11]. The contribution to the gene-tree density from each population 𝑗 consists of two components. The coalescent component, 𝑝( 𝑗 𝑐), is the coalescent rates multiplied by the probability of no coalescent events over the total time period in the population when coalescence was possible, while the migration component, 𝑝( 𝑗 𝑚), is the migration rates multiplied by the probability of no migration events over the total time period available for migration in the population.

Inverse-gamma priors on 𝜽 are conjugate

The gene-tree density (eq. 3) has the form of an inverse-gamma density for any 𝜃 𝑗 for population 𝑗. We can use inverse-gamma priors to integrate out 𝜃s to reduce the space of the Markov chain [6, 12, 21, 30]. Previously Hey and Nielsen [11] integrated out 𝜃s by using uniform priors with user-specified upper bounds; see ref. [10] for a discussion of issues arising from the use of uniform priors in these models [5, ]. Note that from the inverse-gamma density

𝑔˜(𝑦|𝛼,𝛽) = 𝛽 𝛼 Γ(𝛼) 𝑦−𝛼−1e−𝛽/𝑦, (5)

we have

∫ ∞ 0 𝑦−𝛼−1e−𝛽/𝑦 d𝑦 = Γ(𝛼) 𝛽 𝛼 . (6)

Let 𝜃 𝑗 ∼ IG(𝛼 𝑐 ,𝛽 𝑐). If we integrate out 𝜃 𝑗 from the gene-tree density, the coalescent-component of the gene tree density becomes

𝑝¯ ( 𝑗 𝑐) = ∫ 𝑔˜(𝜃 𝑗|𝛼 𝑐 ,𝛽 𝑐)𝑝( 𝑗 𝑐) d𝜃 𝑗

= ∏︁ 𝑖 (cid:16) 2 ℎ 𝑖 (cid:17) 𝑐 𝑗𝑖 · 𝛽 𝑐 𝛼 𝑐 Γ(𝛼 𝑐+𝑐 𝑗) Γ(𝛼 𝑐) (𝛽 𝑐+𝐶∗ 𝑗 )𝛼 𝑐+𝑐 𝑗 (7)

where

𝑐 𝑗 = ∑︁ 𝑖 𝑐 𝑗𝑖 , (8)

𝐶∗ 𝑗 = ∑︁ 𝑖 ℎ 𝑖 𝐶 𝑗𝑖 .

In other words, when all population size parameters are integrated out, the gene-tree density becomes 𝑝(𝑮|𝝉,𝝕) = ∏︁ 𝑗 𝑝¯ ( 𝑗 𝑐)𝑝( 𝑗 𝑚) instead of eq. 3.

However, when population size parameters are integrated out, gene trees at different loci are not independent, and 𝑝¯ (𝑐) for population 𝑗 is not a product over loci (𝑖). Thus MCMC steps updating gene-tree node ages or gene-tree topologies (via subtree pruning and regrafting or SPR) [1] cannot be parallelised. Instead we use the conditional distributions of 𝜃s to propose new values (i.e., the Gibbs sampler) and modify MCMC steps to make use of the knowledge of the conditional.

Given the species tree (𝑆) and the gene trees (𝐺), the conditional distribution of 𝜃 𝑗 is inverse-gamma IG(𝛼(𝑐)∗ 𝑗 ,𝛽(𝑐)∗ 𝑗 ), with

𝛼(𝑐)∗ 𝑗 =𝛼 𝑐+𝑐 𝑗 , 𝛽(𝑐)∗ 𝑗 = 𝛽 𝑐+𝐶∗ 𝑗 . (9)

The Gibbs sampler can thus be used to update 𝜃 𝑗. The algorithm works also under models which assume that several populations share the same 𝜃; one simply counts the coalescent events (𝑐 𝑗) and calculates the total coalescent waiting time (𝐶∗) by summing over all populations that share the same 𝜃.

Gamma priors for 𝜽 are not conjugate but approximate conditionals can be used to achieve similar performance to integration

While the inverse-gamma prior on 𝜃 is mathematically convenient, it is heavy-tailed, which may cause MCMC mixing problems as very large 𝜃 values appear during the MCMC. Heterozygosity of 𝜃 > 10% is rare in extant species of plants and animals and is expected to be rare in extinct ancestral species as well. Thus the inverse-gamma prior may place probability in regions of the parameter space that are not supported by the data and may cause MCMC mixing problems. The light-tailed gamma distribution appears to be a better choice. However, the gamma prior is not conjugate, and the conditional distribution 𝑝(𝜃 𝑗|𝑋,𝜏,𝐺) is intractable, so that the Gibbs sampler is impossible.

We have developed gamma and inverse-gamma approximations to the conditional and implemented Metropolised Gibbs algorithms to update 𝜃s in bpp. We describe the approach in Supplementary Note 2 in the context of Bayesian inference under a Poisson process in which a gamma prior is assigned to the expected waiting time for the event. We evaluated the approximation under different scenarios (e.g., with the prior being consistent or conflicting with the data) and found that the approximation is good even with only 2 or 10 events (fig. S1). Both the gamma and inverse-gamma approximations converge to the true conditional as the data size increases (i.e., the total number of coalescent events increases), and the approximation is good when the prior is not too far from the data.

We also note an interesting difference between the gamma and inverse-gamma approximations when the Markov chain is far out in the tail (fig. S2). In such a case, the light tail of the gamma may lead to rejection of the proposal, causing convergence issues (fig. S2), and the inverse-gamma should be preferable. The situation is similar to rejection sampling or importance sampling, in which the sampling density should have a heavier tail than the target [e.g., 31, p. 60–63, p. 122–3]. When the chain has reached stationarity, both the gamma and inverse-gamma approximations are effective.

Gamma priors on migration rates (𝝕) are conjugate

The gene-tree density of eq. 3 has the form of the gamma density for any migration rate 𝜛 𝑠𝑗 from population 𝑠 to population 𝑗 [11, 12]. Note that from the gamma density

𝑔(𝑦|𝛼,𝛽) = 𝛽 𝛼 Γ(𝛼) 𝑦𝛼−1e−𝛽𝑦, 𝑦 > 0, (10)

we have

∫ ∞ 0 𝑦𝛼−1e−𝛽𝑦 d𝑦 = Γ(𝛼) 𝛽 𝛼 . (11)

Let 𝜛 𝑠𝑗 ∼ G(𝛼 𝑚 ,𝛽 𝑚). If we integrate out 𝜛 𝑠𝑗 from the gene-tree density (eq. 3), the migration component of the gene tree density becomes

𝑝¯ ( 𝑗 𝑚) = ∫ 𝑔(𝜛 𝑠𝑗|𝛼 𝑚 ,𝛽 𝑚)𝑝( 𝑗 𝑚) d𝜛 𝑠𝑗

= ∏︁ 𝑠 (cid:20) 𝛽 𝑚 𝛼 𝑚 Γ(𝛼 𝑚+𝑤 𝑠𝑗) Γ(𝛼 𝑚) (𝛽 𝑚+𝑊 𝑠𝑗)𝛼 𝑚+𝑤 𝑠𝑗 (cid:21) . (12)

In other words, when all migration rates are integrated out, the gene-tree density becomes 𝑝(𝑮|𝝉,𝜽) = ∏︁ 𝑗 𝑝( 𝑗 𝑐)𝑝¯ ( 𝑗 𝑚) instead of eq. 3.

Furthermore, if we assign the gamma prior 𝜛 𝑠𝑗 ∼ G(𝛼 𝑚 ,𝛽 𝑚), the conditional distribution of 𝜛 𝑠𝑗 will be gamma G(𝛼 𝑠 (𝑚 𝑗 )∗,𝛽 𝑠 (𝑚 𝑗 )∗), with

𝛼 𝑠 (𝑚 𝑗 )∗ =𝛼 𝑚+𝑤 𝑠𝑗 , 𝛽 𝑠 (𝑚 𝑗 )∗ = 𝛽 𝑚+𝑊 𝑠𝑗 . (13)

Thus one can sample 𝜛 𝑠𝑗 from its conditional distribution, given the species split times (𝜏s) and the gene trees.

Integrating out both 𝜽 and 𝝕

Under inverse-gamma priors on population sizes, 𝜃 ∼ IG(𝛼 𝑐 ,𝛽 𝑐), and gamma priors on migration rates, 𝜛 ∼ G(𝛼 𝑚 ,𝛽 𝑚), we see from eq. 2 that 𝜽 and 𝝕 are conditionally independent given the coalescent and migration histories in the gene trees. We can integrate out both 𝝕 and 𝜽 from eq. 3, to give the gene-tree density

𝑝(𝑮|𝝉) = ∏︁ 𝑗 𝑝¯ ( 𝑗 𝑐)𝑝¯ ( 𝑗 𝑚), (14)

with 𝑝¯ (𝑐) and 𝑝¯ (𝑚) given in eqs. 7 & 12.

To allow parallelisation, we keep 𝜽 and 𝝕 in the MCMC and instead sample from their conditional distributions given the species split times 𝝉 and the gene trees, 𝜃 𝑗 ∼ IG(𝛼∗ 𝑗 ,𝛽∗ 𝑗 ) and 𝜛 𝑠𝑗 ∼ G(𝛼 𝑠 ∗ 𝑗 ,𝛽∗ 𝑠𝑗 ). In the case of gamma priors on 𝜃s, either approximate gamma or inverse-gamma conditionals are used instead.

Improved MCMC algorithms

While our theory above concerning the treatment of the population size parameters (𝜃s) is developed under the MSC-M model, it applies to all MSC models including the simple MSC model with no gene flow and the discrete MSC-I model [15]. We note that in the MSC-I model the beta distribution is a conjugate prior for the introgression probability (𝜑), and 𝜑 can be treated similarly to the migration rate 𝜛 in the MSC-M model.

Modifications to the MCMC steps in bpp [1, 7, 15, 32] are described in Supplementary Note 1 (fig. 3). Here we provide an overview. If there are 𝐿 loci, with 𝑛 sequences per locus, there will be 𝐿 · (𝑛 − 1) coalescent times, while the number of migration times on gene trees is arbitrary depending on the migration rates. Step 2 applies gene-tree SPR or simulation to update gene tree topologies (and coalescent and migration times), looping over loci and for each locus over branches on the gene tree. With thousands of loci and dozens of sequences per locus, the gene trees constitute a complex space of huge dimensions. A major effort has been to parallelise those two steps which update gene trees at all loci.

Step 4 updates species split times (𝝉s) using the rubber-band proposal under the MSC, MSC-I, and MSC-M models, with coordinated deterministic changes to coalescent times and migration times on gene trees [1, 7]. Step 5 is a scaling or mixing move that multiplies all species split times on the species tree (𝝉) and coalescent and migration times on all gene trees by the same scale factor. In both steps, we sample 𝜽, 𝝋 (in MSC-I) and 𝝕 (in MSC-M) from their (approximate) conditionals given the newly proposed species split times and the gene trees.

Other steps are used to update 𝜽 and 𝝕 under MSC-M or 𝜽 and 𝝋 under MSC-I.

Simple example for a bivariate Gaussian target to illustrate the algorithms

In this section, we define the mixing efficiency (𝐸) for an MCMC algorithm and use a simple problem involving a bivariate Gaussian target to illustrate the new algorithms implemented in bpp and tested in the Results section using three empirical datasets (table S1).

Let (𝜙 1 ,𝜙 2 ,··· ,𝜙 𝑁) be a sample for any parameter 𝜙 from the MCMC, and we estimate the expectation of any function of the parameter, ℎ(𝜙), by the sample mean ℎ˜ = 1/𝑁 ∑︁ 𝑖 ℎ(𝜙 𝑖). The mixing efficiency of the MCMC is then defined as the variance ratio

𝐸 = lim 𝑁→∞ V(ℎ¯) V(ℎ˜) = 1 1+2(𝜌 1 +𝜌 2 +···) , (15)

where ℎ˜ is the estimate based on the MCMC sample and ℎ¯ is the estimate based on an independent sample of the same size, and where 𝜌 𝑘 = corr(ℎ(𝜙 𝑖),ℎ(𝜙 𝑖+𝑘)) is the lag-𝑘 autocorrelation [e.g., 33–35]. Thus efficiency is measured by the ratio of the variance for an estimator based on an independent sample to the variance based on the MCMC sample of the same size. The effective sample size (ESS) is given as 𝑁 · 𝐸. We also use the lag-1 autocorrelation 𝜌 1 , which is simpler to compute.

Besides the mixing efficiency 𝐸 (eq. 15), the computational cost for each MCMC iteration may differ among algorithms. However unless stated otherwise explicitly, the algorithms tested in this study have similar running time per iteration. Also the mixing efficiency may differ widely for different parameters in the model, and as a rule of thumb, the MCMC should be run long enough so that ESS exceeds 100 or 200 for every parameter. All runs in this study exceed this expectation.

Consider MCMC algorithms sampling from a two-dimensional target, 𝜋(𝑥,𝑦) = 𝜋(𝑥)𝜋(𝑦|𝑥), (16) to estimate 𝜇 𝑥; that is, the function ℎ(𝜇 𝑥) = 𝜇 𝑥 in eq. 15. We use the bivariate Gaussian (𝑁 2 ) target so that 𝜋(𝑥) = 𝜙(𝑥), 𝜋(𝑦|𝑥) = 𝜙(𝑦;𝜌𝑥,1−𝜌2), and

𝜋(𝑥,𝑦) = 𝜋(𝑥)𝜋(𝑦|𝑥) = 𝜙 2 (cid:16)(cid:2) 𝑦 𝑥(cid:3) ; (cid:2)0 0 (cid:3),(cid:2) 𝜌 1 𝜌 1 (cid:3)(cid:17) , (17)

where 𝜙(𝑥;𝜇,𝜎2) is the probability density function (PDF) for 𝑁(𝜇,𝜎2), with 𝜙(𝑥;0,1) ≡ 𝜙(𝑥), and where 𝜙 𝑘(𝒙;𝝁,𝚺) is the 𝑘-variate Gaussian density with the mean vector 𝝁 and variance-covariance matrix 𝚺. Note that for the algorithms discussed here, mixing efficiency and acceptance rate depend on the correlation 𝜌 but not the means or variances of the 𝑁 2 target (fig. 4).

In algorithm A1, we integrate out 𝑦, and sample from the marginal distribution 𝜋(𝑥) = 𝜙(𝑥) using a Gaussian sliding-window proposal, with 𝑞(𝑥′|𝑥) = 𝜙(𝑥′;𝑥,𝑠2). It is well-known that the optimal step length is 𝑠 = 2.5, at which efficiency (for estimating the target mean 𝜇 𝑥) is 𝐸 = 0.228 (acceptance rate 43%) [35, 36].

In algorithm A2, we generate 𝑥′ from the sliding window, 𝑞(𝑥′|𝑥) = 𝜙(𝑥′;𝑥,2.52) as in A1, and also sample 𝑦′ from its conditional, with 𝑞(𝑦′|𝑥′,𝑥) = 𝜋(𝑦′|𝑥′), with the (𝑥,𝑦) → (𝑥′,𝑦′) proposal accepted with probability

𝛼 =min 1, 𝑞(𝑥|𝑥′)𝑞(𝑦|𝑥,𝑥′) 𝜋(𝑥′,𝑦′) 𝑞(𝑥′|𝑥)𝑞(𝑦′|𝑥′,𝑥) 𝜋(𝑥,𝑦) (18)

=min 1, 𝑞(𝑥|𝑥′) 𝜋(𝑥′) 𝑞(𝑥′|𝑥) 𝜋(𝑥) . (19)

Efficiency 𝐸 (for estimating 𝜇 𝑥) and the acceptance rate are then independent of 𝜌 (fig. 4, A1 or A2), and indeed algorithms A1 and A2 are equivalent as far as sampling of 𝑥 is concerned.

Algorithms A1 and A2 here are analogous to algorithms A1 and A2 of table S1 implemented in bpp; the variable 𝑥 represents 𝜏s and gene trees in the MSC models while 𝑦 represents 𝜽, 𝝕 or 𝝋. As the inverse-gamma priors are conjugate for 𝜃s, we can integrate 𝜽 out (in A1) or sample them while updating the species split times and gene trees (in A2), and both algorithms have the same efficiency. As 𝜋(𝑦|𝑥) is intractable, we may use a proposal density that is close to the conditional, 𝑞(𝑦′|𝑥′,𝑥) = 𝑞(𝑦′|𝑥′) ≈ 𝜋(𝑦′|𝑥′). Then eq. 18 holds but eq. 19 does not. This mimics algorithms A4 and A5 under MSC-I, and B4 and B5 under MSC-M (table S1), in which the gamma priors on 𝜃s are not conjugate, and the conditionals for 𝜃 are approximated using the gamma (A4 and B4) or inverse-gamma (A5 and B5) distributions. If the approximation is good, the algorithm may achieve nearly the same efficiency and acceptance rate as if the exact conditionals were used.

In algorithm 3, we use two 1-D Gaussian sliding-window proposals, updating 𝑥 and 𝑦 separately. With the step lengths adjusted to have the acceptance rate of 40%, the efficiency (for estimating 𝜇 𝑥) decreases with the increase of correlation, dropping to 0.024 at 𝜌 = 0.9 and 0.012 at 𝜌 = 0.95 (fig. 4). This mimics algorithms A3, A6 and B6 of table S1, algorithms in earlier versions of bpp.

More efficient estimation of posterior means for parameters by use of conditionals

Population sizes (𝜃s), introgression probabilities (𝜑), and migration rates (𝜛) are of biological interest, so it is useful to generate their posterior distributions, in particular, posterior means and highest probability density (HPD) credibility intervals (CIs), even when they are integrated out in the MCMC to improve efficiency. Furthermore even if those parameters are sampled during the MCMC (as in the algorithms of this paper), their conditionals contain more information than the sample values and can be used to produce more precise posterior summaries than the MCMC sample.

Suppose we have an MCMC algorithm operating on 𝑋 with the variable 𝑌 integrated out. Let (𝑥 1 ,𝑥 2 ,··· ,𝑥 𝑁) be a sample from this algorithm. Given each 𝑥 𝑖,𝑖 = 1,··· ,𝑁, the conditional distribution of 𝑌 is available. Thus one can sample from the conditional distribution 𝑦 𝑖 ∼ 𝑝(𝑌|𝑥 𝑖) and use the sample mean 𝑦˜ = 1/𝑁 ∑︁ 𝑦 𝑖 to estimate 𝜇 𝑦. However, we can also calculate the conditional means, 𝑢 𝑖 = 𝐸(𝑌|𝑥 𝑖), and use the sample mean 𝑢˜ = 1/𝑁 ∑︁ 𝑢 𝑖 to estimate 𝜇 𝑦, with a reduced variance relative to 𝑦˜. The efficiency of 𝑢˜ relative to 𝑦˜ (or the ratio of the variances for estimating 𝜇 𝑦) is

𝐸 𝑢˜,𝑦˜ ≡ lim 𝑁→∞ V(𝑦˜) V(𝑢˜) = 1 1−(1−𝑐)𝐸 𝑦˜ >1, (20)

where 𝐸 𝑦˜ is the efficiency of 𝑦˜ relative to an independent sample, and where 𝑐 = 𝜎 𝑢 2 /𝜎 𝑦 2 is the ratio of the variance of the conditional mean 𝑢(𝑋) to the variance of 𝑌. The theory is developed in Supplementary Note 3. Note that here we are treating 𝜇 𝑦 as the expectation of a deterministic function of 𝑋: 𝜇 𝑦 = E 𝑋(𝑢) where 𝑢 = 𝑢(𝑋) = E(𝑌|𝑋) is a function of 𝑋.

In the statistics literature, the conditional expectation is commonly used to improve an estimator, a technique known as Rao-Blackwellisation [37]. For an estimate based on an independent sample the conditional expectation always lead to reduced variance, but this may not be true for a dependent sample [38, 39]. Here our use of the conditional, 𝑢 𝑖 = E(𝑌|𝑥 𝑖), is guaranteed to reduce the variance in the estimate, with V(𝑢˜) < V(𝑦˜). See Supplementary Note 3. Furthermore, the improvement in mixing efficiency requires virtually no extra computational cost after the MCMC sample (𝑥 𝑖) has been generated.

The improvement of 𝑢˜ over 𝑦˜, measured by 𝐸 𝑢˜,𝑦˜ , is a monotonically increasing function of 𝑐 and 𝐸 𝑦˜ (eq. 20). Here 𝑐 is determined by the inference problem or the posterior distribution, and it is an interesting question how 𝑐 is influenced by the model, the parameters and features of the data (such as the number of loci, the number of sequences per locus, &c.). In contrast, 𝐸 𝑦˜ is affected by the MCMC algorithm operating on 𝑋. In this regard, the more efficient 𝑦˜ already is, the greater will be the improvement achieved by switching to 𝑢˜ (eq. 20).

Consider for instance two algorithms for sampling from 𝑋 ∼ G(10,100) and 𝑌|𝑋 ∼ G(10,0.02+0.01𝑥). With the sliding-window algorithm, 𝐸 𝑦˜ = 0.014, while 𝐸 𝑢˜ = 0.58, with an improvement of 𝐸 𝑢˜,𝑦˜ = 1.35 folds (eq. 21). In the other algorithm [40], we use the 𝑦 sampled from 𝑝(𝑦|𝑥). Then 𝐸 𝑦˜ = 1.13, while 𝐸 𝑢˜ = 3.58, with an improvement of 𝐸 𝑢˜,𝑦˜ = 3.58 folds (eq. 21). See Supplementary Note 3 for details.

The theory also works for a function of the variables, such as E(𝑌 1 𝑌 2), where both 𝑌 1 and 𝑌 2 may be integrated out in the MCMC. This is useful for generating the posterior summaries for the population migration rate 𝑀 𝐴𝐵 = 𝜛 𝐴𝐵 𝜃 𝐵/4 in the MSC-M model, when both 𝜛 𝐴𝐵 and 𝜃 𝐵 are integrated out or sampled from their conditionals.

Besides the posterior mean, 𝜇 𝑦 = E(𝑌) or E(𝑌 1 𝑌 2), the whole posterior distribution, or the probability density function (PDF) for 𝑌 can be recovered from the conditionals sampled throughout the MCMC, 𝑝(𝑌|𝑥 𝑖). We also generated the 95% HPD CI as well as the 95% equal-tail CI for 𝑌. The detailed analyses and proofs are presented in Supplementary Note 3 (fig. S3, table S2).

Three datasets for evaluating the mixing efficiency of the new algorithms

We used three real datasets with different features to examine the mixing efficiency of our new algorithms. The first dataset consists of 344 target-enrichment loci from eight baobab species in the genus Adansonia (fig. 5a). Previous analyses of genetic data found considerable uncertainty in the species phylogeny [41, 42]. We inferred the species phylogeny under the MSC model using bpp [43, 44]. The maximum a posteriori (MAP) tree (fig. 5a) has a posterior of ∼ 100%, possibly because the full likelihood method uses the information in the data more efficiently than summary methods used in previous studies. This species tree differs from those inferred in ref. [41] using subsampled data and summary methods such as astral, and is similar to the maximum likelihood tree from plastid data [41, fig. 6]. The phylogeny, with the monophyly of the six Malagasy species and strong support for the two Malagasy sections (Longitubae and Brevitubae), appears to simplify the interpretation of the evolution of morphological characters such as flower colours and pollination [45]. We introduced potential gene-flow events on to the species phylogeny based on previous tests of ref. [41] using the D-statistic [46] and SNaQ [47]. Only one of the two events of figure 5a was supported in our analysis. Below we run MSC-I and MSC-M models with the 𝑥 → 𝑦 gene flow (fig. 5a).

The second dataset consists of 4,133 noncoding loci from chromosome 2L1 from six species of African mosquitoes in the Anopheles gambiae species complex: A. gambiae (G), A. coluzzii (C), A. arabiensis (A), A. melas (L), A. merus (R), and A. quadriannulatus (Q) [48, 49]. The species tree is shown in figure 6a&b, with two gene-flow events (A→GC and R→Q) accounted for using the MSC-I and MSC-M models.

The third dataset consists of 5341 noncoding loci from chromosome 1 for three species of Heliconius butterflies: Heliconius hecale (H), H. cydno (C), and H. Melpomene (M) [17, 50]. The species tree is shown in figure 7a&b, with C→M gene flow modeled using either the discrete MSC-I model or the continuous MSC-M model.

Mixing efficiency of different algorithms applied to empirical data

We analysed the three empirical datasets using six MCMC algorithms under the MSC-I model (A1–A6) and three algorithms under the MSC-M model (B4–B6) (table S1). The mixing efficiency 𝐸 (and 𝜌 1 ) were used to compare algorithms. The results are summarised in figures 5–7. Computational tests were conducted on a Lenovo ThinkSystem SR850 server with 4x Intel Xeon Gold 6154 18C processors. Apart from algorithm A1, all other algorithms involve similar amount of computation, so the improvements are mostly due to reduced variance or improved mixing efficiency. We discuss the results for the baobabs (fig. 5) in detail, and the results for the Anopheles and Heliconius (figs. 6 & 7) follow the same format.

First in algorithm A1 (I-IG-int), we integrate out 𝜽 using the conjugate inverse-gamma priors, while algorithm A2 (I-IG-gIG) under the same MSC-I model and same inverse-gamma priors uses a Gibbs sampler to sample 𝜃 from the inverse-gamma conditional. Here the key ‘I-IG-int’ means the MSC-I model, with inverse-gamma priors on 𝜃s, and with 𝜃s integrated out. In general, our key to algorithms is in the format ‘model-prior-algorithm’, with model to be ‘I’ for MSC-I and ‘M’ for MSC-M; with prior for 𝜃s to be ‘IG’ for inverse-gamma and ‘G’ for gamma (beta for 𝜑 in MSC-I and gamma for 𝜛 in MSC-M are always used); and with algorithm to be ‘slide’ for sliding-windows, ‘it’ for integrating out 𝜃s, ‘gIG’ for (Metropolised) Gibbs based on inverse-gamma conditional for 𝜃s, ‘gG’ for Metropolised Gibbs based on the gamma conditional for 𝜃s.

Figure 5b shows the posterior means and 95% HPD CIs for parameters produced by the nine algorithms. The posterior should be the same under the same model and prior, that is, among algorithms A1 (I-IG-int), A2 (I-IG-gIG), and A3 (I-IG-slide) for the MSC-I model under the inverse-gamma priors on 𝜃s; among algorithms A4 (I-G-gG), A5 (I-G-gIG), and A6 (I-G-slide) under the MSC-I model with gamma priors on 𝜃s; and among B4 (M-G-gG), B5 (M-G-gIG), B6 (M-G-slide) under the MSC-M model with gamma priors on 𝜃s. The gamma and inverse-gamma priors on 𝜃s produced virtually identical posteriors (fig. 5b), as the dataset is large. Also species split times &c. are very similar between the two models (MSC-I and MSC-M). In particular, the posterior mean (and 95% HPD CI) for the rate of gene flow is 𝜑 𝑥𝑦 = 0.428(0.352,0.503) under the MSC-I model, and 𝜛 𝑥𝑦 = 20.0 (4.0, 37.0) under MSC-M (fig. 5b, table S3).

Next we examine the mixing efficiency of different algorithms. Algorithm A2 (I-IG-gIG) is expected to have the same mixing efficiency as A1 (I-IG-int). This is reflected by the efficiency 𝐸 and autocorrelation 𝜌 1 (eq. 15) for the two algorithms falling on the diagonal line for all parameters in figure 5c (A2 = A1). Note that good mixing is indicated by high 𝐸 and low 𝜌 1 , and there is large variation in mixing efficiency among parameters. While 𝜽 are integrated out in A1, in A2 we keep 𝜽 and 𝝕 as part of the Markov-chain state and sample them from their conditional distributions given the species tree and the gene trees. This allows parallelisation of MCMC proposals that update the gene trees at all loci. Thus A2 has a large computational advantage on multi-core computers over A1, achieving a 8.2-fold reduction in running time (fig. 8).

The same idea works even under the gamma priors on 𝜃s, for which analytical integration over 𝜃s (as in algorithm A1) is impossible. By using approximate gamma or inverse-gamma conditionals, algorithms A4 (I-G-gG) and A5 (I-G-gIG) under the gamma priors on 𝜽 achieve similar efficiency as A2 (I-IG-gIG) for the inverse-gamma priors (fig. 5c, A4 ≈ A5, A5 ≈ A2). See also similar results in figure 5c, B4 ≈ B5 under the MSC-M model.

In algorithm A6 (I-G-slide) for the MSC-I model with gamma priors on 𝜽, parameters 𝜽 are updated using sliding windows, and no change of 𝜽 is made in the rubber-band or mixing (scaler) steps which updates the species split times with coordinated changes to the gene tree node ages. This is effectively the algorithm in the previous version of bpp (version 4.7). We expect the new algorithm A5 (I-G-gIG) to have higher mixing efficiency than the old algorithm A6 (I-G-slide), and this is indeed the case (fig. 5c, A5 ≫ A6). Similarly algorithm B5 (M-G-gIG) is much better than B6 (M-G-slide) under the MSC-M model (fig. 5c, B5 ≫ B6). Note that the gamma priors on 𝜽 are not conjugate, so it is impossible to integrate out 𝜽 analytically. Our use of the approximate conditionals thus make it possible to achieve the impossible.

Similarly in analyses of the Anopheles and Heliconius datasets, all algorithms produced the same posterior under the same model and priors (figs. 6 & 7), and the gamma and inverse-gamma priors on 𝜽 had minimal impact on the posteriors of parameters. For Anopheles (fig. 6d, table S4) the posterior means and 95% HPD CIs for the rate of gene flow are 𝜑 𝐴→𝐺𝐶 = 0.963(0.950,0.975) and 𝜑 𝑅→𝑄 = 0.016 (0.008,0.024) under the MSC-I model, and 𝜛 𝑅→𝑄 = 0.137 (0.007,0.291) and 𝜛 𝐴→𝐺𝐶 = 203(191,215) under MSC-M. For Heliconius the posterior means and 95% HPD CIs are 𝜑 𝑐𝑚 = 0.465 (0.430,0.500) under MSC-I and 𝜛 = 2.25 (0.05,5.38) under MSC-M (fig. 7d, table S5).

Similarly in the Anopheles and Heliconius datasets our expectations for the mixing efficiency of different algorithms are also confirmed (figs. 6 & 7). The mixing efficiency of algorithms under the MSC-I model is expected to be in the order: A2 (I-IG-gIG) = A1 (I-IG-int), A2 (I-IG-gIG) ≫ A3 (I-IG-slide), A2 (I-IG-gIG) ≈ A4 (I-G-gG) ≈ A5 (I-G-gIG), A5 (I-G-gIG) ≫ A6 (I-G-slide). For example A2 (I-IG-gIG) has the same mixing efficiency (𝐸 and 𝜌 1 ) as A1 (I-IG-int) (figs. 6d & 7d), but a large computational advantage over A1 on multithreaded computers, achieving a reduction in running time of 3.1× and 3.4× in the two datasets (table 1). Similarly under the MSC-M model, with gamma prior on 𝜽s, the mixing efficiency of the algorithms is expected to be in the order B4 (M-G-gG) ≈ B5 (M-G-gIG), and B5 (M-G-gIG) ≫ B6 (M-G-slide). These expectations are confirmed in all datasets (fig. 6d, and fig. 7d). Relative to version 4.7 of bpp [7, 15] which uses the sliding window moves to update 𝜽s and 𝑀 (i.e., algorithms A3 and A6 for MSC-I and B6 for MSC-M), the improvement in mixing efficiency of the new algorithms (A2, A4 and A5 for MSC-I and B4 and B5 for MSC-M) is up to 4-fold for baobabs (table S3), 51-fold for Anopheles (table S4) and 19-fold for Heliconius (table S5). Note that the improvement in mixing efficiency varies among parameters.

Methods
We use three empirical datasets to examine the mixing efficiency of MCMC algorithms under the MSC-I and MSC-M models.

Baobab species in the genus Adansonia

We analyse the dataset of 344 target-enrichment loci (i.e., putative single-copy nuclear genes) from eight baobab species in the genus Adansonia generated in ref. [41]. We removed the two outgroups Bombax ceiba and Pseudobombax croizatii as they are very distant [54] and kept the close outgroup Scleronema micrantha. Multiple samples for the same species are kept in the same alignment at each locus without subsampling. There are 28–38 sequences per locus, and the sequence length ranges over 759–9574 sites (median 1679). We note that Wan et al. [42] generated genome assemblies for the eight species of baobabs, but only one sample is sequenced per species so that the data may not be ideal for inferring gene flow. There is considerable uncertainty in the species phylogeny and gene flow. In particular, the relationships among A. digitata, A. gregorii, and the clade of six Malagasy species is highly uncertain, as is the position of A. rubrostipa (fig. 5a) [41, 42].

We conducted species tree estimation under the MSC model with no gene flow (the A01 analysis, 55) using different starting trees. The inferred species tree is in figure 5a. We then included gene-flow events on the phylogeny based on previous analyses [41]. Only one gene-flow event in figure 5a (𝑥 → 𝑦) was supported in our analysis. Thus we fit MSC-I and MSC-M models with the 𝑥 → 𝑦 gene flow (fig. 5a).

We used the gamma prior 𝜃 ∼ G(2,200) or inverse-gamma prior IG(3, 0.02) with mean 0.01. The age of the root was assigned the gamma prior 𝜏 𝑟 ∼ G(2,100) with mean 0.02. Under the MSC-I, we assigned the prior 𝜑 ∼ U(0,1) on the introgression probability. We used the linked-msci option so that a branch before and after an introgression event are assigned the same 𝜃. Under the MSC-M model, we assigned the gamma prior 𝜛 ∼ G(2,1) with mean 2. The MSC-I model was implemented using MCMC algorithms A1–A6, and the MSC-M model was implemented using algorithms B4–B6 (table S1). We used a burn-in of 4×104 iterations and took 5×105 samples, sampling every iteration.

African mosquitoes in the Anopheles gambiae group

We analyse a dataset of noncoding loci from chromosome arm L1 from six species of African mosquitoes in the Anopheles gambiae species complex: A. gambiae (G), A. coluzzii (C), A. arabiensis (A), A. melas (L), A. merus (R), and A. quadriannulatus (Q) (fig. 6a&b). The data consist of 4133 noncoding loci, with 12 sequences per locus (two sequences per species). The data consist of “haploid consensus sequences”, with heterozygotes resolved into the majority nucleotide [56]. The data were originally published and analysed in ref. [48], and recompiled by Thawornwattana et al. [49], who analysed coding and noncoding data from all chromosome arms (see also 7, 15). Here we used noncoding loci from L1 only.

We fitted the MSC-I and MSC-M models (fig. 6a&b), using the algorithms of table S1 with similar settings to the above. We used the gamma priors 𝜃 ∼ G(2,100) or IG(3,0.04) with prior mean 0.02, and 𝜏 𝑟 ∼ G(2,20) with prior mean 0.1. Under the MSC-I, we assigned the prior 𝜑 ∼ U(0,1) on the introgression probability. Under the MSC-M model, we assigned the gamma prior 𝜛 ∼ G(2,1). We used a burn-in of 4×104 iterations and took 2 × 105 samples, sampling every iteration. The results under both the MSC-I and MSC-M models match those published in table 1 of ref. [7]. Here we focus on the mixing or computational efficiency of the algorithms.

Heliconius butterflies

We analysed a dataset of 5341 noncoding loci from chromosome 1 from three species of Heliconius butterflies: Heliconius hecale (H), H. cydno (C), and H. Melpomene (M) (fig. 7a&b), from ref. [17]. The species tree is (H,(C,M)), with gene flow between C and M. Thawornwattana et al. ([17], tables 2 & 3) detected ongoing gene flow from C→M but rejected gene flow in the opposite direction. See tables S4–S6 and figure S8 in ref. [17] for analyses of both coding and noncoding data from all 21 chromosomes. Here we use the MSC-I and MSC-M models to account for the C → M gene flow to analyse the noncoding data from chromosome 1 only.

We use the gamma prior for the age of the root, 𝜏 𝑟 ∼ G(2,200). We used the priors G(2,200) and IG(3,0.04) for 𝜃, both with prior mean 0.01. For the rate of gene flow, we assigned the uniform prior 𝜑 ∼ U(0,1) under MSC-I and the gamma prior 𝜛 ∼ G(2,1) under MSC-M. We used a burn-in of 2×104 iterations and took 5×105 samples, sampling every iteration.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors describe an interesting approach to achieve improved performance for Bayesian inference under the\nmultispecies coalescent model, and are clearly experts in this matter. The manuscript is well structured and sufficiently easy\nto read through; I list my comments below.\nLine 31 : can this sentence be better connected to the previous one, i.e. explicitly state which of (i) and / or (ii) you’re tackling\nwith your new method? If you are tackling (i) but not (ii) then you should clearly mention this in the abstract (I don’t mean that\nin any bad way).\nLine 36: this leaves doubt as to whether mixing efficiency is expressed as effective sample size, or as effective sample size\nper time unit (the latter is the one that matters); please be explicit about this in the abstract.\nLine 42: please avoid using terms such as ‘Bayesian MCMC’, as there’s nothing Bayesian about MCMC (in this case: simply\nremove ‘Bayesian’).\nLine 63 (left column): this could do with a citation.\nLine 63 (left column): this does not change the serial nature of MCMC however, and you should point that out to the readers.\nWhat you are doing is more work in a single iteration / generation, which can indeed be achieved through parallelization and\n/ or more complex transition kernels. But neither break / avoid / bypass / change the serial nature of MCMC, and it would be\ngood to point this out to the readers to avoid giving the wrong impression.\nLine 91: ‘referred to’, do you mean ‘which you will refer to in this manuscript’? If not – and it’s a term that has seen use before\n– this needs a citation.\nLine 91 (right column): please number your figures in the order they are referred to in the text.\nLines 130-137 (right column): I appreciate the authors pointing out and clarifying this!\nLine 174 (right column): but those uniform priors were improper, right? Perhaps best to explicitly state this?\nLines 242-243 (left column): While the authors clearly point to supplementary text, I do think that a bit more information in the\nmain text is justified, i.e. one or two sentences – or a statement between parentheses – of which data dimensions are\nneeded for the approximation to be good would benefit this section.\nLine 192 (right column): same question as before, these are essentially improper priors that were used previously? Please\nclarify in the text.\nLine 295: ‘Furthermore,’\nLine 281: ‘is influenced’\nLine 332: my first concern here is using only 8 baobab species, but I hope to see larger data sets (in terms of number of\nspecies) in the subsequent examples\nLine 336: ‘with a posterior of 100%’: what does this mean? Or what are you trying to say here? I consider this to be very\nvague statement.\nLine 354: even fewer species in this example?\nLine 362: and only three species in the final example? If the authors do not analyse data sets that contain more species than\nthe ones they have analysed here, they will need to make corresponding changes in some of the statements in their\nmanuscript. Don’t get me wrong: yes, these are big data sets but ‘only’ in 1 dimension (the number of loci, but not the number\nof sequences). For example, lines 45-46 should be rephrased to something like: ‘the massive number of loci in modern-day\ndata sets’; similar changes should be made to the manuscript where necessary.\nFurther, I would like the authors to add a paragraph in the introduction of the manuscript that states that they explicitly focus\non data sets with a limited number of species that contain a massive number of loci, and explain why that’s often a key focus\nor use case (rather than data sets that contain dozens of species).\nLine 331 (right column): how many replicate analyses were run per MCMC algorithm tested?\nLine 385 (left column): ‘multi-core machines’, I have not yet seen any machine specifications (processor, memory, number of\ncores, …) so far; it’s important to indicate that in the main text, and to also mention if you used different machines with\ndifferent core counts and how the performance scales with available cores.\nLine 386: so the most impressive reduction for the data set with the most species? I would urge the authors to include a data\nset of at least a few dozen species!\nLine 404 (right column): but you only have a few sequences per locus, right? As I mentioned before, I think the matter of how\ngood the approximation is deserves more attention in the main text.\nLine 664: 5h18m for A2, then certainly a data anlalysis with more species should be feasible (you don’t necessarily have to\ncompare computational performance with the slower methods, just show that it’s now possible to tackle data sets that have\nmore species)!\n(Remarks on code availability)\nI lacked the time to thoroughly test the software myself, apologies.", "reviewer_2": "This paper implements novel improvements to the ambitious and useful BPP package. Perhaps it would help to make\nclearer to general readers where this updated BPP sits in relation to other packages. It is difficult to have a feel for how\nsignificant this advance is. Having not followed the gene-flow side of the MSC for a while, I can see that there have been\nsome recent developments. Papers by G. Jones are cited, and looking at those papers I have an impression that these have\nmotivated the broader approach of analytically integrating out thetas and migration parameters through the use of conjugate\npriors (while recognising the earlier work by Hey and Nielsen in this area). Jones's factorisation in equation 3 seems crucial\nto the approach. It looks as though novelty is in scale, and the tricks outlined in ll 187-248. The examples are impressive in\nterms of loci and numbers of taxa.\n159-160 So this is because the lineage no longer exists because of branching? (determined by \\taus?).\nAlthough there is a lot of description of mixing and mixing efficiency generally throughout the text, there seems no advice on\nappropriate ESS (and whether this actually conforms to what one might expect from independent runs...). Also, judging from\nFig 4c the mixing efficiencies for different parameters are at quite different scales (e.g. \\tau \\varphi in all algorithms), so jointly\none might assume (unless effectively independent in the posterior) the \\taus and \\varphi are a limiting step here.\nIn fact if one uses the approaches described from l. 289 presumably ESS is not necessarily a useful criterion (and even in\nsimple MCMC it can be misleading if burn-in not well chosen) and better would be comparison from independent runs?\nSomething along these lines might be convincing in the Supp section.\nIt may be helpful to have a table with the main parameters defined. For example, having missed the definition of \\varphi in\nthe context of MSC-I (line 268) it took a lot of searching to find out what this parameter is.\nIt looks as though the stopping criteria in Table 1 were just from fixed-length runs. Is there a comparison of the quality of\ninferences? I assume in principle one could obtain posteriors for the \\thetas in A1 from the description starting at l. 289\nIn fact generally the plethora of models in Table S1 lead to a bit of confusion, which is not really handled in the Discussion.\nSo, 4c (also issues below) suggests mixing performance of A1 and A2 are similar, with Table 1 suggesting A2 is faster than\nA1 - and this is the figure highlighted in the paper. But what are the times of all the other variants? and which one should a\nresearcher use? In the last two species B5 seems to have quite high efficiencies for some parameters.\nPossibly more detail is needed to explain what is going on in Fig 4b. It took quite a bit of re-reading to work out that the\nindividual point and interval estimates are for different models.\nSimilarly 4c is quite hard. Remind readers that the text at the top of each figure is saying what the authors expect to be\nobserved. Perhaps some guidance on which model is in which axis (not always obvious).\nPossibly in a more general journal Fig 5 could go in the Supp. It does not do much for me, unless explained better.\n(Remarks on code availability)", "reviewer_3": "Review of Flouri et al: Scaling up Bayesian inference…\nFredrik Ronquist 2025-10-22\nThis paper introduces a set of improvements in the efficiency of MCMC inference for the Multispecies Coalescent (MSC)\nmodel with migration or introgression/hybridization (MSC-M and MSC-I models, respectively, with the acronyms used in the\npaper). Together, these improvements represent a huge leap forward in inference efficiency on these biologically very\nimportant problems. Therefore, I am supportive of publication of this manuscript in a journal that would reach a wide\nauthorship. Landmark computational papers tend to be underrepresented in such journals, which potentially means that\nthere will be unnecessary delays in the relevant research community discovering and taking advantage of these\ndevelopments.\nIn general, I think the paper is solid and well written but it is rather rich in technical detail and light on accessible\nexplanations of the nature of the improvements introduced. I think it would enhance the impact of the paper considerably if\nthis could be addressed by the authors. Most importantly, I think the authors should introduce the major tricks they use\nthrough simple examples that are explained in detail and visualized, so that readers who are not familiar with the details of\nMCMC inference can get an intuitive understanding of why the tricks are so effective compared to the state of the art. The\ncurrent introduction of the trick based on conditional distributions (Left column, L91-96) is too short. I misunderstood it the\nfirst time I read it (see below), despite being fairly experienced with respect to MCMC algorithms. The crucial difference\nbetween sampling from the conditional of one parameter, and keeping that sample when sampling other parameters\n(conditional Gibbs sampling), and the idea of resampling from (or simply updating) the appropriate conditional distributions\n_while_ sampling other parameters, needs to be emphasized and explained. Similarly, I think the idea of approximate\nconditionals could be explained better, potentially with reference to Fig. 8 (which needs an improved caption). I could not\nfind a single reference in the text to this figure. Finally, I think it would be helpful to add an illustration demonstrating the\nefficiency gain of sampling from the expectation of the conditionals rather than from samples of them.\nTo make room for this expansion, explaining the inference improvements in a more intuitive way to readers unfamiliar with\nMCMC techniques or statistical inference in general, I think it is possible to move more of the technical details into the online\nsupplemental information, without losing any of the necessary rigor. Aside from these general comments, I do have some\nminor suggestions for improvement, as follows.\nLeft L91-96. “[VDROP] proposes certain parameters from their conditional distribution (achieving the same mixing efficiency\nas analytical integration)”. I first understood this incorrectly, so I think this needs rephrasing or more explanation. Consider,\ne.g., a bivariate normal with strong correlation structure. Updating each one of the two variables separately from their\nconditional posteriors, or one with an MH step and the other one separately from its conditional posterior will clearly not have\nthe same mixing efficiency as integrating one or both of them out. What does have the same mixing efficiency is sampling\none of them using some MH step, while also resampling the other one from its conditional posterior (or simply updating the\nconditional posterior of the other variable) given the proposed new state of the first variable. Please clarify this. See general\ncomments for expansion on this and related issues.\nRight L86. “differences of the” should probably be “differences compared to the”\nRight L104 (appr). “\\tau and \\theta are measured in the expected number of mutations per site”. This is a very unusual\nmeasure of population size (\\theta). How does this work? I assume this requires some notion of the sequence variation in the\npopulation?\nRigh L111 (appr) “for two species” should probably be “for the two species”\nRight L118-119. Maybe it is unnecessary to point out this but the history also includes the location of the migration events on\nthe tree.\nRight L141: “Let G_j be part of” should probably be “Let G_j be the part of”\nRight L238-243: “To allow parallelization…” is a critical paragraph and the approach could be described better.\nLeft L249ff: “A1(I-IG-int)” and “A2 (I-IG-g)” seem overly complicated names for the two algorithms, and it is difficult to\nunderstand how they were derived. Please simplify the names, or explain what the acronyms mean. The explanation follows\non Right L335ff but it needs to be here, at least in a simplified form and a reference to the more detailed explanation.\nLeft L264ff: “the mixing steps…”. Difficult to follow, please rewrite.\nRight L281: “c in influenced” should be “c is influenced”\nRight L288: “Consider for instance…” I find the example relatively uninformative with respect to the gains one might expect\nin the types of problems considered here. I suggest deleting this paragraph.\nLeft L334. “The BPP analysis of the data under the MSC model inferred the species tree…”. In Figure 2, it is not clear how\nthe species tree is inferred; it appears that only species split times are changed, not the species tree topology. Please clarify\nhow the algorithm samples over species tree topologies.\nRight L313ff: This paragraph is entirely focused on mixing efficiency per generation or iteration of the MCMC algorithm. I\nthink it would be good to indicate already here that the clock time and the computational resources needed are also\nimportant factors to consider in comparing the performance of different algorithms (as is also done in the discussion of\nresults later in this section).\nRight L340: “it” should be “int”\nRight L363ff: It is not quite clear what algorithms are identical to the ones used in previous versions of BPP, and how the\nnew algorithms compare to them. Please clarify (Is A6 the same as that used in BPP previously, or not? Similarly, which B\nalgorithm if any is identical to the previous BPP algorithm).\nLeft L475: “events” should be “event”\nFig. 1. I suggest using the standard conventions of graphical models rather than this ad-hoc format. If so, then \\tau, \\theta and\n\\omega should not be surrounded by a square; in typical graphical models, a square with rounded corners or other features\nindicating replication is used for a template or plate of replicated iid variables, but these three random variables are not\nreplicated. Instead, these three variables should be free circles (independent random variables), each pointing\nindependently to G_i. Furthermore, G_i and X_i should be on the same plate, not on separate plates, and X_i should be\nshaded to indicated that it is observed / clamped. Finally, the squares with priors (Gamma-Dirichlet prior etc) should be\nreplaced by the parameters of the gamma-dirichlet etc priors inside simple squares, as they are constant nodes in the model.\nThe distributions themselves are typically not included in the graphical model.\nFig. 2. It would be better if this figure used a different drawing style from Fig. 1, so that it is clear that it is an\nMCMC algorithm illustration and not a graphical model. It is unclear how species tree topologies are sampled; please clarify.\nAlso, it is unclear here and elsewhere in this text whether hybridization/introgression events are assumed a priori, or whether\nthe MSC-I algorithms sample over hybridization/introgression events. The Results section seems to indicate that those events\nare specified in the prior and not inferred, and the fact that sampling over them is not considered in this figure points in the same\ndirection, but there needs to be more clarity on this point both in the text and in the figure.\nFigs 4, 6", "reviewer_7": "I think it would be better to break the panels in the c) part of these figures into one set for the\nMSC-M model and one for the MSC-I model. Alternatively, introduce a vertical line separating the estimates\nfrom each of the two models. These are different models and it is not surprising if the estimates are different.\nAlso, please comment on the few cases where estimates do seem to differ between algorithms generating\nestimates from the same model. There seems to be a few of those cases in Fig. 7, for instance. Could they be\ndue to some of the approximations used or are there other likely explanations?\nWe have followed the reviewer’s suggestion to insert some space and a thin vertical line separating results for\nthe two models in figures 5, 6, 7. We have also edited the legend to figure 5 explain things more properly.\nSome of the differences in parameter estimates are due to the use of different priors (gamma versus inverse-\ngamma for s) under the same model.\nFig. 8. Please explain in the caption what this figure is meant to illustrate; it is not clear from the current\ncaption. I am not convinced that the figure is needed, I could not find a single reference to it from the text.\nWith a better caption, there might be reasons to keep it (see general comments).\nWe have now deleted this figure. We have included a brief summary of figure S1 (p.4 right column), which is\nan expanded version of the old figure 8.\n12", "reviewer_author_response": "RESPONSE TO REVIEWERS' COMMENTS\nReviewer #1 (Remarks to the Author):\nThe authors describe an interesting approach to achieve improved performance for Bayesian inference under\nthe multispecies coalescent model, and are clearly experts in this matter. The manuscript is well structured\nand sufficiently easy to read through; I list my comments below.\nLine 31 : can this sentence be better connected to the previous one, i.e. explicitly state which of (i) and / or (ii)\nyou’re tackling with your new method? If you are tackling (i) but not (ii) then you should clearly mention this\nin the abstract (I don’t mean that in any bad way).\nWe are making improvements in both (i) MCMC mixing and (ii) parallelisation, so we have added (p.1)\n“To deal with both issues,…”\nLine 36: this leaves doubt as to whether mixing efficiency is expressed as effective sample size, or as effective\nsample size per time unit (the latter is the one that matters); please be explicit about this in the abstract.\nThe algorithms discussed in our study have similar running times, so the use of ESS and ESS per minute to\ncompare algorithms will be equivalent. The more advanced algorithms involve more effort in coding but the\nrunning time is similar. We have included the following text (p.6 left column):\n“unless stated otherwise explicitly, the algorithms tested in this study have similar running time per iteration.”\nLine 42: please avoid using terms such as ‘Bayesian MCMC’, as there’s nothing Bayesian about MCMC (in this\ncase: simply remove ‘Bayesian’).\nChanged as suggested.\nLine 63 (left column): this could do with a citation.\nWe have added a few citations.\nLine 63 (left column): this does not change the serial nature of MCMC however, and you should point that out\nto the readers. What you are doing is more work in a single iteration / generation, which can indeed be\nachieved through parallelization and / or more complex transition kernels. But neither break / avoid / bypass /\nchange the serial nature of MCMC, and it would be good to point this out to the readers to avoid giving the\nwrong impression.\nIt is rare for a major high-dimensional MCMC algorithm to have only one proposal changing all dimensions.\nMore commonly multiple proposals each changing a subset of the parameters are used in each MCMC\niteration/generation. Typically all those proposals have to be done sequentially. By relying on conditional\nindependence, we are doing some of those proposals in parallel. If we consider the full MCMC iteration, it is\ntrue that the next iteration has to wait until the current iteration is completed, but within each iteration, we\nare bypassing or changing the serial nature of MCMC.\nIn response to this comment, we have added the following text (p.2 left column):\n“Most high-dimensional MCMC algorithms use multiple proposal steps, each updating a subset of parameters,\nwhile collectively the steps ensure that every dimension of the state space is updated in each MCMC iteration.\niteration can be applied in parallel by relying on conditional independence.”\nLine 91: ‘referred to’, do you mean ‘which you will refer to in this manuscript’? If not – and it’s a term that has\nseen use before – this needs a citation.\nWe have changed “referred to” to “which we refer to”.\nLine 91 (right column): please number your figures in the order they are referred to in the text.\nChanged as suggested.\nLines 130-137 (right column): I appreciate the authors pointing out and clarifying this!\nLine 174 (right column): but those uniform priors were improper, right? Perhaps best to explicitly state this?\nProper uniform priors have been used in different versions of the IMa program. These uniform priors appear\nto create problems in comparisons of models. We have inserted the following text (p.4 left column):\n“using uniform priors with user-specified upper bounds; see Thawornwattana et al. (2026) for a discussion of\nissues arising from the use of uniform priors in the IMa3 program (Hey et al., 2018).”\nLines 242-243 (left column): While the authors clearly point to supplementary text, I do think that a bit more\ninformation in the main text is justified, i.e. one or two sentences – or a statement between parentheses – of\nwhich data dimensions are needed for the approximation to be good would benefit this section.\nThe idea of approximate conditionals appears to be new and is potentially useful in other major MCMC\napplications. We have now expanded our summary in the main text, mentioning figures S1 and S2 (p.4 right\ncolumn). The data size here is the number of coalescent events and increases with the number of loci and the\nnumber of sequences. We have the following (p.4):\n“We evaluated the approximation under different scenarios (e.g., with the prior being consistent or conflicting\nwith the data) and found that the approximation is good even with only 2 or 10 events (fig. S1). Both the\ngamma and inverse-gamma approximations converge to the correct conditional when the amount of data\n(i.e., the total number of coalescent events over all loci in that population) approaches infinity, so that the\napproximation is more accurate in larger datasets with more loci and more sequences.”\nLine 192 (right column): same question as before, these are essentially improper priors that were used\npreviously? Please clarify in the text.\nAs before, proper uniform priors have been used in the IMa program. We have inserted the words “with user-\nspecified upper bounds”.\nLine 295: ‘Furthermore,’\nLine 281: ‘is influenced’\nFixed.\nLine 332: my first concern here is using only 8 baobab species, but I hope to see larger data sets (in terms of\nnumber of species) in the subsequent examples\nWhile it is worthwhile to improve the algorithms to handle large phylogenies, we note that the biological\nproblem of gene flow exists on small phylogenies with only two species. The number of loci is the sample size\nin the inference problem, and it is important to implement algorithms that can deal with thousands of loci so\nthat inference of gene flow becomes meaningful. We have added the following text to emphasized the\nimportance of handing lots of loci (p.2):\n5\n“Indeed under the MSC model (with and without gene flow), sequence data at different loci are independently\nand identically distributed (i.i.d.), so that the number of loci is the sample size in the inference problem\n(Rannala and Yang, 2003; Jiao et al., 2021). Among many factors that affect the information content in the\ndata (such as the number of loci, the number of sequences sampled per species, the sequence length, and the\nmutation rate), the number of loci is often the most important (e.g., Huang et al., 2020; Thawornwattana et\nal., 2025, 2026). For example, to estimate the rates of gene flow reliably, thousands of loci are often needed\nand it is important to develop algorithms that can accommodate many loci even if for a relatively small\nphylogeny.”\nLine 336: ‘with a posterior of 100%’: what does this mean? Or what are you trying to say here? I consider this\nto be very vague statement.\nPrevious studies found only weak support for estimated phylogenies. Here we found 100% support. We have\nedited the text to offer some possible explanations for the discrepancy (p.7):\n“Previous analyses of genetic data found considerable uncertainty in the species phylogeny (Karimi et al.,\n2020; Wan et al., 2024). We inferred the species phylogeny under the MSC model using bpp (Yang and\nRannala, 2014; Rannala and Yang, 2017). The maximum a posteriori (MAP) tree (fig. 5a) has a posterior of ∼\n100%, possibly because the full likelihood method uses the information in the data more efficiently than\nsummary methods used in previous studies.”\nLine 354: even fewer species in this example?\nPlease see our major comment #2 to the Editor.\nLine 362: and only three species in the final example? If the authors do not analyse data sets that contain\nmore species than the ones they have analysed here, they will need to make corresponding changes in some\nof the statements in their manuscript. Don’t get me wrong: yes, these are big data sets but ‘only’ in 1\ndimension (the number of loci, but not the number of sequences). For example, lines 45-46 should be\nrephrased to something like: ‘the massive number of loci in modern-day data sets’; similar changes should be\nmade to the manuscript where necessary.\nAs far as we know, the three datasets analysed in our paper are beyond the power of competing software\nprograms including *beast, phylonet/MCMC-seq, G-PhoCS, etc. Bpp has been used to analyze larger datasets\nthan those used here (e.g., 17,000 loci, each with 12 sequences on a phylogeny for 6 species in Flouri et al.\n2023, table 1). In this study our objective is to evaluate MCMC algorithms rather than testing the limits of\nbpp. From the viewpoint of statistical inference, the number of loci is the sample size, so it is important to\naccommodate many loci. We have added some text to emphasize this. Please see our major comment #1 to\nthe Editor above.\nFurther, I would like the authors to add a paragraph in the introduction of the manuscript that states that they\nexplicitly focus on data sets with a limited number of species that contain a massive number of loci, and\nexplain why that’s often a key focus or use case (rather than data sets that contain dozens of species).\nWe have followed this suggestion and added a paragraph in the Introduction.\nLine 331 (right column): how many replicate analyses were run per MCMC algorithm tested?\nWe did two replicate runs. The reliability of the runs is also clear from the nearly identical results (posterior\nmeans and CIs) between algorithms applied to the same data under the same model and priors in figures 5, 6,\n7.\nLine 385 (left column): ‘multi-core machines’, I have not yet seen any machine specifications (processor,\nmemory, number of cores, …) so far; it’s important to indicate that in the main text, and to also mention if you\nused different machines with different core counts and how the performance scales with available cores.\nIntel Xeon Gold 6154 18C processors” (p.7). We used the same computer in our tests.\nLine 386: so the most impressive reduction for the data set with the most species? I would urge the authors to\ninclude a data set of at least a few dozen species!\nWe have already included three empirical datasets, with realistic sizes. As explained above, the number of\nspecies is not the most important factor to influence the utility of our algorithms. The number of species in\nthe three empirical datasets is 8, 6, and 3. We have seen consistent improvements of the new algorithms in\nthe three datasets, so we see no reason to believe that similar improvements are not possible when data of 20\nor 50 species are analyzed. The objective of our study is to develop MCMC algorithms with improved mixing,\nrather than finding out the computational limits of the current version of bpp in terms of the number of\nspecies.\nLine 404 (right column): but you only have a few sequences per locus, right? As I mentioned before, I think the\nmatter of how good the approximation is deserves more attention in the main text.\nWe have now added extra information for the 3 empirical datasets, such as the number of sequences per\nlocus and the sequence length. For the baobab dataset, there are 9 species, 344 loci, 28-38 sequences per\nlocus, and 759-9574 sites per sequence. The dataset is not small.\nHere the data size for the approximation is the number of coalescent events and the total lineage-pair\ncoalescent waiting time. We have added some explanatory text (p.8):\n“(i.e., datasets with more loci and/or more sequences per species per locus leading to more coalescent\nevents)”\nLine 664: 5h18m for A2, then certainly a data anlalysis with more species should be feasible (you don’t\nnecessarily have to compare computational performance with the slower methods, just show that it’s now\npossible to tackle data sets that have more species)!\nWe estimate that for simple models it should be feasible to run bpp on a species tree with 100 species.\nHowever, we have already included three datasets, and the paper includes many highly technical sections, so\nthat a fourth dataset analyzed in differently from the other three would take space to describe. As we\nexplained above, we have now added text to explain the importance of being able to handle many loci. Please\nsee our comment above.\nReviewer #1 (Remarks on code availability):\nI lacked the time to thoroughly test the software myself, apologies.\nReviewer #2 (Remarks to the Author):\nThis paper implements novel improvements to the ambitious and useful BPP package. Perhaps it would help\nto make clearer to general readers where this updated BPP sits in relation to other packages. It is difficult to\nhave a feel for how significant this advance is. Having not followed the gene-flow side of the MSC for a while, I\ncan see that there have been some recent developments. Papers by G. Jones are cited, and looking at those\npapers I have an impression that these have motivated the broader approach of analytically integrating out\nthetas and migration parameters through the use of conjugate\npriors (while recognising the earlier work by\nHey and Nielsen in this area). Jones's factorisation in equation 3 seems crucial to the approach. It looks as\nthough novelty is in scale, and the tricks outlined in ll 187-248. The examples are impressive in terms of loci\nand numbers of taxa.\nOur Discussion section has a summary of advancements made in this paper. The factorisation is mostly due to\nHey and Nielsen (2017), and Jones (2019) is independent or duplicated work with no correct software\nimplementation.\nwith\n“Currently bpp is the only program that appears to implement the MSC-M model correctly …”\n159-160 So this is because the lineage no longer exists because of branching? (determined by \\taus?).\nYes. We have added the following text to explain when the indicator I = 1 (p.3):\nsjk\n“(that is, if both species s and j exist during time segment k and are permitted to exchange migrants in the\nmodel)”\nChanged as suggested.\nAlthough there is a lot of description of mixing and mixing efficiency generally throughout the text, there\nseems no advice on appropriate ESS (and whether this actually conforms to what one might expect from\nindependent runs...). Also, judging from Fig 4c the mixing efficiencies for different parameters are at quite\ndifferent scales (e.g. \\tau \\varphi in all algorithms), so jointly one might assume (unless effectively\nindependent in the posterior) the \\taus and \\varphi are a limiting step here.\nCheck.\nIn fact if one uses the approaches described from l. 289 presumably ESS is not necessarily a useful criterion\n(and even in simple MCMC it can be misleading if burn-in not well chosen) and better would be comparison\nfrom independent runs? Something along these lines might be convincing in the Supp section.\nCheck.\nIt may be helpful to have a table with the main parameters defined. For example, having missed the definition\nof \\varphi in the context of MSC-I (line 268) it took a lot of searching to find out what this parameter is.\nWe have redrawn figure 1 to include the MSC-I model as well as the MSC-M model. We define parameters in\nboth models in the figure legend.\nIt looks as though the stopping criteria in Table 1 were just from fixed-length runs. Is there a comparison of\nthe quality of inferences? I assume in principle one could obtain posteriors for the \\thetas in A1 from the\ndescription starting at l. 289\nYes, our MCMC runs are of fixed length. The different algorithms produce the same posterior so the ‘quality’\nof inference should be the same. Their differences concern the computational cost; it takes longer for an\ninefficient algorithm to produce estimates at the same precision than an algorithm with improved mixing.\nIn fact generally the plethora of models in Table S1 lead to a bit of confusion, which is not really handled in the\nDiscussion. So, 4c (also issues below) suggests mixing performance of A1 and A2 are similar, with Table 1\nsuggesting A2 is faster than A1 - and this is the figure highlighted in the paper. But what are the times of all\nthe other variants? and which one should a researcher use? In the last two species B5 seems to have quite\nhigh efficiencies for some parameters.\nPossibly more detail is needed to explain what is going on in Fig 4b. It took quite a bit of re-reading to work\nout that the individual point and interval estimates are for different models.\nSimilarly 4c is quite hard. Remind readers that the text at the top of each figure is saying what the authors\nexpect to be observed. Perhaps some guidance on which model is in which axis (not always obvious).\nis now expanded to have a detailed description of one set of results (for the baobabs of figure 5).\nPossibly in a more general journal Fig 5 could go in the Supp. It does not do much for me, unless explained\nbetter.\nThe old figure 5 is now figure 8. This contrasts the computational cost of algorithms A1 and A2, which is a\nmajor result of our paper. We have edited the figure legend to explain why the two algorithms have very\ndifferent CPU usage when run on multi-threaded computers and why their running times are so different.\nReviewer #3 (Remarks to the Author):\nReview of Flouri et al: Scaling up Bayesian inference…\nFredrik Ronquist 2025-10-22\nThis paper introduces a set of improvements in the efficiency of MCMC inference for the Multispecies\nCoalescent (MSC) model with migration or introgression/hybridization (MSC-M and MSC-I models,\nrespectively, with the acronyms used in the paper). Together, these improvements represent a huge leap\nforward in inference efficiency on these biologically very important problems. Therefore, I am supportive of\npublication of this manuscript in a journal that would reach a wide authorship. Landmark computational\npapers tend to be underrepresented in such journals, which potentially means that there will be unnecessary\ndelays in the relevant research community discovering and taking advantage of these developments.\nIn general, I think the paper is solid and well written but it is rather rich in technical detail and light on\naccessible explanations of the nature of the improvements introduced. I think it would enhance the impact of\nthe paper considerably if this could be addressed by the authors. Most importantly, I think the authors should\nintroduce the major tricks they use through simple examples that are explained in detail and visualized, so\nthat readers who are not familiar with the details of MCMC inference can get an intuitive understanding of\nwhy the tricks are so effective compared to the state of the art. The current introduction of the trick based on\nconditional distributions (Left column, L91-96) is too short. I misunderstood it the first time I read it (see\nbelow), despite being fairly experienced with respect to MCMC algorithms. The crucial difference between\nsampling from the conditional of one parameter, and keeping that sample when sampling other parameters\n(conditional Gibbs sampling), and the idea of resampling from (or simply updating) the appropriate conditional\ndistributions _while_ sampling other parameters, needs to be emphasized and explained. Similarly, I think the\nidea of approximate conditionals could be explained better, potentially with reference to Fig. 8 (which needs\nan improved caption). I could not find a single reference in the text to this figure. Finally, I think it would be\nhelpful to add an illustration demonstrating the efficiency gain of sampling from the expectation of the\nconditionals rather than from samples of them.\nWe thank the reviewer for the encouraging comments and for the summary of our results.\nWe have followed the reviewer’s suggestion and added a simple example with the bivariate Gaussian target to\nillustrate the algorithms. Please see our major comment #1 to the Editor above.\nTo make room for this expansion, explaining the inference improvements in a more intuitive way to readers\nunfamiliar with MCMC techniques or statistical inference in general, I think it is possible to move more of the\ntechnical details into the online supplemental information, without losing any of the necessary rigor. Aside\nfrom these general comments, I do have some minor suggestions for improvement, as follows.\nLeft L91-96. “[VDROP] proposes certain parameters from their conditional distribution (achieving the same\nmixing efficiency as analytical integration)”. I first understood this incorrectly, so I think this needs rephrasing\nor more explanation. Consider, e.g., a bivariate normal with strong correlation structure. Updating each one of\nthe two variables separately from their conditional posteriors, or one with an MH step and the other one\nseparately from its conditional posterior will clearly not have the same mixing efficiency as integrating one or\nboth of them out. What does have the same mixing efficiency is sampling one of them using some MH step,\nwhile also resampling the other one from its conditional posterior (or simply updating the conditional\nposterior of the other variable) given the proposed new state of the first variable. Please clarify this. See\n9\ngeneral comments for expansion on this and related issues.\nAll of the reviewer’s interpretations here are correct. We have now followed the reviewer’s suggestion and\nadded a section titled “Simple example for a bivariate Gaussian target to illustrate the algorithms”, as well as a\nnew figure 4, to illustrate the algorithms (p.5-6). Please see our major comment #1 to the Editor above.\nRight L86. “differences of the” should probably be “differences compared to the”\nChanged as suggested.\nRight L104 (appr). “\\tau and \\theta are measured in the expected number of mutations per site”. This is a very\nunusual measure of population size (\\theta). How does this work? I assume this requires some notion of the\nsequence variation in the population?\nThe average coalescent waiting time is 2N so that the average sequence distance or heterozygosity is  = 4N.\nViewing the mutation rate  as a scalar, the heterozygosity  is also known as the population size parameter.\nThis is a fundamental parameter in population genetics models as it among other things determines the\ncoalescent rate.\nRigh L111 (appr) “for two species” should probably be “for the two species”\nWe think “for two species” should be ok.\nRight L118-119. Maybe it is unnecessary to point out this but the history also includes the location of the\nmigration events on the tree.\nYes, the migration history at the locus includes the number, directions and timings of migration events. We\nhave now mentioned this.\nRight L141: “Let G_j be part of” should probably be “Let G_j be the part of”\nChanged as suggested.\nRight L238-243: “To allow parallelization…” is a critical paragraph and the approach could be described better.\nWe hope that the newly added example of the bivariate Gaussian target should make all our algorithms easy\nto understand.\nLeft L249ff: “A1(I-IG-int)” and “A2 (I-IG-g)” seem overly complicated names for the two algorithms, and it is\ndifficult to understand how they were derived. Please simplify the names, or explain what the acronyms\nmean. The explanation follows on Right L335ff but it needs to be here, at least in a simplified form and a\nreference to the more detailed explanation.\nThis is lost during editing.\nWe have made sure that the abbreviations are properly explained when they are used for the first time.\nLeft L264ff: “the mixing steps…”. Difficult to follow, please rewrite.\nWe have changed this to (p.5 right column):\n“Step 5 is a scaling or mixing move that multiplies all species split times on the species tree () and coalescent\nand migration times on all gene trees by the same scale factor.”.\nRight L281: “c in influenced” should be “c is influenced”\nFixed.\nmight expect in the types of problems considered here. I suggest deleting this paragraph.\nWe considered this suggestion and decided to keep the brief summary in the main text. Our theory not only\nproposes a better estimator (u versus y) but also allows calculation of its performance gain in the reduced\nvariance. The example is partly used to illustrate the calculation.\nLeft L334. “The BPP analysis of the data under the MSC model inferred the species tree…”. In Figure 2, it is not\nclear how the species tree is inferred; it appears that only species split times are changed, not the species tree\ntopology. Please clarify how the algorithm samples over species tree topologies.\nWe have edited the text to read as follows:\n“We inferred the species phylogeny under the MSC model using bpp (Yang and Rannala, 2014; Rannala and\nYang, 2017). The maximum a posteriori (MAP) tree (fig. 4a) has a posterior of ∼ 100%, possibly because the\nfull likelihood method uses the information in the data more efficiently than summary methods used in\nprevious studies.”\nThe cited papers (Yang and Rannala, 2014; Rannala and Yang, 2017) describe algorithms that search in the\nspace of species phylogenies.\nRight L313ff: This paragraph is entirely focused on mixing efficiency per generation or iteration of the MCMC\nalgorithm. I think it would be good to indicate already here that the clock time and the computational\nresources needed are also important factors to consider in comparing the performance of different algorithms\n(as is also done in the discussion of results later in this section).\nWe have inserted the following text (p.6 left column):\n“…unless stated otherwise explicitly, the algorithms tested in this study have similar running time per\niteration.”\nRight L340: “it” should be “int”\nRight L363ff: It is not quite clear what algorithms are identical to the ones used in previous versions of BPP,\nand how the new algorithms compare to them. Please clarify (Is A6 the same as that used in BPP previously, or\nnot? Similarly, which B algorithm if any is identical to the previous BPP algorithm).\nYes, we now mention that A6 (and A3 and B6) are equivalent to algorithms in the previous version of bpp\n(ver4.7). In tables S3-S5 we calculate the performance gain relative to the old algorithms.\nLeft L475: “events” should be “event”\nFixed.\nFig. 1. I suggest using the standard conventions of graphical models rather than this ad-hoc format. If so, then\n\\tau, \\theta and \\omega should not be surrounded by a square; in typical graphical models, a square with\nrounded corners or other features indicating replication is used for a template or plate of replicated iid\nvariables, but these three random variables are not replicated. Instead, these three variables should be free\ncircles (independent random variables), each pointing independently to G_i. Furthermore, G_i and X_i should\nbe on the same plate, not on separate plates, and X_i should be shaded to indicated that it is observed /\nclamped. Finally, the squares with priors (Gamma-Dirichlet prior etc) should be replaced by the parameters of\nthe gamma-dirichlet etc priors inside simple squares, as they are constant nodes in the model. The\ndistributions themselves are typically not included in the graphical model.\nWe have redrawn figure 2 (the old figure 1) to conform with the DAG syntax.\nFig. 2. It would be better if this figure used a different drawing style from Fig. 1, so that it is clear that it is an\nMCMC algorithm illustration and not a graphical model. It is unclear how species tree topologies are sampled;\nplease clarify. Also, it is unclear here and elsewhere in this text whether hybridization/introgression events are\n11\nassumed a priori, or whether the MSC-I algorithms sample over hybridization/introgression events. The\nResults section seems to indicate that those events are specified in the prior and not inferred, and the fact that\nsampling over them is not considered in this figure points in the same direction, but there needs to be more\nclarity on this point both in the text and in the figure.\nWe have redrawn figure 3 (the old figure 2) so that it does not look like a DAG diagram.\nYes, the models of gene flow (MSC-I and MSC-M) are implemented as fully specified models.\nWe have algorithms that change the species tree topology under MSC with no gene flow, but don’t have\ncross-model algorithms that change the species phylogeny or gene-flow events under models of gene flow.\nIn response to this comment and to reviewer #2’s comment “where this updated BPP sits in relation to other\npackages”), We have included a brief discussion of the strengths and weaknesses of bpp, starting with (p.9)\n“Currently bpp is the only program that appears to implement the MSC-M model correctly... However, only\nwithin-model algorithms for estimating parameters such as the rates of gene flow are implemented in bpp,\nwhile cross-model MCMC moves for searching in the space of gene-flow models are lacking…”\nFigs 4, 6"}, "author_response": ""}, {"round": 2, "reviewer_comments": {"reviewer_1": "I have carefully reviewed the changes the authors have made to their manuscript, along with their rebuttal letter. In their\nresponses, the authors refer a few times to their replies to the Editor but I was unable to locate these, which was quite\nunfortunate. I must hence leave it to the Editor to carefully check those responses; I support publication of this revised\nmanuscript in Nature Communications.\n(Remarks on code availability)", "reviewer_2": "I have read through the authors' response letter and revised manuscript and I can see that my comments have been well\naddressed. I think it is a significant and methodologically sound contribution to the field.\n(Remarks on code availability)"}, "author_response": ""}]

Decision letter