Inference of Gene Flow between Species under Misspecified Models

Abstract Genomic sequence data provide a rich source of information about the history of species divergence and interspecific hybridization or introgression. Despite recent advances in genomics and statistical methods, it remains challenging to infer gene flow, and as a result, one may have to estimate introgression rates and times under misspecified models. Here we use mathematical analysis and computer simulation to examine estimation bias and issues of interpretation when the model of gene flow is misspecified in analysis of genomic datasets, for example, if introgression is assigned to the wrong lineages. In the case of two species, we establish a correspondence between the migration rate in the continuous migration model and the introgression probability in the introgression model. When gene flow occurs continuously through time but in the analysis is assumed to occur at a fixed time point, common evolutionary parameters such as species divergence times are surprisingly well estimated. However, the time of introgression tends to be estimated towards the recent end of the period of continuous gene flow. When introgression events are assigned incorrectly to the parental or daughter lineages, introgression times tend to collapse onto species divergence times, with introgression probabilities underestimated. Overall, our analyses suggest that the simple introgression model is useful for extracting information concerning between-specific gene flow and divergence even when the model may be misspecified. However, for reliable inference of gene flow it is important to include multiple samples per species, in particular, from hybridizing species.


Introduction
Hybridization can enhance variation in recipient species, and has long been recognized as an important process in plants that can stimulate the origin of new species (e.g., Anderson 1949;Mallet 2007). Analyses of genomic data in the past decade have highlighted the prevalence of hybridization or introgression in animals as well, including bears (Liu et al. 2014;Kumar et al. 2017), birds (Ellegren et al. 2012), and butterflies (Martin et al. 2013). Between-species gene flow may involve either sister or non-sister species and may play an important role in ecological adaptation (Mallet et al. 2016;Martin and Jiggins 2017). Gene flow can be a major contributor of genealogical variation across the genome and gene tree-species tree discordance, in addition to ancestral polymorphism or delayed coalescence (Maddison 1997;Nichols 2001).
There is a long history of studies in population genetics of models of population subdivision and migration (Wright 1943;Malecot 1948;Slatkin 1987), and a number of methods have been developed to estimate the migration rate between populations (Beerli andFelsenstein 1999, 2001;Bahlo and Griffiths 2000). An important limitation of models of population subdivision, when applied to data from different species or subspecies, is that they do not account for the divergence history of the populations or species. Introducing a population/ species phylogeny into models of population subdivision not only improves the realism of the model but also opens up opportunities for addressing a number of interesting questions in evolutionary biology, such as estimating species divergence times and ancestral population sizes, delineating species boundaries, and inferring the direction, rate, and timing of gene flow .
Two classes of models of gene flow have been developed that accommodate the phylogeny of the species, both of which are extensions of the multispecies coalescent (MSC) model (Rannala and Yang 2003). The first is the MSC-with-migration model (MSC-M, or isolation-with-migration or IM model, Hey and Nielsen 2004;Hey 2010;Zhu and Yang 2012;Dalquen et al. 2017;Hey et al. 2018), which assumes that two MBE species exchange migrants at a certain rate over an extended time period. The rate of gene flow from species A to B is measured by the proportion of migrants (m AB ) in the receiving population B or by the population migration rate, M AB = N B m AB , the expected number of immigrants from A to B per generation, where N B is the (effective) population size of species B. We note that the isolation-with-initial-migration (IIM) model of Costa and Wilkinson-Herbots (2017), which assumes that gene flow occurs initially after species divergence but stops after a period of time when reproductive isolation has been fully established, is an instance of the MSC-M or IM model (see below). The second class of models of gene flow is the MSC-with-introgression (MSci) model , also known as multispecies network coalescent model (MSNC; Wen and Nakhleh 2018;Zhang et al. 2018), which assumes that gene flow occurs at fixed time points in the past. The rate of gene flow is measured by the introgression probability (φ or γ), which is the proportion of successful immigrants in the population at the time of introgression.
In the real world, introgressed alleles may be removed by natural selection because they are involved in hybrid incompatibility and are deleterious in the genetic background of the recipient population (Dobzhansky 1937;Muller 1942) or because they are linked to such loci (Petry 1983; Barton and Bengtsson 1986;Uecker et al. 2015). Thus the rate of gene flow (M in MSC-M or φ in MSci), when those models are used to analyze genomic sequence data, reflect the long-term effects of selection and drift as well as hybridization or introgression (Martin and Jiggins 2017). Such an effective rate of gene flow may be expected to vary across the genome, influenced by the presence of loci in the genomic region important in ecological adaptation and by the local recombination rate (Bürger and Akerman 2011;Aeschbacher and Bürger 2014;Akerman and Bürger 2014;Schumer et al. 2018;Edelman et al. 2019;Martin et al. 2019). The rate may also vary over time, depending on geological or ecological events that cause changes in the ecology and distribution of the species and in the chance for two species to exchange genes. One can envisage models of gene flow in which the rate varies over time and across genomic regions. For the present, such extended models are not yet implemented in the MSC framework, and the feasibility of fitting such parameter-rich models to genomic datasets is unexplored. MSC-M and MSci models implemented to date (Dalquen et al. 2017;Hey et al. 2018;Wen and Nakhleh 2018;Zhang et al. 2018;Flouri et al. 2020) assume constant rates, and should be considered first approximations when applied to genomic sequence data.
In this paper, we use mathematical analysis and computer simulation to examine the impact of model misspecification on estimation of parameters under the MSci model, such as species divergence and introgression times, population sizes, and introgression probabilities. We use the Bayesian program Bayesian phylogenetics and phylogeography (BPP) (Flouri et al. 2018 to analyze multilocus sequence data simulated under various MSci and MSC-M models. Although BPP is our own implementation of the MSci model, our results should apply to similar exact or likelihood methods (Wen and Nakhleh 2018;Zhang et al. 2018). Our results may also apply to approximate methods, which use summaries of the data such as the genome-wide sitepattern counts (as in the D-statistic, Green et al. 2010 andHYDE, Meng andKubatko 2009;Blischak et al. 2018), reconstructed gene trees (as in SNAQ, Solis-Lemus and Ane 2016), or other summary statistics used in Approximate Bayesian Computation (Dittberner et al. 2022). However, approximate methods do not make a full use of information in the data and may not identify all parameters in the model. For example, the D-statistic is agnostic of the mode of gene flow (migration versus introgression) and cannot be applied to data sampled from only two species or populations. The computational strengths and statistical weaknesses of approximate methods have been discussed by a number of authors (Degnan 2018;Elworth et al. 2019;Zhu and Yang 2021;Hibbins and Hahn 2022;Ji et al. 2022). In contrast, likelihood methods integrate over all possible gene trees underlying the sequence alignments, making use of all information about the model and parameters in the sequence data. They typically involve a heavy computational load. However, recent algorithmic improvements have made it possible to apply the MSci model to genome-scale datasets with more than 10,000 loci ). Inferring introgression events or constructing an introgression model using genomic sequence data, however, remains a challenging task, even when a binary species tree is specified, onto which introgression events can be added (Ji et al. 2022;Thawornwattana et al. 2022); see Discussion for an overview of currently available methods for inferring gene flow on a species phylogeny. For these and many other reasons, the model of gene flow assumed in our data analysis may often be incorrect. An important question is to what extent inference of gene flow and estimation of the timing and rate of gene flow can still be achieved when the model of gene flow is misspecified. The impact of model misspecification on estimation of other evolutionary parameters such as species divergence times is also of major concern.
Although there are many ways in which the assumed model is wrong, we are particularly interested in a few types that are likely in real data analyses (Finger et al. 2022;Thawornwattana et al. 2022). First, gene flow may be occurring continuously during a time period but an MSci model is fitted to the genomic data, which assumes that gene flow occurred at a particular time point (e.g., Wen and Nakhleh 2018;Jiao et al. 2020). We are here interested in whether species divergence times and ancestral population sizes are affected by the misspecification, and how the migration rate in the migration model (M) corresponds to the introgression probability in the MSci model (φ). The case of two species is analytically tractable. We study the limit of the maximum-likelihood estimates (MLEs) of introgression probability and introgression time when the data size (the number of loci) approaches infinity when the data are generated under the MSC-M MBE model. We use computer simulation to verify and extend the analytical calculation.
Second, the introgression event may be assigned to a wrong branch on the species tree, for example, to a parental or daughter branch of the genuine introgression lineage. Alternatively, introgression may involve species that have since gone extinct or are not included in the data sample. The presence of such ghost species is known to mislead inference of the history of gene flow for the sampled species (Beerli 2004;Tricou et al. 2022). Thus we conducted simulation to examine the impact of unsampled species on the inference of gene flow. In general, our results demonstrate the usefulness of the simple introgression model in inferring gene flow using genomic sequence data. Following Jiao et al. (2020), we study the asymptotic behavior of Bayesian parameter estimation under the introgression (MSci) model when the data are generated under the migration (IM) model in the case of two species, with one sequence per species per locus ( fig. 1). Here we focus on this simple case because it is analytically tractable. Note that our Bayesian implementation in BPP  accommodates an arbitrary number of species and an arbitrary number of sequences per species per locus, and the likelihood calculation averages over the gene genealogy for the sequences at each locus. We assume an infinite number of loci, and the data at each locus consist of a pair of sequences (a, b) from the two species, with x differences at n sites. The coalescent time t for the locus is unknown and underlies the observed difference. Jiao et al. (2020) analyzed the MSC-M model ( fig. 1a) assuming infinite sequence length (n = ∞) so that the true coalescent time between the two sequences (t) is known. Here we accommodate random fluctuations in the number of mutations due to finite sequence length and consider three variants of the migration model.

