StarBeast3: Adaptive Parallelized Bayesian Inference under the Multispecies Coalescent

Abstract As genomic sequence data become increasingly available, inferring the phylogeny of the species as that of concatenated genomic data can be enticing. However, this approach makes for a biased estimator of branch lengths and substitution rates and an inconsistent estimator of tree topology. Bayesian multispecies coalescent (MSC) methods address these issues. This is achieved by constraining a set of gene trees within a species tree and jointly inferring both under a Bayesian framework. However, this approach comes at the cost of increased computational demand. Here, we introduce StarBeast3—a software package for efficient Bayesian inference under the MSC model via Markov chain Monte Carlo. We gain efficiency by introducing cutting-edge proposal kernels and adaptive operators, and StarBeast3 is particularly efficient when a relaxed clock model is applied. Furthermore, gene-tree inference is parallelized, allowing the software to scale with the size of the problem. We validated our software and benchmarked its performance using three real and two synthetic data sets. Our results indicate that StarBeast3 is up to one-and-a-half orders of magnitude faster than StarBeast2, and therefore more than two orders faster than *BEAST, depending on the data set and on the parameter, and can achieve convergence on large data sets with hundreds of genes. StarBeast3 is open-source and is easy to set up with a friendly graphical user interface. [Adaptive; Bayesian inference; BEAST 2; effective population sizes; high performance; multispecies coalescent; parallelization; phylogenetics.]