Notation and Definition of Parameters
In the basic IM model ( fig. 1a), species A and B diverged at time τ R and there has since been gene flow from A to B at the rate of M AB migrants per generation. The IIM model ( fig. 1b) assumes that migration occurred initially after species divergence but stopped at time τ T > 0 (Costa and Wilkinson-Herbots 2017), and is represented by an MSC-M model for three species including a ghost species. Here the ghost does not necessarily represent a real species but is a mathematical device for specifying the IIM model. The IIM model becomes identical to the IM model when τ T = 0. We also consider a secondary contact (SC) model ( fig. 1c), in which two species initially had complete isolation but came into contact at a certain time point (τ T ) with ongoing gene flow at the rate of M AB ever since (Costa and Wilkinson-Herbots 2021). This is similarly specified using a ghost species at time point τ T (fig. 1c). The migration model involves three types of parameters: species divergence times (τ R , τ T ), population sizes for extant, and extinct species (θ A , θ B , θ T , θ R ), and the (population) migration rate M AB . The population size parameter for any species with (effective) population size N is defined as θ = 4Nμ, where μ is the mutation rate per site per generation. We refer to a branch on the species tree by its daughter node so that branch RA is also branch A, with population size parameter θ A . Both divergence times (τ) and population sizes (θ) are measured by the expected number of mutations per site.

Asymptotic Theory
We first consider the IIM model of figure 1b, of which the IM model of figure 1a is a special case with τ T = 0. The backwards-in-time process of coalescent and migration in time interval (τ T , τ R ) is described by a Markov chain with

MBE
three states: AB, AA, and A (Notohara 1990). Here AB is the initial state, with two sequences in the sample, one in A and another in B; AA means both sequences are in A (in other words, sequence b is traced back into A); and A means one sequence in A (in other words, sequence b is traced back into A and has coalesced with sequence a). Note that in the Markov chain, time runs backwards, so the transition from AB to AA means migration of a sequence from A to B in the real world. The generator matrix for the Markov chain is (see, e.g., Notohara 1990;Jiao et al. 2020) where w = m AB /μ = 4M AB /θ B is the mutation-scaled migration rate, and 2/θ A is the coalescent rate in population A, with one time unit being the expected time taken to accumulate one mutation per site. Q has eigenvalues λ 1 = 0, λ 2 = −2/θ A , and λ 3 = −w.
Let the transition probability matrix over time t be P(t) = {p ij (t)} = e Qt , where p ij (t) is the probability that the Markov chain will be in state j time t later given that it is in state i at time 0. This is The probability density of coalescent time t is thus This is a function of w = 4M AB /θ B but not of M AB and θ B individually. The parameters specifying the density are thus Θ iim = (w, θ A , θ R , τ R , τ T ). Note that the density under the IM model f im (t) is given by f iim (t) with τ T = 0 ( fig. 1b  and c).
Similarly under the secondary-contact (SC) model ( fig.  1c), the coalescent-with-migration process over the time interval (0, τ T ) is described by the Markov chain of equation (1). Given the parameters Θ m , the probability density of coalescent time t is Under the MSci model ( fig. 1d), we have (e.g., Jiao et al. 2020) This is a function of parameters Given the coalescent time t for a locus, the probability of observing x differences at n sites under the JC mutation model (Jukes and Cantor 1969) is given by the binomial probability The marginal probability of observing x differences at n sites, under both the migration (IM, IIM, SC) and introgression (MSci) models, is where f(t|Θ) is given by equations (3), (4), or (5). For analytical tractability of the likelihood (eq. 7), we assume the infinite-sites mutation model instead of JC, and replace the binomial likelihood by a Poisson approximation where the subscript "m" stands for any of the three MSC-M models ("im" for IM, "iim" for IIM, or "sc" for SC, fig. 1a-c). The KL divergence is a measure of distance from the fitting introgression model to the true migration model: here Θ m are fixed, whereas Θ i are being estimated.
The limiting values Θ * i as L → ∞ are also known as the pseudo-true parameter values for the misspecified MSci model. The BFGS optimization routine in PAML (Yang 2007) is used to minimize equation (9) to obtain the MLEs.
We are in particular interested in the introgression probability φ and the introgression time τ S . Note that under the migration model, the probability that any lineage from species B traces back to A is where Δτ is the time period of gene flow ( fig. 1a-c). Equation (10) gives the expected proportion of migrants under the true migration model. When M AB is small, Δτ, which is also given by equating the expected total number of migrants under the two models: is the expected number of migrants per generation and Δτ/μ is the number of generations with gene flow. It may be noted that the theory of equation (9) can be used to study the limiting parameter estimates (when L → ∞) in the migration model when the true model is the introgression model. One has only to flip the roles of f m (x | Θ m ) and f i (x | Θ i ) in equation (9). This is not pursued in this paper.

Asymptotic Results under the IM Model
We used the asymptotic theory (eq. 9) to obtain the MLEs (Θ We use five methods (a-e) to fit the MSci model, with method d estimating all five parameters, whereas the others have some parameters fixed ( fig. 2). We examined the effects of the sequence length (n) and the migration rate (M AB ). Note that five parameters are identifiable under the MSci model: . 1d), and θ A , θ B , θ H are unidentifiable as no coalescent events can occur in those populations given one sequence per species per locus. Population size θ S is identifiable as it is possible for both sequences a and b to be traced back to population S. Nevertheless, one expects the information concerning θ S to be weak in datasets of two sequences per locus. In methods c and d, θ S and φ are estimated as free parameters. Application of the misspecified MSci model (to data generated under the IM model) led to unreasonably large estimates of θ S (as large as 0.5 mutations per site), and the poor estimates of θ S caused φ to be poorly estimated as well. This is due partly to our use of one sequence per species per locus and partly to the confounding effects between θ S and φ. We discuss both effects later when we describe the simulation results. Here we focus on methods a, b, and e, in which θ S is fixed at the true value θ 0 (in methods a and b) or constrained to be equal to θ R (in method e).
In the IM model, migration occurs throughout the time interval (0, τ R ), at the rate of M AB migrants per generation ( fig. 1a). When such data are analyzed under the introgression model, a simple expectation might be that the introgression time τ S should be the average τ R /2, whereas the introgression probability might be given by the expected proportion φ 0 of equation (10). However, as we show below, this expectation is too simplistic.
First, we discuss the introgression time τ S , assuming the true coalescent time (or n = ∞). Given the datagenerating IM model, there is a strictly positive probability for 0 < t < ϵ for any small constant ϵ > 0 (supplementary fig. S1, Supplementary Material online). In other words, there must exist loci at which t is arbitrarily close to 0. In the MSci model, sequences a and b cannot coalesce until they are in the same population S, so that τ S < t. When the MSci model is fitted to data generated under the IM model, τ S is dominated by the minimum rather than the average coalescent time, and τ S → τ * S = 0 when the number of loci L → ∞ (and when the true coalescent time t is known). Even though migration occurs throughout the time interval (0, τ R ), the MSci model has to lump all migration events to one time point, τ * S = 0 ( fig. 2b and e). With finite sequences (n < ∞), t is not observed and is reflected in the number of mutations (x). Whatever the true t, there is a positive probability of observing no mutations between the two sequences, so that an absence of mutations (x = 0) is not strong evidence for t = 0. The MLE τ * S reflects not only the minimum coalescent time, but also the whole distribution (supplementary fig. S1, Supplementary Material online). Thus τ * S > 0, different from the case where the coalescent time is known without error (n = ∞). Nevertheless, one expects τ * S to be closer to 0 than to τ R , especially if the number of sites is large. Indeed in our calculations, τ * S ≪ 1 2 τ R ( fig. 2b and e).
Next we consider the introgression probability φ and again focus on methods a, b, and e ( fig. 2). The estimate φ * increases nearly linearly when M AB is small (< 1 4 , say) but tails off at large M AB . All estimates are smaller than φ 0 of equation (10)  . We defer to a later section a detailed discussion of the estimation of φ, contrasting the IM, IIM, and SC models.
Finally the estimated divergence time between the two species (τ R ) matched the true values at low migration rates but was underestimated at high migration rates, with the ancestral population size θ R overestimated ( fig. 2). It may be tempting to interpret the underestimation of τ R (and overestimation of θ R ) by the MSci model as being due to the difficulty of distinguishing complete isolation with recent species divergence from introgression or of distinguishing migration and coalescent events close to species divergence from ancestral polymorphism. However, this does not appear to be a correct interpretation.
We examined the true and fitted distributions of the coalescent time (supplementary fig. S1, Supplementary Material online). If there is no migration (M AB = 0), the MSci model (with φ = 0) will be correct, and the parameter estimates will converge to the true values, with a perfect fit to the density f m (t). At low migration rates (M AB ≤ 0.1, say), the MSci model fits the density f m (t) very well, with the discontinuity point in the true and fitting distributions coinciding: τ * R = τ m R . At the intermediate rate of M AB = 1, the species divergence time τ R is still correctly estimated even though the fit to the density is poor (supplementary fig. S1, Supplementary Material online). At high rates (with M AB ≥ 1.4, say), the true density has a mode in the interval (0, τ m R ), dropping off at τ m R . The best fitting density starts from 0, with an exponential decay, and has a discontinuity point at τ R with again an exponential decay. This best-fitting density is a poor fit, and the discontinuity point τ * R is moved to smaller values as an attempt to accommodate the migration and coalescent events in the middle of the interval (0, τ m R ) to improve the fit (judged by the KL divergence). Thus τ R is underestimated (τ * R < τ m R ). As a result, the population size parameter θ R is overestimated, as those two parameters tend to be strongly negatively correlated (e.g., Burgess and Yang 2008). In other words, the intermediate coalescent times in the interval (0, τ R ), which occur at a large proportion of loci, are accommodated or misinterpreted by the MSci model using a smaller τ R and larger θ R . Coalescent times in the range τ * R < t < τ m R , which represent true migration events, are misinterpreted as coalescent events in species R, and φ * is much less than φ 0 (eq. 10). show similar patterns to those under the IM model discussed above. Similarly, θ S is difficult to estimate using two sequences per locus in methods c and d, and the poor estimates of θ S affects the estimation of φ. Thus we focus on methods a, b, and e, in which θ S is fixed or constrained, and on the introgression time and introgression probability. In the IIM model, migration events occur throughout the time interval (τ T , τ R ) ( fig. 1b), but the estimate of the introgression time is dominated or influenced by the minimum coalescent time, so that τ * S = τ T when n = ∞, and τ * S > τ T when n is finite. In the latter case, τ * S is much closer to τ T than to τ R (supplementary fig. S2, Supplementary Material online).
The introgression probability φ * grew almost linearly with M AB when M AB was small (with M AB ≤ 0.2, say), and this estimate was close to the expectation φ 0 of equation (10)  The results show patterns similar to those under the IM and IIM models discussed above. The species divergence time under the MSci model τ * R = τ (sc) R when the migration rate M AB is small but drops at very high rates (with M AB > 2). The introgression time is dominated by the minimum coalescent time, so that τ * S = 0 when n = ∞, and τ * S is much closer to 0 than to τ T when n is finite (supplementary fig. S4, Supplementary Material online). Note that in the true model migration occurs throughout the time interval (0, τ T ).
The introgression probability φ * grew almost linearly with the migration rate M AB when M AB was small (with M AB ≤ 1 4 , say), and was close to the expectation φ 0 (eq. 10) when M AB < 2 (supplementary fig. S4a, b and e, Supplementary Material online). At very high rates (M AB > 2), φ * was much smaller than φ 0 , and this 'bias' was accompanied by an underestimation of τ R and overestimation of θ R . Similarly to the IM and IIM models discussed above, this is due to the attempt of the MSci model to accommodate the coalescent times in the middle of the interval (0, τ T ) (supplementary fig. S5, Supplementary Material online).