Existing methods for testing macroevolutionary and macroecological questions have not kept pace with the explosion of next-generation sequence data now available (Blom et al. 2016b;Bragg et al. 2017;Stenson et al. 2017). Despite burgeoning databases of withinand between-species genomic diversity (Blom et al. 2016b;Bragg et al. 2017;Stenson et al. 2017), it is still common practice to ignore the gene-tree discordance that underlies any species phylogeny inferred from multilocus sequences and instead infer species ancestry based on concatenated sequence data taken to represent all underlying gene histories (Degnan and Rosenberg 2009;Heled and Drummond 2010;Jones 2017;Ogilvie et al. 2017;Rannala and Yang 2017). While this approach can perform well for inferring topologies when branches are long and incomplete lineage sorting (ILS) is absent, these conditions are rarely met.
Species trees inferred from concatenated sequences are often topologically incorrect (Degnan and Rosenberg 2009;Heled and Drummond 2010;Ogilvie et al. 2017), provide biased estimates for branch lengths and substitution rates (Kubatko et al. 2011;Ogilvie et al. 2016;Mendes and Hahn 2016), and underestimate uncertainty in tree topology, resulting in an unjustified degree of confidence in the wrong tree (Heled and Drummond 2010;Ogilvie et al. 2017). Such biases are exacerbated by subsampling of incongruent genes (Edwards et al. 2016;Mendes and Hahn 2016) and hold even for deep splits in the tree (Oliver 2013). These are crucial concerns in themselves and, more generally, can lead to biased estimates and erroneous inferences about fundamental evolutionary and ecological processes that require accurate phylogenetic trees, such as rates of speciation and extinction (Cadena et al. 2011;Rowe et al. 2011;Pepper et al. 2013), rates of substitution in DNA sequences (Bouckaert et al. 2013) and morphological characters (Pepper et al. 2013), species ancestry and ancestral age estimation (Mitchell et al. 2014), geographical history and origins (Lemey et al. 2009;Bouckaert 2016), and species delimitation (Yang and Rannala 2010;Grummer et al. 2013;Leaché et al. 2014;Yang and Rannala 2014).
The multispecies coalescent (MSC; Maddison 1997;Edwards 2009;Liu et al. 2009) is an approach designed to minimize these potential biases by modeling macroevolution as a distribution of gene trees constrained by a species tree (Degnan and Rosenberg 2009;Heled and Drummond 2010;Jones 2017;Ogilvie et al. 2017;Rannala and Yang 2017). In doing so, the MSC provides a more biologically realistic framework for phylogenetic inference that captures the process of ILS underlying most multilocus phylogenies. Furthermore, by explicitly modeling both species and gene trees, the MSC can address questions that cannot be addressed under a concatenation approach-such as automatic species delimitation (Fujita et al. 2012), with important implications for biodiversity assessment and conservation (Bickford et al. 2007).
A number of software packages have implemented the MSC in various ways (see review by ). Our work at the Centre for Computational Evolution at the University of Auckland has led the development of *BEAST (STARBeast; Heled and Drummond 2010) and StarBeast2 (Ogilvie et al. 2017)-full Bayesian MSC frameworks for species-tree estimation from multilocus sequence data-and UglyTrees for visualizing these models (Douglas 2020). By explicitly modeling the MSC and avoiding the biases associated with concatenation methods (Heled and Drummond 2010; Ogilvie et al.
VOL. 71 FIGURE 1. Depiction of the multispecies coalescent model, with k = 2 gene trees constrained within a single species tree S with n = 8 species. In this depiction, node heights (age) run along the y-axis and species-tree node widths are proportional to effective population sizes (arbitrary units). The relative molecular substitution rate of each species-tree branch is proportional to line thickness. Tree was built from a Gopher data set (Belfiore et al. 2008) and visualized using UglyTrees (Douglas 2020(Douglas ). 2016(Douglas , 2017, an analysis using either of these software packages can significantly improve the conclusions drawn from data. However, despite some advances in computational efficiency of the full Bayesian MSC (Jones 2017;Ogilvie et al. 2017;Rannala and Yang 2017), these complex models remain computationally intractable for large next-generation sequence data sets of 100's of sequenced loci across hundreds of individuals (i.e., 10 4 -10 6 samples×loci). As a result, existing applications of the approach have tended to consider smaller data sets (Kang et al. 2014;Blom et al. 2016a) or to ignore much of the available data (Blom et al. 2016b;Bragg et al. 2017;Stenson et al. 2017), which reduces accuracy and increases uncertainty in species-tree estimates (Song et al. 2012;Ogilvie et al. 2017). One approach to this problem has been the development of much simpler summary coalescent methods that utilize distributions of estimated gene-tree topologies as input to rapidly process large data sets . These include the rooted triplet method MP-EST (Liu et al. 2010) and the quartet method ASTRAL (Mirarab et al. 2014). However, summary coalescent methods are sensitive to gene-tree errors (Mirarab and Warnow 2015;Xi et al. 2015) and produce trees in coalescent units, and thus time and population size estimates used by downstream analyses are confounded.
Here, we aim to perform Bayesian inference on large data sets using the Markov chain Monte Carlo (MCMC) algorithm as our workhorse. As illustrated in Figure 1, the number of parameters involved is quite large, as is the accompanying state space. We develop a set of new MCMC proposals to explore state space in a much more efficient way than previous implementations and demonstrate we can handle data sets several times faster than *BEAST and StarBeast2. The resulting software package StarBeast3 is available as an open-source BEAST 2 package (Bouckaert et al. 2019).

The MSC
Our objective is to develop efficient methods in a Bayesian framework for analyzing models where there is a phylogeny, S, such as a species or language tree, that forms a constraint on a set of k trees G ={g 1 ,...,g k }, such as gene trees. Each taxon within G is assigned to a single taxon within S, from some fixed individual-tospecies mapping function (Fig. 1). Species tree S = (T S ,t S ) consists of a topology T S and divergence times t S , as does the set of gene trees G = (T G ,t G ).
All trees are assumed to be binary rooted time trees, where branch lengths describe the passing of time from the root of the tree down to the tips. Taxon node heights are assumed to be fixed and are typically extant (with height 0). Each gene tree g i consists of 2n g i −1 nodes and 2n g i −2 branches for taxon count n g i , while S consists of 2n S −1 nodes and 2n S −1 branches, including a root branch, for species count n S . Genetree taxa are associated with data D ={D 1 ,...,D k }, for example, nucleotide sequences or cognate data. Let be a set of model parameters, for instance, those related to the speciation or nucleotide substitution processes. Consider the posterior density function p(S,G,|D): The MSC model is therefore hierarchical. S can follow a range of tree prior distributions p(S|), such as the Yule (Yule 1925) or birth-death models (Nee et al. 1994). Whereas, each gene tree g i is assumed to follow the MSC process (Degnan and Rosenberg 2009;Heled and Drummond 2010;Jones 2017;Ogilvie et al. 2017;Rannala and Yang 2017), under which speciestree branches are associated with independently and identically distributed (effective) population sizes N e which govern the coalescent process of G, where |N e |= 2n S −1. Gene trees are thus assumed to be contained within S (Fig. 1).
Site evolution is assumed to follow a continuoustime Markov process (Felsenstein 1981) under some substitution model and clock model: (2) i can adopt a range of molecular substitution models, such as the HKY nucleotide evolution model (Hasegawa et al. 1985) or the WAG amino acid evolution model (Whelan and Goldman 2001). Tree g i has relative molecular substitution rate i ∈ i . Branches in S are associated with substitution rates r, which govern the rate of site evolution of G along the respective branch, where |r|=2n S −1 (Fig. 1). Branch rates r are assumed to be independently and identically distributed under a log-normal distribution with standard deviation (i.e., the MSC relaxed clock model; Drummond et al. 2006;Ogilvie et al. 2017). Lastly, the clock rate can be estimated when accompanied by time-calibration data, such as ancient fossil records (Sauquet et al. 2011;Heled and Drummond 2012;Ballesteros and Sharma 2019), or left fixed when no such data are available. Overall, the total substitution rate of any given branch in g i is the product of i , , and a subset of the elements in r (weighted by their coverage of the gene-tree branch; Ogilvie et al. 2017).
In this article, we develop tools that allow the MSC to be applied to large data sets using complex models of evolution. Although we focus on MSC models, we anticipate that in the future other models of the form expressed in Eq. (1) will be developed, for example, models that allow some lateral gene transfer and therefore allow some gene-tree branches to cross species boundaries in the species tree. We design a number of MCMC operators which generate proposals that explore the state space more efficiently-using a Gibbs sampler for population sizes, a combination of Bactrian (Yang and Rodríguez 2013;Thawornwattana et al. 2018) and adaptable variance multivariate normal (Baele et al. 2017) proposal kernels, a parallel operator for sampling gene trees and substitution model parameters, and an MCMC operator which selects other operators based on their exploration efficiency (Douglas et al. 2021b). Moreover, in the special case of the multispecies relaxed clock model (Ogilvie et al. 2017), we introduce methods for operating on the species tree, the gene trees, and the clock model simultaneously (Zhang and Drummond 2020;Douglas et al. 2021b).

Effective Population Size Gibbs Operator
The StarBeast2 (Ogilvie et al. 2017) and DISSECT (Jones et al. 2015) packages have the capability of integrating effective population sizes N e when using an inverse gamma distributed prior on N e , based on a technique introduced by  and detailed out by (Jones, 2017). This approach greatly reduces the state space. However, consequently the posterior Eq. 1 can no longer be broken down in a product over components over individual gene trees: Thus, the technique is not suitable for gene-tree operator parallelization, and therefore, we estimate N e instead.
Suppose that N e b ∈ N e , for species-tree branch b, follows an inverse gamma prior distribution Inv-( N , N ), where the shape N is fixed at 2 and therefore the scale N is the expected value (because E(N e b ) = −1 ). Following the results by (Jones, 2017), the posterior of N e b follows an inverse gamma Inv-( N , N ), such that N = N +a and N = N +c where a is the total number of coalescent events of all gene trees in branch b and Here, p j is the ploidy of gene g j , c jbi the size of the ith coalescent interval for gene g j in branch b, and n jb the number of lineages of gene tree g j at the tip-side of branch b (so that n jb −i is the number of lineages at the start of the ith coalescent interval for g j ).
Instead of integrating N e , our GibbsPopulation operator samples from the posterior. All 2n S −1 elements in N e are proposed simultaneously. As demonstrated later, this turns out to be more efficient than standard N e random walk operators, with the added advantage of sampling effective population sizes-which may be a parameter of interest-as well as the ability to parallelize gene-tree proposals. This technique is readily applicable for periodically sampling and logging N e to implementations that do integrate this term out.

Bactrian Operators for Trees
The step size of a proposal kernel should be such that the proposed state x is sufficiently far from the current state x to explore vast areas of parameter space, but not so far that the proposal is rejected too often (Gelman et al. 1997). The Bactrian distribution (Yang and Rodríguez 2013;Thawornwattana et al. 2018) has minimal probability mass around the center, and a higher concentration flanking the center, akin to the humps of a Bactrian camel ( Fig. 2; left). This distribution is a preferred alternative to standard uniform-or normaldistributed random walk kernels, as it places minimal probability on step sizes that are too large or too small, and has successfully improved phylogenetic inference in VOL. 71 FIGURE 2. Depiction of random walks → under varying proposal kernels. Left: The random walk occurs from the origin between the two modes, where the vertical axis shows the probability density function of the kernel p( |) (Yang and Rodríguez 2013). Right: A 2D random walk on inversely correlated parameters = ( 1 , 2 ) with different domains (Baele et al. 2017). Contours describe the joint probability density function p( 1 , 2 ) under a transformed multivariate normal distribution learned during MCMC.
In this article, we apply Bactrian proposals to trees. The standard set of tree node height proposals in BEAST 2 consists of a Scale operator which embarks all nodes in the tree on a random walk (in log-space), a RootScale operator which does so for only the root of a tree, an UpDown operator which changes species/gene node heights and various continuous parameters simultaneously (Drummond et al. 2002), a SubtreeSlide operator which slides a node up or down branches (Hohna et al. 2008), and constant distance operators when a relaxed clock model is applied (Zhang and Drummond 2020). Each operator would normally draw a random variable from a uniform distribution, but here we instead use a Bactrian distribution and apply appropriate transformations. We also introduce the Interval operator, which transforms parameters with lower-and upper-bounds (such as tree node heights) by applying a Bactrian random walk in their real-space transformations.

Adaptive Variance Multivariate Normal Operator
An adaptive variance multivariate normal (AVMN) operator (Baele et al. 2017) provides proposals for a set of real-space parameters by learning the posterior throughout the run of the MCMC algorithm and approximating it as a multivariate normal distribution to capture correlations between parameters ( Fig. 2; right). The space spanned by such a set of continuous parameters may need to be transformed (in order to satisfy the assumption that all parameters lie in realspace), by applying a log-transformation to parameters with positive domains (such as substitution rates), or a log-constrained sum transformation to multivariate parameters with unit sums (such as nucleotide frequencies), for instance. AVMN has been demonstrated to be more efficient in estimating phylogenetic parameters than standard random walk or scale operators (Baele et al. 2017;Bouckaert 2020;Douglas et al. 2021b).
Consider a single gene tree g i and its substitution model i , consisting of substitution rates and nucleotide frequencies for instance. Performing a single proposal for any single parameter would require a full recalculation of the tree likelihood p(D i |g i , i ,r,) (see peeling algorithm by (Felsenstein, 1981)). Therefore, proposing all site model parameters i simultaneously can reduce the number of likelihood calculations required and thus lower the computational runtime.

Parallel Gene-Tree Operator
During MCMC, operators are typically sampled proportionally to fixed weights (or proposal probabilities), to ensure the chain is ergodic. Here, we present an alternative method, where a single gene tree g i and its substitution model i is selected, and N p operators are sequentially sampled and applied to g i and i , before returning to the full parameter space. This is equivalent to running a small MCMC chain of N p steps-applying only gene tree and substitution model operators on g i and i -and then accepting the resulting g i and i afterwards with probability 1, as if it were a single Gibbs sampling operation (Geman and Geman 1984).
Observe that because only g i and its associated parameters change, part of Eq. (1) can be rewritten as: Thus, the posterior distribution can be decomposed into the product of contributions of individual gene trees and their substitution models. Assuming that substitution model parameters i are distinct for each gene tree g i , an N p -step MCMC chain could be run for each of g i and g j for (i = j) in parallel, and the resulting g i and g j each accepted with probability 1, as if two Gibbs operators were sequentially applied. Because the posterior density for g i is proportional to p(D i |g i ,) p(g i |S,) and that of g j proportional to p(D j |g j ,) p(g j |S,) then provided that any shared parameters (such as r, S, and N e ) are not being operated on, these two N p -step MCMC chains can run in parallel.
Where there are N t threads available, the k gene trees are split into N t groups (assuming k ≥ N t ). The N t sets of N p -step MCMC chains are run in parallel and the resulting gene trees g are accepted into the main MCMC chain. Here, we introduce a parallel operator ParallelGeneTreeOperator(G, ). This operator partitions gene trees into N t threads and operates on their topologies, node heights, and substitution models. Tree node height proposals employ the Bactrian kernel where applicable (Fig. 2), and substitution model proposals invoke the AVMN kernel (Fig. 2). The chain length N p of each thread is learned during MCMC (Fig. 3).
Since each small MCMC chain for a thread can be considered a single Gibbs proposal, for N t threads in principle N t steps should be added to the main chain. If the operator is selected just before logging a state, in principle some threads may need to be disregarded before logging in order to maintain exactly equal intervals in the trace log. Due to the low frequency at which the operator is selected, and the logging intervals being orders of magnitude larger than the number of threads, this does not appear to be a problem in practice.

Species Tree Relaxed Clock Model Operators
The constant distance operator family exploits the negative correlations between divergence times and branch substitution rates by proposing both terms simultaneously (Zhang and Drummond 2020). This technique has yielded a parameter convergence rate of one to two orders of magnitude faster, particularly for large data sets that come with peaked posterior distributions (Douglas et al. 2021b). Under the MSC relaxed clock model used by StarBeast2, the branch rate of gene-tree branch b is the length-weighted branch rate r of all species-tree branches that contain b (Ogilvie et al. 2017). Moreover, effective population sizes N e are positively correlated with divergence times, so this correlation could also be readily exploited.
Extending the work by (Zhang and Drummond, 2020), we introduce the ConstantDistanceMSC operator. This operator proposes a node height t X for speciestree internal node X, the three branch rates (elements of r) and population sizes (elements of N e ) incident to X, and heights for all gene-tree non-leaf nodes that are contained within these three incident branches (Fig. 4). t X is embarked on a Bactrian random walk (Yang and Rodríguez 2013) to give t X , then r and the node heights in G are proposed such that all genetic distances are conserved following the change in t X , and N e is proposed such that the positive correlation between itself and the branch lengths incident to X is respected (see Algorithm S1).
Previously, we introduced the narrow exchange rate (NER) operator (Douglas et al. 2021b). This operator combined the simple NarrowExchange operator (i.e., a proposal which swaps a subtree with its uncle subtree; (Drummond et al., 2002)) with the ConstantDistance operator (Zhang and Drummond 2020), by applying a small topological change to the tree and then recomputing branch substitution rates such that evolutionary distances are preserved. We demonstrated that this operator assisted the traversal of tree topology space on longer alignments compared with shorter ones.
Here, we combine this work with the CoordinatedExchange operator implemented by (Ogilvie et al., 2017)-based on work by (Jones, 2017) and (Rannala and Yang, 2017)-and introduce the coordinated narrow exchange rate (CNER) operator. This operator exchanges a species-tree node with its uncle node adjusts gene-tree topologies g to preserve compatibility with S, and proposes three nearby branch rates in r to preserve genetic distances (Algorithm S2).

Adaptive Operator Weighing
Previously, we developed the AdaptableOperator Sampler(x) operator (Douglas et al. 2021b). This operator learns the weights (or proposal probabilities) behind a set of suboperators during MCMC, by rewarding operators which bring about large changes to parameter x in short computational runtime, with respect to some distance function: Euclidean distance when x is real, and RNNI distance (Collienne and Gavryushkin 2021) when x is tree topology. This approach can account for the scenario when an operator's performance is conditional on the data set. When a data set contains very little signal with respect to a certain parameter x and its prior distribution, then resampling that parameter from its prior distribution using the SampleFromPrior(x) operator may be more efficient than embarking x on a random walk, for instance (Douglas et al. 2021b). In contrast, data sets with more signals are likely to prefer smarter operators which account for correlations in the posterior distribution, such as the constant distance or NER operators (Zhang and Drummond 2020;Douglas et al. 2021b).
Here, we have applied the AdaptableOperator Sampler to seven areas of parameter space: the species and gene-tree node heights (t s and t G ), the relaxed clock model rates r and standard deviation , the mean effective population size N , the species-tree birth rate (assuming a Yule speciation model; Yule 1925), and FIGURE 3. Optimization of gene-tree parallel operator chain lengths. Top: The time limit of each parallel MCMC chain is randomized on each call so that the overhead (intercept) and time-per-proposal (slope) can be learned as a linear regression model. Bottom: The linear regression model is applied, and parallel MCMC chain lengths are set such that the slowest thread attains the user-specified target overhead (i.e., the bottom thread has attained 20% overhead in the example above).
the species-tree topology T S . These operator schemes are explicated in Tables 1 and 2.

RESULTS
In this section, we first validate the correctness of StarBeast3 through a well-calibrated simulation study. Then, we demonstrate that StarBeast3 is efficient at doing Bayesian inference on large data sets compared with StarBeast2. We did not compare to *BEAST directly, since it does not provide relaxed clock models on species trees, but note that (Ogilvie et al., 2017) benchmarked StarBeast2 against *BEAST for strict clocks and found StarBeast2 to be an order faster than *BEAST, so any gain over StarBeast2 will be more so over *BEAST.

Validation
In order to validate the correctness of StarBeast3, we performed two well-calibrated simulation studies. These were achieved by simulating nucleotide alignments (of two varying sizes) using parameters directly sampled from the prior distribution, and then recovering the posterior estimates of these parameters by doing Bayesian inference on the simulated alignments using StarBeast3. For each study, the 95%-coverage of each parameter was approximately 95% (meaning that the true parameter estimate was within the 95% highest posterior density interval approximately 95% of the time). Therefore, these experiments provide confidence in StarBeast3's correctness and are presented in Figure 5 and Section S4 of Supplementary material available on Dryad at http://dx.doi.org/10.5061/dryad.f1vhhmgzk.

Performance Benchmarking
We evaluated the performance of StarBeast3 for its ability to achieve multispecies coalescent parameter convergence in a Bayesian framework, compared with FIGURE 4. An example of an ConstantDistanceMSC proposal, acting on species nodes X and its two children L and R. First, the height of X (t X ) is increased to t X . Then, the relative substitution rates of branches L (r L ) and R (r R ) are decreased to r L and r R , and r X is increased to r X . These compensations in branch length ensure that the genetic distance of each branch (X, L, and R) is maintained. The thicknesses of the species node lines are proportional to these substitution rates. Finally, the effective population size of L and R are increased to N e L and N e R , while that of X is decreased to N e X . These compensations in node height ensure that the ratio between branch length and branch population size are maintained. Species node widths are proportional to their effective population size. During this operation, gene-tree nodes always remain constrained by the species tree. Figure was generated by UglyTrees (Douglas 2020). that of StarBeast2. Although it is a nontrivial problem to determine if an MCMC chain has converged, the effective sample size (ESS) can serve as a useful metric. Thus, we computed the number of effective samples generated per hour (ESS/h) across multiple replicates of MCMC, using three real and two simulated data sets ( Table 3). The ESS of any parameter should be over 200 in order to estimate its posterior distribution (Tracer; Rambaut et al. 2018). To allow both software packages to perform at their best, effective population sizes were integrated by StarBeast2, but were estimated by StarBeast3. This section provides a general comparison of StarBeast3 and StarBeast2; however, the performances of individual operators can be found in Sections S5 and S6 of Supplementary material available on Dryad.
The ESS/h was evaluated in five distinct areas of parameter space. First, we considered generic summaries of convergence: the ESS/h of the posterior density p(|D), the likelihood p(D|), and the prior density p(). Second, species tree S convergence was evaluated in terms of its height h S , its length l S , and hyperparametersthe Yule model birth rate (Yule 1925)-and N -the mean effective population size. In the case of StarBeast3, where effective population sizes are estimated, we also measured the mean ESS/h associated with species-tree leaf nodes of N e . Third, gene-tree convergences were evaluated by their heights h G , their lengths l G , and the RNNI distances (Collienne and Gavryushkin 2021) to their UPGMA D UPGMA (Sokal 1958) and neighborjoining D NJ trees (Saitou and Nei 1987). As there are multiple gene trees, we only considered the mean ESS/h of each term. Fourth, substitution model convergence (HKY substitution model; Hasegawa et al. 1985) was measured from the transition-transversion ratio , nucleotide frequencies f , and gene-tree substitution rates , where the ESS/h of each term was averaged across all k substitution models. Lastly, relaxed clock model convergence was evaluated by considering the mixing of branch rate empirical meanr and variance var(r), as well as the relaxed clock standard deviation parameter .
These results showed that, depending on the data set, the "slowest" parameter generally converged considerably faster for StarBeast3 than it did for StarBeast2 (see the min term in Figs. 6 and 7). On the smallest data set considered (Frog), StarBeast2 and 3 performed comparably well overall (and no significant difference in min). However, StarBeast3 performed better on all of the other data sets, with the "slowest" parameter converging between 4 and 37× as fast, and the posterior density p(|D) converging between 2 and 36× as fast, often at a statistically significant level. For StarBeast3, the absolute time needed to converge varied a lot across the data sets, and even across multiple replicates of the same data set (see final column of Table 3). The fastest data sets -Frog and Simulated(12)-required 1-2 days to converge, while the Spider data set required over a month.
Notably, relaxed clock model parameters converged up to 35× as fast under StarBeast3. This was credited to the use of a real-space branch rate parameterization (where branch rates are real numbers as opposed to discrete bins, as implemented in StarBeast2) as well as constant distance operators, which adjust branch rates and divergence times simultaneously (Zhang and Drummond 2020;Douglas et al. 2021b). The disparity between StarBeast3 and StarBeast2 was less extreme for the smaller k = 26 gene tree Frog data set (Barrow et al. 2014), consistent with previous experiments (Douglas et al. 2021b).
Substitution model parameters generally converged faster for StarBeast2 than they did for StarBeast3. Note, Bactrian Operators for Trees

Species Tree Relaxed Clock Model Operators
Ogilvie et al. (2017) Gene trees/site models ParallelMCMC(T G ,t G , ) 3.42 Table 2 Tree hyperparameters GibbsPopulation(N e ) 50 Notes: The ParallelMCMC operator weight was set such that it is sampled 1% of the time. Further operator details can be found in Drummond and Bouckaert (2015). a Bactrian kernel applied to random walk (Yang and Rodríguez 2013). however, that this is by design. The total operator weight assigned to parameters was 50% smaller in StarBeast3, in order to ensure balanced convergence across all areas of parameter space. In all data sets considered, substitution models converged significantly faster than any other area of parameter space, despite receiving relatively little operator weight, and therefore computational resources that were being spent on the substitution model were better off spent in "slower" areas of parameter space, such as gene-tree node heights.

Effective Population Size Gibbs Operator
The AdapableOperatorSampler operators (Table 1) confirmed the value in the NER and ConstantDistanceMSC operators for operating on their respective areas of parameter space. The ConstantDistanceMSC operator almost always outperformed other operators at proposing species node heights t S (Table 4). The exception to this was the Skink data set, for which the UpDown operator was superior at proposing branch lengths, and the Frog data set, for which ConstantDistanceMSC, CoordinatedExponential, and UpDown were all on a par. In general, very little operator weight was rewarded to the Uniform, Interval, TreeScale CoordinatedUniform, and CoordinatedExponential operators for their abilities to propose species node heights. Similarly, among NarrowExchange variants evaluated by AdaptableOperatorSampler(T S ), the NER operator was marginally favored by all data sets (Table 5). This was due to the operator making larger or more StarBeast3 parallel operator scheme for gene trees and their associated site models (assumed to be an HKY model with transition-transversion ratio and nucleotide frequencies f )

Operator
Weight Reference Bactrian Operators for Trees ∀ i∈1,...,k AdaptableOperatorSampler(t g i ) 100 Notes: Each operator is applicable to a single gene tree g i or its site model i . AVMN( i ,t g i ) generated proposals for the site model and complete set of tree node heights simultaneously. Operator weights are normalized into proposal probabilities within a single MCMC chain called by ParallelMCMC. Further operator details can be found in (Drummond and Bouckaert, 2015). a Bactrian kernel applied to random walk (Yang and Rodríguez 2013). frequent topological changes to the species tree, in faster computational runtime, especially compared with CoordinatedNarrowExchange and CNER.
Overall, this experiment reinforced the value of learning operator weights on a problem-by-problem basis. A full breakdown of the remaining four adaptive operators can be found in Section S6 of Supplementary material available on Dryad. Lastly, we evaluated the effect of threading on StarBeast3, by comparing its performance under 1, 2, 4, 8, and 16 threads allotted to the ParallelMCMC gene-tree operator (Fig. 8). There was a positive-but-modest correlation between the number of threads and the overall rate of convergence among the terms considered, with an overall log-linear slope coefficient of 0.19. This can be interpreted as follows: across the range of threads and data sets considered, doubling the number of threads was associated with an increase in mixing by 14%. Multithreading provided the strongest boost for the Skink and Spider data sets and made little difference to the simulated data set (48 taxa). This is an unexpected result, because the Skink and Spider data sets have fewer genes (k = 50 compared with k = 100), and may be due to the former data sets having more taxa and thus larger trees.

Benchmarking on Large Data Sets
We benchmarked the performance of StarBeast3 on simulated data sets with increasingly large numbers of gene trees k, ranging from 250 to 1000 genes. Each gene was 200 nucleotides in length. In order to achieve convergence in a timely manner, we performed inference under a strict clock model (i.e., all branch rates fixed r = 1) and with a small sample size (n s = 4 species n G = 12 taxa). These experiments showed that StarBeast3 required more time to produce one sample for larger k, and therefore more time to produce one effective sample, as expected (Fig. 9). The k = 250 gene data set would require ≈ 10 h for the average ESS to exceed 200 in all areas of parameter space, while the k = 1000 gene data set would require ≈ 40 h. Furthermore, we confirmed that gene-tree parallelization gave a noticeable-but-modest improvement to runtime (Fig. 9). Although the trees were small, this experiment showed that StarBeast3 is indeed capable of running on large data sets with several hundred genes.

The Next Generation of Bayesian MCMC Operators
In recent years, Bayesian MCMC proposals have advanced significantly beyond that of the unidimensional random walk. The use of adaptive algorithms and advanced proposal kernels have become increasingly prevalent (Haario et al. 2001;Vihola 2012;Yang and Rodríguez 2013;Benson and Friel 2018). In phylogenetic inference in particular, tree proposals have been guided by conditional clade probabilities and parsimony scores (Höhna and Drummond 2012;, and mirror kernels learn target distributions which act as "mirror images" (Thawornwattana et al. 2018), for instance.
Here, we introduced a range of recently developed MCMC operators to the MSC, including Bactrian proposal kernels (Yang and Rodríguez 2013), which have been successfully applied to bird phylogeny (Maliet et al. 2019), and tree "flex" operators (BICEPS; Bouckaert 2021), which have been applied to coronavirus disease-2019 genomic data (Douglas et al. 2021a). We also FIGURE 5. Well-calibrated simulation study analyzing n S = 16 species, n G = 48 taxa, and k = 50 genes. One-hundred simulations were performed to recover the coverage between "true" simulated values and their estimates under the posterior distribution. 95% highest posterior density (HPD) intervals of parameters are represented by vertical lines. Each line represents a single simulation, and is colored blue when the true value was contained within the 95% interval, or red otherwise. The top of each plot shows the coverage of each parameter (i.e., the number of MCMC simulations for which the "true" parameter value was contained within the 95% HPD).
invoked a series of more meticulous operators which account for known correlations, such as the AVMN kernel (Baele et al. 2017), constant distance operators (Zhang and Drummond 2020), and the NER operator (Douglas et al. 2021b), as well as adaptive operators that improve over the course of MCMC, such as the adaptable operator sampler (Douglas et al. 2021b), parallel gene-tree operators, and the AVMN kernel Notes: Fifty gene trees were subsampled from the Skink and Spider data sets. The simulated data sets were directly sampled from the model specification used during Bayesian inference (described in Section S3 of Supplementary material available on Dryad). In the final column, we estimate the time required for the MCMC chain to converge using StarBeast3 with 16 threads (min-max across 5 MCMC replicates). These terms were approximated as the time to achieve an effective sample size of 200 for the posterior density p(|D), with a 20% burn-in.
FIGURE 6. Performance benchmarking the two simulated data sets. Each point is the geometric-mean ESS/h across five replicates, for either StarBeast2, or StarBeast3 with 16 threads. The geometric-mean relative performance of StarBeast3, compared with StarBeast2, is indicated above each term, and a * is present if the difference across five replicates is significant according to a Student's t-test. Note that the y-axis is in log-space. VOL. 71 FIGURE 7. Performance benchmarking the two biological data sets. See Figure 6 caption for figure notation. (Baele et al. 2017). Indeed, these operators have yielded a software package which outperforms StarBeast2 by up to one-and-a-half orders of magnitude, depending on the data set and the parameter.
While StarBeast3 provides a clear advancement to the problem, Bayesian MCMC is still lagging behind the volumes of next-generation genomic data. Therefore, the continued development of efficient, meticulous, and adaptive MCMC operators is essential.

Efficient Parallelized Bayesian Inference under the MSC
As genomic data becomes increasingly available, concatenating genomic sequences and inferring the phylogeny of the species as that of the genes can Notes: Note that the timer starts at the beginning of the proposal and ends when the proposal has accepted or rejected. NE = narrow exchange; NER = narrow exchange rates; CNE = coordinated narrow exchange; CNER = coordinated narrow exchange rates. The operator which was rewarded the highest proposal probability for each data set is indicated by a .
FIGURE 8. Effect of threading on StarBeast3 performance. Each point represents the ESS/h of the posterior density P(|D) (averaged across five replicates), for the indicated thread count and data set. These terms are normalized to enable comparison across data sets, by dividing it by that of one thread. A linear model was fitted to the ESS/h and number of threads, each in log 2 space, and is reported at the top of the plot. The positive coefficient of the slope indicates that performance increased with the number of threads, across the range of threads considered. Parallel MCMC chain lengths were optimized using the adaptive scheme presented in Figure 3. become enticing. However, this approach makes for an inconsistent estimator of topology when divergence times are small (Pamilo and Nei 1988), and a biased estimator of species divergence times and substitution rates when ILS is present (Arbogast et al. 2002;Mendes and Hahn 2016;Ogilvie et al. 2016). MSC methods address these issues, but at the drawback of their demanding computational runtimes.
Therefore, as multithreading technologies become increasingly affordable, the appeal in parallelizing multispecies inference becomes clear. StarBeast3 exploits the assumption of conditional independence between gene trees, by doing Bayesian inference on gene trees in parallel, and therefore it scales with the size of the problem. StarBeast3 can handle large data sets (with hundreds of genes) and achieve convergence several times faster than its predecessors.

A Balanced Traversal Through Parameter Space
All areas of parameter space should be explored approximately evenly during MCMC. If one area of parameter space is being explored more rapidly than another, then computational resources allotted to the former should be diverted to the latter. This is best exemplified by the phylogenetic substitution model which, VOL. 71 despite requiring relatively little attention to converge, still requires full recalculation of the tree likelihood upon every proposal (Felsenstein 1981). Conversely, tree topologies often converge rather poorly and can require significant attention to be rescued from local optima. By fine tuning our MCMC operator proposal probabilities, we have achieved a balanced traversal through all areas of the MSC parameter space. Although some parameters converge slower for StarBeast3 than they do for StarBeast2 (such as those in the substitution model), the slowest parameters converge significantly faster for the former; up to 37× as fast (see the min term in Figs. 6 and 7).
For StarBeast3, we employed adaptable operators which are able to learn the proposal probabilities of other operators based on their ability to explore a single area of parameter space (Douglas et al. 2021b). However, there would be a great benefit in an adaptable operator scheme which learns and applies a balanced exploration across different areas of parameter space on a problem-by-problem basis.

CONCLUSION
Here we introduce StarBeast3-a software package for performing efficient Bayesian inference on genomic data under the MSC model. We verified StarBeast3's correctness and we benchmarked its performance against StarBeast2, which is an order of magnitude faster than its still popular predecessor *BEAST. We showed that StarBeast3 is significantly faster than StarBeast2. Notably, relaxed clock parameters converged between 3 and 30× faster, but most importantly even the "slowest" parameters converged up to 36× faster. Our adaptive operator scheme allows proposal probabilities to be learned on a problem-by-problem basis, making StarBeast3 suitable for a range of data sets. By estimating effective population sizes (instead of analytically integrating the term out), we were able to parallelize gene-tree proposals and demonstrated that doubling the number of allotted threads was associated with an increase in performance by around 14%. StarBeast3 is highly effective at performing fast Bayesian inference on large data sets with over 100 genes.