The Amount of Gene Flow under the IM, IIM, and SC Models
Although the expected total amount of gene flow measured by φ 0 (eq. 10) is the same under the IM, IIM, and SC models of figure 1a-c, the estimates under the MSci model differ, as summarized in figure 1e.
At low migration rates, τ R , θ S , and θ R in the MSci model are nearly accurately estimated to match those in the true model (figs. 2, supplementary figures S2 and S4, Supplementary Material online). Consider the case of infinitely long sequences with known coalescent time. Let and let the introgression time be τ * S = 0 for the IM and SC model, and τ * S = τ T for the IIM model. We also match the probability density of coalescent time t, with f i (t) = f m (t), for t > τ R . With those simplifying assumptions, φ * that minimizes the KL divergence (eq. 9) can be derived as At low migration rates, equation (11) (11), we have In other words, recent gene flow (as in SC) is easier to recover by the MSci model than ancient gene flow (as in IM or IIM). Note that φ * (IIM) = φ * (IM) holds only when one sequence is sampled per species; as there is no coalescent over (0, τ T ), IIM is essentially the same model as IM with a time shift ( fig. 1). This will not be the case when multiple sequences per species are sampled or when the sequence length is finite.

Simulation Results
As our asymptotic theory was limited to a single sequence per species per locus, we used simulation to verify and augment our analytical calculations above. We simulated data under the IM, IIM, or SC models of figure 1a-c, using the same parameter values as above, and analyzed them using BPP under the MSci model ( fig. 1d). The JC mutation model (Jukes and Cantor 1969) was assumed. In the basic setting, we used S = 4 sequences per species per locus, n = 1, 000 sites per sequence, and L = 4, 000 loci in each dataset, with the migration rate M AB = 0.2. We varied n, S, L, M AB MBE to examine their effects. With multiple sequences per species (S > 1), all eight parameters of the MSci model, 1d) are identifiable . The results are summarized in figure 3. We note a few common features first. In nearly all cases, population sizes for extant species (θ A , θ B ) were very well estimated, with posterior means close to the true values and with very narrow highest-probability-density (HPD) credibility intervals (CIs). The exception was parameter θ B under the IM model (note that B is the species receiving immigrants), which was less well estimated when the dataset was small and had either short sequences (n = 250) or few loci (L ≤ 500), or when the migration rate was very high. The poorer estimation of θ B appeared to be related to the underestimation of φ and τ R ; see below. The population size for the common ancestor θ R was mostly well estimated, although overestimated at very high migration rates. Population sizes for the ancestral species (θ S , θ H ) are harder to estimate; indeed they had larger CIs and were influenced by misspecification of the model of gene Inference of gene flow · https://doi.org/10.1093/molbev/msac237 MBE flow. As expected from the asymptotic results, τ R was very well estimated, except at very high migration rates, in which case τ R was underestimated (and θ R overestimated).
Next we examine the effects of n, S, L, M AB in turn. First, the number of sites (n) had a relatively small impact on MSci parameters, when other factors were fixed (at the basic setting of S = 4, L = 4,000, and M AB = 0.2). When n = 250, CIs for parameters such as the introgression time and probability (τ S and φ) were wide. When n ≥ 1, 000, the CIs were much smaller for all parameters. The introgression time τ S decreased slightly as the sequence became longer. This is consistent with the asymptotic analysis, which suggests that τ * S is dominated by the smallest coalescent time or sequence divergence and should be 0 for the IM and SC models or τ T for the IIM model (when n → ∞) (figs. 2, supplementary figures S2 and S4, Supplementary Material online). Similarly, for the IM and SC models, φ increased with the increase of n when M AB was low, as observed in the asymptotic analysis. Under the IIM model, small datasets with short sequences (n = 250) produced very uncertain estimates of φ and θ H (and, to a lesser extent, τ S and θ S ). The two parameters are nearly confounded; this is discussed below when we examine the impact of the migration rate (M AB ).
Second, we varied the number of sequences per species (S). When one sequence per species is in the data (S = 1), only five parameters in the MSci model ( fig. 1d) are identifiable: θ R , θ S , τ R , τ S , φ). When multiple sequences were sampled per species, all eight parameters are identifiable. They were well estimated when the dataset was large (say, with S ≥ 2 for IM and SC or S ≥ 4 for IIM). Even with S = 4 sequences per species, estimates of φ from data generated under the IIM model involved wide CIs, with τ S being close to τ R , and θ S and θ H being very imprecise as well ( fig. 3). This is due to the semi-unidentifiability or the confounding effects of the parameters, and will be discussed below. Here we note that the problem disappeared and all parameter estimates were well-behaved in large datasets when many sequences were sampled (S ≥ 2 for IM and SC or S ≥ 4 for IIM; fig. 3).
Third, we examined the impact of the number of loci (L). The IIM model was hard to fit in small datasets with a small number of loci (L ≤ 1, 000), generating large CIs for parameters φ and θ H . This is the same pattern as in the case of short loci (n = 250) or few sequences (S ≤ 2), discussed above. In large datasets, the parameters were well estimated. Note that the number of loci L is the sample size in the statistical model as data at different loci are independently and identically distributed. Theory predicts that in large datasets the variance should be proportional to 1/L (see O'Hagan and Forster 2004 for the case of correctly specified models and Yang and Zhu 2018 for the case of misspecified models), and thus the CI should decrease at the rate of L −1/2 . This prediction held for parameters that were well estimated ( fig. 3). As discussed earlier, the introgression time τ S is dominated by the smallest coalescent time or smallest sequence divergence. Thus increasing the number of loci led to a decrease in the estimated introgression time, and the trend was in particular apparent for the IIM model (under which τ S → τ T when L → ∞ if n = ∞). In all cases, the estimated introgression time (τ S ) was closer to the more recent end of the time interval for gene flow than to the midpoint, that is, τ S < τ R /2 for IM, τ S < (τ R + τ T )/2 for IIM, and τ S < τ T /2 for SC (see fig. 1a-c).
Finally, we evaluated the impact of the migration rate (M AB ) ( fig. 3). Under the IM model, there is a near linear relationship between the introgression probability φ and M AB at low rates. The amount of gene flow estimated under the MSci model is less than the true amount expected under the IM model (φ 0 of eq. 10) but the two were close at low rates (with M AB < 0.1, say). At very high rates (with M AB > 1.0, say), divergence time τ R was increasingly underestimated and the population size θ R was overestimated. These patterns are the same as observed in the asymptotic analysis of infinite data (L = ∞), and are due to the attempt of the Under the IIM model, φ involved very large uncertainties at low rates (M AB < 0.04, say), with θ S , τ S , and θ H affected as well. Given the small M AB , why did φ not converge to ∼ 0 with narrow CIs? Note that if M AB = 0 in the IIM model, the MSC model with no gene flow or MSci with φ = 0 will be the correct model. Similarly in figure 3 where M AB = 0.2 was fixed, wide CIs for those parameters were observed in small datasets with short loci (n ≤ 250), few sequences (S ≤ 2), or few loci (L ≤ 1, 000), as noted above. Also in the asymptotic analysis (with L = ∞), we noted that θ * S and φ * were grossly wrong but had no sampling errors because the data size was L = ∞ ( fig. 2; supplementary figures S2 and S4, Supplementary Material online, methods c, d). We suggest that all those results are due to the near unidentifiability of parameters in the MSci model (in particular, θ S and φ); in other words, the parameters are confounded.
If M AB = 0 in the MSC-M model, MSci with φ = 0 will be the correct model, but φ > 0 with a large appropriately adjusted θ S may provide a very good fit to the data (of two sequences per locus). When M AB is small but nonzero, the MSci model will never achieve a perfect fit, and a large φ with appropriately adjusted θ S may provide a better fit than a small φ. Thus in infinite data (L = ∞), we may get grossly wrong estimates with no uncertainty (supplementary figure S2, Supplementary Material online, methods c, d). In finite datasets (L < ∞), there will be a ridge in the posterior surface involving φ and θ S , leading to wide CIs for those parameters, influenced by both model misspecification and the prior ( fig. 3). Including multiple samples from the same species (S > 1) is useful for improving the information content in the data, but strong correlation between φ and θ S may be expected nevertheless. In this regard, the large uncertainties in posterior estimates of parameters may be useful as they help MBE the investigator avoid incorrect inferences of a large φ when gene flow is minimal.

Introgression Events Assigned to Wrong Branches
We conducted simulations to examine the bias in parameter estimates when the introgression event is assigned on either the parental or daughter branch of the lineage genuinely involved in introgression. The data were simulated under model trees A or B and analyzed under models A or B of figure 4a and b.
In the A-A and B-B settings ( fig. 4e), the correct MSci model was assumed, and the performance of the method serves as a reference for comparison. Most parameters, Introgression probability is φ = 0.2. In the IM model C, τ R = 4θ 0 , τ S = 3θ 0 , and τ T = 2θ 0 , with migration occurring from species A to B over time period (0, τ T ) at the rate M = 0.1 migrants per generation. In the IM model D, τ R = 4θ 0 , τ S = 3θ 0 , and τ T = θ 0 , with migration from species A to ST over time period (τ T , τ S ) at the rate M = 0.1. (e) The 95% HPD CIs for parameters in 100 replicate datasets of L = 250, 1,000, and 4,000 loci. Column labels refer to the simulation model followed by the analysis model; for example, "B-A" means that data were simulated under model B and analyzed under model A. The values of θ and τ parameters are multiplied by 10 3 . Black solid line indicates the true value.

MBE
including the species divergence times (τ R , τ S , τ T , and τ X = τ Y ) and population sizes for extant species (θ A , θ B , θ C , θ D ), were well estimated. For well-estimated parameters, the CI width reduced by a half as the number of loci (L) quadrupled, as predicted by theory. Population sizes for ancestral species (θ R , θ S , θ T , θ X , and θ Y ) were less well estimated, although performance improved with sample size: with L = 4, 000 loci, these parameters were well estimated. Introgression probability (φ) was well estimated, but thousands of loci were necessary to obtain precise estimates with narrow CIs under the standard settings used here (four sequences per species per locus and 500 sites per sequence).
In the other settings ( fig. 4e), there was mismatch between the models used to simulate and to analyze data. We note that population sizes for extant species (θ A , θ B , θ C , θ D ) were well estimated, as was the age of the root (τ R ). Performance for estimation of those parameters was very similar whether or not there was model misspecification (e.g., the A-B setting versus the B-B setting and C-A versus A-A). Below we focus on estimation of the other parameters.
In the A-B setting ( fig. 4e), data were simulated under model A with A → B introgression ( fig. 4a) but analyzed under model B with introgression incorrectly assigned to the parental branch ST. Ancestral population sizes θ R and θ S were well estimated, similar to the B-B setting. Divergence times τ R and τ S were well estimated, but τ T and τ X were stuck together. We expect τ (B) T to be mostly determined by the smallest sequence divergence (t bc ) between B and C, which should be close to τ (A) T = 2θ 0 = 0.004. Here, we use the superscript to indicate the model in which the parameter is defined. In the fitting model B, the introgression time τ (B) X (which is >τ (B) T ) should reflect the smallest sequence divergence t ab , whereas in the true model A, t ab is mostly determined by τ X (which is <τ T ). Thus misidentification of the introgression lineage caused τ (B) X to be stuck at τ (B) T ( fig. 5a). There was virtually no information for θ T as the population was estimated to have near-zero time duration with no chance for coalescence. The introgression probability was seriously underestimated, converging to φ * A−B ≈ 0.12 when the number of loci L increases (table 1), whereas the true value was 0.2. This smaller estimate of introgression probability is explained by the distribution of coalescent times between species in the true and fitting models (supplementary fig. S6, Supplementary Material online, true model A). Under the true model A, sequences from A and B are more similar than those between A and C due to the A → B introgression, with an excess of small coalescence time t ab . Under the analysis model B, t ab and t ac have the same distribution. Thus the true model predicts an excess of small t ab , whereas the fitting model predicts an excess of small t ac , and having a smaller φ in the fitting model helps to reduce the discrepancy.
In the B-A setting ( fig. 4e), the simulation model (MSci model B of fig. 4b) assumes introgression involving the ancestral branch ST but the analysis model (model A) assigned introgression to the daughter branch TB. Posterior means and CIs for divergence times τ R and τ T were similar to those in the A-A setting. Note that τ (A) T should be mostly determined by the smallest sequence divergence (t bc ) between B and C, and given that this is τ (B) T = 2θ 0 = 0.002, τ (A) T was well estimated, unaffected by mis-assigned introgression event. Although the true introgression time τ X was 0.003, it was forced to be less than τ T by the analysis model A. As the number of loci increases, τ (A) X became stuck at τ (A) T ( fig. 5b). However, τ (A) S was seriously underestimated. This may be explained as follows. In the analysis model A, τ (A) S was mostly determined by the shortest sequence distance between A and C. In the true model B, this should be close to τ (B) X = 1.5θ 0 = 0.003, due to introgression. With mutational fluctuations in the sequences, one can expect τ (A) S to lie between (τ (B) X , τ (B) S ) = (1.5θ 0 , 3θ 0 ), but closer to τ (B) X in large datasets with many sites and/or many loci. Population sizes θ (A) S and θ (A) Y were affected by the mis-assigned introgression events as well, as those populations are close to the introgression branches. In particular, θ (A) Y was very imprecise as branch YT was very short, and θ (A) S was overestimated because τ (A) S was seriously underestimated (as those two parameters are negatively correlated). Finally, the introgression probability (φ) was underestimated, apparently converging to φ * B−A ≈ 0.02 when the number of loci L increased (table 1), whereas the true value was 0.2. This greatly reduced introgression probability appeared to reflect the very poor fit of the misspecified model A to data generated under model B (see the large differences between the true and fitting distributions of coalescent times in supplementary fig. S6, Supplementary Material online, second row). As τ (A) X and τ (A) S are seriously underestimated, an excess of small coalescent times (t ab , t ac ) is expected in the fitting model A but does not appear in the data, so that having a smaller φ improves the fit.
In summary, assigning introgression events to a wrong parental or daughter branch led to biased estimates of introgression times (causing the introgression events to collapse onto speciation events) and to seriously underestimated introgression probabilities.

Continuous Migration versus Episodic Introgression
In this set of simulations, we generated data under the IM models C and D of figure 4c and d and analyzed them under the MSci models A and B, with the mode of gene flow misspecified and with gene flow assigned to either the correct branch or a wrong branch on the species tree.
In the C-A and D-B settings ( fig. 4e), gene flow occurred continuously but the data were analyzed under the MSci model assuming introgression at a time point. The mode of gene flow was misspecified, but the lineages involved were correctly identified. In the C-A setting, gene flow was between non-sister species, whereas in the D-B setting it was between sister species. Speciation times (τ R , τ S , τ T ) and population sizes (θ) were well estimated, similar to the A-A setting. Surprisingly ancestral population sizes MBE θ T , θ X , θ Y appeared to be even better estimated, with narrower CIs, in the C-A setting than in A-A. Speciation times and population sizes were extremely similar between settings D-B and B-B. Those results were consistent with the results for the case of two species ( fig. 3), which showed that at low migration rates, species divergence times and population sizes were well estimated under the MSci model when the data were generated under the IM model.
In the C-A setting, the estimated introgression time τ (A) X appeared to converge (when L increased) to 0.0011, much more recent than the average time of gene flow (τ (C) T /2 = 0.002), and the introgression probability φ C−A appeared to converge to φ * C−A = 0.12 (table 1), smaller than the expected proportion of total migrants: T /θ (C) B = 0.148. As discussed earlier for the case of two species, the limiting value for τ (A) X was nonzero, as the sequence length is finite, and the MLE φ C−A slightly underestimated the true amount of gene flow. In the D-B setting, the introgression time τ (B) X appeared to converge to 0.0027, larger than τ (D) T = 0.002 but much smaller than the average time of gene flow, 1 2 (τ (D) S + τ (D) T ) = 0.004, and the introgression probability φ D−B appeared to converge to φ * D−B = 0.08 (table 1), much smaller than φ 0 = 0.148 from equation (10). In both the C-A and D-B settings, the estimated introgression time was within the time interval of gene flow, but closer to the time when gene flow stopped, whereas the amount of gene flow was underestimated These patterns are consistent with our analysis of the two-species case at low migration rates (eq. 12, fig. 3), which suggested that gene flow after a period of isolation (the SC model) is easier to recover than gene flow that starts at speciation but stops some time afterwards (the IIM model).
In the C-B and D-A settings ( fig. 4e), the mode of gene flow was misspecified and furthermore gene flow was assigned onto the wrong branch of the species tree. In the C-B setting, divergence time τ T was underestimated slightly, due to gene flow assigned to the wrong branch, as observed in the A-B setting. Ancestral population sizes θ T and θ Y were affected by gene flow, similar to the A-B setting. Model B forces τ X > τ T . Thus we expect τ (B) X and τ (B) T to get stuck together, with both being smaller than   Inference of gene flow · https://doi.org/10.1093/molbev/msac237 MBE τ (C) T = 2θ 0 = 0.004; as the number of loci L increased, τ (B) X appeared to converge to 0.0029, and φ C−B to φ * C−B = 0.10 (table 1).
In the D-A setting, the divergence time τ S was underestimated, due to gene flow assigned to the wrong branch, similarly to the B-A setting. The ancestral population sizes θ R and θ X were well estimated as in the A-A setting, but θ T had a slight positive bias. The ancestral population sizes θ S and θ Y were affected by the gene flow, similar to the B-A setting. The introgression time and probability (τ X and φ) do not exist in the simulation model D. Model A forces τ X < τ T , so we expect τ (A) X to be close to τ (D) T = θ 0 = 0.002; when the number of loci L increased, τ (A) X appeared to converge to 0.00186, and Those results are consistent with our early results for fitting the MSci model to data generated under the migration model in the two-species case (eq. 12, fig. 3), and with the results for the A-B and B-A settings that assignment of gene flow to a wrong branch reduces the estimate of φ.
In summary, the estimated introgression probabilities, at 0.12, 0.08, 0.10, and 0.02 for the C-A, D-B, C-B, and D-A settings, respectively, even though the total amount of gene flow was the same in models C and D (table 1), suggest the following general patterns. First, the MSci model underestimates the total amount of gene flow if gene flow occurs continuously in every generation (i.e., φ C−A < φ 0 ,φ D−B < φ 0 ), as discussed in our analysis of the two-species case. Second, assigning gene-flow events to wrong lineages led to serious underestimation of the amount of gene flow (i.e., φ C−B <φ C−A ,φ D−A <φ D−B ). Third, recent gene flow in the data is more easily recovered (i.e., φ C−A >φ D−B ,φ C−B >φ D−A ).

Isolation with Initial Migration (IIM) Model
Next, we assessed the effects of taxon sampling when the mode of gene flow is misspecified. We used the IIM model for three species of figure 6a to simulate data and analyzed them under the MSci model of figure 6b. Species divergence times (τ R , τ S ) and population sizes (θ A , θ B , θ C , θ R , θ S , and even θ X and θ Y ) were well estimated. We expect the estimated introgression time τ X to converge to τ T = θ 0 = 0.002 if the sequence length is infinite and to a higher limit for finite sequence length. In our simulation τ X ≈ 0.00283 at L = 4, 000 (table 1). The estimated introgression probability (φ) converged to a nonzero limit, ∼ 0.08 (table 1), compared with φ 0 = 0.148 by equation (10).
The IIM model of figure 6a is very similar to the two-species model of figure 1b except that here the tree is larger with more species, and serves to highlight the fact that the impact of the misspecification of the model of gene flow is local. The case is also similar to the D-B setting of figure 4, with the only difference that here the hybridizing species T had only one descendent species sampled in the data, whereas in figure 4 (D-B) it had two descendent species sampled. Thus estimates of parameters such as the introgression probability and introgression time were similar to those in the D-B setting of figure 4 but with wider CIs (table 1). Unlike approximate methods designed to work with species triplets or quartets only, the Bayesian approach accommodates an arbitrary number of species in the data (with arbitrary data configurations at each locus), so that the difference in taxon sampling has only the effect of affecting the information content in the data.

Ghost Species
We considered two scenarios in which a species that contributed migrants to extant species has gone extinct or is otherwise unsampled in the data. Note that existence of extinct or unsampled species that received genetic materials from ancestors of extant species in the sample is not relevant to the analysis of the sampled data and does not constitute a model misspecification. In the first scenario, model A ′ of figure 7a ′ is used to simulate data, which assumes that species XUV contributed migrants to species B but is not included in the sample. Note that this model is equivalent to model A of figure 7a. When we fit model B ( fig. 7b), the only incorrect assumption is the constraint that τ X = τ Y . This is a minor misspecification. Indeed all parameters shared between the simulation model and the analysis model were well estimated ( fig. 7c). The estimates of introgression time, τ X =τ Y ≈ 0.0028 (table 1), were close to the average of the two parameters in the true model (0.0025). Introgression probability φ ≈ 0.21 (table 1) was also close to the true value (0.2). The existence of the ghost species (XUV) had very little effect on the inference.
In the second scenario ( fig. 8a), the true model assumes continuous migration involving intermediate ancestral species that have gone extinct, and the MSci model ( fig. 8b) was fitted to data sampled from extant species. Divergence times τ R and τ T were very well estimated, as were the population sizes shared between the simulation and analysis models (θ A , θ B , θ C , θ R ). We expect τ T in model B to be dominated by the minimum coalescent time t ab between sequences from A and B, and this is given by T . Gene flow from branches RC to SU over the time interval (τ U , τ S ) and then from SU to TB during (τ U , τ T ) was interpreted as introgression in the MSci model. The effective rate for this migration may be close to M CU M UB = 0.04, fig. 8c, table 1). The introgression time τ X should be between τ U = θ 0 = 0.002 and τ T = 2θ 0 = 0.004 and the estimate was ≈ 0.0030 ( fig. 8c, table 1). Note that both θ T and θ Y were overestimated ( fig. 8c). Branch T of figure 8b corresponds to branches RS and ST of figure 8a, with population size θ 0 = 0.002. Branch Y corresponds to a segment of branch TB over the time interval (τ U , τ T ), with θ 1 = 0.01. Overestimation of θ Y (and θ T ) may be because there is a deficit of t bb over the interval (τ U , τ S ) due to gene flow, and the fitting MSci model, with the amount of gene flow underestimated (φ < φ 0 ), used large θ Y and θ T to compensate.

Discussion
The Mode of Gene Flow and the Utility of Misspecified Introgression Models The asymptotic theory, even though based on only two species with one sequence sampled per species per locus, has been very useful. It generated a number of insights that were confirmed and extended in our simulation. Together the theory and simulation suggest the following correspondence between the MSC-M and MSci models. When gene flow occurs continuously over an extended time period after divergence of two species and we fit the introgression model, the estimated introgression time tends to be closer to the more recent end of the time period of gene flow, because the introgression time is dominated by the most recent coalescent time or the minimum sequence divergence between species. If the true coalescent time is known and used as data, the introgression time will converge to the time when gene flow stopped. At low migration rates (M < 1 4 , say), the species divergence time is correctly estimated by the MSci model, and the introgression probability φ is lower than but close to the expected proportion of migrants (φ * < φ 0 ). The estimate is particularly close under the secondary-contact model (supplementary fig. S4, Supplementary Material online). At very high migration rates, the estimated introgression probability φ * may be much less than φ 0 , and furthermore the species divergence time is underestimated to account for intermediate coalescent times generated under the MSC-M model. Recent gene flow (as in the SC model) is easier to recover (with φ * closer to φ 0 ) than ancient gene flow (as in the IIM model).
The accurate estimation of species divergence times under the MSci model despite the misspecification, at least at lower migration rates (e.g., τ R for M ≤ 0.3 in fig. 3), may be worth emphasizing. It is well known that ignoring gene flow between two species may lead to serious underestimation of the species divergence time. Here our results suggest that if gene flow is continuous, the MSci model assuming introgression at a fixed time point still gives reliable estimates of the species divergence time. The estimated introgression probability (φ) may also serve as a useful guide even though it reflects both the migration rate per generation (m or M) and the time duration of the period of gene flow (eq. 10). Even if gene flow occurs continuously over time (so that the migration model is a more realistic model), the MSci model is effective in extracting historical information about species divergence times and

MBE
population sizes. Note that on the evolutionary time scale, a few hundred or thousand generations may count as a fixed time point, in which case the MSci model may provide an adequate approximation. Both the asymptotic theory and simulation have highlighted the semi-unidentifiability or confounding effects between the introgression probability (φ) and the population size of the donor species (θ S in fig. 1d) (e.g., fig. 2, methods c and d). The problem is particularly acute under the IIM model applied to small datasets (with short loci, few sequences per species, or few loci), where high estimates of φ with wide CIs are produced even though migration occurs at very low rates ( fig. 3). One such case has been observed in a recent analysis of genomic data from the erato group of Heliconius butterflies (Thawornwattana et al. 2022). The estimated H. sara → H. demeter introgression probability was high with wide CIs for some chromosomal regions with a small number of loci (e.g., chromosome 21 with 4350 noncoding and 3628 coding loci, and an inversion on chromosome 15 with 149 noncoding and 167 coding loci), with the introgression time close to the species divergence time, whereas for the other large chromosomes, the estimates were nearly zero (φ < 0.01). The true rate in this case appeared to be φ ≈ 0, but the limited data from small chromosomal segments led to poorly supported large introgression rates, as in our simulations ( fig. 3).
We demonstrated that including multiple samples from the same species (in particular, from recipient species) is important to resolving unidentifiability issues or confounding effects, as well as boosting up the information content concerning the rate of gene flow in the data. In this regard, it may be noted that many approximate methods are designed to use only one sample per species, and it has been claimed that "adding more samples provides little new information with respect to introgression" (Hibbins and Hahn 2022). We suggest that this may not be a generally correct statement.
Overall, our simulations using larger species trees with more than two species suggest that misspecification of the mode of gene flow (continuous migration versus episodic hybridization/introgression) has relatively small and localized effects, restricted to divergence times and population sizes around the lineages involved in gene flow, while species divergence times, population sizes for extant species and for ancestral species not involved in gene flow are largely unaffected. If gene flow occurs between species A and B but more distantly related species are included in the data sample, parameters outside the AB clade are largely unaffected (e.g., compare results for the IIM model for two species of fig. 3 with those for three species of fig. 6). Similarly, if A represents a clade rather than one species, divergence times and population sizes inside the A clade are not affected by gene flow involving the branch ancestral to the A clade (e.g., compare the D-B setting of fig. 4 with the IIM model of fig. 3).
Assigning gene flow to parental or daughter branches causes the introgression probability to be underestimated, and the introgression time to collapse onto the species divergence time. This result may be used to diagnose the mis-assignment of introgression lineages in real data analysis (Ji et al. 2022). A number of authors have discussed the impact of ghost species on detection of betweenspecies gene flow (Beerli 2004;Ottenburghs 2020). Tricou et al. (2022) used simulations to demonstrate that D-statistics can be misled to detect false signals of introgression when the model involved an unsampled (ghost) species. In our simulations, the impact of ghost species on Bayesian estimation of introgression rate and time was minor provided we considered the rate of gene flow in the migration and introgression models to reflect both indirect gene flow via intermediate species and direct gene flow.

Testing Models of Gene Flow
In this study, we fixed the model of introgression in our analyses, with all introgression events pre-identified, to examine the effects of model misspecification. One may ask what happens if different introgression models (which for example assign introgression events onto different branches of the species tree) are compared using genomic data. Currently, both *BEAST and PHYLONET have implemented cross-model MCMC algorithms under the MSci model, which insert and delete introgression events on the species tree, allowing the Markov chain to move between models. Those algorithms are computationally expensive and currently the two programs can handle only very small datasets (with <100 loci, say). In the BPP program, one may use the Bayes factor to compare two MSci models, using thermodynamic integration (Gelman and Meng 1998;Lartillot and Philippe 2006) combined with Gaussian quadrature to calculate the marginal likelihood values (Rannala and Yang 2017). In the case where the compared models are nested (e.g., one with introgression and another without), the Bayes factor may also be calculated through the Savage-Dickey density ratio (Dickey 1971), which uses only a within-model MCMC run under the more general model (Ji et al. 2022). This has a computational advantage over reversible jump MCMC (Green 1995), and has recently been applied to formulate and compare introgression models in an analysis of genomic data from the Tamias quadrivittatus group of North American chipmunks (Ji et al. 2022). Calculation of marginal likelihood values or Bayes factors may be feasible if we have only a small number of well-specified models but may not be feasible for searching in the space of MSci models for a given set of species.
Approximate methods have also been developed to infer introgression events or the so-called phylogenetic networks using summaries of the multilocus sequence data. For example, estimated gene tree topologies may be MBE treated as data, as in PHYLONET/GT (Wen et al. 2016). Some methods are designed to detect gene flow in a small tree with three or four species, including summary methods based on genome-wide site-pattern counts (such as D and HYDE discussed earlier) or on estimated gene trees (e.g., SNAQ) and maximum likelihood applied to multilocus sequence alignments (e.g., 3S, Zhu and Yang 2012;Dalquen et al. 2017). Results for species subsets may then be combined to formulate an introgression model on the large tree for all species, which is a challenging task (Edelman et al. 2019;Thawornwattana et al. 2022). In summary, there is currently an acute need for improving the computational efficiency of Bayesian MCMC algorithms for inference under the MSC model with gene flow and the statistical efficiency of approximate methods.
It will also be interesting to use the same genomic data to compare the MSC-M and MSci models. The two classes of models often predict very different distributions of gene trees and coalescent times (e.g., supplementary figs. S1, S3, S5, Supplementary Material online; see also . Thus, genomic data may be informative to distinguish them. A stochastic search in the combined space of MSC-M and MSci models may be infeasible, as the two types of models are very different. However, they can be compared using Bayes factors.

Simulation to Establish a Correspondence between the Migration and Introgression Models in the Case of Two Species
We analyzed the relationships between parameters when data are generated under the continuous migration model (IM, IIM, and SC; fig. 1a-c) and analyzed under the episodic introgression (MSci) model ( fig. 1d). Our theory assumed an infinite number of loci (L = ∞), a finite number of sites per sequence (n), with only one sequence per species per locus. We conducted computer simulations to augment the theoretical analysis. Data of multilocus sequence alignments were simulated under the IM, IIM, and SC models of figure 1a-c, and analyzed under the MSci model ( fig. 1d). Population sizes on the species tree ( fig. 1) were θ 0 = 0.002 for the thin branches and θ 1 = 0.01 for the thick branches. Migration occurred from species A to B after their divergence at τ R = θ 0 in the IM model, between τ R = 2θ 0 and τ T = θ 0 in the IIM model, and between τ T = θ 0 and the  1A in Flouri et al. 2020) assumes that τ X > τ Y and τ T > τ Y and can represent scenario (a ′ ) in which species X split into two species (A and U), and species XUV contributed migrants into species TYB at time τ Y but has since become extinct. This model was used to simulate data, with θ 0 = 0.002 for the thin branches and θ 1 = 0.01 for the thick branches, τ R = 4θ 0 , τ S = 3θ 0 , τ T = 2θ 0 , τ X = 1.5θ 0 , and τ Y = θ 0 . The introgression probability is φ = 0.2. The number of sequences is S = 4, and the sequence length is n = 500. (b) MSci model B ( fig. 1B in Flouri et al. 2020) used to analyze the data, which incorrectly assumes τ X = τ Y . (c) The 95% HPD CIs for parameters, with θs and τs multiplied by 10 3 and black solid line indicating the true value.

MBE
continuous in the true model but the MSci model assumes episodic introgression at a particular time point, so that the mode of gene flow is misspecified. In settings C-B and D-A, the mode of gene flow was similarly misspecified but we had in addition mis-assignment of gene flow to wrong branches on the species tree. Other parameter settings were the same as above. With two trees, three numbers of loci (L), a total of 600 datasets were generated, each analyzed twice (under models A and B).

Isolation with Initial Migration (IIM) Model
Data were simulated under the IIM model A of figure 6a, with A → B migration over the time period (τ T , τ S ), and analyzed under the MSci model of figure 6b, assuming introgression at time τ X = τ Y . The IIM model was specified using a ghost species (U) from which no sequences were available. We generated 100 replicate datasets, each of L = 250, 1,000, or 4,000 loci, with a total of 300 datasets simulated. MCMC settings were the same as above.

Ghost Species
To assess the effects of unsampled ghost population, we simulated data under MSci model A ′ (see fig. 1A in Flouri et al. 2020) of figure 7a ′ and analyzed them under the MSci model B of figure 7b, with τ X = τ Y incorrectly assumed. Here introgression involved a ghost species XUV which went extinct or was otherwise unsampled in the data. This scenario is equivalent to model A of figure 7a. With the three values for L (250, 1,000, 4,000), 300 datasets were generated, all analyzed under the MSci model ( fig. 7b).
We also used the IIM model of figure 8a to generate data, with migration from species RC to SU and from SU to TB, and with V and W to be unsampled ghost species. Data (i.e., sequences from A, B and C) were analyzed under the MSci model of figure 8b. We used three values for L (250, 1,000, 4,000) and 100 replicates, with 300 datasets simulated in total. Other settings were the same as above.

Supplementary Material
Supplementary data are available at Molecular Biology and Evolution online.